Surface tension and wetting at free surfaces in smoothed particle hydrodynamics

Abstract Surface tension and wetting are dominating physical effects in microscale and nanoscale flows. We present an efficient and reliable model of surface tension and equilibrium contact angles in smoothed particle hydrodynamics for free-surface problems. We demonstrate its robustness and accuracy by simulating several three-dimensional free-surface flow problems driven by interfacial tension.


Introduction
At the micro and nano scale, fluid flows are dominated by interfacial tension that arises at the interface between different phases.When a liquid drop is in contact with a plane solid substrate, the wetting force acting at the triple line leads to an equilibrium contact angle as illustrated in Fig. 1.Surface tension and wetting are essential for many phenomena, solid liquid gas Θ ∞ solid-liquid interface three-phase contact line liquid-gas interface contact point solid-gas interface including technological processes like oil recovery [1,2], two-phase heat transfer [3], and inkjet printing [4].The wetting properties of a solid surface depend on its chemical composition and morphology, determining its repellent or attractive behavior [5,6].
In Smoothed Particle Hydrodynamics (SPH) [7,8], surface tension can be modeled either by pairwise forces between the particles [9,10,11,12] or by a phenomenological Continuum Surface Force (CSF) [13,14,15].The concept of pairwise forces between SPH particles is motivated by the molecular origin of surface tension [16] and readily applies to the simulation of free surfaces.Here, the attractive forces between particles of different phases must be calibrated to obtain a desired equilibrium contact angle.This approach is widely used for solving problems in droplet dynamics, drop interaction with textured surfaces [17], and contact angle hysteresis [18].When used in simulations, however, such models require a large support domain for particles, affecting the computational performance significantly [12].CSF can reliably predict the dynamics due to surface tension gradients at phase boundaries, known as the Marangoni effect.However, its application to free surfaces is problematic since the truncation of the kernel near the free surface leads to unacceptable errors in the local curvature.Incompressible Gas-liquid two phase flow problems can be reliably modelled as free surface flows whenever the shear stress due to the gas phase is negligible.Thus the limitation in modeling surface tension at the free surface has severely restricted the scope of 3D free surface SPH simulations to high Weber number flows.This kernel truncation error at the free surface in CSF has been addressed using several strategies.Mirror particles were employed by Ordoubadi et al. [19] to produce robust twodimensional flow simulations in weakly compressible SPH.Hirschler et al. [20] used CSF to simulate two-dimensional droplet collisions at the free surface using kernel correction parameters.Unfortunately, only two-dimensional systems can be described using any of these approaches.For three dimensional cases, [21] introduced an analytical coefficient to account for truncation of kernels near free surface in SPH operators.The curvature estimate at the free surface was recently improved by Fürstenau et al. [22], while Blank et al. [23] developed a method to describe surface tension using a Young-Laplace pressure boundary condition in SPH.When using CSF, the wetting force in the triple line region can be modeled using the smoothed normal correction scheme introduced by Breinlinger et al. [24].Here, the normal vectors (pointing from the liquid into the gas phase), computed at the positions of SPH particles that are located in the vicinity of the three-phase contact line, are modified.Doing so, the curvature is computed at the positions of those SPH particles, and thus, wetting is incorporated into CSF.This approach can, however, not be applied to free surface flow problems due to the inaccurate force computation caused by insufficient support of particles near the triple contact line region.
The current paper proposes a CSF model that relies on an improved smoothed normal correction scheme.We will show that this model reliably describes three-dimensional free surface problems.We simulate several problems using Incompressible Smoothed Particle Hydrodynamics (ISPH) to demonstrate the model's accuracy and robustness.These problems comprise notoriously difficult systems like droplet oscillation, droplet spreading on a solid substrate, and the pinning of a drop at the contact line of two inclined planes.

Smoothed Particle Hydrodynamics
Smoothed Particle Hydrodynamics (SPH) is a numerical method for solving continuum problems [7,8,25].Unlike in mesh-based methods, the domain is discretized by SPH particles with a certain mass at which properties such as velocity, u, pressure, p, or density ρ, are defined.
The integral interpolant, Φ I (r), of a physical property, Φ(r), at position r, can be obtained from a spatial convolution with a kernel function, W (r − r ′ , h), having a compact support: Here, h is the smoothing length, and Ω c denotes the whole continuous domain.We approximate the integral in Eq. ( 1) by a sum over the set Ω of all SPH particles, where r b , m b , and ρ b are the position, mass, and density of particle b.In particular, the summation interpolant at the position r a of an SPH particle reads In the same approximation, the spatial derivative of a physical quantity, ∇Φ (r), at position r a , can be expressed in terms of the kernel gradient ∇W (r a − r b , h) by We introduce shorthand notations and rewrite Eqs. ( 3) and (4): Here we use the Wendland C 2 kernel function normalized for three spatial dimensions [26], We describe the dynamics of a Newtonian incompressible fluid using the Navier-Stokes equation and the continuity equation where u, ρ, p, ν, are the velocity, density, pressure, and kinematic viscosity, f b is an acceleration due to a body force such as gravity, and t is time.The acceleration caused by surface tension, f s , is incorporated into Eq.( 7) as a body force according to the Continuum Surface Force (CSF) model [13].
The total acceleration of a particle a given by Eq. ( 7) is, thus, a sum of accelerations due to pressure gradient, f p a , viscosity, f v a , surface tension, f s a , and body forces per unit mass, f b a .The first contribution, the acceleration due to pressure gradient, can be approximated by [25] where 10) ensures conservation of linear momentum, and is appropriate for SPH particles, a, that are fully embedded by other SPH particles, b.This is the case with particles that are in the bulk of the material, far from surfaces.This approximation cannot be applied to particles, a, near a free surface.This case will be discussed in Sec. 4.
The second contribution to the acceleration in Eq. ( 9), the viscous acceleration of particle a, is given by [27] where u a and u b are the velocities of particle a and b, η a and η b are the dynamic viscosity coefficients assigned to particles a and b, and F ab is given by Here the term (0.01h) 2 in the denominator avoids divergence of F ab when particles a and b come very close to each other.The third contribution in Eq. ( 9), the acceleration of a particle a due to surface tension, f s , will be described in Sec. 5.

