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.

1. Introduction

At the microscale and nanoscale, 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 figure 1. Surface tension and wetting are essential for many phenomena, including technological processes such as oil recovery (Babadagli 2002; Babadagli & Boluk 2005), two-phase heat transfer (Ebadian & Lin 2011) and ink-jet printing (Calvert 2001). The wetting properties of a solid surface depend on its chemical composition and morphology, determining its repellent or attractive behaviour (Kam, Bhattacharya & Mazumder 2012; Lai et al. 2013).

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

In smoothed particle hydrodynamics (SPH) (Gingold & Monaghan 1977; Lucy 1977), surface tension can be modelled either by pairwise forces between the particles (Nugent & Posch 2000; Tartakovsky & Meakin 2005; Tartakovsky & Panchenko 2016; Nair & Pöschel 2018) or by a phenomenological continuum surface force (CSF) (Brackbill, Kothe & Zemach 1992; Morris 2000; Sirotkin & Yoh 2012). The concept of pairwise forces between SPH particles is motivated by the molecular origin of surface tension (Rowlinson & Widom 2013) 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 (Shigorina, Kordilla & Tartakovsky 2017) and contact angle hysteresis (Bao et al. 2019). When used in simulations, however, such models require a large support domain for particles, affecting the computational performance significantly (Nair & Pöschel 2018).

A CSF can reliably predict the dynamics due to surface tension gradients at phase boundaries, known as the Marangoni effect. A CSF is also widely used in Eulerian, mesh-based two phase flow simulations and in the lattice Boltzmann method (Bagheri, Siavashi & Mousavi 2023), where the dynamics of the contact line is only empirically modelled (see e.g. Bogdanov et al. (2022)), relating the contact angle to the velocity, based on experimental observations (Fries & Dreyer 2008; Saha & Mitra 2009). Eulerian mesh-based methods can also employ phase fields for simulating wetting phenomena, which are limited to low density ratio. In the SPH method, known for its Lagrangian nature and mass conservation at free surfaces, the application of CSF to the free surface 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 modelling surface tension at the free surface has severely restricted the scope of three-dimensional 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, Yaghoubi & Yeganehdoust (2017) to produce robust two-dimensional flow simulations in weakly compressible SPH. Hirschler et al. (2017) 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, Geara et al. (2022) introduced an analytical coefficient to account for truncation of kernels near a free surface in SPH operators. The curvature estimate at the free surface was recently improved by Fürstenau, Weißenfels & Wriggers (2020), while Blank, Nair & Pöschel (2023) 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 modelled using the smoothed normal correction scheme introduced by Breinlinger et al. (2013). 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 plane Poiseuille flow, the Laplace jump of a droplet, the oscillation and relaxation of a perturbed droplet, the equilibrium shape of a droplet in contact with a wetting/non-wetting substrate and its relaxation to this equilibrium, the flattening of a droplet under the action of gravity and the interaction of a droplet with a barrier under the action of gravity.

2. Smoothed particle hydrodynamics

Smoothed particle hydrodynamics is a numerical method for solving continuum problems (Gingold & Monaghan 1977; Lucy 1977; Monaghan 2005). Unlike in mesh-based methods, the domain is discretized by SPH particles with a certain mass at which properties such as velocity, $\boldsymbol {u}$, pressure, $p$, or density, $\rho$, are defined.

The integral interpolant, $\varPhi ^{I}(\boldsymbol {r})$, of a physical property, $\varPhi (\boldsymbol {r})$, at position, $\boldsymbol {r}$, can be obtained from a spatial convolution with a kernel function, $W(\boldsymbol {r}-\boldsymbol {r}',h),$ having a compact support,

(2.1)\begin{equation} \varPhi^{I} (\boldsymbol{r})= \int_{\varOmega_c} \varPhi(\boldsymbol{r}') W(\boldsymbol{r}-\boldsymbol{r}',h) \,{\rm d}\boldsymbol{r}^\prime. \end{equation}

Here, $h$ is the smoothing length, and $\varOmega _c$ denotes the whole continuous domain. We approximate the integral in (2.1) by a sum over the set $\varOmega$ of all SPH particles,

(2.2)\begin{equation} \varPhi^{s}\left(\boldsymbol{r}\right) = \sum_b \frac{m_b}{\rho_b} \varPhi\left(\boldsymbol{r}_b\right) W\left(\boldsymbol{r}-\boldsymbol{r}_b,h \right), \end{equation}

where $\boldsymbol {r}_b$, $m_b$ and $\rho _b$ are the position, mass and density of particle $b$. In particular, the summation interpolant at the position $\boldsymbol {r}_a$ of an SPH particle reads

(2.3)\begin{equation} \varPhi^{s}\left(\boldsymbol{r}_a\right) = \sum_b \frac{m_b}{\rho_b} \varPhi\left(\boldsymbol{r}_b\right) W\left(\boldsymbol{r}_a-\boldsymbol{r}_b,h \right). \end{equation}

In the same approximation, the spatial derivative of a physical quantity, $\boldsymbol {\nabla }\varPhi (\boldsymbol {r})$, at position $\boldsymbol {r}_a$, can be expressed in terms of the kernel gradient $\boldsymbol {\nabla } W(\boldsymbol {r}_a-\boldsymbol {r}_b,h )$ by

(2.4)\begin{equation} \boldsymbol{\nabla}\varPhi^{s}\left(\boldsymbol{r}_a\right) = \sum_b \frac{m_b}{\rho_b} \varPhi\left(\boldsymbol{r}_b\right) \boldsymbol{\nabla} W\left(\boldsymbol{r}_a-\boldsymbol{r}_b,h \right). \end{equation}

We introduce shorthand notations and rewrite (2.3) and (2.4),

(2.5a,b)\begin{equation} \varPhi^{s}_{a} = \sum_b \frac{m_b}{\rho_b} \varPhi_b W_{ab},\quad \boldsymbol{\nabla} \varPhi^{{s}}_{a} = \sum_b \frac{m_b}{\rho_b} \varPhi_b \boldsymbol{\nabla} W_{ab} . \end{equation}

Here we use the $C^2$ kernel function by Wendland (1995) normalized for three spatial dimensions,

(2.6a,b)\begin{align} W_{ab}\equiv W\left(\boldsymbol{ r}_a - \boldsymbol{ r}_b, h\right) = \begin{cases} \dfrac{21}{16{\rm \pi}} \dfrac{1}{h^3}\left(1-\dfrac{r_{ab}}{2h}\right)^4\left(\dfrac{2r_{ab}}{h}+1\right) & r_{ab} \leq 2h \\ 0 & r_{ab} > 2h \end{cases},\quad r_{ab} \equiv \left\lVert\boldsymbol{ r}_a - \boldsymbol{ r}_b\right\rVert. \end{align}

We describe the dynamics of a Newtonian incompressible fluid using the Navier–Stokes equation

(2.7)\begin{equation} \frac{\mathrm{D}\boldsymbol{u}}{\mathrm{D} t} =-\frac{1}{\rho}\boldsymbol{\nabla} p + \nu \nabla^{2} \boldsymbol{u} + \boldsymbol{f}^{{s}} + \boldsymbol{f}^{{b}}, \end{equation}

and the continuity equation

(2.8)\begin{equation} \frac{\partial \rho}{ \partial t} =- \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho \boldsymbol{u}\right), \end{equation}

where $\boldsymbol {u}$, $\rho$, $p$, $\nu$, are the velocity, density, pressure and kinematic viscosity, $\boldsymbol {f}^{{b}}$ is an acceleration due to a body force such as gravity and $t$ is time. The acceleration caused by surface tension, $\boldsymbol {f}^{{s}}$, is incorporated into (2.7) as a body force according to the CSF model by Brackbill et al. (1992).

The total acceleration of a particle $a$ given by (2.7) is, thus, a sum of accelerations

(2.9)\begin{equation} \frac{\mathrm{D} \boldsymbol{u}_a}{\mathrm{D} t} = \boldsymbol{f}^{{p}}_a + \boldsymbol{f}^{{v}}_a + \boldsymbol{f}^{{s}}_a + \boldsymbol{f}^{{b}}_a \end{equation}

due to pressure gradient, $\boldsymbol {f}^{p}_{a}$, viscosity, $\boldsymbol {f}^{v}_{a}$, surface tension, $\boldsymbol {f}^{{s}}_a$ and body forces per unit mass, $\boldsymbol {f}^{{b}}_a$.

The first contribution, the acceleration due to pressure gradient, can be approximated by (Monaghan 2005)

