A sharp interface approach for wetting dynamics of hydrophobic coated droplets and soft particles

,


I. INTRODUCTION
The wetting of a solid surface by a liquid coincides with its ability to preserve the contact with the liquid [1][2][3][4].The wettability of a solid substrate by a pure droplet is quantified by the droplet's equilibrium contact angle θ eq , which, in turn, is determined by the balance between adhesive and cohesive forces of the three phases involved (solid, liquid, vapor).At the macroscopic scale, Young's equation [5] describes this balance as where σ sl , σ sg and σ are the solid-liquid, solid-gas and liquid-gas surface , respectively.Eq. ( 1) also estimates the degree of wettability, making the distinction between poor (θ eq > 90 • ) and good (θ eq < 90 • ) wetting regimes.Out of equilibrium, the additional complexities arising from time dependence and viscous dissipation make dynamic wetting critical to a wide range of phenomena including droplet spreading, capillary rise, imbibition, and more complex situations like fluid displacement in porous media or multiphase flow in oil recovery [6][7][8][9].The recent development of new catalytic devices, for example, requires the usage of liquid metals and metal alloys in the form of catalytic liquid droplets adsorbed on a porous solid support [10,11].However, several liquid metals such as gallium and galliumbased alloys oxidize when exposed to air and an inherent oxide layer appears on the top of the surface.This oxide layer acts as a solid-like "skin", encapsulating a liquid metal core [12,13] and changes the wetting properties of the droplet [14][15][16][17].Another example concerns the so called liquid marbles, realised by rolling a small liquid droplet in an poorly wetting powder.Because of the layer of powder grains at the liquid-air interface, the wetting of these droplets is inhibited [18,19], as required in some recent technological and microfluidic applications [20,21].The lattice e X P z L 6 2 c 6 q v s 5 i 9 N M 0 5 g s F 0 U Z h z q B 8 x j g g E l K N J 8 a g o l k 5 l Z I R l h Y y 6 9 F 6 s V 5 X r T l r P X M K f s B 6 + w T E J 5 P I < / l a t e x i t > ' ' ' 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " A r i 6 p S m z m 1 8 N i 8 q Z P m j I

Elastic energy:
Volume conservation: Wall-mesh interaction: < l a t e x i t s h a 1 _ b a s e 6 4 = " T + i g 8 i M Z x k q E R 9 k c W G e a P y k a T / g = " > A e r H e r Y 9 5 a 8 k q Z g 7 A H 1 i f P x Q N l H s = < / l a t e x i t > ' ' ' S i < l a t e x i t s h a 1 _ b a s e 6 4 = " 2 X J S d i o S 9 L i + V r p z y   The interface is resolved with a 3D triangular mesh.On each triangular face j, some force contributions ϕ ϕ ϕ are computed and distributed to the vertices i ∈ j, with the aim to consider (i) the interface elasticity/rigidity, (ii) the volume conservation, and (iii) the wall-particle interaction.We also report the corresponding involved parameters.Non-pre-stressed capsule 0.7 1.6 0.2 0.75 TABLE I. List of system models that can be explored by tuning the parameter α1, α2, and α3 in the interface model for a generic particle reported in Eq. ( 2).The left and right tables refer to the pre-stressed and non-pre-stressed particle classes, respectively.For each model we display the corresponding fitting parameters a, b, c and d appearing in Eq. (25).
Boltzmann (LB) method has been used since decades to address problems in wettability [22][23][24][25][26][27][28][29][30][31][32][33].A typical strategy used to simulate droplets within the framework of LB includes introducing non-ideal interface force models, such as the Shan-Chen [34] and the Free-Energy ones [35].However, these approaches model a diffuse interface.Simulating droplets with a complex rheology, specifically coated droplets including for example liquid metal ones with an oxide layer or liquid marbles, require the use of a constitutive law for the interface.In this case it is more convenient to use a method that reproduces the sharp-interface limit of hydrodynamics.In addition, pseudopotential or free-energy approaches do not easily allow to model a behaviour that, as it is typical for coated droplets, is intermediate between the case of a pure droplet and that of a capsule, as is the case for coated droplets.For these reasons, here we have opted for combining the LB model with an immersed boundary (IB) method, which naturally preserves the hydrodynamic sharp-interface limit, to simulate the complex droplet's wetting dynamics (see Fig. 1).
Our goal is to introduce a comprehensive numerical approach that allows modelling droplets with complex interfacial properties in a consistent way.Here, we model the interface as a 3D triangular mesh and employ a constitutive law based on the theory of Barthès-Biesel and Rallison [36] to explore the case of coated droplets.This approach allows us to describe in a continuous way the transition from pure droplet to capsule-like models by minimising the number of involved parameters.We provide a validation of our IBLB numerical simulations against the theoretical prediction in case of a simple shear flow experiment.Then, we perform wetting dynamics simulations, which show a good agreement with experimental observations in the case of a pure droplet, and we explore the range of accessible contact angles in terms of the involved parameters and the intensity of the interaction with the wall.With this approach we aim at providing a qualitative approximation of the mechanical behaviour of droplets with a complete and wide range of interfacial properties.Nevertheless, the model can be further refined to include additional interface properties, enhancing its accuracy and applicability, such as an extended model to mimic the oxidized layer thickness when dealing with liquid metal droplets.The paper is organised as follows: in Sec.II we describe the interface model introduced in the IBLB framework.Then, in Sec.III we summarise the main features of the IBLB model employed.The benchmark of a droplet in a simple shear flow is shown in Sec.IV.In the context of the wetting dynamics, Sec.V A report a model validation, while all wetting dynamics facets are analysed and discussed in Sec.V B. Results are summarised in Sec.VI.