Incompressible Smoothed Particle Hydrodynamics (ISPH)
For incompressible fluids, the velocity field is divergence-free, and the continuity equation, Eq. ( 8), turns into The idea of Incompressible Smoothed Particle Hydrodynamics (ISPH) is to enforce Eq. ( 13) by imposing a pressure field in such a way that the velocity field becomes divergence-free [28].In ISPH, we use the predictor-corrector integration scheme by Cummins and Rudman [28] to advance particles in space.At time step n, the prediction step computes the positions of the SPH particles from their current positions and velocities, The corresponding velocities are The pressure, p * a , corresponding to the predicted velocity can be obtained by solving the pressure Poisson equation (PPE) From the predicted velocities, u * a , and the corresponding pressure, p * a , the correction step yields the divergence-free velocities at time step n + 1, The corresponding particle positions at time step n + 1 are The left-hand side of the PPE, Eq. ( 16), can be discretized to obtain [28] ∇ and its right-hand side results in [25] ∇ Substituting Eqs.(19) and (20) into Eq.( 16) Equation ( 21) can be uniquely solved using any iterative linear solver such as the BiCGSTAB method [29], provided the linear system is not singular, e.g., if a Dirichlet boundary condition is applied.The discretized PPE, Eq. ( 21), is an appropriate approximation of Eq. ( 16) for well-embedded SPH particles, that is, for particles in the bulk of the material.The corresponding approximation for SPH particles near a free surface is described in the subsequent section.

Kernel correction for particles near a free surface
In the vicinity of boundaries, the summations in Eqs. ( 3) and (4) underestimate Φ I (r) and ∇Φ I (r), due to lacking SPH particles in the neighborhood of r.We call such SPH particles insufficiently supported by neighboring SPH particles.To account for such particles, semianalytical [30] or numerical normalization techniques [31,32] are used near free boundaries.
In the current paper, akin to [30] we do not model the gas phase by SPH particles, but we describe the influence of the ambient gas on the liquid through boundary conditions for the liquid phase at the free surface.
To identify SPH particles with insufficient support we use the Shepard filter [33], However, unlike the approaches described in [34] where a thin layer of particles comprises the pressure boundary condition (which deteriorates the accuracy and the stability, see [30,35]), in the present approach, the pressure at these particles is also solved for.
Particles with insufficient support are close to interfaces, therefore, this method defines the fluid-gas boundary region.Whether a particle a belongs to the subset of particles representing the bulk, Ω b ⊂ Ω, or to the subset representing the boundary, For well supported particles, a ∈ Ω b , we evaluate Eq. ( 10) using the pressure gradient approximation [25]; for a ∈ Ω fs we use the approximation given in [23]: Here, p o a is the external pressure representing the Dirichlet boundary condition of a particle a.Since pressure appears in the Navier-Stokes equation only as a gradient, we set p o a = 0.In this work, we assume that the motion of the gas phase follows the liquid phase, that is, the relative velocity between the gas and liquid phase ceases.As a consequence, there is no contribution of the gas phase to the viscous acceleration of SPH particles Eq. (11).Analogously, using the same assumption the contribution of the gas phase to the divergence of velocity approximation on the right-hand side of the PPE in Eq. ( 21) is zero.Therefore, Eqs. ( 11) and ( 21) can be used for both, SPH particles located in the vicinity of a free surface, and SPH particles located in the bulk of the simulated material.To approximate the left-hand side of Eq. ( 21) for SPH Particles that are located in the vicinity of a free surface, we use [30,23] Note that the term marked by β is constant if the mass and volume assigned to the SPH particles is uniform throughout the domain.
In conclusion, we approximate the PPE by Eq.( 21) if S a > 0.95 Eq.( 25) if S a ≤ 0.95 .
This ambient pressure application can be used to model surface tension by computing the Young-Laplace pressure boundary condition, as shown in [23].