(2.10)\begin{equation} \boldsymbol{f}^{p}_{a} =- \sum_b m_b\left(\frac{p_a }{\rho_a^2} + \frac{p_b }{\rho_b^2}\right) \boldsymbol{\nabla} W_{ab}, \end{equation}

where $p_a \equiv p(\boldsymbol { r}_a)$, $p_b \equiv p(\boldsymbol { r}_b)$, $\rho _a \equiv \rho (\boldsymbol { r}_a)$, $\rho _b \equiv \rho (\boldsymbol { r}_b)$. Equation (2.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 § 4.

The second contribution to the acceleration in (2.9), the viscous acceleration of particle $a$, is given by (Morris, Fox & Zhu 1997)

(2.11a,b)\begin{equation} \boldsymbol{f}^{v}_{a} = \sum_b \frac{m_b\left(\eta_a + \eta_b\right) \boldsymbol{u}_{ab}}{\rho_a\rho_b} F_{ab},\quad \boldsymbol{u}_{ab} \equiv \boldsymbol{u}_a - \boldsymbol{u}_b, \end{equation}

where $\boldsymbol {u}_a$ and $\boldsymbol {u}_b$ are the velocities of particle $a$ and $b$, $\eta _a$ and $\eta _b$ are the dynamic viscosity coefficients assigned to particles $a$ and $b$, and $F_{ab}$ is given by

(2.12a,b)\begin{equation} F_{ab} \equiv F\left(\boldsymbol{r}_{a}-\boldsymbol{r}_{b}\right) = \frac{\boldsymbol{r}_{ab} \boldsymbol{\cdot} \boldsymbol{\nabla} W_{ab} }{r_{ab}^2 + \left(0.01h\right)^2}, \quad \boldsymbol{r}_{ab}\equiv \boldsymbol{r}_a-\boldsymbol{r}_b. \end{equation}

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 (2.9), the acceleration of a particle $a$ due to surface tension, $\boldsymbol {f}^{{s}}$, will be described in § 5.

3. Incompressible smoothed particle hydrodynamics

For incompressible fluids, the velocity field is divergence-free, and the continuity equation, (2.8), turns into

(3.1)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0. \end{equation}

The idea of ISPH is to enforce (3.1) by imposing a pressure field in such a way that the velocity field becomes divergence-free (Cummins & Rudman 1999). In ISPH, we use the predictor–corrector integration scheme by Cummins & Rudman (1999) 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,

(3.2)\begin{equation} \boldsymbol{r}_{a}^\ast= \boldsymbol{r}_{a}^n + \boldsymbol{u}_{a}^n \Delta t;\quad a=1,2,3,\dots \, . \end{equation}

The corresponding velocities are

(3.3)\begin{equation} \boldsymbol{u}_{a}^\ast= \boldsymbol{u}_{a}^n + \left(\boldsymbol{f}_{a^\ast}^{{v}} + \boldsymbol{f}_{a^\ast}^{{s}} +\boldsymbol{f}_{a^\ast}^{{b}}\right)\Delta t \quad \text{with}\ \boldsymbol{f}_{a^\ast}^{v} \equiv \boldsymbol{f}^{v}\left(\boldsymbol{r}^\ast_a\right), \quad \text{etc.} \end{equation}

The pressure, $p^\ast _a$, corresponding to the predicted velocity can be obtained by solving the pressure Poisson equation (PPE)

(3.4)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \frac{\boldsymbol{\nabla} p_{a}^{{\ast}}}{\rho_{a}}\right) = \frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}_{a}^\ast}{\Delta t}. \end{equation}

From the predicted velocities, $\boldsymbol {u}_a^\ast$, and the corresponding pressure, $p^\ast _a$, the correction step yields the divergence-free velocities at time step $n+1$,

(3.5)\begin{equation} \boldsymbol{u}_{a}^{n + 1} = \boldsymbol{u}_{a}^\ast- \frac{1}{\rho_{a}}\boldsymbol{\nabla} p_{a}^{{\ast}}\Delta t. \end{equation}

The corresponding particle positions at time step $n+1$ are

(3.6)\begin{equation} \boldsymbol{r}_{a}^{n + 1} = \boldsymbol{r}_{a}^n + \left(\frac{\boldsymbol{u}_{a}^{n} + \boldsymbol{u}_{a}^{n + 1}}{2}\right) \Delta t. \end{equation}

The left-hand side of the PPE, (3.4), can be discretized to obtain (Cummins & Rudman 1999)

(3.7)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\left( \frac{\boldsymbol{\nabla} p_{a}^{{\ast}}}{\rho_{a}}\right) = \sum_b \frac{m_b}{\rho_b} \frac{4}{\rho_a + \rho_b} \left(\,p_a^{{\ast}} - p_b^{{\ast}}\right) F_{ab}, \end{equation}

and its right-hand side results in (Monaghan 2005)

(3.8a,b)\begin{equation} \frac{\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{a}^{{\ast}}}{\Delta t} = \sum_b -\frac{m_b}{\rho_b} \frac{\boldsymbol{u}^{{\ast}}_{ab} \boldsymbol{\cdot}\boldsymbol{\nabla} W_{ab}}{\Delta t}, \quad \boldsymbol{u}^{{\ast}}_{ab} \equiv \boldsymbol{u}^{{\ast}}_{a} - \boldsymbol{u}^{{\ast}}_{b}. \end{equation}

Substituting (3.7) and (3.8a,b) into (3.4) yields the discretized PPE

(3.9)\begin{equation} \sum_b \frac{m_b}{\rho_b} \frac{4}{\rho_a + \rho_b} \left(\,p_a^{{\ast}} - p_b^{{\ast}}\right) F_{ab} = \sum_b -\frac{m_b}{\rho_b} \frac{\boldsymbol{u}^{{\ast}}_{ab}\boldsymbol{\cdot} \boldsymbol{\nabla} W_{ab}}{\Delta t}. \end{equation}

Equation (3.9) can be uniquely solved using any iterative linear solver such as the BiCGSTAB (biconjugate gradient stabilized) method by Sleijpen, Van der Vorst & Fokkema (1994), provided the linear system is not singular, e.g. if a Dirichlet boundary condition is applied. The discretized PPE, (3.9), is an appropriate approximation of (3.4) 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.

4. Kernel correction for particles near a free surface

In the vicinity of boundaries, the summations in (2.3) and (2.4) underestimate $\varPhi ^{I}(\boldsymbol {r})$ and $\boldsymbol {\nabla }\varPhi ^{I}(\boldsymbol {r})$, due to lacking SPH particles in the neighbourhood of $\boldsymbol {r}$. We call such SPH particles insufficiently supported by neighbouring SPH particles. To account for such particles, semianalytical (Nair & Tomar 2014) or numerical normalization techniques (Bonet & Lok 1999; Oger et al. 2007) are used near free boundaries. In the current paper, akin to Nair & Tomar (2014) 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 (Shepard 1968),

(4.1)\begin{equation} S_a = \sum_b \frac{m_b}{\rho_b} W_{ab}. \end{equation}

