Elasto-viscoplastic Spreading: from Plastocapillarity to Elastocapillarity

We study the spreading of elastoviscoplastic (EVP) droplets under surface tension effects. The non- Newtonian material flows like a viscoelastic liquid above the yield stress and behaves like a viscoelastic solid below it. Hence, the droplet initially flows under surface tension forces but eventually reaches a final equilibrium shape when the stress everywhere inside the droplet falls below the resisting rheological stresses. We use numerical simulations and combine Volume-of-Fluid (VOF) method and an EVP constitutive model to systematically study the dynamics of spreading and the final shape of the droplets. The spreading process examined in this study finds applications in coating, droplet-based inkjet printing, and 3D printing, where complex fluids such as paints, thermoplastic filaments, or bio-inks are deposited onto surfaces. Additionally, the computational framework enables the study of a wide range of multiphase interfacial phenomena, from elastocapillarity to plastocapillarity.


I. INTRODUCTION
The deposition and spreading of fluids over surfaces occur in a wide range of industrial applications, such as spray coating and printing [1][2][3][4][5].The fluids used in these applications often contain microscopic constituents like polymers or colloids, resulting in highly nonlinear macroscopic rheological features like elasticity, plasticity, and shear-dependent viscosity.Hence, understanding the rheological effects on the fluid mechanics of spreading is essential for further optimizing current systems or designing new materials for specific applications [6].In the past few years, many studies have investigated the impact and spreading of droplets on surfaces for fluids with various rheological properties, including viscous fluids [7][8][9][10], viscoelastic fluids [11][12][13][14][15], and yield stress materials [16][17][18][19][20][21][22][23][24][25].However, systematic experiments are generally complicated due to the high-dimensional interconnected parameter space.To this end, computer simulations can be employed to investigate the significance of rheological parameters independently from one another.Our focus will be on a general class of nonlinear materials named elastoviscoplastic (EVP) fluids which can exhibit both elastic and plastic properties, as well as viscous behaviour.Many models of EVP fluids have been proposed previously.Saramito [26,27] presented a description of EVP constitutive equations which combine classical viscoplastic models (Bingham or Herschel-Bulkley) with viscoelastic models (Oldroyd-B or PTT).Adopting a different perspective, de Souza Mendes [28] used a generalized viscoelastic model where material properties are functions of the strain rate in order to incorporate the yielding behaviour as well as the thixotropic effects.Dimitriou and McKinley [29] included isotropic and kinematic hardening in their EVP formulation to capture various steady and unsteady flow responses.Recent reviews of the development and thorough comparison of different EVP models can be found in Fraggedakis et al. [30] and Saramito and Wachs [31], where they comment both on the mathematical/physical properties of the models, as well as the numerical characteristics.Moreover, many recent studies have incorporated the models above into computational frameworks to study various fluid dynamics problems such as flow around particles [32][33][34], flow in channels [35][36][37], flow in porous media [38,39], and flow around bubbles [40,41].More related to the present investigations are the computational studies on droplets with EVP properties.Oishi et al. [20] studied the flow of materials on an inclined plane, considering elastic properties, where they captured the so-called avalanche effect where a decrease in viscosity (triggered by a stress field) induces a motion that successfully creates another decrease in viscosity.In follow-up studies, Oishi et al. [22,42] have also included surface tension and thixotropic effects on the impact of EVP droplets on normal and inclined solid surfaces.Izbassarov and Tamissola [43] investigated how an EVP droplet deforms inside a Newtonian medium under simple shear and observe that the droplet deformation can present a non-monotonic behaviour as elasticity is changed.
In the present work, our main goal is to study the spreading of EVP droplets on a surface.To this end, we aim to extend the previous numerical analysis [44] to consider not only viscoplastic rheology, but also the addition of elasticity.By using Saramito's EVP model [26], a parametric study will be carried over in order to understand how the addition of elasticity can affect the transient spreading dynamics and also the final shapes of elastoviscoplastic droplets.
The paper is organized as follows.Section 2 presents a description of the problem, the governing equations used to model the flow and the rheology of the EVP material, and also an overview of the numerical method used in this work.In section 3, results are presented and various limits of the problem are visited.Section 4 concludes the results and presents future perspectives.Additional numerical details can be found in the appendices.