Surface tension at free boundaries
The Continuum Surface Force (CSF) model [13] evaluates the acceleration of a particle a due to surface tension by where F s a is a force per unit volume [14] Here, the index a to F s a and its constituents must be interpreted such that the quantity refers to a local property of the surface, but is assigned to SPH particle a.In this sense, σ a is the local coefficient of surface tension of the phase boundary, assigned to particle a.In the SPH literature, we speak briefly of "surface tension of particle a".Similarly, κ a , is the local value of the mean curvature of the surface, assigned to particle a.It is termed "curvature of particle a".Analogously, the local values of the unit normal vector to the interface pointing from phase I to phase II are assigned to the SPH particles.The contribution na assigned to particle a is called "unit normal vector of particle a" and can be computed by The same applies to the surface gradient of the surface tension, ∇ s σ a .For the translation between surface tension into a volumetric force at the position of an SPH particle a we employ the surface delta function, δ s a , that peaks at the interface.The second term in Eq. ( 28) drives fluid flow tangential to the interface due to the surface tension gradient.In the current paper, we assume constant surface tension.Therefore, the surface tension gradient and, thus, the second term in Eq. ( 28) vanish.The first term in Eq. ( 28) describes a force (per unit volume) directed perpendicular to the interface.This force counteracts the curvature of the interface and, thus, minimizes the surface area.We represent the surface delta function of an SPH particle in Eq. ( 28) by [22] where Ω l ⊂ Ω is the subset of SPH particles representing the liquid phase.
Both the absolute value and the direction of n a obtained from Eq. ( 29) depend sensitively on the positions of neighboring particles b, therefore, naïve computation of the normal unit vector vector, na ≡ n a / ∥n a ∥, is subject to large fluctuations leading eventually to inaccurate F s a when used in Eq. ( 28).Therefore, instead of using the naïve value, we use a smoothed normal vector [27,22]: The smoothing increases the region where particles contribute to surface tension.Following our previous notation, n b is a normal vector at the position of particle b.In the bulk of the fluid, the smoothed normal vectors have small magnitudes and their orientation is of low significance.Therefore, we discard the normal unit vectors of such SPH particles a ∈ Ω l , whose normal vector's magnitude is smaller than a threshold, ε n ∈ 0.1 h , 0.2 h .In the simulations presented in this paper, we have used ε n = 0.01 throughout.Subsequently, the smoothed and filtered normal vectors from Eq. ( 31) are normalized by The set of unit vectors provided by Eq. ( 32), is used to compute the mean curvature which is needed to evaluate the first term in Eq. ( 28) [36,14]: To increase the accuracy of the curvature approximation in Eq. ( 33), we substitute the kernel gradient with ∇W ab ≡ C a ∇W ab [31,32], where is the correction matrix computed from all neighboring SPH particles b.This yields the following curvature approximation with O(h 2 ) convergence [31,32] Using Eqs. ( 28), ( 27),( 32), (35), and assuming a constant surface tension coefficient in the simulated material, yields the acceleration of a particle a in the normal direction Similar to [23], where the Shepard filter is used to increase the robustness and accuracy of the obtained Young-Laplace pressure jump, we modify Eq. (36) to Here, S n a is the Shepard filter computed by where Ω n ⊂ Ω is the subset of SPH particles satisfying ∥ñ a ∥ ≥ ε n .The use of the factor (1 + 1/S n a ) which replaces the theoretical factor of 2 is empirical and is obtained from numerical experimentation.

Wetting forces at free boundaries
The spreading force per unit length at the three-phase contact line can be described as a function of the contact angle, Θ, and the equilibrium contact angle, Θ ∞ , by [37] According to Breinlinger et al. [24], Eq. ( 39) can be imposed on SPH particles located in the vicinity of the three-phase contact line by modifying their normal vector orientations.As a result, the curvature is modified, and Eq. ( 37) includes the acceleration due to wetting phenomena.The scheme of [24] cannot be directly applied to free surface problems since the gas phase is not represented by SPH particles.Therefore, in the following subsections, we describe two alternative approaches for modeling wetting and non-wetting contact for free surface problems.