However, unlike the approaches described by Lee et al. (2008) where a thin layer of particles comprises the pressure boundary condition (which deteriorates the accuracy and the stability (see Asai et al. (2012) and Nair & Tomar (2014)), 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, $\varOmega ^{{b}}\subset \varOmega$, or to the subset representing the boundary, $\varOmega ^{{fs}}\subset \varOmega$ is determined by

(4.2)\begin{equation} \begin{cases} a \in \varOmega^{{fs}} & \text{if}\ S_a \leq 0.95 \\ a \in \varOmega^{{b}} & \text{else} \end{cases}. \end{equation}

For well supported particles, $a\in \varOmega ^{{b}}$, we evaluate (2.10) using the pressure gradient approximation (Monaghan 2005); for $a \in \varOmega ^{{fs}}$ we use the approximation given by (Blank et al. 2023),

(4.3)\begin{equation} \boldsymbol{f}^{p}_{a} = \begin{cases} - \sum\limits_b m_b\left(\dfrac{p_a }{\rho_a^2} + \dfrac{p_b }{\rho_b^2}\right) \boldsymbol{\nabla} W_{ab} & a\in\varOmega^{{b}}\\ -\sum\limits_{b} m_{b} \left( \dfrac{p_b- p^{o}_{a} }{\rho_{b}^2} \right) \boldsymbol{\nabla} W_{a{b}} & a\in\varOmega^{{fs}} \end{cases}. \end{equation}

Here, $p_a^{{o}}$ 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 (2.11a,b). 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 (3.9) is zero. Therefore, (2.11a,b) and (3.9) 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 (3.9) for SPH particles that are located in the vicinity of a free surface, we use (Nair & Tomar 2014; Blank et al. 2023)

(4.4)\begin{align} &\left(\,p_a - p_{a}^{o}\right) \underbrace{\sum_{b}\frac{m_b}{\rho_b}\frac{4}{\rho_a + \rho_b} F_{ab}}_{\beta} - \sum_{b} \frac{m_{b}}{\rho_{b}} \frac{4p_{b}}{\rho_a + \rho_{b}} F_{a{b}} \nonumber\\ &\quad = \sum_{b} \frac{m_{b}}{\rho_{b}} \left( -\frac{\boldsymbol{u}^{{\ast}}_{a{b}}\boldsymbol{\cdot}\boldsymbol{\nabla} W_{a{b}}}{\Delta t} - \frac{4p_{a}^{o}}{\rho_a + \rho_{b}}F_{a{b}}\right). \end{align}

Note that the term marked by $\beta$ is constant if the mass and volume assigned to the SPH particles are uniform throughout the domain.

In conclusion, we approximate the PPE by

(4.5)\begin{equation} \begin{cases} (3.9) & \text{if}\ S_a > 0.95\\ (4.4) & \text{if}\ S_a \leq 0.95 \end{cases}. \end{equation}

This ambient pressure application can be used to model surface tension by computing the Young–Laplace pressure boundary condition, as shown by Blank et al. (2023).

5. Surface tension at free boundaries

The CSF model by Brackbill et al. (1992) evaluates the acceleration of a particle $a$ due to surface tension by

(5.1)\begin{equation} \boldsymbol{f}^{{s}}_a = \frac{1}{\rho_a} \boldsymbol{F}^{{s}}_a, \end{equation}

where $\boldsymbol {F}^{{s}}_a$ is a force per unit volume (Morris 2000),

(5.2)\begin{equation} \boldsymbol{F}^{{s}}_a = \left(2\sigma_a \kappa_a \hat{\boldsymbol{n}}_a + \boldsymbol{\nabla}^{{s}} \sigma_a\right) \delta^{s}_a. \end{equation}

Here, the index $a$ to $\boldsymbol {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, $\sigma _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, $\kappa _a$, is the local value of the mean curvature of the surface, assigned to particle $a$. It is termed the ‘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 $\hat {\boldsymbol {n}}_a$ assigned to particle $a$ is called the ‘unit normal vector of particle $a$’ and can be computed by

(5.3a,b)\begin{equation} \boldsymbol{n}_a =-\sum_{b}\frac{m_b}{\rho_b} \boldsymbol{\nabla} W_{ab},\quad \hat{\boldsymbol{n}}_a = \frac{\boldsymbol{n}_a}{\left\lVert\boldsymbol{n}_a\right\rVert}. \end{equation}

The same applies to the surface gradient of the surface tension, $\boldsymbol {\nabla }^{{s}} \sigma _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, $\delta ^{{s}}_a$, that peaks at the interface.

The second term in (5.2) 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 (5.2) vanish. The first term in (5.2) 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 (5.2) by (Fürstenau et al. 2020)

(5.4)\begin{equation} \delta^{{s}}_a = \left\lVert\boldsymbol{n}_a\right\rVert,\quad a\in \varOmega^{l}, \end{equation}

where $\varOmega ^{l}\subset \varOmega$ is the subset of SPH particles representing the liquid phase.

Both the absolute value and the direction of $\boldsymbol {n}_a$ obtained from (5.3a,b) depend sensitively on the positions of neighbouring particles $b$, therefore, naïve computation of the normal unit vector, $\hat {\boldsymbol {n}}_a \equiv \boldsymbol {n}_a/\left \lVert \boldsymbol {n}_a\right \rVert$, is subject to large fluctuations leading eventually to inaccurate $\boldsymbol {F}^{{s}}_a$ when used in (5.2). Therefore, instead of using the naïve value, we use a smoothed normal vector (Morris et al. 1997; Fürstenau et al. 2020),

(5.5)\begin{equation} \tilde{\boldsymbol{n}}_a = \frac{1}{S_a}\sum_{b} \frac{m_b}{\rho_b} \boldsymbol{n}_b W_{ab}. \end{equation}

The smoothing increases the region where particles contribute to surface tension. Following our previous notation, $\boldsymbol {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\in \varOmega ^{l}$, whose normal vector's magnitude is smaller than a threshold, $\varepsilon ^{{n}}\in [{0.1}/{h}, {0.2}/{h}]$. In the simulations presented in this paper, we have used $\varepsilon ^{{n}}=0.01$ throughout. Subsequently, the smoothed and filtered normal vectors from (5.5) are normalized by

(5.6)\begin{equation} \hat{\tilde{\boldsymbol{n}}}_a = \frac{\tilde{\boldsymbol{n}}_a}{\left\lVert \tilde{\boldsymbol{n}}_a\right\rVert}. \end{equation}

The set of unit vectors provided by (5.6), is used to compute the mean curvature which is needed to evaluate the first term in (5.2) (Monaghan 1992; Morris 2000),

(5.7)\begin{equation} {\kappa}_a =-\frac{1}{2} \sum_b\frac{m_b}{\rho_b} \left(\hat{\tilde{\boldsymbol{n}}}_a - \hat{\tilde{\boldsymbol{n}}}_b \right) \boldsymbol{\cdot}\boldsymbol{\nabla} W_{ab}. \end{equation}

To increase the accuracy of the curvature approximation in (5.7), we substitute the kernel gradient with $\hat {\boldsymbol {\nabla }}W_{ab} \equiv \boldsymbol {C}_a \boldsymbol {\nabla } W_{ab}$ (Bonet & Lok 1999; Oger et al. 2007), where

(5.8)\begin{equation} \boldsymbol{C}_a = \left(\sum_{b}\frac{m_b}{\rho_b}\boldsymbol{ r}_{ab} \otimes \boldsymbol{\nabla} W_{ab} \right)^{-1} \end{equation}

is the correction matrix computed from all neighbouring SPH particles $b$. This yields the following curvature approximation with $O(h^2)$ convergence (Bonet & Lok 1999; Oger et al. 2007):

(5.9)\begin{equation} {\kappa}_a =-\frac{1}{2} \sum_b\frac{m_b}{\rho_b} \left(\hat{\tilde{\boldsymbol{n}}}_a - \hat{\tilde{\boldsymbol{n}}}_b \right) \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} W_{ab}. \end{equation}

Using (5.2), (5.1),(5.6), (5.9) and assuming a constant surface tension coefficient in the simulated material, yields the acceleration of a particle $a$ in the normal direction

(5.10)\begin{equation} \boldsymbol{f}_{a}^{{s}} = 2 \frac{\sigma_{a}}{\rho_a} {\kappa}_{a} \hat{\tilde{\boldsymbol{n}}}_a \delta_a^{s}. \end{equation}

Similar to Blank et al. (2023), where the Shepard filter is used to increase the robustness and accuracy of the obtained Young–Laplace pressure jump, we modify (5.10) to

(5.11)\begin{equation} \boldsymbol{f}_{a}^{{s}} = \left(1+\frac{1}{S^{{n}}_a}\right) \frac{\sigma_{a}}{\rho_a} {\kappa}_{a} \hat{\tilde{\boldsymbol{n}}}_a \delta_a^{s}. \end{equation}

Here, $S^{{n}}_a$ is the Shepard filter computed by

(5.12)\begin{equation} S^{{n}}_a = \sum_{b\in\varOmega^{n}} \frac{m_b}{\rho_b} W_{ab}, \end{equation}

where $\varOmega ^{n} \subset \varOmega$ is the subset of SPH particles satisfying $\left \lVert \tilde {\boldsymbol {n}}_a\right \rVert \geq \varepsilon ^{n}$. The use of the factor $( 1+1/S_a^n)$ which replaces the theoretical factor of $2$ is empirical and is obtained from numerical experimentation.

6. 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, $\varTheta$, and the equilibrium contact angle, $\varTheta _{\infty }$, by (de Gennes 1985)

(6.1)\begin{equation} {S}^{{w}} = \sigma \left(\cos \varTheta_{\infty} - \cos\varTheta \right). \end{equation}

According to Breinlinger et al. (2013), (6.1) 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 (5.11) includes the acceleration due to wetting phenomena. The scheme by Breinlinger et al. (2013) 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 modelling wetting and non-wetting contact for free-surface problems.

6.1. Wetting contact angle, $\varTheta _{\infty } \leq 90^\circ$

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. Following Brackbill et al. (1992) the normal vector associated with a desired equilibrium contact angle, $\varTheta _{\infty }$, can be computed for a particle $a\in \varOmega ^{l}$ by

(6.2)\begin{equation} \hat{\tilde{\boldsymbol{n}}}^{\infty}_{a} = \hat{\tilde{\boldsymbol{t}}}_a^{{sf}} \sin \varTheta_{\infty} - \hat{\tilde{\boldsymbol{n}}}_a^{{sf}} \cos \varTheta_{\infty}. \end{equation}

Here, $\hat {\tilde {\boldsymbol {t}}}_a^{{sf}}$ is the smoothed normalized tangent vector between the solid and fluid phases, and $\hat {\tilde {\boldsymbol {n}}}_a^{{sf}}$ is the smoothed normalized normal vector pointing from the solid to the fluid phase computed at the position of particle $a$. The normal vector, $\hat {\tilde {\boldsymbol {n}}}_a^{{sf}}$, at the position of particle $a$, is given by

(6.3a–c)\begin{equation} \hat{\tilde{\boldsymbol{n}}}_a^{{sf}} = \frac{\tilde{\boldsymbol{n}}_a^{{sf}}}{\lVert\tilde{\boldsymbol{n}}_a^{{sf}}\rVert}, \quad \tilde{\boldsymbol{n}}_a^{{sf}} = \frac{1}{S_a^{{s}}} \sum_{b \in \varOmega^{{s}}} \frac{m_b}{\rho_b} \boldsymbol{n}_a^{{sf}} W_{ab}, \quad \boldsymbol{n}_a^{{sf}} =-\sum_{b \in \varOmega^{{s}}} \frac{m_b}{\rho_b} \boldsymbol{\nabla} W_{ab}, \end{equation}

where $\varOmega ^{{s}}\subset \varOmega$ is the subset of all SPH particles representing the solid phase, and $S_a^{{s}}$ is the Shepard filter computed by

(6.4)\begin{equation} S_a^{{s}} = \sum_{b \in \varOmega^{{s}}} \frac{m_b}{\rho_b} W_{ab}. \end{equation}

Analogous to the computation of $\hat {\tilde {\boldsymbol {n}}}_a$, we discard the unit normal vectors of SPH particles if $\left \lVert \tilde {\boldsymbol {n}}\right \rVert _a < \varepsilon ^{n}$. The unit tangent vector in (6.2), $\hat {\tilde {\boldsymbol {t}}}_a^{{sf}}$, is computed by

(6.5a,b)\begin{equation} \tilde{\boldsymbol{t}}^{{sf}}_{a} = \hat{\tilde{\boldsymbol{n}}}_a - \left(\hat{\tilde{\boldsymbol{n}}}_a \boldsymbol{\cdot} \hat{\tilde{\boldsymbol{n}}}^{{sf}}_a \right)\hat{\tilde{\boldsymbol{n}}}^{{sf}}_a ,\quad \hat{\tilde{\boldsymbol{t}}}^{{sf}}_{a} = \frac{\tilde{\boldsymbol{t}}_{a}^{{sf}}}{\lVert\tilde{\boldsymbol{t}}_{a}^{{sf}}\rVert}. \end{equation}

The normal vectors computed in (6.2) could be used to replace the normal vectors given by (5.6) to model wetting.

However, as shown by Breinlinger et al. (2013), 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 $\hat {\tilde {\boldsymbol {n}}}_a$ to $\hat {\tilde {\boldsymbol {n}}}^{\infty }_a$ (6.2) depending on a particle's distance to the solid phase (wall) should be applied:

(6.6)\begin{equation} \hat{\tilde{\boldsymbol{n}}}_a^{{lg}} = \frac{f_a \hat{\tilde{\boldsymbol{n}}}_a + \left(1 - f_a\right)\hat{\tilde{\boldsymbol{n}}}^{{\infty}}_{a}}{\left\lVert f_a \hat{\tilde{\boldsymbol{n}}}_a + \left(1 - f_a\right)\hat{\tilde{\boldsymbol{n}}}^{\infty}_{a}\right\rVert}. \end{equation}

Here, $f_a$ is given by

(6.7)\begin{equation} f_a =\begin{cases} 0 & \forall \ d_{a}^{{s}} < 0,\\ d_{a}^{{s}}/r_{max} & \forall \ 0 \leq d_{a}^{{s}} \leq r_{max}\\ 1 & \forall \ d_{a}^{{s}} > r_{max} \end{cases}, \end{equation}

where $r_{max}$ is the kernel radius (${r}_{max} = 2h$ when using the Wendland kernel in (2.5a,b)), and $d_{a}^{{s}}$ is the shortest distance of a particle $a\in \varOmega ^{l}$ to the particles $b \in \varOmega ^{s}$:

(6.8)\begin{equation} d_{a}^{{s}} = \mathrm{min} \left(\left(\boldsymbol{r}_a - \boldsymbol{r}_b\right) \boldsymbol{\cdot} \hat{\tilde{\boldsymbol{n}}}_{ a}^{{sf}}\right) - \Delta x. \end{equation}

Here, $\Delta x$ is the spacing between two SPH particles when arranged on a square lattice, and $\Delta 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 $\Delta x$ are prescribed with $\hat {\tilde {\boldsymbol {n}}}_a^\infty$. 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\in \varOmega ^{l}$ can be computed by

(6.9)\begin{equation} {\kappa}_a =-\frac{1}{2}\sum_{b} \frac{m_b}{\rho_b}\left(\hat{\tilde{\boldsymbol{n}}}_a^{{lg}} - \hat{\tilde{\boldsymbol{n}}}_b^{{lg}} \right)\hat{\boldsymbol{\nabla}} W_{ab}. \end{equation}

In the course of this work, it was found that computing the curvature in (6.9) is not accurate enough to obtain stable equilibrium droplet shapes for $\varTheta _{\infty }\leq 90^{\circ }$. 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\in \varOmega ^{l}$ by

(6.10)\begin{equation} {\kappa}_a =-\frac{1}{2}\sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b}\left(\hat{\tilde{\boldsymbol{n}}}_a^{{lg}} - \hat{\tilde{\boldsymbol{n}}}_b^{{lg}} \right)\hat{\boldsymbol{\nabla}} W_{ab} -\frac{1}{2}\sum_{b \in \varOmega^{{s}}} \frac{m_b}{\rho_b}\left(\hat{\tilde{\boldsymbol{n}}}_a^{{lg}} - \hat{\tilde{\boldsymbol{n}}}^{{s}}_{b}\right)\hat{\boldsymbol{\nabla}} W_{ab}, \end{equation}

where $\hat {\boldsymbol {n}}^{{s}}_b$ is the normal vector of a neighbouring particle $b\in \varOmega ^{s}$ (particles representing the solid phase) given by

(6.11)\begin{equation} \hat{\tilde{\boldsymbol{n}}}^{{s}}_{b} = \begin{cases} \hat{\tilde{\boldsymbol{t}}}_a^{{sf}} \sin \varTheta_{\infty}^{{s}} - \hat{\tilde{\boldsymbol{n}}}_a^{{sf}} \cos\varTheta_{\infty}^{{s}} & \text{if}\ \boldsymbol{r}_{ab} \boldsymbol{\cdot} \hat{\tilde{\boldsymbol{n}}}_a^{{lg}} \geq 0\\ \hat{\tilde{\boldsymbol{n}}}_a^{{lg}} & \text{else} \end{cases}. \end{equation}

Here, $\varTheta _{\infty }^{{s}}$ is a calibration parameter used to obtain corrected normal vectors of neighbouring SPH particles $b\in \varOmega ^{s}$. By setting $\varTheta _{\infty }^{{s}} > \varTheta _{\infty }$, the curvatures at the position of SPH particles $a \in \varOmega ^{l}$ are shifted to larger values, which prevents the continuous spreading of the drop on the solid substrate, which otherwise happens if $\varTheta _{\infty }$ 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 the table in § 8.4. The condition in (6.11) measures the distance of a neighbouring particle $b\in \varOmega ^{s}$ to particle $a$ in the normal direction. The sign of the distance allows us to distinguish between neighbouring particles $b\in \varOmega ^{s}$ located on the liquid or on the gas side of particle $a$, as illustrated in figure 2. Neighbouring SPH particles that are located on the gas side of particle $a$, $\boldsymbol {r}_{ab} \boldsymbol{\cdot} \hat {\tilde {\boldsymbol {n}}}_a^{{lg}}>0$, do not contribute to the curvature estimate obtained in (6.10).

Figure 2. Identification of neighbouring particles $b\in \varOmega ^{s}$ which contribute to the curvature computation of a particle $a\in \varOmega ^{l}$. This schematic shows SPH particles near the three-phase contact of a droplet resting on a solid boundary. Grey and blue spheres represent the solid boundary and the liquid phase, respectively. The normal vector $\hat {\tilde {\boldsymbol {n}}}_a^{{lg}}$ is used to compute the distance of a neighbouring particle $b\in \varOmega ^{s}$ to the tangent plane (here shown as a dashed black line) by $\boldsymbol {r}_{ab}\boldsymbol{\cdot}\hat {\tilde {\boldsymbol {n}}}_a^{{lg}}$. Particles $b\in \varOmega ^{s}$ which satisfy $\boldsymbol {r}_{ab}\boldsymbol{\cdot}\hat {\tilde {\boldsymbol {n}}}_a^{{lg}} \geq 0$ contribute to the curvature computation and are shown by black spheres.

Figure 3 shows cuts through the axis of symmetry of two drops resting on a plane surface. For a desired equilibrium contact angle of $\varTheta _{\infty } = 30^{\circ }$ (figure 3a) and $\varTheta _{\infty } = 60^{\circ }$ (figure 3b), the proposed approach modifies the normal vectors of the red-coloured SPH particles during the curvature computation. Note that each SPH particle $a \in \varOmega ^{l}$ in the vicinity of the three-phase contact line has its own set of SPH particle neighbours $b \in \varOmega ^{s}$.

Figure 3. Identified solid SPH particles (red) representing an extension of the liquid–gas interface: (a) $\varTheta _{\infty } = 30^{\circ }$; (b) $\varTheta _{\infty } = 60^{\circ }$.

The acceleration due to surface tension and wetting experienced by a particle $a\in \varOmega ^{l}$ is now given by

(6.12)\begin{equation} \boldsymbol{f}_{a}^{{s}} = \left(1+\frac{1}{S^{{n}}_a}\right) \frac{\sigma_{a}}{\rho_a} {\kappa}_{a} \hat{\tilde{\boldsymbol{n}}}_a^{{lg}} \delta_a^{{s}}. \end{equation}

Here, $S^{{n}}_a$ is the Shepard filter computed by

(6.13)\begin{equation} S^{{n}}_a = \sum_{b\in\varOmega^{{n}}_{l}} \frac{m_b}{\rho_b} W_{ab} + \sum_{b\in\varOmega^{{n}}_{s}} \frac{m_b}{\rho_b} W_{ab}, \end{equation}

where $\varOmega ^{{n}}_{l} \subset \varOmega ^{l}$ is the subset of SPH particles satisfying $\left \lVert \tilde {\boldsymbol {n}}_a\right \rVert > \varepsilon ^{n}$ and $\varOmega ^{{n}}_{s} \subset \varOmega ^{s}$ is the subset of neighbouring SPH particles satisfying $\boldsymbol {r}_{ab} \boldsymbol{\cdot} \hat {\tilde {\boldsymbol {n}}}_a^{{lg}} \geq 0$ in (6.11).

Using (6.12) and (6.13), the acceleration of an SPH particle $a$ due to surface tension and wetting is now given by

(6.14)\begin{equation} \boldsymbol{f}_{a}^{{s}} = \left(1+\frac{1}{S^{{n}}_a}\right) \frac{\sigma_{a}}{\rho_a} {\kappa}_{a} \hat{\tilde{\boldsymbol{n}}}_a^{{lg}} \delta_a^{{s}}. \end{equation}

6.2. Non-wetting contact angle, $\varTheta _{\infty } > 90$

Instead of computing the curvature using the approach presented in § 6.1, we model non-wetting contact angles by computing the curvature of an SPH particle $a \in \varOmega ^{l}$ using only SPH particle neighbours $b \in \varOmega ^{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 § 6.1, the normal vector assigned to an SPH particle $a\in \varOmega ^{l}$ is

(6.15)\begin{equation} \hat{\boldsymbol{n}}^{{l},\infty}_{a} = \hat{\boldsymbol{t}}_a^{{sf},{l}} \sin \varTheta_{\infty} - \hat{\tilde{\boldsymbol{n}}}_a^{{sf}} \cos \varTheta_{\infty}. \end{equation}

Here, $\boldsymbol {t}^{{sf},{l}}_{a}$ is the tangent vector at the position of a liquid SPH particle given by

(6.16a,b)\begin{equation} \boldsymbol{t}^{{sf},{l}}_{a} = \hat{\tilde{\boldsymbol{n}}}^{{lg},{l}}_a - \left(\hat{\tilde{\boldsymbol{n}}}_a ^{{lg},{l}}\boldsymbol{\cdot} \hat{\tilde{\boldsymbol{n}}}^{{sf}}_a \right)\hat{\tilde{\boldsymbol{n}}}^{{sf}}_a,\quad \hat{\boldsymbol{t}}^{{sf},{l}}_{a} = \frac{\boldsymbol{t}_{a}^{{sf},{l}}}{\lVert\boldsymbol{t}_{a}^{{sf},{l}}\rVert}, \end{equation}

where the normal vector, $\hat {\tilde {\boldsymbol {n}}}_a^{{l}}$, is computed from liquid SPH particles by

(6.17a–c)\begin{align} \hat{\tilde{\boldsymbol{n}}}_a^{{l}} = \begin{cases} 0 & \text{if}\ \left\lVert\tilde{\boldsymbol{n}}_a^{{l}}\right\rVert< \varepsilon^{n}\\ \dfrac{\tilde{\boldsymbol{n}}_a^{{l}} }{ \lVert \tilde{\boldsymbol{n}}^{{l}}_a\rVert} & \text{otherwise.} \end{cases}, \quad \tilde{\boldsymbol{n}}_a^{{l}} = \frac{1}{S_a^{{l}}} \sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b} \boldsymbol{n}_a^{{l}} W_{ab}, \quad \boldsymbol{n}_a^{{l}} =-\sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b} \boldsymbol{\nabla} W_{ab}. \end{align}