II. INTERFACE MODEL
In this section, we describe the theoretical model employed in this work to simulate a generic soft particle.Following Barthès-Biesel & Rallison [36], we consider a two-dimensional, isotropic and homogeneous elastic interface with no bending resistance.Its mechanical response is characterised by an interface strain energy w S = w S (I 1 , I 2 , α 1 , α 2 , α 3 ) which is written in terms of the principal strain invariants I 1,2 and three parameters α 1,2,3 as [36] w where w S,0 is a reference value.I 1 and I 2 quantify the strain and dilation state of the membrane, respectively.The parameters α 1,2,3 , instead, characterize the material properties: The pre-stress α 1 is an isotropic tension without an applied load; α 2 is the resistance against area dilatation and α 3 the resistance against shear deformation (i.e., the strain modulus).For the sake of simplicity, hereafter we will refer to these three parameters as α = (α 1 , α 2 , α 3 ).
Concerning their choice, we distinguish the two main classes of pre-stressed (α 1 > 0) and non-pre-stressed (α 1 = 0) particles.For the latter class, an appropriate combination of α 2 and α 3 leads to the well-known Skalak [37] and Neo-Hookean [38] models.By conveniently tuning these parameters, one can switch from a pure droplet (α = (σ, 0, 0), where σ is the surface tension) to a pure elastic capsule (α = (0, 0, α 3 > 0)) and describe intermediate and more complex situations.In particular, in this work we consider also the classes of particles with α = (α 1 > 0, 0, α 3 > 0) and α = (α 1 > 0, α 2 > 0, α 3 > 0).We call the first type of particle a "softly coated droplet" because it has the characteristic surface tension term of a droplet, but also a strain modulus because of α 3 > 0. Because of the presence of a dilatational term α 2 > 0, we call the second type "rigidly coated droplet".Here we argue that the latter case could be used to describe the interfacial properties of particles like liquid metal droplets with oxidized surface.In Table I we summarise the type of particles investigated in this work.
In order to check that Eq. ( 2) leads to the correct particle dynamics, we consider the case of a shear flow experiment, where a particle with initial radius R and dynamic viscosity µ is placed between two distant moving walls which generate a shear rate γ (see top panel of Fig. 2).In this setup, time-dependent motion of the particle shape may be distinguished into two contributions: a solid body rotation and a stretching.Notice that the interface may rotate via a tank-treading motion despite of the particle reaching its steady state.This means that at this stage the interface deformation is constant in time at an Eulerian point x, but it is not constant by looking at a point X on the interface.Thus, following concepts and notation of Ref. [36], the time-evolution of the dimensionless position x (i.e., the position divided by the initial radius R), which is representative of the deformation field, reads where β 1 is the expansion coefficient around the initial spherical position (we truncate the equation at the leading order in β), while J and K are two symmetric and traceless second-rank tensors which depend only on time.It follows that the instantaneous external shape of the particle r can be computed in terms of the norm of Eq. ( 3), together with its normal n [36], Since in Eqs.(4) only tensor J appears, it means that J describes the overall deformation (i.e., the stretching contribution), while K describes the motion on the interface (i.e., the solid body rotation) [36].We remark that J is traceless because of the volume conservation constraint, whereas this property for K can be checked a posteriori.In the limit of small deformations, the evolution equations for J and K are given by [36] where is the Jaumann derivative [36] applied to a generic tensor A, which takes into account the rotation of the particle with the vorticity of the external fluid.E and Ω are the symmetric and asymmetric part of the velocity gradient, respectively, λ is the viscosity ratio between inside and outside fluids, and By numerically integrating Eq. ( 5), it is possible to obtain information on the transient deformation dynamics since the tensor J is directly related to the particle-deformation as [39] where Ca= µR γ/α is the capillary number.In the latter definition, one can consider α = α 1 = σ for a pure droplet or α = α 3 for a pure capsule.For the sake of simplicity, we fix the values of parameters α 1 , α 2 , and α 3 to be equal to the same value ᾱ and we will refer to this triad of values simply as α = ᾱ(0/1, 0/1, 0/1), with the vector elements turned on (1) and off (0) with the corresponding model (see Table I).