A. Problem description: capillary spreading of a droplet
We consider the spreading of an axisymmetric EVP droplet on a wetted surface.The choice of geometry has two main reasons.Firstly, in several applications, including 3D printing and spray coating, it is typical for droplets to spread and deposit on an existing layer of the same fluid [9,10].Secondly, opting for a geometry that eliminates the presence of a triple contact line provides theoretical and computational advantages as it simplifies the associated complex physics related to boundary conditions and stress singularities.
As illustrated in Figure 1a, a droplet is initially placed over a thin film of the same material and allowed to spread due to capillary forces.Viscous forces oppose this spreading by slowing it down, and yield-stress is capable of stopping it completely as shown previously in [44] for a purely viscoplastic scenario.The role of elasticity in this process, however, is less clear and a numerical framework will be implemented in order to understand such effects.To this end, we will model the rheology of the EVP material by the Saramito model [26,27], which generalizes both the Bingham viscoplastic and the Oldroyd-B viscoelastic models.Figure 1b shows a mechanical analog of this model.Below the yield stress the material behaves like a Kelvin-Voigt solid, while after yielding it flows as an Oldroyd-like fluid.Fig. 1: a) Sketch of the geometry of a spreading elastoviscoplastic droplet.We assume the droplet to be axisymmetric around the z-axis and spread on a thin pre-wetted film of the same material with thickness h∞.The radius and the height of the droplet are denoted respectively as R(t) and H(t).Gravity is neglected, hence, the spreading is purely due to surface tension effects.The droplet includes microscopic constituents, resulting in macroscopic EVP models.b) Mechanical analog of Saramito's EVP model.

B. Governing equations
The governing equations for the isothermal incompressible bi-phase flow, are the continuity and momentum conservation given by where u and p are the velocity and pressure fields.The gravitational force is defined as f g = ρg where ρ is the fluid mass density.In our numerical method, the surface tension force is also defined as a body force f σ = σκδ s n, where κ is the curvature of the interface, σ the constant surface tension coefficient, n is the unit vector normal to the interface, and δ s is the Dirac delta function centered on the interface [45].
The deviatoric stress tensor τ is the sum of the solvent τ s and polymeric τ p contributions, The solvent stress contribution presents a Newtonian behaviour and the polymeric stress includes the memory effects, hence, where µ s is the solvent viscosity and D = 1 2 ∇u + (∇u) T is the strain rate tensor.The polymer stress τ p includes the elasto-viscoplastic contribution, emerging from the microstructures.For viscoelastic fluids, it is common to model the contribution of the polymeric stress as τ p = µp λ f (A), where λ is the relaxation time, and f (A) is a strain function of the conformation tensor A [46][47][48].We assume f (A) = (A−I) and A follows a liner relaxation law , where ▽ (•) is the upper-convected time derivative.Adding a critical yield stress at which the material switches between liquid and solid phase, we arrive at a constitute law, first formulated by Saramito [26,27].The model combines the Bingham (viscoplastic) [49,50] and the upper convected Maxwell (viscoelastic) models [46-48, 51, 52].The polymeric stress then is governed by, where, τ 0 is the yield stress, ∥τ p ∥ = tr(τ p2 ) , and A dimensionless version of these equations can be obtained by scaling the variables as follows where x is the position vector, t is time, U = σ ρ d L is the characteristic velocity, ρ d is the droplet density, and the length scale is L = [3V/(4π)] 1/3 with V being the volume of the droplet (L can be seen as the radius of a corresponding spherical droplet with same volume V).
Removing, for convenience, the bars in (7), the dimensionless governing equations for the droplet phase are given by with the dimensionless groups: Ohnesorge numbers (Oh s and Oh p ), Bond number (Bo), plastocapillary number (J ), and Deborah number (De) defined respectively as The Ohnesorge numbers compare the characteristic Rayleigh timescale and the visco-capillary time scale and, for this problem, could be seen as normalized solvent or polymeric viscosities.The Bond number compares the gravitational stresses and capillary pressure.In the present study, we focus on pure capillary spreading, where gravity is negligible, hence Bo = 0 in all simulations.Although we admit that the limit of large Bo and negligible surface tension is interesting and important for geophysical problems like lava flows, mudslides, and snow avalanches where large-scale EVP fluids flow under gravity [53,54].The plastocapillary number compares the yield stress and the capillary stresses, and the Deborah number is the ratio of the polymeric relaxation time to the characteristic time-scale of the problem.Note that many authors choose to replace parameters Oh s and Oh p by Oh = Oh s +Oh p and β = Oh s /(Oh s +Oh p ) [55], where β can be seen as the normalized ratio of solvent to apparent viscosity.We also note that, using the present multiphase method, we also have fluid motion in the air phase.Therefore, a set of equations similar to equations ( 8)-( 10), but with Newtonian rheology, is also solved for the air flow.Consequently, two other nondimensional parameters are also relevant: the density ratio between droplet and air (ρ d /ρ a ) and viscosity ratio (µ s /µ a ), where subscript a denotes air properties.Both of these ratios are kept constant with a value of 100.
[Case 2: De = 0 -Plastocapillarity Regime] Excluding viscoelastic (memory) effects results in a purely viscoplastic (Bingham) free surface flow under capillary stresses [61].Note that, substituting De = 0 in equation (10), from equation ( 9), we arrive at the deviatoric stress of for yielded regions.For a Bingham fluid, D = 0 when the material is unyielded.Hence, numerical regularization is required.In practice, for this limit, instead of the EVP constitutive model, we will solve a regularized version of the generalized Newtonian model above [44,62].
[Case 3: J = 0 and De = 0 -Newtonian Regime I] Without elasto-plastic rheology, we simply have a Newtonian fluid spreading due to surface tension, and the deviatoric stress tensor is In other words, the polymers behave like passive scalars and do not contribute to the dynamics of the problem.This regime might not be physical because, in reality, the viscous and plastic dissipations might also depend on the relaxation time.Note, if Oh p /De (∼ elastic module) remained finite, then equation (10) converged to the constitutive model of an elastic solid [48].
[Case 5: J → ∞ -Elastocapillarity Regime] For a finite De, when the plastocapillary number is large, equation ( 10) reduces to i.e., the stress inside the material remains below the yield stress and the polymeric response is elastic.In fact, Oh p /De = µL/σλ is the elastocapillary number (or the inverse of it [56]).In this limit, the droplet behaves like a Kelvin-Voigt solid that spreads under capillary stress [63][64][65].This regime is highly related to other soft wetting phenomena, where capillarity forces deform soft solids [66][67][68].
In all of the cases above, the system of equations is closed with appropriate boundary conditions.We apply a no-slip condition on the solid wall.Rotational symmetry is applied at the center of the droplet and outflow boundary conditions are used at the two boundaries distant from the droplet.