Here, $S_a^{{l}}$ is the Shepard filter applied to liquid SPH particles

(6.18)\begin{equation} S_a^{{l}} = \sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b} W_{ab}. \end{equation}

Finally, the normal vector of an SPH particle $a\in \varOmega ^{l}$ located in the vicinity of the three-phase contact line region is given by

(6.19)\begin{equation} \hat{\tilde{\boldsymbol{n}}}_a^{{slg},{l}} = \frac{f_a \hat{\tilde{\boldsymbol{n}}}_a^{{l}} + \left(1 - f_a\right)\hat{\boldsymbol{n}}^{{l},\infty}_{a}} {\left\lVert f_a \hat{\tilde{\boldsymbol{n}}}_a^{{l}} + \left(1 - f_a\right)\hat{\boldsymbol{n}}^{{l}, \infty}_{a}\right\rVert}. \end{equation}

The modified normal vectors, $\hat {\tilde {\boldsymbol {n}}}_a^{{slg},{l}}$, replace the normal vectors computed in (6.16a,b) if the SPH particle is located in the vicinity of the three-phase contact line

(6.20)\begin{equation} \hat{\tilde{\boldsymbol{n}}}^{{lg}}_a = \begin{cases} \hat{\tilde{\boldsymbol{n}}}^{{slg},{l}}_a & \text{if}\ a \in \varOmega^{{slg}}\\ \hat{\tilde{\boldsymbol{n}}}^{{lg},{l}}_a & \text{else} \end{cases}. \end{equation}