III. NUMERICAL IMPLEMENTATION
The dynamics of the inner and outer fluid is simulated using a single-component lattice Boltzmann (LB) method in terms of the fluid particle populations f i (x, t).The latter represents the probability distribution function of finding a fluid particle in a discrete lattice (Eulerian) node x at a discrete time t.The corresponding macroscopic behaviour is recovered in the long-wavelength limit, which allows the link with the Navier-Stokes equations.Indeed, the solutions of the Navier-Stokes equation for the total density and momentum are easily accessible from the populations as ρ(x, t) = i f i (x, t) and ρ(x, t)u(x, t) = i c i f i (x, t), respectively, with c i representing a set of 19 discrete velocities (i = 0, . . ., 18) living on a three-dimensional lattice (i.e., we employ a D3Q19 LB model).The dynamics of f i is ruled by a continuous succession of propagation and collision steps, as highlighted by the discretized Boltzmann equation [40,41] where ∆t is the time step.The propagation of f i on the lattice is described by the l.h.s. of Eq. ( 9) with the help of c i , while the single-relaxation-time BGK approximation of the collision operator appears as the first term in the r.h.s.The latter has the aim of modelling the relaxation of f i towards the equilibrium distribution f (eq) i (x, t), represented as the local Maxwellian distribution The relaxation process lasts for a relaxation time τ .In Eq. ( 10), f (eq) i is weighted by the lattice-dependent weights w i1 and depends on the speed of sound c s = ∆x/( √ 3∆t), where ∆x is the lattice spacing.The last term of Eq. ( 9) refers to the forcing implementation following the Guo scheme [42], where F is the force acting on the fluid.Notice that, to guarantee the second-order space-time accuracy, this forcing scheme modifies the fluid velocity as ρ(x, t)u(x, t) = i c i f i (x, t) + F∆t/2.In our simulations, we keep fixed to unity both ∆x and ∆t.Furthermore, the fluid dynamic viscosity µ in LB models is related to the relaxation time τ as µ = c2 s ρ(τ − 1/2).Here, we keep the viscosity ratio λ fixed to unity, since the investigation on the role played by λ goes beyond the purpose of this work.Then, to simulate the interface of a coated droplet or soft particle immersed in the surrounding LB fluid, we model the spherical particle interface using a 3D triangular mesh generated from a recursive refining of an icosahedron.Thus, the mesh resolution is defined in terms of the total number of triangular faces N f (see Fig. 6 for a pictorial view of particles with different resolutions).To couple the soft particle dynamics with that of the surrounding fluid, we use the immersed boundary (IB) method, i.e., a fluid-mesh interaction method developed for the first time by Peskin [43] and based on the distinction between interface (Lagrangian) nodes q(t) and fluid (Eulerian) nodes x.The resulting coupling is distinct in two operations, i.e., interpolation and spreading.The interpolation operation consists of the computation of the i−th interface-node velocity qi (t) from the fluid one (u(x, t)) as 2 [41] qi (t) = x u(x, t)δ D (x − q i (t))∆x 3 . ( This operation allows to update the node position q i (t) as: Then, the spreading operation is an interpolation of the interface nodal force to the fluid one which allows to make the latter aware of the presence of the interface: at this step, the total force (volume-)density the particle exerts on the fluid at the Eulerian node x is given by where ϕ ϕ ϕ i is the total force on the Lagrangian node i and the sum runs over all Lagrangian nodes.Both operations involve the so-called discrete delta function δ D , which is used to approximate the Dirac delta function on our lattice and is defined as [41,43,44] where φ 4 (r) is the "interpolation stencil" involving four Eulerian nodes along each coordinate axis [45] and defined as follows: The resulting IBLB method has been largely used to simulate the dynamics of capsules [45][46][47][48][49] and red blood cells [50][51][52][53][54].However, only a few works employed this method also for simulating droplet dynamics [55][56][57].A detailed step-by-step description of the IBLB algorithm implementation can be found in Ref. [41].
In our implementation, the total nodal force ϕ ϕ ϕ i , appearing in Eq. ( 13) and acting on the i-th node at position r i at time t, is given by the sum of several contributions, i.e., Each contribution plays a distinct role.First of all, ϕ ϕ ϕ S i incorporates the information on the elastic properties of the interface.Thus, we compute this nodal force term as where w S is the generalised strain energy defined in Eq. ( 2).Eq. ( 17) is calculated using a first order finite element method as described in Ref. [45].Then, because we are dealing with incompressible fluids, we need to consider a volume conservation constraint.With this aim, we follow Ref. [45] and we write the nodal volume force contribution ϕ ϕ ϕ V i as [58] ϕ ϕ ϕ where w V = k V (V − V 0 ) 2 /2V 0 is the volume energy.In this definition of the volume energy, k V refers to the volumeforce coefficient and it is kept fixed to 1. V = j V j is the instantaneous total particle volume, with the index j running over the number of faces N f3 , while V 0 is the initial total particle volume.Further details on how to compute nodal force contributions in Eqs. ( 17) and ( 18) can be found in Ref. [58].The last contribution in Eq. ( 16) corresponds to the wall-particle interaction, the key element for wetting dynamics simulations.The IBLB approach used in this work involves only one single fluid component and it is not possible to control the wall-fluid surface tensions σ sl and σ sg by introducing two different interactions as, for example, Huang and coworkers did in the case of multi-component pseudopotential LB models [24].However, many implementation of fluid-wall interactions in the case of single-component LB models [59][60][61] use a pseudopotential-like fluid-wall interaction which does not control σ sl and σ sg separately, but only their overall effect.These models work well in describing the wetting dynamics of droplets, despite of this limitation.In this work, we follow the same approach and we introduce a Lennard-Jones interaction on behalf of the wall-particle interaction: where d i is the shortest displacement vector between the centroid of the triangle to which node i belongs and the wall surface, and It means that the force is computed once for each triangle, and then it is distributed on its vertices.The choice of employing a Lennard-Jones interaction potential results from the necessity to model adhesive and repulsive forces between the droplet or particle and the surface.It follows the large amount of works based on molecular dynamics simulations which successfully studied the behaviour of nanodroplets or ridges on chemically patterned substrates [62][63][64][65][66]. Notice that in this work, we set ξ = 0.5∆x to have an interface-wall interaction that decays to zero after one lattice spacing, thus respecting as much as possible the microscopic range, and tune the interaction by changing only .The nodal force contributions along with the corresponding parameters are summarised in Fig. 1.By summarizing, the IBLB algorithm implemented in this work matches the following steps [41]: 1. Compute the nodal force ϕ i on each node i (Eq.( 16)); 2. Spread the nodal force to obtain the force acting on the fluid F(x, t) via Eq.( 13); 3. Perform the LB integration step: compute equilibrium distributions (Eq.( 10)), then apply the collision and perform the propagation.At this stage, F(x, t) enters in r.h.s. of Eq. ( 9); 4. Compute the fluid velocity u(x, t) from LB populations; 5. Interpolate the fluid velocity to compute the Lagrangian node velocity (Eq.( 11)); 6. Update the position of each node q(t) via Eq.( 12); 7. Iterate from step 1.
All simulations have been performed in a periodic domain in the x-and y-directions, while two walls are placed along the (vertical) z-direction.A half-node bounce-back rule implements second-order no-slip boundary conditions at the walls [41].Dimensional quantities are shown in lattice Boltzmann units (lbu).Note that, although this method is not able to capture particle breakup and coalescence, it provides an easy way to model different systems by simply tuning the α i parameters, as detailed in Section II.The introduction of the wall-particle interaction (19) induces an accumulation of interface nodes on the wall-particle contact area.Such an aggregation is more prominent for high values of that are required for observing small contact angles, and is responsible for a numerical instability in that regime.Re-meshing technique may help mitigating this problem, but since we are interested in large contact angles, this accumulation does not affect the results presented in this paper.

IV. BENCHMARK: SHEAR FLOW DYNAMICS
To showcase the versatility of the interface model proposed in Sec.II, we perform a double analysis by measuring the deformation of coated droplets and soft particles undergoing a shear flow.Indeed, on the one hand, we benchmark our model with what is known in the literature for the case of a pure droplet, while, on the other hand, we explore the different scenarios associated to each particle case listed in Table I.In this setup, we run simulations for particles with an initial radius R = 19 lbu, placed in a channel with a distance between the two walls H=128 lbu.The system has the same size along the other two directions, x and y.In order to vary the capillary number Ca, we systematically FIG. 2. Simulated experiment of a single particle under shear flow for the pre-stressed (panels (a)-(c)) and non-pre-stressed (panels (d)-(f)) particle models.In all panels, different symbols/colours refer to different models.Top panels: a sketch of the shear experiment and final shape of the particle for Ca= 0.3.Panels (a) and (d): time evolution of the deformation index D(t) (Eq.( 20)) as a function of time for capillary number Ca=0.05.Time is shown normalised with the shear rate γ, and dashed lines refer to the analytical solutions of Eqs.(5).Panels (b) and (e): the steady-state value of the deformation D as a function of Ca.Dashed lines draw the theroretical predictions: Eq. ( 22) (salmon line) for the pure droplet case and Eq. ( 21) (other colors lines) for all the other models.Panels (c) and (f): the steady-state value of the inclination angle Θ as a function of Ca.To validate the model in the case of a pure droplet, we report black crosses from Ref. [67], and we draw dotted lines for Eq.(23).
tune the values of ᾱ, keeping fixed the shear rate γ by the constraint of low Reynolds number (Re=10 −2 ).Without loss of generality, we set the fluid density ρ = 1 lbu and the relaxation time τ = 1 lbu, resulting in a dynamic viscosity of the particle µ=1/6 lbu.In Fig. 2(a) and (d), we report simulation data for the time-evolution of the deformation index defined as where r 1 and r 3 are the main particle semi-axes in the shear plane (see the top of Fig. 2 for a sketch).The simulation time is normalised with γ. Results show a very good agreement between simulations and the time-evolution of the deformation D defined in Eq. ( 8), obtained from the analytical solutions of Eqs.(5) (dashed lines).In addition, the steady-state value of the deformation D can be analytically estimated as a function of the triad of α as r/R (a) µ eq =88 ± µ eq =110 ± µ eq =128 ± µ eq =138 ± µ eq =157 ± where Ca=µR γ/α 1 .Since Eq. ( 21) has been computed with α 2 , α 3 = 0, it does not hold for a pure droplet, for which α = ᾱ(1, 0, 0).In the latter case we have [36] with λ = 1 in the present work.Fig. 2(b) and (e) confirm the agreement between simulation data and Eqs. ( 21) and ( 22) in the limit of small deformations (i.e., small Ca), while it diverges for larger values of D. Note that in Fig. 2(b) the theoretical prediction for cases α = ᾱ(1, 0, 1) and α = ᾱ(1, 1, 1) are so close to not being distinguishable.
In addition, to complete the picture of soft particle dynamics under shear flow, we report in Fig. 2(c) and (f) the inclination angle Θ (see top panel) as a function of the capillary number Ca.In the limit of small deformations and the case of a pure droplet with λ = 1, these results are again in agreement with what is expected from simulations [67] (black crosses) and the theory of Chaffey and Brenner [68] (dotted black line) which reads For non-pre-stressed particle models, we observe a stronger dependency on Ca, probably due to the higher rigidity.
To summarize, we find a good agreement for the time-evolution of the particle deformation D(t) and its steady-state value D between the analytical solution of the model equations ( 5) and our numerical model for both coated droplets and soft particles.This is true in the limit of small deformation, which is the basic assumption behind the theory [36].Furthermore, in the case of a pure droplet, both D and the inclination angle Θ follow the analytical predictions.It is worth noting that this benchmark contribute also to the validation of our generalised interface model against the interface response to an external flow.

V. WETTING DYNAMICS A. Model validation
We now analyse the wetting dynamics of coated droplets and soft particles simulated using the interface model discussed in Sec.II.In this kind of experiment, we consider a single particle with initial radius R and placed close to a flat wall, i.e., its initial position is such that the z-coordinate of its centre-of-mass Z CM is at a distance R from the wall to let it feel the action of an attractive wall-interface interaction with intensity (see Eq. ( 19) and Fig. 1 for a pictorial view).In our implementation, since we do not have direct control on σ sl , σ sg appearing in Eq. ( 1), we consider playing the role of an effective solid surface tension as ∝ σ sl − σ sg .I).Data refer to simulations with a number of triangular faces equal to N f = 16820.
Before entering into the details of the wetting dynamics of pre-stressed and non-pre-stressed particles, we validate our implementation by quantitatively investigating the spreading dynamics of a pure droplet, α = ᾱ(1, 0, 0), by comparing the time evolution of the radius r of the contact area with the literature.Indeed, it has been observed that this observable scales in time as where both the prefactor C and the exponent δ can vary.When capillary forces drive the droplet spreading and inertial effects are negligible, Eq. ( 24) coincides with the Tanner's law [69], predicting an exponent δ = 1/10.Contrariwise, when capillary and inertial forces are balanced, it has been observed that the value of the exponent can vary with some factors, such as viscosity [70], surface tension [70], droplet initial shape [71] and wettability [70,[72][73][74].In particular, a value of δ = 1/2 has been observed in the case of very small contact angles.In Fig. 3(a) we report the time evolution of r, normalised to the initial radius R at varying equilibrium contact angles θ eq .A scaling law following Eq.( 24) is observed, with the exponent δ slightly decreases at increasing θ eq , in agreement with Ref. [72] (see Fig. 3(b)).This implies that in our simulations of wetting dynamics, the inertia is not negligible and it plays a role in resisting the deformation.Note that, as later highlighted in Fig. 4(a), our model can capture only cases of large contact angles (88 • ≤ θ eq ≤ 180 • ), indeed δ approaches but does not reach a value of close to 1/2 which is characteristic for small contact angles [70,[72][73][74].This is comparable, for example, to the case of a liquid metal droplet, which has been observed to never assume values of θ eq < 100 • also in the case of oxidization, suggesting that our model can be used to study such kind of system.Concerning the prefactor C (Fig. 3(c)), it decreases as θ eq increases, once again in agreement with Ref. [72].Note that the "jumps" in r that are visible in Fig. 3(a) for θ eq = 157 • originate from the numerical error in measuring very small variations of the contact area. .All data refer to the case with all components of α = (α1, α2, α3) turned on, but different symbols/colours refer to a different "extreme" cases, where one of the three parameters is increased by an order of magnitude.
The number of triangular faces N f is the same as for the data in Fig. 4.