Wetting contact angle, Θ ∞ ≤ 90 •
To model wetting contact angles (hydrophilic substrates), we compute the normal vectors of those SPH particles located near the three-phase contact line as a function of the required equilibrium contact angle.According to [13] the normal vector associated with a desired equilibrium contact angle, Θ ∞ , can be computed for a particle a ∈ Ω l by Here, tsf a is the smoothed normalized tangent vector between the solid and fluid phases, and nsf a is the smoothed normalized normal vector pointing from the solid to the fluid phase computed at the position of particle a.The normal vector, nsf a , at the position of particle a, is given by where Ω s ⊂ Ω is the subset of all SPH particles representing the solid phase, and S s a is the Shepard filter computed by Analogous to the computation of na , we discard the unit normal vectors of SPH particles if ∥ñ∥ a < ε n .
The unit tangent vector in Eq. ( 40), tsf a , is computed by The normal vectors computed in Eq. ( 40) could be used to replace the normal vectors given by Eq. ( 32) to model wetting.However, as shown in [24], it is preferable to avoid an instantaneous change of the normal vectors of those SPH particles located near the three-phase contact line region.Instead, a smooth transition from na to n∞ a (Eq.( 40)) depending on a particle's distance to the solid phase (wall) should be applied Here, f a is given by where r max is the kernel radius (r max = 2h when using the Wendland kernel in Eq. ( 6)), and Here, ∆x is the spacing between two SPH particles when arranged on a square lattice, and ∆x 3 is the volume assigned to an SPH particle.The lower limit of f a ensures that particles that are located closer to the wall than ∆x are prescribed with n∞ a .The upper limit of f a restricts the normal correction to particles in a range of r max to the solid phase.The curvature at the position of a particle a ∈ Ω l can be computed by In the course of this work, it was found that computing the curvature in Eq. ( 47) is not accurate enough to obtain stable equilibrium droplet shapes for Θ ∞ ≤ 90 • ).In particular, underestimated curvatures of SPH particles located near the three-phase contact line lead to unlimited spreading of drops on a plane solid substrate.For this reason, we propose to compute the curvature at the position of an SPH particle a ∈ Ω l by where ns b is the normal vector of a neighboring particle b ∈ Ω s (particles representing the solid phase) given by Here, Θ s ∞ is a calibration parameter used to obtain corrected normal vectors of neighboring SPH particles b ∈ Ω s .By setting Θ s ∞ > Θ ∞ , the curvatures at the position of SPH particles a ∈ Ω l are shifted to larger values, which prevents the continuous spreading of the drop on the solid substrate, which otherwise happens if Θ ∞ is used.This parameter may require calibration for different discretization parameters such as the kernel and the smoothing length.The value of this parameter is listed where relevant, for example, in Table 3.The condition in Eq. ( 49) measures the distance of a neighboring particle b ∈ Ω s to particle a in the normal direction.The sign of the distance allows to distinguish between neighboring particles b ∈ Ω s located on the liquid or on the gas side of particle a, as illustrated in Fig. 2. Neighboring SPH particles that are located on the gas side of particle a, r ab • nlg a > 0, do not contribute to the curvature estimate obtained in Eq. (48).The acceleration due to surface tension and wetting experienced by a particle a ∈ Ω l is now given by Here, S n a is the Shepard filter computed by where Ω n l ⊂ Ω l is the subset of SPH particles satisfying ∥ñ a ∥ > ε n and Ω n s ⊂ Ω s is the subset of neighboring SPH particles satisfying r ab • nlg a ≥ 0 in Eq. (49).Using Eqs.(50) and (51), the acceleration of an SPH particle a due to surface tension and wetting for is now given by 6.2.Non-wetting contact angle, Θ ∞ > 90 Instead of computing the curvature using the approach presented in Sec.6.1, we model non-wetting contact angles by computing the curvature of an SPH particle a ∈ Ω l using only SPH particle neighbors b ∈ Ω l representing the liquid phase.In the following, we use the superscript l to denote that a property is computed using only the subset of SPH particles representing the liquid phase.Analogous to the correction scheme in Sec.6.1, the normal vector assigned to an SPH particle a Here, t sf,l a is the tangent vector at the position of a liquid SPH particle given by where the normal vector, nl a , is computed from liquid SPH particles by Here, S l a is the Shepard filter applied to liquid SPH particles Finally, the normal vector of an SPH particle a ∈ Ω l located in the vicinity of the three-phase contact line region is given by The modified normal vectors, nslg,l a , replace the normal vectors computed in Eq. (55) if the SPH particle is located in the vicinity of the three-phase contact line Finally, the mean curvature of a liquid SPH particle is given by where the renormalized kernel gradient, ∇l W ab , is computed from liquid SPH particles b ∈ Ω l by [31,32] ∇l Here, C l a is the correction matrix computed for SPH particles b ∈ Ω l .The acceleration of an SPH particle a ∈ Ω l due to surface tension and wetting forces is then given by with For desired equilibrium contact angles approaching 90 • , a slight discrepancy between the wetting and the non-wetting approach is observed, with the wetting approach yielding greater accuracy.

Extreme wetting and thin films
Obtaining stable drops on a plane surface with Θ ∞ ≤ 30 • is challenging because the particles near the three-phase contact line have much fewer neighbors within their support radius (b ∈ Ω l ) than for larger contact angles.Further, for low contact angles, the resulting sharp curvatures cannot be estimated accurately due to the smoothing of the normals (see Eqs. (58) and ( 49)).When the liquid phase spreads thinner than the SPH kernel radius on the substrate, the aforementioned approach will result in all particles being identified as near the contact line because all particles will be near the free surface and the substrate.Consequently, the normal correction in Eqs. ( 32) and (55) will be applied to all particles a ∈ Ω l in the liquid layer, resulting in unphysical surface forces.
For the case of small contact angle, as simulated in Sec.8.5, we substitute the curvature computation in Eq. ( 48) by Eq. ( 59) (where a renormalized kernel gradient is used) for an SPH particle a ∈ Ω slg if there is an insufficient neighbor particle support which we identify by S a < 0.6.
For a thin film of liquid, the thickness of the film may be less than the kernel support radius.This leads to incorrect identification of particles.However, these particles are not necessarily near the triple contact line, but their surface normals can be parallel to the substrate normals.Therefore, we use the threshold nlg a • nsf a ≤ 0.995, for the normal correction (Eq.44) of particles near the triple contact line to avoid the particles in a thin film being classified as located near the triple line.