Finally, the mean curvature of a liquid SPH particle is given by

(6.21)\begin{equation} {\kappa}_a =- \frac12\sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b}\left(\hat{\tilde{\boldsymbol{n}}}_a^{{lg}} - \hat{\tilde{\boldsymbol{n}}}_b^{{lg}} \right)\hat{\boldsymbol{\nabla}}^{{l}} W_{ab}, \end{equation}

where the renormalized kernel gradient, $\hat {\boldsymbol {\nabla }}^{{l}}W_{ab}$, is computed from liquid SPH particles $b\in \varOmega ^{l}$ by (Bonet & Lok 1999; Oger et al. 2007)

(6.22a,b)\begin{equation} \hat{\boldsymbol{\nabla}}^{{l}} W_{ab} \equiv \boldsymbol{C}^{{l}}_a \boldsymbol{\nabla} W_{ab},\quad \boldsymbol{C}^{{l}}_a = \bigg(\sum_{b \in \varOmega^{{l}}} \frac{m_b}{\rho_b}\boldsymbol{ r}_{ab} \otimes \boldsymbol{\nabla} W_{ab} \bigg)^{-1}. \end{equation}

Here, $\boldsymbol {C}^{{l}}_a$ is the correction matrix computed for SPH particles $b \in \varOmega ^{l}$. The acceleration of an SPH particle $a\in \varOmega ^{l}$ due to surface tension and wetting forces is then given by