B. Results
With the aim of simulating the wetting dynamics of droplets with complex interface properties, we explore both pre-stressed (α 1 > 0) and non-pre-stressed (α 1 = 0) particles.For each system, characterized by α = (α 1 , α 2 , α 3 ) = ᾱ(0/1, 0/1, 0/1), we apply the same strategy used for the benchmark, which we summarize again for the sake of clarity.First, we fix the value of ᾱ to 10 −4 lbu.This choice fixes both the surface tension for pre-stressed particles and the strain modulus for non-pre-stressed particles.Then we measure θ eq as a function of the wall-particle interaction energy .In Appendix A we discuss also a resolution test for the wetting dynamics.In Fig. 4(a) and (b) we report the measured θ eq as a function of the ratio /ᾱ for pre-stressed and non-pre-stressed particle models, respectively.In Fig. 4(c) and (d) we report for convenience the corresponding values of cos θ eq .After the largest value of /ᾱ reported in each plot, numerical instabilities appear in the contact area region; these set the limit of applicability of our approach in terms of contact angles that can be modeled.Concerning pre-stressed particles, pure droplet (circles) appear to be marginally more stable with respect to the choice of than the other particles (softly coated particles, α = ᾱ(1, 0, 1), upward triangles; rigidly coated particles, α = ᾱ(1, 1, 1), pentagons), but can reach only a slightly higher contact angle (θ eq = 88 vs. θ eq = 79 • ).Notice that the cases α = ᾱ(1, 0, 1) and α = ᾱ(1, 1, 1) are very similar, meaning that when the system is very rigid, the dilatational contribution given by α 2 is not relevant for the equilibrium contact angle θ eq .This result is in contrast with what we observed in the shear flow.Non-pre-stressed particle models (Fig. 4(b) and (d)) are stable for a more limited range of values of /ᾱ.After around the value of /ᾱ = 1.5 the contact angle does not drop significantly anymore.Similarly to the pre-stressed case, α 2 does not seem to have any influence on θ eq .The behaviour of cos θ eq as a function of /ᾱ follows very well the empirical behavior cos θ eq a tan −1 [b( /ᾱ − c)] − d where a, b, c and d are fitting parameters depending on the type of particle (see Table I and dashed lines in Fig. 4(c) and (d)).Eq. ( 25) differs from Eq. (1) because, as mentioned above, the model we present in this work misses the direct control of the wall surface tensions but rather drives the mechanical interaction between the particle and the solid surface.Obviously, the fitting constants represent (unknown) functions of the parameters α 1 , α 2 , α 3 and .By increasing separately by a factor of ten each of the components of α, as reported in Fig. 5(a), we can understand that the leading order behavior is dictated by α 1 , while α 2 and α 3 provide relatively minor changes in θ eq .In addition, from Fig. 5(b) one can see that α 1 must enter in Eq. ( 25) in the ratio with , because data for the increased α 1 (downward triangles) collapses onto the original case (pentagons) once plotted as a function /α 1 .The remaining parameters α 2 and α 3 , instead, do not appear to provide a similar scaling, thus meaning that these two parameters may functionally enter in the other fitting parameters of Eq. (25).