Wall boundary
To model wetting phenomena on an impermeable substrate using SPH, the wall must satisfy the no penetration and the no-slip boundary conditions [38].Several techniques were proposed to ensure no penetration.For instance, wall boundaries can be modeled using repulsive forces [39], dummy particles [40], mirror particles [27], immersed boundary methods [41], or semi-analytical approaches [42].To simultaneously satisfy no-penetration and no-slip boundary conditions at walls, we modify the PPE in Eq. ( 19) following [43].We define the set of wall particles Ω w as a subset of solid particles, Ω w ⊂ Ω s .Wall particles, a ∈ Ω w , are thus near fluid particles.In mathematical terms, an SPH particle represents the solid wall boundary, a ∈ Ω w , if a ∈ Ω s and S l a > 10 −3 .
The threshold S l a > 10 −3 must be chosen large enough to allow the convergence of the iterative solver solving the PPE.The PPE for particles representing the wall (a Thus, the neighborhood of wall particles includes exclusively fluid particles.In summary, the pressure of a particle a follows from different PPEs depending on whether it is a liquid particle in the bulk of the fluid, a liquid particle near the free surface, or a wall particle: Eq.( 21) if a / ∈ Ω w and S a > 0.95 Eq.( 25) if a / ∈ Ω w and S a ≤ 0.95 Eq.(64) if a ∈ Ω w .

(65)
To enforce the no-slip boundary condition at the wall, the wall's velocity is calculated via an SPH interpolation of the liquid neighborhood of a ∈ ω w as Subsequently, the wall particle velocity, u w a , that replaces the velocity of neighboring particles b ∈ Ω w in Eq. (11), is given by u w a = 2u a − ũa . (67) 8. Validation of the proposed CSF and wetting model 8.1.Numerical parameters Table 1 summarizes the numerical parameters used in all tests unless otherwise stated.The parameter chosen, albeit arbitrary, belong to the realm of realistic fluids.Stable simulations with properties of water may also be achieved with a higher spatial resolution than provided.The integration time step was chosen due to the Courant-Friedrich-Lewy condition [14] from the maximum acceleration a max , the maximum velocity u max , the viscous diffusion, and the surface tension: The solid SPH particles (including the wall particles) are kept at zero velocity throughout the simulations.

Test case: Plane Poiseuille flow
We simulate plane Poiseuille flow between two stationary, infinite planes, located at x = ±L with L = 0.5 m.The fluid is accelerated in the y-direction by the body force g y .The analytical solution for the y-component of the time-dependent velocity, u y , as a function of the x-position reads [27] In the simulation, we apply periodic boundary conditions in y-and z-directions, thus, in this example, there is no free boundary.Consequently, to solve the PPE, we apply a Dirichlet boundary condition for the pressure to the solid SPH particles, which are not identified as wall SPH particles.The pressure of the solid phase is computed implicitly by solving Eq.( 64) if a ∈ Ω w Eq.( 25) if a ∈ Ω fs ∧ a / ∈ Ω w Eq.( 21) else . (70) In this validation example, we disregard surface tension effects, hence the momentum balance of the fluid reads The flow is characterized by the Reynolds number, Re = 2Lu max y /ν = 125.The SPH particles representing the liquid phase are initially placed on a rectangular lattice in the three-dimensional interval x ∈ (−L, L) 3 .At x = ±L, the boundary is modeled by SPH wall particles, placed at L ≤ x ≤ L + 2r max and L − 2r max ≤ x ≤ L. The thickness of the walls is twice the Wendland kernel radius, r max = 2h.This is required to apply the zero-pressure Dirichlet boundary condition in Eq. ( 25) to solid SPH particles, in order to solve the PPE.The liquid phase is represented by L/∆x ∈ {10, 20, 40} SPH particles in each spatial dimension leading to a total of {8 000, 64 000, 512 000} liquid SPH particles in the simulation.Figure 4 compares the analytical solution, Eq. ( 69), of the transient velocity field with the numerical result for L/∆x = 40.To evaluate the precision of the numerical model, we  compute the relative deviation of the numerical SPH particles' velocities from the corresponding analytical values for SPH particles located in the interval −∆x ≤ y ≤ ∆x and −∆x ≤ z ≤ ∆x, as a function of their x-position.For tν/L 2 = 4 for L/∆x = 40 we obtain the mean relative error 1.6 %.Reducing the spatial discretization to L/∆x = 10 we obtain the mean relative error 3.8 %.

Test case: Laplace pressure jump
At equilibrium, the pressure inside a drop exceeds the ambient pressure due to the surface tension.This effect is termed the Laplace pressure jump.In this simulation, we arrange N SPH particles on a rectangular lattice inside a sphere of radius r 0 = 10 −3 m, using the discretization given in Tab. 2. All smoothed normal vectors of magnitude smaller than ε n = 0.3 h are discarded from the computation of the curvature.According to the Young-Laplace equation, the pressure drop across a liquid-gas interface at equilibrium is given by [44] ∆p YL = 2σκ .
For the parameters given above, we find that the pressure inside the drop, p i , exceeds the ambient pressure by 20 Pa according to Figure 5 shows the equilibrium pressure along a cut through the symmetry axis of the drop after a sufficient relaxation time (t = 1 s) when equilibrium was achieved.For increasing number of SPH particles (increasing resolution), the equilibrium pressure inside the droplet converges to the analytical solution, given by Eq. ( 73).The relative error is 18.0 % when using 4 224 SPH particles, 0.8 % using 113 104 SPH particles, and 1.0 % when 268 096 SPH particles are used to discretize the droplet.

Test case: Droplet oscillations
In this validation case, we study the damped oscillation of a viscous drop.The analytical solution for the shape of a drop as a function of time reads [45] r Here, θ is the polar angle, r 0 is the radius of the (spherical) droplet in equilibrium, and is the second-order Legendre polynomial.The governing parameter  contains the oscillation frequency and the damping constant, For the initial condition, we assume [45] r (θ, which describes a weakly deformed sphere.For the simulation, we assume η = 5 × 10 −3 Pa s and the values given in Tab. 1.For initialization, the shape described in Eq. ( 77) is represented by 33 240 SPH particles placed on a square lattice, see Fig. 6. Figure 7 shows the major and minor axis oscillation caused by the initial deformation of the drop from its equilibrium shape, as obtained from the numerical simulation For a quantitative comparison with the analytical result, we fit the parameters of a damped harmonic oscillator, to the simulation data.We take the amplitude A osc from the initialization and obtain the best fit for the damping constant λ osc = 22.44 s −1 and the oscillation frequency ω osc = 261.5 s −1 .
We compare these numerical values with the analytical quantities.For the given material and system parameters, according to Eq. (76) , they read ω 2,0 = 283, 03 s −1 and λ 2 = 25 s −1 , resulting in the relative deviation |ω 2,0 − ω osc | /ω 2,0 ≈ 7% and |λ 2 − λ osc | /λ 2 ≈ 5.1% An appropriate value of the normal threshold, ε n , was found to be crucial for stable oscillations, see Eq. (55).Choosing ε n too small results in incorrect curvature computation due to contributions from SPH particles located far away from the interface.Choosing ε n too large results in large statistical errors since only a small number of SPH particles contributes to the surface tension forces.An appropriate choice of ε n compromises between these limits.

Test case: Equilibrium shape of a drop in contact with a plane
The numerical simulation of a drop in contact with a plane is a challenging problem.Depending on the choice of the equilibrium contact angle, Θ ∞ , which is a parameter of the simulation (see Eq. ( 39)), we describe wetting or de-wetting contact of the liquid drop on the solid substrate.Here, we demonstrate the stability of the numerical method by considering the relaxation of a particular initial state towards equilibrium for varying values of Θ ∞ .The initial condition is described by a hemispherical liquid drop resting on a solid substrate.We place liquid SPH particles on a square lattice with spacing ∆x inside a hemisphere of radius r 0 = 10 −3 m.The flat side of the hemisphere is in contact with the plane, modeled as a circular disk of radius 3 × 10 −3 m and thickness 5∆x.In all cases studied, the disk's radius was much larger than the equilibrium radius of the drop.The disk is represented by solid SPH particles arranged on a square lattice of spacing ∆x.At time t = 0, we start the simulation, and the initially hemispherical drop relaxes to its asymptotic equilibrium shape.
For the simulation, we use the parameters in Tab.(1), the smoothing length h = 2.5∆x, and ε n = 0.2/h.The equilibrium contact angles used to modify the normal vectors are shown in Tab. 3.
Table 3: Equilibrium contact angles used to correct the normal vectors of liquid and solid SPH particles in the vicinity of the three-phase contact line.For Θ∞ ≤ 95 • we employ Θ s ∞ ≥ Θ∞, see Sec.  decays by several orders of magnitude, see Fig. 9.We note that the asymptotic value of the kinetic energy of drops with non-wetting exceeds the value for wetting contacts by at least two orders of magnitude: At t = 1 s we find for wetting contact E kin ≈ 1.Similarly, the fluctuations of the kinetic energy are much larger for non-wetting contact than for wetting contact.This behavior agrees with physical reality: Drops of ultrahydrophobic surfaces are much less bound than wetting drops and can reveal interesting dynamics, e.g. the lotus effect [46].In terms of the numerical model, this effect can be understood from the influence of the normal threshold, ε n , used to identify SPH particles with valid normal vectors.The vectors of SPH particles that are located near the three-phase contact line have magnitude, below the threshold ε n and do, therefore, not experience surface tension.This leads to a vortex flow of particles in this region, which also promotes the translational motion of the liquid drop across the solid substrate which in turn increases the total kinetic energy.
Figure 10 shows the position of the free surface of the drop in equilibrium at position y = 0, that is, a cut through the drop, for wetting and non-wetting contact.The legend shows in brackets the desired equilibrium contact angle (Θ ∞ ) of the drop (input parameter), and the momentary value where H and B are the drop's height and base radius, and is the radius of the drop.The precision of the numerical model can be evaluated by comparing the simulation results for the drop's height and base radius with the analytical results,   10: Cut through a drop at y = 0 at large time, t = 1 s , when it assumed its equilibrium shape.The lines show the free liquid-gas interface for different contact angles.The legend shows the momentary value Θ of the contact angle defined in Eq. ( 79) and the equilibrium contact angle, Θ∞, in brackets with the analytical value of the drop's radius, Figure 11 shows this comparison for various contact angle values, Θ ∞ .
In the range 15 • ≤ Θ ∞ ≤ 135 • , the deviation of the simulated contact angle from the analytical value is below 4.5 %, see Fig. 11b, and increases to 14.3 % for Θ ∞ ≤ 150 • .This deviation is due to the unprecise contributions to the force from SPH particles located near the three-phase contact line.The precision could be improved by decreasing ε n , however, this would deteriorate the accuracy of the computed curvature.As in the preceding example, we have to compromise between large and small values of ε n .

Test case: Droplet deformation due to gravity
In the test cases presented so far, we did not consider the action of gravity.In the current test case, we study the deformation of a drop resting on a horizontal plate as a function of the value of the gravity constant, g.Gravity gives rise to a body-force acceleration, f g = [0, 0, −g], acting on the SPH particles.This force leads to flattening the drop, which shall be studied here.We use the properties given in Tab. 1 and the setup geometry introduced in Sec.8.5: A total of 16 776 liquid SPH particles are placed on a square lattice inside a hemisphere of radius r 0 = 10 −3 m, resting on a plane solid substrate of radius 4 × 10 −3 m.The equilibrium contact angle is Θ ∞ = 50 • (Θ s ∞ = 55 • ). Figure 12 shows the height of the drop, H, as a function of the acceleration due to gravity, g.Here, gravity is expressed by the Bond number that quantifies the ratio of gravitational to surface tension forces: The drop height, H, is scaled by the drop height in the absence of gravity H 0 where Bo = 0.In Fig. 12, the Bond number covers the interval Bo ∈ [10 −3 , 10 2 ], and the function H(Bo) decays monotonously.For small Bo, the situation is surface-tension dominated, thus, H → H 0 for Bo → 0.
For large gravity, Bo → ∞, the regime is gravity-dominated, and he height of the drop converges to the capillary length [47] H Figure 13 shows the drop profile for several values of Bo, that is, the position of the free surface at y = 0.For real-life applications, a dynamic contact angle model (such as [48]) may be implemented for taking into account chemical heterogeneity of the substrate.

Test case: Droplet pinning
The contact line of two planes with different inclinations represents a barrier for liquid droplets.To pass the barrier, the droplet's contact angle must exceed Θ ∞ + Ψ [24], where Ψ is the difference between the inclinations of the planes.To overcome this threshold, sufficient force is needed.In the absence of such a force, the droplet remains attached at the contact line of the planes, see Fig.    smooth counteracting force, enabling some SPH particles to pass the contact line.The last column of Fig. 15 shows the droplets at t = 150 ms.For Ψ ∈ {22.5 • , 45 • , 67.5 • } the drop passes the discontinuity.For Ψ = 90, the external acceleration is insufficient to generate the threshold contact angle Θ ∞ + Ψ = 165 • .Figure 16 shows the total kinetic energy of the liquid SPH particles as a function of time.
During the relaxation, t ∈ [0, 50 ms], the drop finds its equilibrium shape, therefore, energy decays.For t > 50 ms, the decline is accelerated by f b acting on the liquid SPH particles leading to a sharp increase in the total kinetic energy.At t ≈ 55 ms, the drop approaches the contact line between the two planes, indicated by a decrease of the kinetic energy for t > 55 ms.As expected, the smallest deceleration of the drop is experienced for Ψ = 22.5 • and the most significant for Ψ = 90 • .The subsequent increase in kinetic energy for t ≳ 85 ms, for ψ < 90 • corresponds to the drop's passing of the discontinuity.For Ψ = 90 • , the kinetic energy remains about three orders of magnitude below, indicating the pinning of the drop at the contact line.
The kinetic energy curves in Fig. 16 show sharp peaks.These peaks are caused by SPH particles on the rear side of the moving droplet.Due to the small number of liquid SPH particles attached to the solid phase behind the moving droplet, the calculated curvatures of these particles can adopt large, inaccurate values.As a result, these particles experience either high accelerations toward the solid substrate or into the surrounding environment.The wall boundary condition prevents liquid SPH particles from entering the solid plane due to the large pressure of the solid SPH particles.However, this also leads to a strong repulsion of the liquid particles into the environment.This problem could be solved by calculating the mean curvature of the liquid SPH particles only when there is a sufficient number of neighbors of the liquid phase.Otherwise, the particles should remain unaffected by the surface tension.

Conclusion
We introduce a surface tension model based on the CSF approach and applicable to threedimensional free surface flows.In contact with a substrate, the model can handle wetting and non-wetting behavior in a wide range of parameters.This is made possible with the use of primarily 3 empirical parameters, namely the Shephard sum correction for free surface introduced in Eq. 37, the threshold for magnitude of surface normal vectors (Eq.31), ϵ n and the desired equilibrium contact angle associated with the substrate Θ s ∞ listed in Table 2.We have shown that the model performs numerically stable in several test cases.By directly comparing analytical results and literature results, we could show that the model delivers quantitatively reliable results.To demonstrate the model's performance, we applied the new simulation method to a set of free surface problems, some of which are notoriously difficult to simulate.These test examples include (a) plane Poiseuille flow, (b) the Laplace jump of a droplet, (c) the oscillation and relaxation of a perturbed droplet, (d) the equilibrium shape of a droplet in contact with a wetting/non-wetting substrate and its relaxation to this equilibrium, (e) the flattening of a droplet under the action of gravity, and (f) the interaction of a droplet with a barrier under the action of gravity.For all of the test cases, we provide detailed quantitative comparisons with analytical results and results from the literature and discuss the reasons for deviations.

Figure 1 :
Figure 1: Sketch of a droplet in contact with a solid substrate in equilibrium

r ab • nlg a ≥ 0 Figure 2 :
Figure 2: Identification of neighboring particles b ∈ Ω s which contribute to the curvature computation of a particle a ∈ Ω l .This schematic shows SPH particles near the three-phase contact of a droplet resting on a solid boundary.Gray and blue spheres represent the solid boundary and the liquid phase, respectively.The normal vector nlg a is used to compute the distance of a neighboring particle b ∈ Ω s to the tangent plane (here shown as a dashed black line) by r ab • nlg a .Particles b ∈ Ω s which satisfy r ab • nlg a ≥ 0 contribute to the curvature computation and are shown by black spheres.

Figure 3 Figure 3 :
Figure3shows cuts through the axis of symmetry of two drops resting on a plane surface.For a desired equilibrium contact angle of Θ ∞ = 30 • (left) and Θ ∞ = 60 • (right), the proposed approach modifies the normal vectors of the red-colored SPH particles during the curvature computation.Note that each SPH particle a ∈ Ω l in the vicinity of the three-phase contact line has its own set of SPH particle neighbors b ∈ Ω s .

Figure 4 :
Figure 4: Instantaneous velocity profiles for viscous fluid flow between two parallel plates driven by the pressure gradient: the fluid flows in the y-direction and its non-dimensional velocity is plotted along the normal direction to the plates.

Figure 5 :
Figure 5: Comparison of the computed equilibrium pressure inside a droplet with the analytical solution of the Young-Laplace equation, Eq. (72).The simulation result is averaged over the time interval 0.8 ≤ t ≤ 1 s.

Figure 7 :
Figure 7: Damped oscillation of a viscous droplet.The curves show the drop extension along the major and minor axes as functions of time, normalized by the oscillation period, Tosc = 2π/ω 2,0 = 2.22 × 10 −2 s.

Figure 8 :
Figure 8: At large time, t = 1 s, the drops have assumed their equilibrium shape.Liquid SPH particles are shown in red, and solid particles are shown in blue.For Θ ∞ ∈ {120 • , 150 • }, some particles near the three-phase contact line disintegrated from the body of the liquid phase.Refer to Fig. 10 for the radial profile of the free surface.

Figure 9 :
Figure9: Total kinetic energy of relaxing droplets as a function of time.In agreement with physical reality, drops with wetting contact angles relax to a lower value of kinetic energy than drops with non-wetting contact angles.Similarly, the fluctuations are smaller for wetting contact

Figure
Figure10: Cut through a drop at y = 0 at large time, t = 1 s , when it assumed its equilibrium shape.The lines show the free liquid-gas interface for different contact angles.The legend shows the momentary value Θ of the contact angle defined in Eq. (79) and the equilibrium contact angle, Θ∞, in brackets

6 ΘFigure 11 :
Figure 11: Comparison between simulation and theory for the drop's base radius, height, and simulated contact angle, as a function of the equilibrium contact angle

14 .
We describe the numerical experiment similar to Secs.8.5 and 8.6 with the parameters specified in Tab. 1, viscosity η = 5 × 10 −3 Pa s, and the equilibrium contact angle Θ ∞ = 75 • (Θ s ∞ = 80 • ).Similar to the preceding test cases, we fill a hemisphere with radius r 0 = 10 −3 m with 16 776 liquid SPH particles arranged on a square lattice.This half-sphere rests on a solid horizontal plane modeled as a block of size (2.6 × 3 × 0.25) mm 3 filled by SPH particles arranged in a square lattice.The initial condition is sketched in the left part of Fig. 14.The second plane 10 −3 is modeled in the same way but inclined by Ψ ∈ {22.5 • , 45 • , 67.5 • }, see right part of Fig.14.To prepare our initial conditions, we simulate relaxing the drop on the horizontal substrate to equilibrium.After t = 50 ms of relaxation, the equilibrium shape was assumed.For t > 50 ms a constant body acceleration, f b = [7.5, 0, −7.5], drives the drop towards the contact line of the planes.Figure15shows snapshots of the simulation.Column (a) shows the equilibrated drops at t = 50 ms.Column (b) shows the drops at t = 100 ms under the influence of the acceleration f b .Due to the weighting concept of the SPH method, pinning starts to occur when the contact line of the two planes enters the smoothing length of an SPH particle.Instead of a sharp force, the particles experience a

Figure 13 :Figure 14 :
Figure 13: Cut through the drop at y = 0 (liquid-gas interface position) as a function of the Bond number for Θ∞ = 50 • .

Table 1 :
Numerical parameters used in the test simulations

Table 2 :
Number of SPH particles and discretization spacing used to validate the Young-Laplace pressure boundary condition.
16 × 10 −14 J for Drop height as a function of the Bond number.For small gravity, Bo → 0, the regime is dominated by surface tension, and H approaches H 0 .For large gravity, Bo → ∞, the regime is gravity-dominated, and H approaches H∞ given by Eq. (85).The symbols show the drop height obtained from SPH simulations, and the solid lines show the limits H 0 and H∞.