(6.23)\begin{equation} \boldsymbol{f}_{a}^{{s}} = \left(1+\frac{1}{S^{{n}}_a}\right) \frac{\sigma_{a}}{\rho_a} {\kappa}_{a} \hat{\tilde{\boldsymbol{n}}}^{{lg}}_a \delta^{{s},{l}}_a, \end{equation}

with

(6.24a,b)\begin{equation} S^{{n}}_a = \sum_{b \in \varOmega^{n}_{l}} \frac{m_b}{\rho_b} W_{ab},\quad \delta^{{s},{l}}_a= \left\lVert\boldsymbol{n}_a^{l}\right\rVert. \end{equation}

For desired equilibrium contact angles approaching $90^\circ$, a slight discrepancy between the wetting and the non-wetting approach is observed, with the wetting approach yielding greater accuracy.

6.3. Extreme wetting and thin films

Obtaining stable drops on a plane surface with $\varTheta _{\infty } \leq 30^{\circ }$ is challenging because the particles near the three-phase contact line have many fewer neighbours within their support radius ($b\in \varOmega ^{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 (6.20) and (6.11)). 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 (5.6) and (6.16a,b) will be applied to all particles $a \in \varOmega ^{l}$ in the liquid layer, resulting in unphysical surface forces.

For the case of small contact angle, as simulated in § 8.5, we substitute the curvature computation in (6.10) by (6.21) (where a renormalized kernel gradient is used) for an SPH particle $a\in \varOmega ^{slg}$ if there is an insufficient neighbour 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 $\hat {\tilde {\boldsymbol {n}}}_a^{{lg}} \boldsymbol{\cdot}\hat {\tilde {\boldsymbol {n}}}_a^{{sf}} \leq 0.995$, for the normal correction (6.6) of particles near the triple contact line to avoid the particles in a thin film being classified as located near the triple line.

7. 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 (Violeau 2015). Several techniques were proposed to ensure no penetration. For instance, wall boundaries can be modelled using repulsive forces (Monaghan 1994), dummy particles (Crespo, Gómez-Gesteira & Dalrymple 2007), mirror particles (Morris et al. 1997), immersed boundary methods (Nasar et al. 2019) or semianalytical approaches (Leroy et al. 2014). To simultaneously satisfy no-penetration and no-slip boundary conditions at walls, we modify the PPE in (3.7) following Adami, Hu & Adams (2012). We define the set of wall particles $\varOmega ^{w}$ as a subset of solid particles, $\varOmega ^{w}\subset \varOmega ^{s}$. Wall particles, $a\in \varOmega ^{w}$, are thus near fluid particles. In mathematical terms, an SPH particle represents the solid wall boundary, $a\in \varOmega ^{w}$, if

(7.1a,b)\begin{equation} a\in\varOmega^{s}\quad \text{and}\quad S^{{l}}_a > 10^{-3}. \end{equation}

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\in \varOmega ^{w}$) is

(7.2)\begin{equation} \sum_{b \in \varOmega^{{l}}} \frac{m_{b}}{\rho_{b}} \frac{4}{\rho_{a} + \rho_{b}} \left(\,p_{a}^{{\ast}} - p_{b}^{{\ast}}\right) F_{ab} = \sum_{b\in \varOmega^{{l}}}-\frac{m_{b}}{\rho_{b}} \frac{\boldsymbol{u}_{ab}^\ast \boldsymbol{\cdot}\boldsymbol{\nabla} W_{ab}}{\Delta t}. \end{equation}

Thus, the neighbourhood 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,

(7.3) \begin{equation} \begin{cases} (3.9) & \text{if}\ a \notin \varOmega^{w} \text{ and } S_a > 0.95\\ (4.4) & \text{if}\ a \notin \varOmega^{w} \text{ and } S_a \leq 0.95\\ (7.2) & \text{if}\ a \in \varOmega^{w} \end{cases}. \end{equation}

To enforce the no-slip boundary condition at the wall, the wall's velocity is calculated via an SPH interpolation of the liquid neighbourhood of $a\in \omega ^{w}$ as

(7.4)\begin{equation} \tilde{\boldsymbol{u}}_{a} = \frac{1}{S_a^{{l}}}\sum_{b \in\varOmega^{{l}}} \frac{m_b}{\rho_b}\boldsymbol{u}_{b} W_{ab}. \end{equation}

Subsequently, the wall particle velocity, ${\boldsymbol {u}}^{{w}}_a$, that replaces the velocity of neighbouring particles $b\in \varOmega ^{w}$ in (2.11a,b), is given by

(7.5)\begin{equation} {\boldsymbol{u}}^{{w}}_a = 2 \boldsymbol{u}_{a} - \tilde{\boldsymbol{u}}_{a}. \end{equation}

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 (Morris 2000) from the maximum acceleration $\boldsymbol {a}_{max}$, the maximum velocity $\boldsymbol {u}_{max}$, the viscous diffusion and the surface tension,

(8.1)\begin{equation} \Delta t = \frac14\mathrm{min}\left(\sqrt{\frac{h}{\left\lVert\boldsymbol{a}_{max}\right\rVert}}, \frac{h}{\left\lVert\boldsymbol{u}_{max}\right\rVert},\ \frac{h^2 \rho}{2\eta},\ \sqrt{\frac{h^3 \rho}{2{\rm \pi} \sigma}} \right). \end{equation}

The solid SPH particles (including the wall particles) are kept at zero velocity throughout the simulations.

Table 1. Numerical parameters used in the test simulations.

8.2. Test case: plane Poiseuille flow