VI. CONCLUSIONS
In this work we introduced a novel numerical framework to accurately characterize coated droplets and soft particles.This approach is based on the thory of Barthès-Biesel and Rallison [36] and enables us to capture the unique behavior t e x i t s h a 1 _ b a s e 6 4 = " a 5 r D s k A k D Y + q 3 p A E t g G O P j R B V F w = " > A A A B 8 H i c b V D L S g N B E J y N r x h f U Y 9 e B o P g K e w G X y c J e v E Y w T w k W Z b Z y W w y Z G Z 2 m e k V w p K v 8 O J B E a 9 + j j f / x k m y B 0 0 s aC i q u u n u C h P B D b j u t 1 N Y W V 1 b 3 y h u l r a 2 d 3 b 3 y v s H L R O n m r I m j U W s O y E x T H D F m s B B s E 6 i G Z G h Y O 1 w d D v 1 2 0 9 M G x 6 r B x g n z J d k o H j E K Q E r P U J Q w 9 c Y A i 8 o V 9 y q O w N e J l 5 O K i h H I y h / 9 f o x T S V T Q A U x p u u 5 C f g Z 0 c C p Y J N S L z U s I X R E B q x r q S K S G T + b H T z B J 1 b p 4 y j W t h T g m f p 7 I i P S m L E M b a c k M D S L 3 l T 8 z + u m E F 3 5 G V d J C k z R + a I o F R h i P P 0 e 9 7 l m F M T Y E k I 1 t 7 d i O i S a U L A Z l W w I 3 u L L y 6 R V q 3 o X 1 f P 7 s 0 r 9 J o + j i I 7 Q M T p F H r p E d X S H G q i J K J L o G b 2 i N 0 c 7 L 8 6 7 8 z F v L T j 5 z C H 6 A + f z B z p 8 j 2 U = < / l a t e x i t > t 2 > t 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " i W O B 2 l F f y t z a J X W X a Q K d f y B 1 c e c = " > A A A B 8 H i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 5 M E v X i M Y B 6 S L M v s Z J I M m Z 1 d Z n q F s O Q r v H h Q x K u f 4 8 2 / c Z L s Q a M F D U V V N 9 1 d Y S K F Q d f9 c g p L y y u r a 8 X 1 0 s b m 1 v Z O e X e v a e J U M 9 5 g s Y x 1 O 6 S G S 6 F 4 A w V K 3 k 4 0 p 1 E o e S s c 3 U z 9 1 i P X R s T q H s c J 9 y M 6 U K I v G E U r P W D g k S u C g R u U K 2 7 V n Y H 8 J V 5 O K p C j H p Q / u 7 2 Y p R F X y C Q 1 p u O 5 C f o Z 1 S i Y 5 J N S N z U 8 o W x E B 7 x j q a I R N 3 4 2 O 3 h C j q z S I / 1 Y 2 1 J I Z u r P i Y x G x o y j 0 H Z G F I d m 0 Z u K / 3 m d F P u X f i Z U k i J X b L 6 o n 0 q C M Z l + T 3 p C c 4 Z y b A l l W t h b C R t S T R n a j E o 2 B G / x 5 b + k e V L 1 z q t n d 6 e V 2 n U e R x E O 4 B C O w Y M L q M E t 1 K E B D C J 4 g h d 4 d b T z 7 L w 5 7 / P W g p P P 7 M M v O B / f N 2 6 P Y w = = < / l a t e x i t > t 1 > t 0 Initial time step Final time step < l a t e x i t s h a 1 _ b a s e 6 4 = " z l D c t 7 B 4 a 2 l y 7 y W F Q D O T y k 0 C n 6 c = " > A A A B 6 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K r 2 P Q i 8 e I 5 g H J E m Y n s 8 m Q 2 d l l p l c I S z 7 B i w d F v P p F 3 v w b J 8 k e N L G g o a j q p r s r S K Q w 6 L r f T m F l d W 1 9 o 7 h Z 2 t r e 2 d 0 r 7 x 8 0 T Z x q x h s s l r F u B 9 R w K R R v o E D J 2 4 n m N A o k b w W j 2 6 n f e u L a i F g 9 4 j j h f k Q H S o S C U b T S A / b c X r n i V t 0 Z y D L x c l K B H P V e + a v b j 1 k a c Y V M U m M 6 n p u g n 1 G N g k k + K X V T w x P K R n T A O 5 Y q G n H j Z 7 N T J + T E K n 0 S x t q W Q j J T f 0 9 k N D J m H A W 2 M 6 I 4 N I v e V P z P 6 6 Q Y X v u Z U E m K X L H 5 o j C V B G M y / Z v 0 h e Y M 5 d g S y r S w t x I 2 p J o y t O m U b A j e 4 s v L p H l W 9 S 6 r F / f n l d p N H k c R j u A Y T s G D K 6 j B H d S h A Q w G 8 A y v 8 O Z I 5 8 V 5 d z 7 m r Q U n n z m E P 3 A + f w A H V I 2 l < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " 7 d m i A k b F 4 q 4 9 U / c + T e m D 0 F F D 1 w 4 n H M v 9 9 4 T J p w p j d C H l d v Y 3 N r e y e 8 W 9 v Y P D o / s 4 5 O O i l N J a J v E P J a 9 E C v K W U T b m m l O e 4 m k W I S c d s P J 9 c L v T q l U L I 7 u 9 S y h v s C j i A 0 Z w d r e n M m z o v z 7 n z M W w t O P n O A / s D 5 / A F B O J P 0 < / l a t e x i t > , ✏ < l a t e x i t s h a _ b a s e = " G i s J L o z Y h k D Y F S l A y E s v w = " > A A A B n i c b V D L S g N B E O y N r x h f U Y e B o P g K e y K r P s B g A M / w C m + O c F c d + d j l p w p l D + A P n w c z N o C < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " r O O q w Z c h p w N c w L G 4 c C a o A k w e A L P 4 B W 8 O b n z 4 r w 7 H / P W F a e Y O Q R / 4 H z + A B i f l H 4 = < / l a t e x i t > ' ' ' V i < l a t e x i t s h a 1 _ b a s e 6 4 = " x q L c v k q 1 i b m o / R s h i D 0 b H f 4 b z y k = " > A A A B + 3 i c b V D L S g M x F M 3 4 r P U 1 1 q W b Y B F c l R n x t S y 6 c V n B P q A z D p k 0 b U O T T E g y x T L M r 7 h x o Y h b f 8 S d f 2 P a z k J b D 1 w 4 n H M v 9 9 4 T S 0 a 1 8 b

