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. (Reference Brackbill, Kothe and Zemach1992) 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. (Reference Breinlinger, Polfer, Hashibon and Kraft2013), 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 Reference Bonet and Lok1999; Oger et al. Reference Oger, Doring, Alessandrini and Ferrant2007)

(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.

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 Reference Morris2000) 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. Reference Morris, Fox and Zhu1997)

(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 Reference Mugele and Heikenfeld2019)

(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 Reference Aalilija, Gandin and Hachem2020)

(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. Reference Aalilija, Gandin and Hachem2020)

(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é Reference Lafuma and 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 Reference Dupont and Legendre2010)

(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. Reference Göhl, Mark, Sasic and Edelvik2018) 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. Reference Breinlinger, Polfer, Hashibon and Kraft2013), 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. Reference Breinlinger, Polfer, Hashibon and Kraft2013; Wang, Lin & Zhao Reference Wang, Lin and Zhao2019)). 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.