C. Initial condition & numerical method
We used the open source code Basilisk C to solve the equations described in the previous section.An overview of the numerical procedure will be given in this section, but detailed descriptions of Basilisk can be viewed in [69] (also see appendices A and B).
The simulation is setup by initializing a droplet according to the following shape which represents a half-parabola in an axisymmetric coordinate system.The droplet is placed in a squared domain with dimensions [0, 5R 0 ] × [0, 5R 0 ] which is fully discretized with a non-uniform quadtree grid [70,71].
To accurately resolve the flow structure inside the droplet and its shape, we apply increased refinement levels for the liquid phase and also at the interface (see appendix A).As the droplet spreads over time, the mesh is also adapted so that the refined region follows the interface.
The interface is represented implicitly by the Volume of Fluid (VOF) scheme [72], in which each mesh cell stores a value representing the fraction of droplet fluid.Density and viscosity are then locally determined based on the volume fraction c(x, t) according to where ρ and µ represent the density and dynamic viscosity of a fluid, respectively.The subscripts d and a represent, the droplet and the air, respectively.
This volume fraction field c is then advected over time by solving the equation The numerical code then solves the governing equations using a projection method and a multilevel Poisson solver (see [71,73] for more details of the VOF implementation).

III. RESULTS AND DISCUSSION
We construct the discussion on EVP spreading by first analyzing the results of pure viscoplastic and viscoelastic droplets.In all simulations in this section, the following parameters are fixed: Oh s = 1/90, Oh p = 8/90, Bo = 0.