R
a o A k w e A L P 4 B W 8 O b n z 4 r w 7 H / P W F a e Y O Q R / 4 H z + A B o l l H 8 = < / l a t e x i t > ' ' ' W i Total nodal force < l a t e x i t s h a 1 _ b a s e 6 4 = " 5 V 1 p k w U 6 4 7 r / 2 / a J 5 k H c 3 P p P b Q 4 = " > A A A B + X i c b V D L S s N A F L 2 p r 1 p f U Z d u B o v g q i T i a 1 n Q h c s K 9 g F N C J P p t B 0 6 M w k z k 0 I J / R M 3 L h R x 6 5 + 4 8 2 + c t l l o 6 4 E L h 3 P u 5 d 5 7 4 p Q z b T z v 2 y m t r W 9 s b p W 3 K z u 7 e / s H 7 u F R S y e Z I r R J E p 6 o T o w 1 5 U z S p m G G 0 0 6 q K

FIG. 1 .
FIG.1.Sketch of the wetting dynamics of a generic particle with initial radius R initially placed in contact with a flat wall.The interface is resolved with a 3D triangular mesh.On each triangular face j, some force contributions ϕ ϕ ϕ are computed and distributed to the vertices i ∈ j, with the aim to consider (i) the interface elasticity/rigidity, (ii) the volume conservation, and (iii) the wall-particle interaction.We also report the corresponding involved parameters.