We simulate plane Poiseuille flow between two stationary, infinite planes, located at ${x = \pm L}$ with ${L=0.5~{\rm 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 (Morris et al. 1997)

(8.2)\begin{align} u_y(x,t) = \frac{g_y}{2\nu} x\left(x-L\right) + \sum_{n=0}^{\infty} \frac{4g_y L^2}{\nu {\rm \pi}^3\left(2n+1\right)^3} \sin\left(\frac{{\rm \pi} x}{L}\left(2n +1\right)\right)\exp\left(-\frac{\left(2n+1\right)^2 {\rm \pi}^2 \nu}{L^2} t\right). \end{align}

In the simulation, we apply periodic boundary conditions in the $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

(8.3) \begin{equation} \begin{cases} (7.2) & \text{if}\ a \in \varOmega^{{w}} \\ (4.4) & \text{if}\ a \in \varOmega^{{fs}} \land a \notin \varOmega^{{w}}\\ (3.9) & \text{else} \end{cases}. \end{equation}

In this validation example, we disregard surface tension effects, hence the momentum balance of the fluid reads

(8.4)\begin{equation} \frac{\mathrm{D} \boldsymbol{u}_a}{\mathrm{D} t} = \boldsymbol{f}^{{p}}_a + \boldsymbol{f}^{{v}}_a + \boldsymbol{f}^{{b}}_a. \end{equation}

The flow is characterized by the Reynolds number, ${Re = 2 L u_{y}^{max}/\nu = 125}$.

The SPH particles representing the liquid phase are initially placed on a rectangular lattice in the three-dimensional interval $\boldsymbol {x} \in (-L,\,L)^3$. At ${x = \pm L}$, the boundary is modelled by SPH wall particles, placed at ${L\leq x \leq L + 2 r_{max}}$ and ${L - 2r_{max} \leq x\leq 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 (4.4) to solid SPH particles, in order to solve the PPE. The liquid phase is represented by ${L/\Delta x\in \{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, (8.2), of the transient velocity field with the numerical result for ${L/\Delta 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 $-\Delta x \leq y \leq \Delta x$ and $-\Delta x\leq z \leq \Delta x$, as a function of their $x$-position. For $t\nu /L^2 = 4\,$ for $L/\Delta x = 40$ we obtain the mean relative error $1.6\,\%$. Reducing the spatial discretization to $L/\Delta x = 10$ we obtain the mean relative error $3.8\,\%$.

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.

8.3. 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}\ {\rm m}$, using the discretization given in table 2. All smoothed normal vectors of magnitude smaller than ${\varepsilon ^{{n}} = {0.3}/{h}}$ are discarded from the computation of the curvature.

Table 2. Number of SPH particles and discretization spacing used to validate the Young–Laplace pressure boundary condition.

According to the Young–Laplace equation, the pressure drop across a liquid–gas interface at equilibrium is given by (Mugele & Heikenfeld 2019)

(8.5)\begin{equation} \Delta p^{{YL}} = 2\sigma\kappa. \end{equation}

For the parameters given above, we find that the pressure inside the drop, $p_{i}$, exceeds the ambient pressure by $20\ {\rm Pa}$ according to

(8.6)\begin{equation} p^{i} = p^{amb} + \frac{2\sigma }{r_0}. \end{equation}

Figure 5 shows the equilibrium pressure along a cut through the symmetry axis of the drop after a sufficient relaxation time $(t = 1\ \mathrm {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 (8.6). 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.

Figure 5. Comparison of the computed equilibrium pressure inside a droplet with the analytical solution of the Young–Laplace equation, (8.5). The simulation result is averaged over the time interval $0.8 \leq t \leq 1 \ {\rm s}$. (a) Pressure profile along the $x$-axis ($y/r_0 = z/r_0 =0$) for a drop centred at $(x,y,z)=(0,0,0)$ for different spatial resolutions. The dashed line shows the analytical solution, (8.5). (b) Relative deviation of the simulated equilibrium pressure inside the droplet from the analytical value, (8.5).

8.4. 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 (Aalilija, Gandin & Hachem 2020)

(8.7)\begin{equation} r\left(\theta, t\right) = r_0 \left(1+\varepsilon_2(t) P_2\left( \cos\left(\theta\right)\right)-\tfrac{1}{5}\varepsilon^{2}_2(t) \right). \end{equation}

Here, $\theta$ is the polar angle, $r_0$ is the radius of the (spherical) droplet in equilibrium and $P_2(x)\equiv \frac 12(3x^2-1)$ is the second-order Legendre polynomial. The governing parameter

(8.8)\begin{equation} \varepsilon_2(t) \approx 0.08 \mathrm{e}^{-\lambda_2 t} \cos\left(\omega_{2,0} t\right) \end{equation}

contains the oscillation frequency and the damping constant,

(8.9a,b)\begin{equation} \omega_{2,0} \equiv \sqrt{\frac{8\sigma}{\rho r_0^{3}}},\qquad \lambda_2 \equiv \frac{5\eta}{\rho r_0^2}. \end{equation}

For the initial condition, we assume (Aalilija et al. 2020)

(8.10)\begin{equation} r\left(\theta, t = 0\right) = r_0 \left(1+0.08 P_2\left(\cos\left(\theta\right)\right)-\tfrac{1}{5}0.08^2\right), \end{equation}

which describes a weakly deformed sphere. For the simulation, we assume ${\eta = 5\times 10^{-3}\ {\rm Pa}\ {\rm s}}$ and the values given in table 1. For initialization, the shape described in (8.10) is represented by 33 240 SPH particles placed on a square lattice, see figure 6.

Figure 6. Ellipsoid according to (8.10), represented by 33 240 SPH particles.

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,

(8.11)\begin{equation} r = A_{{osc}}\exp\left(-\lambda_{{osc}}\right) \cos\left(\omega_{{osc}} t\right), \end{equation}

to the simulation data. We take the amplitude $A_{{osc}}$ from the initialization and obtain the best fit for the damping constant $\lambda _{{osc}} = 22.44\ {\rm s}^{-1}$ and the oscillation frequency $\omega _{{osc}} = 261.5\ {\rm s}^{-1}$. We compare these numerical values with the analytical quantities. For the given material and system parameters, according to (8.9a,b), they read $\omega _{2,0} = 283,03\ {\rm s}^{-1}$ and $\lambda _2 = 25\ {\rm s}^{-1}$, resulting in the relative deviation $|\omega _{2,0} - \omega _{{osc}}|/ \omega _{2,0}\approx 7\,\%$ and $|\lambda _2-\lambda _{osc}|/\lambda _2\approx 5.1\,\%$.

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, $T_{{osc}} = 2{\rm \pi} /\omega _{2,0} = 2.22\times 10^{-2}\ {\rm s}$.

An appropriate value of the normal threshold, $\varepsilon ^n$, was found to be crucial for stable oscillations, see (6.16a,b). Choosing $\varepsilon ^{n}$ too small results in incorrect curvature computation due to contributions from SPH particles located far away from the interface. Choosing $\varepsilon ^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 $\varepsilon ^n$ compromises between these limits.

8.5. 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, $\varTheta _{\infty }$, which is a parameter of the simulation (see (6.1)), we describe wetting or dewetting 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 $\varTheta _{\infty }$.

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 $\Delta x$ inside a hemisphere of radius $r_0 = 10^{-3}\ {\rm m}$. The flat side of the hemisphere is in contact with the plane, modelled as a circular disk of radius $3\times 10^{-3}\ {\rm m}$ and thickness $5\Delta 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 $\Delta 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 table 1, the smoothing length $h=2.5\Delta x$ and $\varepsilon ^{{n}}=0.2/h$. The equilibrium contact angles used to modify the normal vectors are shown in table 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 $\varTheta _{\infty }\leq 95^{\circ }$ we employ $\varTheta _{\infty }^{{s}} \geq \varTheta _{\infty }$, see § 6.1.

Figure 8 shows the drop in equilibrium for contact angles $\varTheta _{\infty } \in \{15^{\circ }, 60^{\circ }, 120^{\circ }, 150^{\circ }\}$. During relaxation, the kinetic energy of the drop and, thus, its constituting SPH particles decays by several orders of magnitude, see figure 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\ {\rm s}$ we find for wetting contact $E_{kin}\approx 1.16\times 10^{-14}\ {\rm J}$ for $\varTheta _{\infty } = 15^{\circ }$ and $E_{kin}\approx 6.94\times 10^{-15}\ {\rm J}$ for $\varTheta _{\infty } = 60^{\circ }$; while for non-wetting contact, we find $E_{kin}\approx 7.43\times 10^{-12}\ {\rm J}$ for $\varTheta _{\infty } = 120^{\circ }$ and $E_{kin}\approx 1.38\times 10^{-12}\ {\rm J}$ for ${\varTheta _{\infty } = 150^{\circ }}$. Similarly, the fluctuations of the kinetic energy are much larger for non-wetting contact than for wetting contact. This behaviour 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 (Lafuma & Quéré 2003). In terms of the numerical model, this effect can be understood from the influence of the normal threshold, $\varepsilon ^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 $\varepsilon ^n$ and, therefore, do 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 8. At large time, $t = 1\ {\rm s}$, the drops have assumed their equilibrium shape. Liquid SPH particles are shown in red, and solid particles are shown in blue. For $\varTheta _{\infty } \in \{120^{\circ }, 150^{\circ }\}$, some particles near the three-phase contact line disintegrated from the body of the liquid phase: (a) $\varTheta _{\infty } = 15^{\circ }$; (b) $\varTheta _{\infty } = 60^{\circ }$; (c) $\varTheta _{\infty } = 120^{\circ }$; (d) $\varTheta _{\infty } = 150^{\circ }$. Refer to figure 10 for the radial profile of the free surface.

Figure 9. 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 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 ($\varTheta _\infty$) of the drop (input parameter), and the momentary value

(8.12)\begin{equation} \varTheta = \frac{\rm \pi}{2} + \arctan\left(\frac{H-r}{B}\right), \end{equation}

where $H$ and $B$ are the drop's height and base radius, and

(8.13)\begin{equation} r = \frac{H^2 + B^2}{2 H} \end{equation}

is the radius of the drop.

Figure 10. A cut-through of a drop at $y=0$ at large time, $t=1\ {\rm 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 $\varTheta$ of the contact angle defined in (8.12) and the equilibrium contact angle, ${\varTheta _\infty}$, in brackets. (a) Wetting contact angles; (b) non-wetting contact angles.

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,

(8.14)$$\begin{gather} B_{{anl}} = \sqrt{2r_{{anl}} H_{{anl}} -H_{{anl}} ^2} \end{gather}$$

(8.15)$$\begin{gather}H_{{anl}} = r_{{anl}} \left(1-\cos\left(\varTheta_{\infty}\right)\right), \end{gather}$$

with the analytical value of the drop's radius,

(8.16)\begin{equation} r_{{anl}} = \left(\frac{2r_0^3}{2-3\cos(\varTheta_{\infty}) + \cos^3\left(\varTheta_{\infty}\right)}\right)^{1/3}. \end{equation}

Figure 11 shows this comparison for various contact angle values, $\varTheta _\infty$.

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. (a) Analytical (solid lines) and simulated (circles) values of the drop base diameter ($2B$) and height ($H$). (b) Simulated contact angle, $\varTheta$, as a function of the equilibrium contact angle, $\varTheta _\infty$.

In the range ${15^{\circ } \leq \varTheta _{\infty } \leq 135^{\circ }}$, the deviation of the simulated contact angle from the analytical value is below $4.5\,\%$, see figure 11(b), and increases to $14.3\,\%$ for ${\varTheta _{\infty } \leq 150^{\circ }}$. 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 $\varepsilon ^{{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 $\varepsilon ^{{n}}$.

8.6. 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, ${\boldsymbol {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 table 1 and the set-up geometry introduced in § 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}\ {\rm m}$, resting on a plane solid substrate of radius $4\times 10^{-3}\ {\rm m}$. The equilibrium contact angle is $\varTheta _{\infty } = 50^{\circ }$ ($\varTheta _{\infty }^{{s}} = 55^{\circ }$).

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,

(8.17)\begin{equation} \textit{Bo} \equiv \frac{\rho g r_0^2}{\sigma}. \end{equation}

The drop height, $H$, is scaled by the drop height in the absence of gravity $H_0$ where $\textit {Bo} =0$. In figure 12, the Bond number covers the interval $\textit {Bo}\in [10^{-3},10^2]$, and the function $H(\textit {Bo})$ decays monotonically. For small $\textit {Bo}$, the situation is surface-tension dominated, thus, $H\to H_0$ for $\textit {Bo}\to 0$.

Figure 12. Drop height as a function of the Bond number. For small gravity, $\textit {Bo}\to 0$, the regime is dominated by surface tension, and $H$ approaches $H_0$. For large gravity, $\textit {Bo}\to \infty$, the regime is gravity-dominated, and $H$ approaches $H_\infty$ given by (8.18). The symbols show the drop height obtained from SPH simulations, and the solid lines show the limits $H_0$ and $H_\infty$.

For large gravity, $\textit {Bo}\to \infty$, the regime is gravity-dominated, and the height of the drop converges to the capillary length (Dupont & Legendre 2010)

(8.18)\begin{equation} H_\infty = 2 \sqrt{\frac{\sigma}{\rho g_z}} \sin{\left(\frac{\varTheta_{\infty}}{2}\right)}. \end{equation}

Figure 13 shows the drop profile for several values of $\textit {Bo}$, that is, the position of the free surface at $y=0$. For real-life applications, a dynamic contact angle model (such as Göhl et al. 2018) may be implemented for taking into account chemical heterogeneity of the substrate.

Figure 13. A cut-through of the drop at $y=0$ (liquid–gas interface position) as a function of the Bond number for $\varTheta _{\infty }= 50^{\circ }$.

8.7. Test case: droplet pinning

The three-dimensional test cases presented so far are axisymmetric. While this symmetry makes the analytical solution considerably easier in many cases, it does not simplify the numerical solution using the method presented here. At no point in the numerical procedure is the symmetry of the solution assumed. In this respect, the significance of the test cases is not reduced by their axial symmetry. On the contrary, the fact that the numerical solution exhibits the required symmetry is a further proof of the stability of our method.

The problem in the current section is not axisymmetric.

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 $\varTheta _\infty + \varPsi$ (Breinlinger et al. 2013), where $\varPsi$ 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 figure 14.

Figure 14. Sketch of a pinned drop at the contact line between two planes of different inclinations. The applied body acceleration $\boldsymbol {f}^{{b}}$ drives the droplet towards the contact line between the planes. It can pass the barrier only if the contact angle exceeds the threshold $\varTheta = \varTheta _{\infty } + \varPsi$.

We describe the numerical experiment similar to §§ 8.5 and 8.6 with the parameters specified in table 1, viscosity $\eta = 5\times 10^{-3}\ {\rm Pa}\ {\rm s}$ and the equilibrium contact angle ${\varTheta _{\infty } = 75^{\circ }}$ ($\varTheta _{\infty }^{{s}} = 80^{\circ }$).

Similar to the preceding test cases, we fill a hemisphere with radius $r_0=10^{-3}\ {\rm m}$ with 16 776 liquid SPH particles arranged on a square lattice. This half-sphere rests on a solid horizontal plane modelled as a block of size $(2.6\times 3 \times 0.25)\ {\rm mm}^{3}$ filled by SPH particles arranged in a square lattice. The initial condition is sketched in the left-hand part of figure 14. The second plane is modelled in the same way but inclined by $\varPsi \in \{22.5^{\circ }, 45^{\circ }, 67.5^{\circ }\}$, see the right-hand part of figure 14.

To prepare our initial conditions, we simulate relaxing the drop on the horizontal substrate to equilibrium. After $t = 50\ {\rm ms}$ of relaxation, the equilibrium shape was assumed. For $t > 50\ {\rm ms}$ a constant body acceleration, $\boldsymbol {f}^{{b}} = [7.5, 0, -7.5]$, drives the drop towards the contact line of the planes. Figure 15 shows snapshots of the simulation.

Figure 15. Simulation snapshots of droplets that are accelerated towards the contact line of two planes with different inclinations. Panel (a i–iv) shows the droplets at time $t = 50\ {\rm ms}$ resting in equilibrium on a horizontal solid plane with the equilibrium contact angle $\varTheta _{\infty } = 50^{\circ }$. Panel (b i–iv) shows the droplets at time $t = 100\ {\rm ms}$ when approaching the contact line under the action of acceleration $\boldsymbol {f}^{{b}}$. Panel (c i–iv) at time $t = 150\ {\rm ms}$: For $\varPsi \in \{22.5^\circ, 45^\circ, 67.5^\circ \}$ the contact angle exceeds the threshold, thus, the drop passed over the surface discontinuity. For $\varPsi = 90^{\circ }$, the contact angle is below the threshold, thus, the drop remains pinned.

Figure 15(a i–iv) shows the equilibrated drops at $t = 50\ {\rm ms}$. Figure 15(b i–iv) shows the drops at $t=100\ {\rm ms}$ under the influence of the acceleration $\boldsymbol {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 smooth counteracting force, enabling some SPH particles to pass the contact line. Figure 15(c i–iv) shows the droplets at $t = 150\ {\rm ms}$. For $\varPsi \in \{22.5^{\circ }, 45^{\circ }, 67.5^{\circ }\}$ the drop passes the discontinuity. For $\varPsi = 90$, the external acceleration is insufficient to generate the threshold contact angle $\varTheta _{\infty } + \varPsi = 165^{\circ }$ (based on the contact-line energy balance (Breinlinger et al. 2013; Wang, Lin & Zhao 2019)). A maximum contact angle of approximately $118^{\circ }$ is achieved at the pinned state.

Figure 16 shows the total kinetic energy of the liquid SPH particles as a function of time. During the relaxation, $t\in [0,50\ {\rm ms}]$, the drop finds its equilibrium shape, therefore, energy decays. For $t>50\ {\rm ms}$, the decline is accelerated by $\boldsymbol {f}^{{b}}$ acting on the liquid SPH particles leading to a sharp increase in the total kinetic energy. At $t\approx 55\ {\rm ms}$, the drop approaches the contact line between the two planes, indicated by a decrease of the kinetic energy for $t > 55\ {\rm ms}$. As expected, the smallest deceleration of the drop is experienced for $\varPsi =22.5^{\circ }$ and the most significant for $\varPsi = 90^{\circ }$. The subsequent increase in kinetic energy for $t \gtrsim 85\ {\rm ms}$, for $\psi <90^{\circ }$ corresponds to the drop's passing of the discontinuity. For $\varPsi = 90^{\circ }$, the kinetic energy remains approximately three orders of magnitude below, indicating the pinning of the drop at the contact line.

Figure 16. Kinetic energies of the drops over time for different inclinations, $\varPsi$, of the plane. The time ${t< 50\ {\rm ms}}$ corresponds to the relaxation of the drop to find its equilibrium shape. For $\varPsi \in \{22.5^{\circ }, 45^{\circ }, 67.5^{\circ }\}$ the drop passes the discontinuity. For $\varPsi = 90$ it remains pinned.

The kinetic energy curves in figure 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 towards 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 neighbours of the liquid phase. Otherwise, the particles should remain unaffected by the surface tension.

9. Conclusion

We introduce a surface tension model based on the CSF approach and applicable to three-dimensional free-surface flows. In contact with a substrate, the model can handle wetting and non-wetting behaviour in a wide range of parameters. This is made possible with the use of primarily three empirical parameters, namely the Shepard sum correction for free surface introduced in (5.12), the threshold for magnitude of surface normal vectors, (5.5), $\epsilon ^{n}$ and the desired equilibrium contact angle associated with the substrate $\varTheta _\infty ^s$ listed in table 3. We have shown that the model performs numerically stable in several test cases. By directly comparing analytical results and the literature results, we were able to show that the model delivers quantitatively reliable results.

To demonstrate the model's performance, we applied the new simulation method to (i) plane Poiseuille flow and a set of free-surface problems, including (ii) the Laplace jump of a droplet, (iii) the oscillation and relaxation of a perturbed droplet, (iv) the equilibrium shape of a droplet in contact with a wetting/non-wetting substrate and its relaxation to this equilibrium, (v) the flattening of a droplet under the action of gravity and (vi) 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.