A. Plastocapillarity
Theoretical, computational, and experimental results for viscoplastic droplet spreading (De = 0) have been previously presented [25,44].We will revisit this limit first to validate the simulation and also to extend the available results for higher plastocapillary numbers.To this end, we will focus on the properties of the final shape of the droplet when t → ∞.Unlike Newtonian droplets, droplets of yield-stress fluids will be arrested at this limit, resulting in a finite final radius R f and height H f .In a pure viscoplastic case, these features can be explained by balancing the capillary stress and yield stress [44], resulting in theoretical scaling laws: . Note, the pre-factors are obtained from asymptotic analysis [44].Also note, to test our numerical results against these laws, we must choose a stoppage criteria in our simulations since we use a regularized model in the viscoplastic limit.We do this based on the nondimensional kinetic energy ( ) of the droplet, such that the simulation is stopped when E k < 10 −6 .
Figure 2 shows the final radius and height of our simulated droplets versus the plastocapillary number J .As expected, due to the higher influence of yield stress on capillary-driven spreading, the final radius decreases with J , while the final height increases.Both final radius and height reach a plateau after a certain value of plastocapillary number, J ≈ 1 (see Video I).At this limit, the droplets practically get stuck at their initial shape, since the capillary stress is not strong enough to overcome the yield-stress at all.There exists a transition regime between the low yield-stress scaling laws (soft) and the high yield-stress plateau (stiff) regimes, where the numerical results smoothly vary between the two.There is a good agreement between the theory and numerical results when the droplets are soft enough to yield.As plastocapillary number increases, the droplet does not entirely yield and hence break the assumption in the theoretical predictions.Note that the exact point of transition from soft to stiff limit depends on the initial conditions.This could have important implications for some applications such as 3D printing [25,[74][75][76][77][78], when depending on the values of J , the history (shape) of droplet/filament at the time of deposition can influence the final geometry of the print.

B. Visco-elastocapillarity
We continue by considering the case of viscoelastic materials without plasticity, that is, J = 0.In this situation, the constitutive model reduces to an Oldroyd-B fluid.Figure 3 shows the droplet radius (top) and height (bottom) over time for different values of De (see Video II).In the Newtonian limit (De = 0), the spreading eventually follows the rate predicted by Tanner's law [7], i.e., R ∝ t 1/10 and H ∝ t −1/5 .For the viscoelastic case, we note that, in the first moments, the droplet spreads considerably more as we increase De.We anticipate this is due to the increased relaxation time of the fluid.As we increase the relaxation time (or De), the stresses take a longer time to develop as the flow field develops inside the droplet, consequently, the droplet spreads more since the internal stress is smaller during this transient period.As a consequence, interestingly, the spreading curves converge to an apparent Newtonian limit when De → ∞.In this limit, the polymeric stress does not have enough time to develop at all within the timescale of the simulation, and we actually recover a Newtonian droplet that only exhibits the solvent stress, i.e., we have a Newtonian fluid with Ohnesorge number Oh ∞ = Oh s (circles in figure 3).For the intermediate values of De, the interface experiences an oscillatory behaviour (see figure 3), where the droplet height first reaches a local minimum and then increases again as elastic stresses build up, eventually reaching a decaying regime.To further analyse the anatomy of viscoelastic spreading, we inspect the flow field inside the droplet using a flow parameter [79]: where D and Ω are the deformation rate and vorticity tensors (the symmetric and antisymmetric components of ∇u), respectively: D = 1/2 ∇u + (∇u) T and Ω = 1/2 ∇u − (∇u) T .The flow parameter can vary in the range ξ ∈ [−1, 1], where ξ = 1 represents a purely extensional flow, ξ = 0 shear flow and ξ = −1 indicates solid-like rotation.Figure 4 shows the value of ξ inside a viscoelastic droplet with De = 0.245 for different time frames.At the beginning of spreading, the flow close to the axi-symmetry axis and in the bulk of the droplet is predominantly extensional (red color).Meanwhile, a shear-dominant boundary layer forms from the contact line and across the substrate (zone I figure 4a).Two rotating regions also develop close to the interface and grow over time (zones II and III in figure 4b) but decay as the droplet further spreads.Eventually, the flow inside the droplet is mainly extensional in the centre (zone IV in figure 4d) and shear-dominant close to the substrate.
In the next step, we will study the additional effects of plasticity on the complex flow structures shown above in an EVP spreading.

C. Elastoviscoplastic spreading
We now consider the general scenario of elastoviscoplastic spreading, i.e, De and J are both finite.In this case, all the nonlinear properties mentioned earlier (including the coexistence of solid and liquid states in plastic fluids and the time-dependent characteristics of viscoelastic fluids) are simultaneously observed.We inspect the nondimensional polymeric stress τ p and observe in which regions this stress is above (yielded) or below (unyielded) the value of J (see [43]). Figure 5 shows the value of scalar S = log(∥τ p ∥) − log (J ) at different times inside a droplet with J = 0.18 and De = 0.816.S > 0 means the material is fluidized and flows like a viscoelastic fluid, and S < 0 means the material is not yielded and behaves like a viscoelastic solid.As expected, the droplet is initially mostly unyielded (large blue region) since we assume there is no internal stress as our initial condition.The stress begins to increase from the contact line, which is the location of the highest curvature (zone I in figure 5a).After some time, most of the droplet is yielded (red colors) with a particularly higher stress region near the wall and droplet edge (zone II in figure 5b).Meanwhile, a moving plug forms near the interface of the droplet (zone III in figure 5b), correlated with the rotating regions shown in figure 4. As time continues to increase, the droplet solidifies once again with a static plug in the centre (zone IV in figure 5c), leading to a full stoppage.Note that, unlike the pure viscoplastic case, no regularization is needed in the case of EVP materials; the viscoelastic solid will reach an equilibrium state as long as stress everywhere inside the droplet is below the yield stress.
We systematically extended the analysis above by changing the control non-dimensional parameters.Figure 6 0.01  demonstrates the value of S for different combinations of J and De (see Video III).For a given De, increasing the value of J leads to larger unyielded regions inside the droplet, which is expected as we increase the material yield-stress.Hence, for a given De, and at a given time, the droplet spreads less as the values of J increase.For a given J , increasing the Deborah number also generates larger unyielded regions (in blue).Similarly to the explanation in section B, increasing De results in smaller values of τ p and more regions below the yield-stress.Intriguingly, at the high J cases (last row in figure 6), it is notable to observe that by increasing De, the droplet spreads significantly even though it almost entirely behaves like a solid.This is because the current EVP formulation models the unyielded case as a viscolastic solid, allowing for a finite deformation rate (as opposed to a Bingham rigid solid where the deformation rate is zero for the solid state).This means that the material can experience deformation even below the yield stress, as seen clearly in the examples when both J and De are high (elastocapillarity regime).
For a better quantitative analysis of the dynamics shown above, we look at the variation of droplet height over time for different values of J and De, as shown in Figure 7.For EVP droplets the spreading stops and a final shape is reached.This is similar to the plastocapillarity regime, however, the elasticity clearly influences the final shape.Similarities to the Visco-elastocapillarity regime can also be seen, particularly in the local overshoot and the oscillation of the height that is observed for most cases and amplifies with De.At the same time, the impact of De on the dynamics is more pronounced for higher values of J , resulting in a larger difference in the final shape (height and radius).Blue regions indicate stress below the plastocapillary number J (unyielded) and red regions are above J (yielded).We note that the colorbar limits do not include the total range of values present in the data, since we are only interested in visualizing regions that are above or below zero .
Finally, we quantify the Deborah number effects on the final radius of the droplet to complete the picture shown in section A. Figure 8 shows the final radius of droplet versus J for different values of De.The limit De = 0 corresponds to a limit, also shown in figure 2 and the dashed lines represent the theoretical scaling for this case.As De increases, the droplet spreads more, particularly for higher values of J .For higher De values, we also see a decrease in the critical value of J from which the spreading stays almost constant.This happens because, as we saw earlier, the elasticity promotes unyielded droplets, and for droplets that are already unyielded, further increasing the value of J does not introduce further changes.For the elastocapillary limit, when De and J are large, balancing the surface tension and elastic forces results in scaling laws for the final radius and height of the droplet.At the stoppage moment, the surface tension force can be estimated as F σ = σH 2 f /R f .Balancing this with the elastic forces, estimated as F e = (µ p /λ)R 2 f , and a for a given volume of the droplet . These formulations are similar to those for the plastocapillarity limit, except that the plastocapillary number (yield stress) is replaced by elastocapillary number (elastic module).In appendix C, we test these scaling laws by looking at the final radius and height as a function of De.

IV. CONCLUSIONS
We numerically investigated the spreading of elastoviscoplastic fluids under surface tension forces.Direct numerical simulations are performed using the EVP model of Saramito to understand how different non-dimensional groups influence the spreading.We focused our study on the importance of two non-dimensional groups: the plastocapillary number, J , and the Deborah number, De.We confirm that increasing the plastocapillary number in EVP spreading, reduces the droplet spreading, as previously explained for the viscoplastic limit [25,44].
The influence of elasticity on spreading can also be significant.Increasing the Deborah number (while other parameters are fixed) promotes more spreading.This effect is associated with the polymeric relaxation time, which delays the development of the internal droplet stresses, allowing the droplet to spread with less resistance.
Overall, for fixed (solvent and polymeric) Ohnesorge numbers, the De − J parameter space covers a range of regimes, including visco-elastocapillarity, elastocapillarity, and plastocapillarity, which can be further studied by the model presented here.Figure 9 presents a (schematic) overview of these limits.The spreading of EVP fluids plays a key role in many industrial applications, such as coating and 3D printing.Our work, based on a continuum model, sheds light on how elastic and plastic rheological properties alter viscous spreading.The study can be further extended in many different ways.Firstly, further experimental studies are required.The present experimental data [44] are limited to viscoplastic limits and are not sufficient to compare with the simulations here.Ideally, experiments in which the values of plastocapillary and the Deborah numbers can be systematically changed should be performed.This, however, could be a difficult task as chemical or physical changes in microstructure often change the yield stress, viscosity, and elastic moduli of the material, simultaneously.
The present article focuses on the problem of spreading.Still, in principle, the computational framework presented here could be used to study and analyze a range of capillary-driven phenomena such as bubble or drop coalescence [80][81][82][83][84], pinch-off [85,86], multi-component systems like drops on liquid-infused surfaces [87,88] or soft wetting [66,68,89,90], and non-axisymmetric shapes [75].Finally, the model can be extended for materials with more complicated rheological properties such as thixotropy [20,91,92] and also for wetting of dry surfaces, when a triple line exists [23].

A. Viscoplasticity and viscoelasticity
We use a generalized Newtonian Bingham model when De = 0, i.e, a viscoplastic material.In this case, the apparent viscosity is replaced by a regularized Bingham-Papanastasiou viscosity [93] where ϵ vp is the regularization parameter (see [94]). Figure 10 shows an example of convergence tests for the viscoplastic regularization parameter when J = 0.18, Oh s + Oh p = 0.1, De = 0, Bo = 0 and the interface is captured at time t = 10.Performing a series of theses tests, we note that all solutions have converged for values ϵ vp ≤ 10 −3 .Therefore, in all viscoplastic simulations reported in this article, we fixed ϵ vp = 10 −3 .0 0.5 1.0 1.5 2.0 0 0.5 1.0 1.5 2.0 Fig. 11: Non-uniform grid generated for the spreading simulation.The right panel shows the magnified view of the white box in the left.Cells are colored by their level of refinement (the computational domain is a square and has a level 0 refinement; see [70]).We note that this figure does not capture the entire domain, which continues both in the r and z directions.
An adaptive quadtree mesh is used in all the simulations in this work.An example of non-uniform grids is shown in figure 11.The maximum level of refinement is initially applied only near the interface, while one level coarser is used everywhere inside the droplet.Over time, the adaptive mesh refinement may also impose the maximum level in regions inside the droplet, based on error estimations on the velocity and volume fractions.In the outer (air) phase, the size of the grids gradually increases from the interface.We checked for the effect of grid size on the numerical results presented here to ensure the reported outcomes on dynamics of spreading and the droplet shapes are independent of the grid size.In Figure 12a, the maximum mesh refinement level is varied for a fixed ϵ vp = 10 We now perform a convergence test for the case of viscoelastic droplets (J = 0 and De ̸ = 0).In this situation, no viscoplastic regularization needs to be used, so we only check for mesh convergence.Figures 12b and 12c shows the radius and height of the droplet over time for different mesh levels with fixed De = 0.41.We note that little difference is observed with mesh levels above 9.For this reason, all the viscoelastic simulations performed in this work will use the mesh level 10.

B. Elastoviscoplasticity
When De ̸ = 0 and J ̸ = 0, the full Saramito constitutive equations (10) are used.The software Basilisk already contains a well-tested implementation of the Oldroyd-B model for solving viscoelastic flows given the parameters λ (relaxation time) and µ p (polymeric viscosity) [95].To make use of this existing module, we re-wrote the Saramito constitutive equation (10) and ϵ evp is a small threshold parameter.We note that equation (22) has the same form as the traditional Oldroyd-B equation but with non-constant relaxation time and polymeric viscosity coefficients.Therefore, we use the standard Oldroyd-B solver in Basilisk by dynamically setting these two coefficients according to the expressions in equation (22).
Figure 13 shows convergence tests for the EVP threshold parameter and for the mesh refinement parameter in a simulation with De = 0.816 and J = 0.18.We observe that ϵ evp has negligible influence on the shape of the droplet, even at relatively large values.Regarding mesh convergence, we see that for refinement levels above 9, the results are already very similar.For all the following EVP simulations in this article, we used the combination ϵ evp = 10 −7 and Level = 10.

C. Elastocapillarity
Figure 14 presents the data in Figure 8 as a function of De.Note that the pre-factors presented here are simply good fit through the data and are not from asymptotic analysis.A more rigorous analysis like those presented for viscoplastic limits [44] is needed to find the exact pre-factors.e l a s t o c a p i l l a r i t y l i m i t e l a s t o c a p i l l a r i t y l i m i t Fig. 14: Final radius and height as a function of De.The dashed lines show the scaling laws for the elastocapillarity limit, where both J and De are high enough such that the droplet spreads and stops like a soft solid.

[Case 4 :
De → ∞ -Newtonian Regime II]For finite Oh p and J , if De → ∞, we arrive at τ p → 0 and ▽ τ p → 0. Hence, we will have another Newtonian spreading with τ ≈ 2 Oh s D.

Fig. 2 :
Fig.2: Droplet final radius and height as a function of the plastocapillary number.Symbols are the present numerical simulations.The thick dashed lines show the theoretical predictions from[44].The gray horizontal line indicates the initial radius and height of the droplets.

Fig. 3 :Fig. 4 :
Fig. 3: Spreading radius (top) and height (bottom) over time for different values of the Deborah number.

Fig. 5 :
Fig. 5: Distribution of S at different timestamps inside a elastoviscoplastic droplet with De = 0.816 and J = 0.18.Blue regions indicate stress below the plastocapillary number J (unyielded) and red regions are above J (yielded).We note that the value of J used to calculate S is different in each row of this figure.

Fig. 6 :
Fig. 6: Distribution of S inside the elastoviscoplastic droplet for different values of De and J .All snapshots are taken at time t = 8.16 (see Video III).Blue regions indicate stress below the plastocapillary number J (unyielded) and red regions are above J (yielded).We note that the colorbar limits do not include the total range of values present in the data, since we are only interested in visualizing regions that are above or below zero .

Fig. 7 :
Fig. 7: The evolution of droplet height as a function of the plastocapillary and Deborah numbers.

Fig. 8 :
Fig. 8: Final radius (top) and height (bottom) of EVP droplets for a range of J and De.The horizontal dashed gray lines show the initial radius and height.The thick blue and red dashed lines correspond to the plastocapillarity limits (see section A).

Fig. 9 :
Fig.9: De − J parameter space and its different limits.The general case of EVP spreading reaches different regimes, depending on which rheological factor (elasticity or plasticity) dominates.Note that the high De regime of Newtonian fluids is for finite Ohp.

Fig. 10 :
Fig. 10: Convergence tests for the Bingham regularization parameter.The plastocapillary number is J = 0.18 and the mesh is fixed at Level = 10.

Fig. 12 :
Fig. 12: Mesh convergence tests for a Bingham material and a viscoelastic fluid.a) Snapshot of a Bingham droplet at time t = 10 with fixed ϵvp = 10 −3 and J − 0.18.b) and c) Radius and height evolution for a viscoelastic droplet with De = 0.41 and different mesh refinements.

Fig. 13 :
Fig.13: Convergence tests for the EVP regularization parameter (left) and for the mesh refinement parameter (right).In the regularization tests, we keep fixed the mesh level as 10 and in the mesh tests we fix the ϵevp = 10 −7 .All simulations are performed with De = 0.816 and J = 0.18. as