< l a t e x i t s h a 1 _
b a s e 6 4 = " 1 T n a w e f 8 B 2 I c X + 5 J D n L m G T d 6 E x 4 = " > A A A B / 3 i c d V D L S s N A F J 3 U V 6 2 v q O D G z W A R 3 B i S v q y 7 o h u X F e w D m l I m 0 2 k 7 d D I J M x M 1 x C 7 8 F T c u F H H r b 7 j z b 5 y m F V T 0 w I X D O f d y 7 z 1 e y K h U t v 1 h Z B Y W l 5 Z X s q u 5 t f W N z S 1 z e 6 c p g 0 h g 0 s A B C 0 T b Q 5 I w y k l D U c V I O x Q E + R 4 j L W 9 8 P v V b 1 0 R I G v A r F Y e k 6 6 M h p w O K k d J S z 9 w 7 j n q J 6 3 v B b e I q y m N 4 g x i b T H p m 3 r b s c r l Q s K F t F U v F Y t V J S a l Y O Y W O Z a f I g z n q P f P d 7 Q c 4 8 g l X m C E p O 4 4 d q m 6 C h K K Y k U n O j S Q J E R 6 j I e l o y p F P Z D d J 7 5 / A Q 6 3 0 4

3 <
B J 6 N e + P R e D F e Z 6 0 Z Y z 6 z C 3 7 A e P s E / o W W x Q = = < / l a t e x i t > u wall < l a t e x i t s h a 1 _ b a s e 6 4 = " D h / I + B V x F t 3 8 c 6 O 6 v 5 Y h 0 3 L B M 0 4= " > A A A B / n i c d V D L S s N A F J 3 U V 6 2 v q r h y M 1 g E V y F p 2 l p 3 R T c u K 9 g H N C F M p p N 2 6 G Q S Z i Z q C Q V / x Y 0 L R d z 6 H e 7 8 G 6 c P Q U U P X D i c c y / 3 3 h M k j E p l W R 9 G b m l 5 Z X U t v 1 7 Y 2 N z a 3 i n u 7 r V l n A p M W j h m s e g G S B J G O W k p q h j p J o K g K G C k E 4 w u p n 7 n h g h J Y 3 6 t x g n x I j T g N K Q Y K S 3 5 x Y P U z 9 w o i O 8 y V 1 E + h r e I s c n E L 5 Y s 0 6 p W y 2 U L W q Z T c Z y 6 P S M V p 3 Y G b d O a o Q Q W a P r F d 7 c f 4 z Q i X G G G p O z Z V q K 8 D A l F M SO T g p t K k i A 8 Q g P S 0 5 S j i E g v m 5 0 / g c d a 6 c M w F r q 4 g j P 1 + 0 S G I i n H U a A 7 I 6 S G 8 r c 3 F f / y e q k K 6 1 5 G e Z I q w v F 8 U Z g y q G I 4z Q L 2 q S B Y s b E m C A u q b 4 V 4 i A T C S i d W 0 C F 8 f Q r / J+ 2 y a d f M 6 l W l 1 D h f x J E H h + A I n A A b n I I G u A R N 0 A I Y Z O A B P I F n 4 9 5 4 N F 6 M 1 3 l r z l j M 7 I M f M N 4 + A Z F 1 l o 4 = < / l a t e x i t > l a t e x i t s h a 1 _ b a s e 6 4 = " J g m E o s O v O E V p x U 9 D A f m e j B T U t P 0 3 e l M p N S 6 z O P L g C B y D U + C B G m i A a 9 A E b U D A H X g A T + D Z k c 6 j 8 + K 8 L l t z T j Z z C H 7 A e f s E t H q P P Q = = < / l a t e x i t >

< 3 < 3 <FIG. 3 .
FIG.3.Spreading experiment for a pure liquid droplet, α = (α1, 0, 0), on a flat surface.Panel (a): radius of the contact area r as a function of time t, where r and t are reported normalised to the initial radius R and the characteristic time t * = (ρR 3 /σ) 1/2 , respectively.Different symbols/colours refer to different values of the equilibrium contact angle θeq.In all cases, we observe a scaling law r/R = C(t/t * ) δ .The solid line indicates the scaling with δ = 3/10, while the dotted line refers to scaling δ = 3/20.Panels (b) and (c) show the value of the dimensionless exponent δ and the dimensionless prefactor C, respectively, as a function of θeq.

FIG. 4 .
FIG. 4. Experiment of wetting dynamics.Panels (a) and (b): Equilibrium contact angle θeq as a function of the wall-particle interaction intensity , normalised to ᾱ. Panels (c) and (d): Corresponding values of cos θeq.Left panels ((a) and (c)) refer to pre-stressed particle models, while right panels ((b) and (d)) refer to non-pre-stressed particle models.In all panels, different symbols/colours refer to different models, while dashed lines indicate fitting curves with Eq. (25) (values of fitting parameters are listed in TableI).Data refer to simulations with a number of triangular faces equal to N f = 16820.

FIG. 5 .
FIG.5.Equilibrium contact angle θeq as a function of the wall-particle interaction intensity , normalised to ᾱ (panel (a)) and max(αi) with i = 1, 2, 3 (panel (b)).All data refer to the case with all components of α = (α1, α2, α3) turned on, but different symbols/colours refer to a different "extreme" cases, where one of the three parameters is increased by an order of magnitude.The number of triangular faces N f is the same as for the data in Fig.4.