2.1. Governing equations

The continuum model we use is the weakly compressible, barotropic Navier Stokes equations

(2.1) \begin{align} \frac {{\textrm{D}} \rho }{{\textrm{D}} t} &= -\rho (\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}), \end{align}

(2.2) \begin{align} \frac {{\textrm{D}} \boldsymbol{v}}{{\textrm{D}} t} &= \frac {-1}{\rho } \boldsymbol{\nabla } P + \frac {\mu }{\rho } \nabla^{2} \boldsymbol{v} + \boldsymbol{g} + \boldsymbol{F}^{(\textit{pf})}, \end{align}

with $\boldsymbol{x}$ the position, $\boldsymbol{v}$ the velocity, $\rho$ the density, $P$ the pressure, $\mu$ the dynamic viscosity and $\boldsymbol{g}$ the gravitational acceleration. The main focus of this work will be the body force (per unit mass) $\boldsymbol{F}^{(\textit{pf})}$ , which we will construct in § 2.4 to reproduce surface tension and contact angle effects. The derivative ${\textrm{D}}/{\textrm{D}}t$ is the Lagrangian derivative

(2.3) \begin{align} \frac {{\textrm{D}} }{{\textrm{D}} t} = \frac {\partial }{\partial t} + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla }, \end{align}

which, as we shall see, lends itself naturally to discretisation using SPH.

2.2. Discretisation

As is standard in SPH methods, we represent the fluid by $N$ particles, each with their own position, velocity, density and label (to distinguish fluid from solid, for example). Using $i$ and $j$ as particle indices, the discretised SPH form of the model (2.1)–(2.2) is

(2.4) \begin{align}\frac {{\textrm{d}}\rho _i}{{\textrm{d}}t} &= \rho _i \sum _{j=1}^N ( \boldsymbol{v}_{\textit{ij}} \boldsymbol{\cdot }\boldsymbol{\nabla } _i W_{\textit{ij}} ) V_j,\end{align}

(2.5) \begin{align}\frac {{\textrm{d}}\boldsymbol{v}_i}{{\textrm{d}}t} &= \frac {-1}{\rho _i} \sum _{j=1}^N \big(P_i + P_j - \varPi {}\big) \boldsymbol{\nabla } _i W_{\textit{ij}} V_j + \boldsymbol{g} + \boldsymbol{F}_i^{(\textit{pf})},\end{align}

(2.6) \begin{align}\frac {{\textrm{d}}\boldsymbol{x}_i}{{\textrm{d}}t} &= \boldsymbol{v}_i, \notag \\ \varPi &= 2(d+2)\mu \frac {\boldsymbol{v}_{\textit{ij}} \boldsymbol{\cdot }\boldsymbol{x}_{\textit{ij}}}{\| \boldsymbol{x}_{\textit{ij}} \|^2}, \notag \\ P_i(\rho _i) &= \frac {c_0^2 \rho _0}{7} \left [ \left ( \frac {\rho _i}{\rho _0} \right )^7 - 1 \right]\!, \\[9pt] \nonumber\end{align}

where $\boldsymbol{x}_{\textit{ij}} = \boldsymbol{x}_i - \boldsymbol{x}_j$ for notational convenience, $W_{\textit{ij}} = W(\boldsymbol{x}_{\textit{ij}}; H)$ is an SPH kernel with compact support radius $H$ , $V_j = m_j / \rho _j$ is the volume of particle $j$ and $d$ is the spatial dimension (always equal to 3 in this work). The pressure $P_i$ is calculated from the density $\rho _i$ with the equation of state (2.6). The constants $c_0$ and $\rho _0$ are the artificial speed of sound and the reference density, respectively. Note that we have switched from using the material derivative ${\textrm{D}}/{\textrm{D}}t$ to the total derivative ${\textrm{d}}/{\textrm{d}}t$ , as the time derivative of a particle’s property follows the flow by definition.

In this work we use the quintic Wendland function (Wendland Reference Wendland1995) as the SPH kernel, which we parameterise by $H$ , its radius of support, also called the kernel cutoff radius

(2.7) \begin{align} W(\boldsymbol{x} - \boldsymbol{x}'; H) &= \frac {21}{2 \pi H^3} \; w(\| \boldsymbol{x} - \boldsymbol{x}' \| / H), \nonumber \\[3pt] w(\tau ) &= \begin{cases} (1 - \tau )^4 (4\tau + 1), & 0 \leqslant \tau \lt 1, \\ 0, & \tau \geqslant 1. \end{cases} \end{align}

Wendland functions were developed independently of SPH for their smoothness properties, but have been found to give accurate and stable SPH results (Dehnen & Aly Reference Dehnen and Aly2012). In the literature, this kernel is sometimes parameterised by its smoothing length $h$ , which is defined as half the kernel cutoff radius. To avoid this confusion, we parameterise the kernel by its cutoff, which we denote $H$ , and set $H / \Delta x = \kappa = 4$ , where $\Delta x$ is the particle width.

With the exception of the pairwise-force term $\boldsymbol{F}^{(\textit{pf})}$ , which will be the main focus of this work, the discretisation summarised above is well established in the literature: see, for example, Monaghan (Reference Monaghan2005). In this work we will implement and validate a new form of the pairwise-force term to model surface tension and contact angle effects, and highlight some particular properties of $\boldsymbol{F}^{(\textit{pf})}$ that yield stable and accurate simulations of droplets.

2.2.1. Continuity equation

The operator we use for the divergence of velocity is from Monaghan (Reference Monaghan2005)

(2.8) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v} &= - \sum _{j=1}^N ( \boldsymbol{v}_{\textit{ij}} \boldsymbol{\cdot }\boldsymbol{\nabla } _i W_{\textit{ij}} ) V_j, \end{align}

and is specifically constructed to vanish when $\boldsymbol{v}$ is constant.

2.2.2. Pressure gradient and equation of state

The discretisation of the gradient of pressure that we use, namely

(2.9) \begin{align} (\boldsymbol{\nabla } P)_i &= \sum _{j=1}^N (P_i + P_j) \boldsymbol{\nabla }_i W_{\textit{ij}} V_j, \end{align}

is due to Bonet & Lok (Reference Bonet and Lok1999) (and recently used by Domínguez et al. (Reference Domínguez2022)), who showed it to be consistent with (2.4) and to conserve linear momentum exactly. Note that since $\boldsymbol{\nabla }_i W_{\textit{ij}} = -\boldsymbol{\nabla }_j W_{\textit{ji}}$ , we have the anti-symmetry $(\boldsymbol{\nabla } P)_i = -(\boldsymbol{\nabla } P)_j$ . This is the property that conserves momentum, since the contribution of particle $j$ to ${{\textrm{d}}\boldsymbol{v}_i}/{{\textrm{d}}t}$ is equal and opposite to the contribution of particle $i$ to ${{\textrm{d}}\boldsymbol{v}_j}/{{\textrm{d}}t}$ .

The pressure is calculated from the density by the equation of state. We follow Monaghan (Reference Monaghan1994, Reference Monaghan2005) in using

(2.10) \begin{align} P_i &= \frac {c_0^2 \rho _0}{7} \left [ \left ( \frac {\rho _i}{\rho _0} \right )^7 - 1 \right ], \end{align}

which was originally reported by Cole (Reference Cole1948, p. 44) in the study of underwater explosions. Cole notes that they chose the exponent $7$ to approximately fit experimental data. In our case, we find that the results are not at all sensitive to the particular equation of state in use – even the first-order approximation $P(\rho ) = c_0^2 (\rho - \rho _0)$ gives almost indistinguishable results. The coefficient $c_0$ is an artificial speed of sound that controls the compressibility of the fluid. The actual speed of sound in water is around $1500\, {\textrm{m s}}^{-1}$ , but such a value would require the use of extremely small time steps to properly resolve pressure waves travelling at that speed. Instead, in the weakly compressible model, we use an artificial value of $c_0$ of the order of $100\;{\textrm{m s}}^{-1}$ such that density fluctuations are kept to within $1\, \%$ of the reference value $\rho _0$ (Monaghan Reference Monaghan2005). The artificial speed of sound required for a particular problem can be estimated as $10 v_{\textit{max}}$ , an order of magnitude greater than the maximum expected fluid speed.

2.2.3. Viscosity

The discrete operator we use to approximate the Laplacian is that of Monaghan & Gingold (Reference Monaghan and Gingold1983), namely

(2.11) \begin{align} (\nabla^{2} \boldsymbol{v})_i &= 2(d+2) \sum _{j=1}^N \frac {\boldsymbol{v}_{\textit{ij}} \boldsymbol{\cdot }\boldsymbol{x}_{\textit{ij}}}{\| \boldsymbol{x}_{\textit{ij}} \|^2} \boldsymbol{\nabla } _i W_{\textit{ij}} V_j, \end{align}

where $d$ is the spatial dimension. More recently proposed operators exist with favourable properties, but we have chosen this particular discretisation based on the analysis by Colagrossi, Antuono & Le Touzé (Reference Colagrossi, Antuono and Le Touzé2009) that suggests it is more appropriate for the simulation of free surfaces.

2.2.4. Solid boundary treatment

At the fluid–solid interface, we have the no-slip condition: $\boldsymbol{v} = \boldsymbol{0}$ . We impose this condition indirectly, representing the solid substrate with fixed dummy particles, initialised on a regular grid with spacing $\Delta x$ . Multiple layers of these particles are used to avoid an SPH neighbourhood deficiency at the boundary. The solid particles have the same physical properties as the fluid at rest, with mass $\rho _0 \Delta x^3$ and density $\rho _0$ , but with zero velocity. This simple approximation of the no-slip condition has been shown to be consistent in the limit as $H\rightarrow 0$ , although does not reproduce it exactly in simulations. Other schemes utilise mirroring approaches (Macia et al. Reference Macia, Antuono, Gonzalez and Colagrossi2011), in which fluid properties are interpolated and extended into the solid boundary using surface normals. In our study of droplets on rough surfaces, however, this approach would be quite complex to implement for a general surface. For example, it is not straightforward to map any point in space to its closest surface point and thus the appropriate surface normal. Furthermore, some boundary points may be assigned incorrect properties because their associated interpolation point above the solid is in the gas phase, which we do not discretise with SPH particles. The subject of more accurate boundary condition treatment is thus left to future work, and the interested reader is referred to Macia et al. (Reference Macia, Antuono, Gonzalez and Colagrossi2011) for more details.

The ‘dummy’ solid particles are included in the summations over $j$ in the discrete continuity and momentum (2.4) and (2.5), as if they were fluid particles. Pressure forces and the repulsion due to pairwise forces ensure that fluid particles do not penetrate the boundary. Figure 1 shows a typical set-up in which a droplet of fluid particles is about to impact a flat solid boundary.

Figure 1. A side view of a three-dimensional particle simulation in which a droplet of fluid particles (blue) is about to impact solid boundary particles (grey). Insets show closeups of each type of particle. Fluid particles are advected with the flow and thus disorganised, while solid particles are fixed on a cubic lattice. The particles are rendered as spheres of volume $m_j / \rho _j$ , but we note that this is for visualisation only: in the SPH scheme, they are better understood as interpolation nodes.

2.2.5. Time integration

The time integration scheme we use is a second-order, symplectic, position-Verlet scheme. It is modified from Leimkuhler & Matthews (Reference Leimkuhler and Matthews2015) by Domínguez et al. (Reference Domínguez2022) to include a velocity half-step, which is necessary to integrate the continuity equation and calculate the viscous term. We will use the following notation for the discretised governing equations (2.1) and (2.2):

(2.12) \begin{align} \frac {{\textrm{d}}\rho _i}{{\textrm{d}}t} := Q_i, \quad \frac {{\textrm{d}}\boldsymbol{v}_i}{{\textrm{d}}t} := \boldsymbol{A}_i. \end{align}

Dropping particle indices for clarity, and letting the value of a quantity $\varphi$ at time step $k$ be denoted $\varphi ^{(k)}$ , the original position-Verlet scheme (Leimkuhler & Matthews Reference Leimkuhler and Matthews2015, p. 107) reads

(2.13) \begin{align} \boldsymbol{x}^{\left (k+\frac {1}{2}\right )} &= \boldsymbol{x}^{\left (k\right )} + \frac {\Delta t}{2} \boldsymbol{v}^{\left (k\right )}, \nonumber \\ \boldsymbol{v}^{\left (k+1\right )} &= \boldsymbol{v}^{\left (k\right )} + \Delta t \, \boldsymbol{A}^{\left (k + \frac {1}{2}\right )}, \nonumber \\ \boldsymbol{x}^{\left (k+1\right )} &= \boldsymbol{x}^{\left (k+\frac {1}{2}\right )} + \frac {\Delta t}{2} \boldsymbol{v}^{\left (k+1\right )}. \end{align}

However, the viscous term in $\boldsymbol{A}^{ (k + (1/2))}$ involves the velocity at the half time step, so Domínguez et al. (Reference Domínguez2022) introduce the intermediate step

(2.14) \begin{align} \boldsymbol{v}^{\left (k+\frac {1}{2}\right )} = \boldsymbol{v}^{\left (k\right )} + \frac {\Delta t}{2} \boldsymbol{A}^{\left (k\right )}. \end{align}

Simplifying the expression for $\boldsymbol{x}^{ (k+1 )}$ , and including the integration of the continuity equation, the full scheme reads

(2.15) \begin{align}\text{Calculate } \boldsymbol{A}^{\left (k\right )} \text{ and } Q^{\left (k\right )}, \nonumber \\ \boldsymbol{x}^{\left (k+\frac {1}{2}\right )} = \boldsymbol{x}^{\left (k\right )} + \frac {\Delta t}{2} \boldsymbol{v}^{\left (k\right )},\nonumber \\ \boldsymbol{v}^{\left (k+\frac {1}{2}\right )} = \boldsymbol{v}^{\left (k\right )} + \frac {\Delta t}{2} \boldsymbol{A}^{\left (k\right )},\nonumber \\ \rho ^{\left (k+\frac {1}{2}\right )} = \rho ^{\left (k\right )} + \frac {\Delta t}{2} Q^{\left (k\right )},\end{align}

(2.16) $$\text{Calculate } \boldsymbol{A}^{\left (k+\frac {1}{2}\right )} \text{ and } Q^{\left (k+\frac {1}{2}\right )},\\\boldsymbol{v}^{\left (k+1\right )} = \boldsymbol{v}^{\left (k\right )} + \Delta t \boldsymbol{A}^{\left (k + \frac {1}{2}\right )},\\\boldsymbol{x}^{\left (k+1\right )} = \boldsymbol{x}^{\left (k\right )} + \frac {\Delta t}{2} \left [ \boldsymbol{v}^{\left (k\right )} + \boldsymbol{v}^{\left (k+1\right )} \right ],\\\rho ^{\left (k+1\right )} = \rho ^{\left (k\right )} + \Delta t Q^{\left (k+\frac {1}{2}\right )}.$$

The time step is chosen adaptively according to viscous, maximum force and acoustic constraints

(2.17) $$\Delta t_v = \frac {\Delta x^2 \rho _0}{\mu }, \quad \Delta t_a = 0.15 \sqrt { \frac {\Delta x}{\max _j \| \boldsymbol{A}_j \|} }, \quad \Delta t_c = 0.15 \frac {\Delta x}{c_0},\\ \Delta t = \min (\Delta t_v, \Delta t_a, \Delta t_c).$$

In the present application, the acoustic constraint determining $\Delta t_c$ is almost always far smaller than the other two in (2.17), due to the small scale of the droplets and the comparatively large artificial speed of sound $c_0$ .

2.3. Implementation details

The Lagrangian nature of SPH, while making it a very flexible method, also makes it challenging to implement efficiently. We follow the work of Ihmsen et al. (Reference Ihmsen, Akinci, Becker and Teschner2011) in the use of some key data structures and parallel methods, which we will briefly summarise here.

Implemented naively, each sum in the discretised equations (2.4) and (2.5) has a time complexity of ${\mathcal{O}}(N^2)$ . This is made more efficient by pre-computing a neighbour list

(2.18) \begin{align} \mathcal{J}(i) = \left \{ j : \| \boldsymbol{x}_i - \boldsymbol{x}_j \| \lt H \right \}, \end{align}

for each fluid particle $i$ . Any sum over $j = 1,\ldots ,N$ then becomes a sum over $j \in \mathcal{J}(i)$ . Since density fluctuations are small in this weakly compressible scheme, the number of neighbours is approximately constant (around $4\pi \kappa ^3/3$ ). Thus, the complexity of calculating the particle interactions becomes ${\mathcal{O}}(N)$ .

We accelerate the neighbour search by using a background grid with cells of width $H$ . Each grid cell is uniquely identified by a tuple of integer coordinates $(c_1, c_2, c_3)$ , and we keep a hash table of lists of particles contained in each cell. To minimise memory allocations, we pre-allocate enough storage for each of these lists to contain $\kappa ^d$ indices. Density fluctuations are low, so the maximum number of fluid particles that one cell may contain is roughly constant. The neighbour search then only considers particles in neighbouring cells (the number of which is roughly constant), therefore reducing the time complexity to ${\mathcal{O}}(N)$ .

Finally, we make use of shared memory parallelism with many CPU cores wherever possible. The particle–particle interactions are calculated in parallel, as are the neighbour lists. For more optimisation details we refer the interested reader to Ihmsen et al. (Reference Ihmsen, Akinci, Becker and Teschner2011) for CPU implementations, or for GPU implementations see Domínguez et al. (Reference Domínguez, Crespo, Gómez-Gesteira and Marongiu2011), Crespo et al. (Reference Crespo, Dominguez, Barreiro, Gómez-Gesteira and Rogers2011, Reference Crespo, Domínguez, Rogers, Gómez-Gesteira, Longshaw, Canelas, Vacondio, Barreiro and García-Feal2015) and Domínguez et al. (Reference Domínguez2022). Simulations were run using up to 64 processors in parallel.

2.4. Pairwise-force model for interfacial phenomena

Surface tension and wetting are both effects of intermolecular forces (Bormashenko Reference Bormashenko2017). A molecule at an interface misses approximately half of its interactions with neighbours when compared with a molecule in the bulk, and the resulting imbalance of forces leads to free-surface energy. For surface tension, the relevant forces are cohesive (fluid–fluid interactions) while, for wetting, the relevant forces are adhesive (fluid–solid interactions). Depending on the relative strengths of these forces, a droplet could be almost spherical, spread to completely wet a solid or any configuration in between, to minimise the free-surface energy.

We aim to reproduce this behaviour on a much larger scale by mimicking intermolecular forces between SPH particles. This is conceptually similar to the Lennard–Jones potentials used in molecular dynamics simulations (Jones Reference Jones1924). The basis for these forces in SPH is empirical; nevertheless, in § 3 we will show that the pairwise-force SPH model reproduces interfacial phenomena consistently and predictably.

The pairwise particle interaction forces are included in the momentum equation (2.5) via the term $\boldsymbol{F}_i^{(\textit{pf})}$

(2.19) \begin{equation} \boldsymbol{F}_i^{(\textit{pf})} = \frac {H}{m_i} \sum _{j=1}^N s_{\textit{ij}} f_{\textit{ij}}( \| \boldsymbol{x}_{\textit{ij}} \| / H ) \frac {\boldsymbol{x}_{\textit{ij}}}{\| \boldsymbol{x}_{\textit{ij}} \|}, \end{equation}

where $s_{\textit{ij}}$ controls the strength of the pairwise force between particles $i$ and $j$ . Following the analogy with intermolecular potentials, we construct the pairwise-force profile $f_{\textit{ij}}( \| \boldsymbol{x}_{\textit{ij}} \| / H )$ to be repulsive at short distances (less than one particle width), attractive at medium distances and vanish outside the SPH kernel support radius, $H$ . In other works, this term has been constructed as a combination of Gaussians (Tartakovsky & Panchenko Reference Tartakovsky and Panchenko2016), SPH kernels (Kordilla et al. Reference Kordilla, Tartakovsky and Geyer2013; Shigorina et al. Reference Shigorina, Kordilla and Tartakovsky2017) or part of a cosine curve (Tartakovsky & Meakin Reference Tartakovsky and Meakin2005; Nair & Pöschel Reference Nair and Pöschel2018). For simplicity, we will instead use polynomials for the force profile, with some intuitive constraints at key points. Those constraints are

(2.20) \begin{align} &f_{\textit{ij}}(0) = 1, &\text{(avoid trivial solution)} \nonumber \\ &f^{\prime}_{\textit{ij}}(0) = f^{\prime}_{\textit{ij}}(1) = 0, &\text{(smoothness at endpoints)} \nonumber \\ &f_{\textit{ij}}(1) = 0, &\text{(continuity at kernel support radius)} \nonumber \\ &f_{\textit{ij}}(\kappa ^* \Delta x / H) = 0. &\text{(repulsion-to-attraction point)} \end{align}

The last constraint, controlling the zero crossing of the force profile, is intentionally different for the fluid–fluid ( ${\textrm{ff}}$ ) and the fluid–solid ( ${\textrm{fs}}$ ) interactions. We choose $\kappa ^*$ to be $1.1$ and $2$ for the fluid–fluid and the fluid–solid interactions, respectively. The fluid–fluid force profile is attractive at a shorter distance than the fluid–solid to prioritise cohesion over adhesion for stability at the contact line.

We have chosen the zero crossing of the fluid–fluid pairwise-force curve to be approximately the ‘rest distance’ of the particles. To do this, we imagine optimally packed spheres in a face-centred cubic layout. The Voronoi diagram of this arrangement is the rhombic dodecahedral honeycomb. If we equate the volume of one of our particles ( $(\Delta x)^3$ ) with the volume of a rhombic dodecahedron, we find that the distance between neighbouring particles in such a packing is approximately $1.1 \Delta x$ . This is the geometric motivation behind the location of the zero crossing of the force profile in (2.19).

The minimum degree of a polynomial satisfying our five constraints is four, and so we have

(2.21) \begin{gather} f_{\textit{ij}}(\tau ) = \begin{cases} 1 + \alpha _{\textit{ij}}\tau ^2 - 2(2 + \alpha _{\textit{ij}})\tau ^3 + (3 + \alpha _{\textit{ij}})\tau ^4, & 0 \leqslant \tau \leqslant 1, \\ 0, & \tau \gt 1, \end{cases} \end{gather}

where we use $\alpha _{\textit{ff}} = -23.5$ for fluid–fluid interactions, and $\alpha _{\textit{fs}} = -11$ for fluid–solid interactions, to satisfy the zero crossing constraints above. Figure 2 shows these force profiles over the range of possible particle distances.

Figure 2. Distance-dependent pairwise-force profiles $f_{\textit{ff}}$ and $f_{\textit{fs}}$ over the kernel support, $\tau \in [0,1]$ or $\|\boldsymbol{x}_{\textit{ij}}\| / \Delta x \in [0,\kappa ]$ . Positive values indicate repulsion and negative values indicate attraction. Note that we have plotted the force profiles with respect to the particle width $\Delta x$ rather than the dimensionless $\tau$ , to highlight the physically motivated choice of zero crossing discussed in § 4.4.

The parameter $s_{\textit{ij}}$ in (2.19) controls the strength of the pairwise force and, like $f_{\textit{ij}}(\tau )$ , depends on the labels of the particles $i$ and $j$ : fluid or solid. We denote the cohesive force strength $s_{\textit{ff}}$ (fluid–fluid) and the adhesive force strength $s_{\textit{fs}}$ (fluid–solid). Note the inclusion of the length scale $H$ in the form of $\boldsymbol{F}_i^{(\textit{pf})}$ in (2.19): this is a departure from established pairwise-force methods (Kordilla et al. Reference Kordilla, Tartakovsky and Geyer2013; Nair & Pöschel Reference Nair and Pöschel2018; Tartakovsky & Meakin Reference Tartakovsky and Meakin2005; Liu, Meakin & Huang Reference Liu, Meakin and Huang2006; Tartakovsky & Panchenko Reference Tartakovsky and Panchenko2016; Shigorina et al. Reference Shigorina, Kordilla and Tartakovsky2017; Arai et al. Reference Arai, Tartakovsky, Holt, Grace and Ryan2020), which have a scale dependency on $H$ . Our approach ensures that the units of $s$ are that of an interfacial tension. Taking $f_{\textit{ij}}(\tau )$ to be dimensionless, the units of (2.19) are

(2.22) \begin{align} \big[{\textrm{ms}}^{-2}\big] &= [{\textrm{m}}] \big[{\textrm{kg}}^{-1}\big] s_{\textit{ij}} [1] [{\textrm{m}}] \big[{\textrm{m}}^{-1}\big], \nonumber \\ s_{\textit{ij}} &= \big[\textrm{Nm}^{-1}\big]. \end{align}

Thus, the units of $s$ are ${\textrm{Nm}}^{-1}$ , analogous to the true surface tension $\sigma$ , which ensures the simulated surface tension is independent of the resolution.

While the PF-SPH model may seem an oversimplification of the complex (and small-scale) behaviour that is molecular cohesion and adhesion, it is best thought of as a coarse approximation of these phenomena, at the resolution of the SPH discretisation. We will show in the next section that the resulting surface tension is linearly related to the pairwise-force strength $s_{\textit{ff}}$ , and that this simple model captures complex behaviour in a predictable way.

3.1. Laplace pressure

The first test we use will validate the surface tension of the fluid in a static scenario, independent of fluid–solid interactions. The Young–Laplace equation describes the difference in pressure $\Delta P$ due to surface tension $\sigma$ in a spherical droplet of radius $R$ (Bormashenko Reference Bormashenko2017)

(3.1) \begin{align} \Delta P = \frac {2\sigma }{R}. \end{align}

By testing the pressure difference at different droplet radii for a fixed inter-particle force strength $s_{\textit{ff}}$ , we can calibrate (and validate) the resulting surface tension $\sigma$ . Given that we only model the fluid phase, and therefore have $P_{\textit{out}} = 0$ , we need only calculate $P_{\textit{in}}$ . We do this by borrowing a technique from molecular dynamics, calculating the total pressure from a many-body simulation (Hoover Reference Hoover1998). When calculated this way, the pressure due to particle–particle interactions is called the virial pressure, and is defined by Allen & Tildesley (Reference Allen and Tildesley1989), Tartakovsky & Meakin (Reference Tartakovsky and Meakin2005) as

(3.2) \begin{gather} P(r) = \frac {1}{2{d}V(r)} \sum _{i\in \mathcal{I}(r)} \sum _{j=1}^N \boldsymbol{F}_{\textit{ij}} \boldsymbol{\cdot }(\boldsymbol{x}_i - \boldsymbol{x}_j), \end{gather}

where $d$ is the spatial dimension (taken here to be $d=3$ ), $V(r)$ is the volume of a sphere of radius $r$ and $\boldsymbol{F}_{\textit{ij}}$ is the sum of the pressure force and pairwise force that particle $i$ experiences due to particle $j$ . The outer summation (over $i$ ) includes only particles in the set $\mathcal{I}(r)$ of particles within a distance $r$ of the centre of the droplet. This so-called virial radius $r \lt R$ is used in place of the actual droplet radius $R$ to exclude the region near the interface (Tartakovsky & Meakin Reference Tartakovsky and Meakin2005). When $r = R$ , we see a divergence of the pressure near the free surface due to the neighbourhood deficiency in the SPH approximations there. Figure 3 shows the virial pressure measured at different radii, with the pressure clearly diverging as $r$ approaches $R$ . For accurate estimates of the virial pressure we take an average of $P(r)$ over the interval $r \in [R-3H, R-2H]$ . If $R \leqslant 4H$ , we take $r = R - 2H$ .

Figure 3. Virial pressure calculated using (3.2) at different radii $r$ , given here in units of the kernel support radius $H$ . The actual radius of this spherical droplet is $1\, \textrm{mm}$ , approximately $8.5H$ . The calculated virial pressure diverges for $r \gtrsim 6.5H$ due to the neighbourhood deficiency of the particles near the free surface. A red line shows the region over which we average $P(r)$ .

For the simulations, we initialise spherical droplets consisting of particles with properties given in table 1, in zero gravity. We randomly perturb each particle, by no more than $0.1\Delta x$ , to speed up their rearrangement due to inter-particle forces $\boldsymbol{F}^{(\textit{pf})}$ . We then allow the droplet to reach equilibrium over $200 \,{\textrm{ms}}$ before measuring the virial pressure. Figure 4 shows that the pairwise-force model reproduces the linear relationship $\Delta P \propto 2/R$ ; we can measure the surface tension at each value of $s_{\textit{ff}}$ as the slope of each of these lines. We also note that the standard deviation of the virial pressure over the range $r \in [ R-3H, R-2H ]$ is relatively small and thus does not introduce significant uncertainty in the calibrated surface tensions. Plotting the measured surface tension against the model parameter $s_{\textit{ff}}$ reveals a simple linear relationship in figure 5, namely

(3.3) \begin{align} \sigma = 30.96 s_{\textit{ff}}. \end{align}

This is consistent with the dimensional analysis in § 2.4, and makes the prescription of surface tension in the model very simple. We have tested this relationship for particle width $\Delta x$ as small as $2 \times 10^{-5} \,{\textrm{m}}$ and as large as $8 \times 10^{-5} \,{\textrm{m}}$ and found it to be independent of the resolution, as expected.

Table 1. Fluid properties in Laplace pressure simulations.

Figure 4. Laplace pressure validation using calculated virial pressures of spherical droplets. For different values of the cohesive force strength $s_{\textit{ff}}$ , the pressure follows $P \propto 2 / R$ , where $R$ is the radius of the droplet. Markers show measurements from simulations, and black lines show linear fits. Error bars show standard deviations of the virial pressure across the virial radii $r \in [R-3H, R-2H]$ . The slope of each line gives the surface tension $\sigma$ for the corresponding $s_{\textit{ff}}$ .

Figure 5. Calibrating the cohesive force strength $s_{\textit{ff}}$ to measured surface tension in two different tests. Circles show surface tensions calculated from the Laplace pressure $P = 2\sigma /R$ (see figure 4). Squares show surface tensions calculated from ellipsoidal droplet oscillations (e.g. figure 6). The fitted line shows a simple linear relationship between $s_{\textit{ff}}$ and the surface tension.

3.2. Oscillating droplets

With the surface tension now calibrated in a static scenario, we next validate the model for surface tension in a dynamic scenario. For this task we choose to study the oscillation of a free droplet that has been perturbed from its spherical equilibrium. The linear frequency of oscillation of an inviscid droplet (in the eigenmode of interest) was found by Rayleigh (Reference Rayleigh1879) (with a more succinct derivation given by Landau & Lifshitz (Reference Landau and Lifshitz1987)) to be

(3.4) \begin{align} f &= \frac {1}{\pi } \sqrt { \frac { 2 \sigma }{ R^3 \rho } }. \end{align}

With material properties given in table 2, we initialise a spherical droplet of radius $0.788$ mm with particles that we randomly perturb by no more than $0.1\Delta x$ , and allow the particle distribution to settle for 1 ms. Then we ‘stretch’ the spherical droplet into an ellipsoid with the coordinate transform

(3.5) \begin{align} \boldsymbol{x} \leftarrow \underbrace {\frac {1}{\sqrt [3]{abc}} \begin{bmatrix} a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{bmatrix}}_T \boldsymbol{x}. \end{align}

Since $\det (T) = 1$ , this transformation preserves volume, and therefore density of the fluid particles. The elements $a,b,c$ are the relative lengths of the ellipsoid’s semi-axes, which we take to be 1.0, 0.7 and 0.7, respectively.

Table 2. Fluid properties in oscillating droplet simulations.

The simulation then proceeds, with the diameter of the droplet oscillating over 3 ms as shown in figure 6 for $s_{\textit{ff}} = 0.00156$ (corresponding to $\sigma = 47 \,{\textrm{mN m}}^{-1}$ ). Despite the fluid having zero viscosity in the simulation, the oscillations are clearly damped – the amplitude decreases with each oscillation. Nair & Pöschel (Reference Nair and Pöschel2018) investigate possible causes of this damping, finding that some of the system’s energy is dissipated as the particles arrange themselves into a crystal-like lattice.

Despite these effects, we can recover the frequency of the oscillations to determine the surface tension of the droplet. A discrete Fourier transform of the oscillations (figure 7) shows a peak at 138 Hz. With (3.4), we can calculate the surface tension as $\sigma = 47 \,{\textrm{mN m}}^{-1}$ when $s_{\textit{ff}} = 0.00156$ . Figure 5 shows the surface tensions calculated this way, plotted against the cohesive force $s_{\textit{ff}} \in [0,0.002]$ , corroborating the linear relationship from the Laplace pressure test in § 3.1. Thus, we have shown in two independent tests – one static and the other dynamic – that the surface tension induced by the pairwise forces scales linearly with the interaction strength parameter $s_{\textit{ff}}$ .

Figure 6. The oscillating diameters of an inviscid, axis-aligned ellipsoidal droplet over 30 ms, for $s_{\textit{ff}} = 0.00156$ . Material properties are given in table 2.

Figure 7. The power spectrum of the oscillating droplet’s diameters (see figure 6), with a peak at 138 Hz, which gives $\sigma = 47 \,{\textrm{mN m}}^{-1}$ when $s_{\textit{ff}} = 0.00156$ from equation (3.4).

3.3. Young–Laplace profile fitting for contact angle

Having calibrated and validated the surface tension in both static and dynamic scenarios, we now turn our attention to wetting phenomena and the adhesive force strength $s_{\textit{fs}}$ . The angle that a droplet’s liquid–gas interface makes with its liquid–solid interface – the contact angle – is commonly used to measure the ability of the liquid to wet the solid (Bormashenko Reference Bormashenko2017). This angle is often measured experimentally by drawing a tangent line to the liquid–gas interface where it meets the solid; however, this would seem only to validate the model in the immediate vicinity of the contact line. Instead, we will use semi-analytical solutions for whole droplet shapes to validate the model and calibrate the adhesive force to the contact angle.

For this test, we initialise a spherical droplet at a distance of $H$ (the kernel support radius) above a flat solid surface, with a thickness of $H$ , to ensure the neighbourhood of a fluid particle near the boundary is not deficient. In all the contact angle simulations, the fluid particles have density $1000\, {\textrm{kg m}}^{-3}$ , viscosity $0.89\,\textrm{mPa}\boldsymbol{\cdot}\textrm{s}$ and artificial speed of sound $100\,{\textrm{m s}}^{-1}$ . We have simulated droplets of various surface tensions, volumes, and resolutions to ensure the model is not scale-dependent. Initially, the fluid particles are perturbed by no more than $0.1 \Delta x$ and the particle distribution is allowed to settle for $1\,{\textrm{ms}}$ under zero gravity, without interacting with the surface. Then, we switch on gravity ( $g = -9.81\, {\textrm{m s}}^{-2}$ ), allowing the droplet to fall and spread across the solid surface. After $50\,{\textrm{ms}}$ , the droplet settles and we are ready to extract the liquid–gas interface.

The liquid–gas interface of an SPH-simulated droplet is an isosurface of the density field defined by the SPH interpolation $\langle \rho (\boldsymbol{x}) \rangle$

(3.6) \begin{equation} \mathcal{S} = \left \{ \boldsymbol{x} \,:\, \sum _{j \,\in \,\mathcal{I}_{{f}}} m_j W(\boldsymbol{x} - \boldsymbol{x}_j; H) = \frac {\rho _0}{2} \right \}, \end{equation}

where $\mathcal{I}_{{f}}$ is the index set of fluid particles. We realise the implicit surface $\mathcal{S}$ using the marching cubes algorithm, with a grid spacing of $0.1\Delta x$ . Then, to determine the contact angle of the droplet, we fit the shape of an axisymmetric sessile droplet determined by the Young–Laplace (Y–L) equation (Danov et al. Reference Danov, Dimova, Ivanov and Novev2016) to the vertices of the SPH droplet isosurface in cylindrical polar coordinates $(R_j,Z_j)$ where the origin is given by the droplet’s centre of mass.

In order to fit the data points to the Y–L surface we minimise the function

(3.7) \begin{align} L(\theta _{\textit{CA}},Z_{\textit{shift}}) = \sum _j \frac {1}{R_j + \epsilon } \, \textrm{mindist}\left(R_j,\!Z_j-Z_{\textit{shift}};\theta _{\textit{CA}}\right)^2, \end{align}

where $\theta _{\textit{CA}}$ is the contact angle, $Z_{\textit{shift}}$ is a vertical shift of the data position and $\textrm{mindist}(R,Z;\theta _{\textit{CA}})$ is the minimum distance from the point $(R,Z)$ to the Y–L surface for the specified volume, surface tension and contact angle. The circumference of the droplet, and therefore the number of expected isosurface vertices, increases linearly with $R$ . Thus, to account for the varying density of the samples, we weight the errors with the factor $1 / (R_j + \epsilon )$ , where $\epsilon \gt 0$ is a small constant to avoid division by zero.

To construct the distance function $d(R,Z;\theta _{\textit{CA}})$ we must first compute the Y–L droplet shape for a given volume $V_0$ , surface tension $\sigma$ and contact angle $\theta _{\textit{CA}}$ . Following Danov et al. (Reference Danov, Dimova, Ivanov and Novev2016) (using equations first reported by Hartland & Hartley (Reference Hartland and Hartley1976, p. 40)), we solve the system of ordinary differential equations

(3.8) \begin{align} \begin{split} \frac {{\textrm{d}}r}{{\textrm{d}}s}&=\cos \theta ,\\ \frac {{\textrm{d}}z}{{\textrm{d}}s}&=\sin \theta ,\\ \frac {{\textrm{d}}\theta }{{\textrm{d}}s}&=\frac {2}{\mathcal{B}}-\frac {\sin \theta }{r}+\frac {\rho g}{\sigma }z,\\ \frac {{\textrm{d}}V}{{\textrm{d}}s}&=r^2\sin \theta , \end{split} \end{align}

from the ‘initial condition’ $r=z=\theta =V=0$ until $\theta =\theta _{\textit{CA}}$ , where $\mathcal{B}$ is the curvature at the apex of the droplet, chosen such that the correct volume is achieved. Numerically, the Y–L drop surface is given as a series of points $(r_i,z_i)$ with their gradient as an angle $\theta _i$ .

Figure 8. Schematic of the transformation from the cylindrical polar world coordinates to the local polar coordinates of a surface segment.

We approximate the normal distance to the surface by transforming the coordinates for a given surface segment $(r_i,z_i)-(r_{i+1},z_{i+1})$ into a local polar coordinate system

(3.9) \begin{align} (\varrho ,\psi )\rightarrow \left (\sqrt {(r_i-r_c)^2+(z_i-z_c)^2},\arctan \frac {z_i-z_c}{r_i-r_c}\right ), \end{align}

where the centre $(r_c,z_c)$ is given by the intersections of two lines passing through each of the end points of the segment perpendicular to their gradients (figure 8). That is

(3.10) \begin{align} r_c &= \frac {\sin \theta _i(z_{i+1}-z_i)+\cos \theta _i(r_{i+1}-r_i)}{\sin \theta _i\cos \theta _{i+1}-\cos \theta _i\sin \theta _{i+1}}\sin \theta _i+r_i, \end{align}

(3.11) \begin{align} \,\,\,\,\,z_c &= -\frac {\sin \theta _i(z_{i+1}-z_i)+\cos \theta _i(r_{i+1}-r_i)}{\sin \theta _i\cos \theta _{i+1}-\cos \theta _i\sin \theta _{i+1}}\cos \theta _i+z_i. \end{align}

The surface location along the segment is then given by $\varrho _i(\psi )=\varrho _i+(\varrho _{i+1}-\varrho _i)(3t^2-2t^3)$ , where $t=(\psi -\psi _i)/(\psi _{i+1}-\psi _i)$ . The difference $P-\varrho _i(\varPsi )$ is taken as the approximation of the normal distance from the point $(R,Z)$ to the segment. Therefore the minimum distance from the point $(R,Z)$ to the surface is the minimum distance over all segments for which the surface interpolation is valid

(3.12) \begin{align} \textrm{mindist}(R,Z;\theta _{\textit{CA}}) = \min _i \min _{\varPsi \in [\psi _i,\psi _{i+1}]}\left |P-\varrho _i(\varPsi )\right |. \end{align}

Figure 9 shows an example of an SPH droplet isosurface and its best fit Y–L shape. Also shown are the Y–L shapes for the contact angles $\theta _{\textit{low}}$ and $\theta _{\textit{high}}$ , such that $\theta _{\textit{CA}} \in [\theta _{\textit{low}}, \theta _{\textit{high}}]$ and $L(\theta _{\textit{low}}, Z_{\textit{shift}}) = L(\theta _{\textit{high}}, Z_{\textit{shift}}) = \texttt {tol}$ , a fitting tolerance. For each simulated droplet, this approach gives a feasible range of contact angles as well as the optimal fit. The results of measuring the contact angles of simulated droplets in this way are shown in figures 10 and 11, grouped by surface tension and volume, respectively, and plotted against the ratio $s_{\textit{fs}} / s_{\textit{ff}}$ . The contact angles for different surface tensions and volumes coincide for equal values of this ratio, suggesting that the contact angle produced by the PF-SPH model is scale-independent and making the specification of the contact angle straightforward in practice. We find that for higher contact angles, where larger changes in $\theta _{\textit{CA}}$ result in only small changes in the droplet profile, the feasible range $\theta _{\textit{high}} - \theta _{\textit{low}}$ is relatively large. For intermediate contact angles between $30^\circ$ and $120^\circ$ , however, the feasible region is smaller and we can be more precise in our specification of $\theta _{\textit{CA}}$ via $s_{\textit{fs}}$ .

Figure 9. An example of a Y–L profile (a solution to the system (3.8)) fitted to an SPH droplet. Isosurface vertices are shown in cylindrical coordinates about the centre of the top of droplet. The solid line shows the profile of an axisymmetric droplet satisfying the Y–L equation, fit to the isosurface vertices with the contact angle as a free parameter. Dashed and dot-dashed lines show the extent of contact angles satisfying the fitting tolerance.

Figure 10. Contact angles, measured by fitting Y–L profiles to sessile droplet density isosurfaces, are plotted against the ratio of pairwise-force strengths, $s_{\textit{fs}}/s_{\textit{ff}}$ , and grouped by surface tension. A droplet is excluded if the volume enclosed by the density isosurface differs from the actual droplet’s volume by more than 10 %. Inset: the curves without error bars, for clarity.

With this validation of the droplet shape and calibration of the contact angle, we now have a simple procedure for specifying the surface tension $\sigma$ and contact angle $\theta _{\textit{CA}}$ in the PF-SPH scheme. We can calculate the required cohesive strength $s_{\textit{ff}} = \sigma / 30.96$ from (3.3). Then we can interpolate the data from figure 10 to find the required ratio $s_{\textit{fs}}/s_{\textit{ff}}$ , and, knowing $s_{\textit{ff}}$ , calculate $s_{\textit{fs}}$ .

4.1. Microstructured surface

Synthetic surfaces with manufactured microscopic roughness attract interest from scientists and engineers alike for their potential commercial applications (e.g. self-cleaning surfaces, reduced drag on marine vessels, collecting freshwater from fog (Chamakos, Sema & Papathanasiou Reference Chamakos, Sema and Papathanasiou2021). Surfaces with sharp geometric features, however, are challenging to incorporate into computational models with explicit boundary conditions on the contact line. Here we will demonstrate the ease with which the calibrated PF-SPH model can be used to calculate the shape of a sessile droplet on a geometrically patterned surface. To do this, we will follow the experimental set-up of Dupuis & Yeomans (Reference Dupuis and Yeomans2005) in simulating a $3 \, \unicode{x03BC} \textrm{l}$ droplet on a surface featuring square pillars, as shown in figure 12. The parameters used in these simulations are given in table 3.

Table 3. Fluid properties in pillared substrate simulation in figure 13.

Figure 13 shows a comparison between two almost identical numerical experiments – the only difference being the initial kinetic energy of the droplet before ‘impact’ with the surface. Figure 13(a) shows a settled droplet that was initialised with zero velocity. This droplet sits atop the pillars on the substrate, without enough energy to infiltrate the gaps between them. The droplet in figure 13(b) was initialised with a velocity in the vertical direction of $-0.1\,{\textrm{m s}}^{-1}$ , allowing it to overcome the energy barrier discussed by Dupuis & Yeomans (Reference Dupuis and Yeomans2005) and transition from a ‘suspended’ to a ‘collapsed’ state, in which the fluid has infiltrated between the pillars. That the present model reproduces this qualitative behaviour suggests it should be applicable to problems with complex surfaces.

Figure 13. A comparison of nearly identical simulations of a droplet settling on the physically patterned surface of square pillars depicted in figure 12. The addition of a small impact velocity changes the wetting behaviour significantly (note that these are side views of three-dimensional simulations). (a) After 10 ms of zero gravity, then 30 ms under 1 g, this droplet has settled, suspended atop the pillars, (b) 30 ms after a 0.1 m s^–1 impact, this droplet has infiltrated the pillars.

4.2. Chemically patterned surface

Another type of surface heterogeneity that is widely studied is a chemically patterned surface. By manufacturing a surface with, for example, alternating hydrophobic and hydrophilic stripes, one can influence the shape or motion of a droplet. Varagnolo et al. (Reference Varagnolo, Schiocchet, Ferraro, Pierno, Mistura, Sbragaglia, Gupta and Amati2014) note that controlling the motion of very small droplets in this way is one of the key problems in the design of reliable microfluidic devices. We will demonstrate that our PF-SPH is also a suitable model for the shape of a droplet settled on such a patterned surface. To include the variable wettability in the PF-SPH model, we simply modify the adhesive force strength $s_{\textit{fs}}$ across the surface, depending on whether a boundary particle is hydrophobic (low wettability, high $\theta _{\textit{CA}}$ ) or hydrophilic (high wettability, low $\theta _{\textit{CA}}$ ). The fluid properties for this simulation are given in table 4. For the hydrophobic surface type, we used $s_{\textit{fs}}/s_{\textit{ff}}=1.4$ for an equilibrium contact angle of $\theta _{\textit{hydrophobic}}\approx 155^\circ$ , and for the hydrophilic surface type, $s_{\textit{fs}}/s_{\textit{ff}}=3.5$ for $\theta _{\textit{hydrophilic}}\approx 45^\circ$ .

Table 4. Fluid properties in patterned substrate simulation in figure 14.

One of the advantages of our PF-SPH scheme is that the relationship between the contact angle and the position and motion of the contact line is not explicitly prescribed; instead, the simulated contact angle is ultimately a function of the adhesive force strength $s_{\textit{fs}}/s_{\textit{ff}}$ and the geometry of the substrate. Thus, we may use the model to study the transition between contact angles on the hydrophobic and hydrophilic zones. To visualise the contact angle on the contact line, we use the surface normal of the density isosurface (3.6) at a distance of $H/2$ above the substrate. Figure 14 shows side views of the droplet and the preferential spreading on the hydrophilic stripe. The contact line is shown, coloured by the contact angle, with smooth transitions between the contact angles $\theta _{\textit{hydrophobic}}$ and $\theta _{\textit{hydrophilic}}$ . The side views reveal that the contact angle the droplet makes is vastly different on the hydrophilic substrate sections as compared with the hydrophobic sections. At its extremities, the measured contact angle reaches the specified hydrophobic value of $155^\circ$ , whereas the specified hydrophilic contact angle of $45^\circ$ is not realised in the simulation, with the lowest measured contact angle being $55^\circ$ . That our model is capable of producing such behaviour suggests it could be used to design and study manufactured, variable-wettability substrates for droplet motion control.

Figure 14. A droplet settles on a flat hydrophobic surface ( $\theta _{\textit{CA}} \approx 155^\circ$ ) with a hydrophilic stripe ( $\theta _{\textit{CA}} \approx 45^\circ$ ). (a) isometric view. (b) the contact angle is plotted around the contact line, showing a transition from one equilibrium contact angle to another. (c) side views of the droplet, showing the significant difference in apparent contact angles on the different surface types. Surface representations of the droplet are obtained from the density isosurface; (3.6).

4.3. Wheat leaf

The broader context and motivation of this work is to enable the study of droplets impacting and settling on plant leaves, for which a key challenge is calculating the shape of a droplet on a plant leaf with microscopic roughness, and chemical heterogeneity (Mayo et al. Reference Mayo, McCue, Moroney, Forster, Kempthorne, Belward and Turner2015). Here, we will demonstrate the model’s applicability to complex surfaces – specifically, a reconstructed wheat leaf surface. This virtual wheat leaf was reconstructed from a micro computed tomography scan using a radial basis function partition of unity method as described in related work: Whebell et al. (Reference Whebell, Moroney, Turner, Pethiyagoda, Wille, Cooper-White, Kumar, Taylor and McCue2025). We discretise this surface for an SPH simulation by taking a block of boundary particles on a regular grid, and discarding any for which the implicit surface indicator function is negative ( $\mathcal{F}(\boldsymbol{x}) \lt 0$ ). Manually selected boundary particles are labelled as hairs and assigned a higher adhesive force strength to reflect the surface chemistry of real wheat leaves. The parameters for this simulation are given in table 5.

Table 5. Fluid properties in wheat leaf simulation in figure 15.

Figure 15. A PF-SPH simulation of a $0.27\,\unicode{x03BC} \textrm{L}$ water droplet settled on the surface of a wheat leaf after $12\,{\textrm{ms}}$ . The droplet’s liquid–gas interface is rendered by polygonising an isosurface of the density field, $\langle \rho (\boldsymbol{x})\rangle = \rho _0 / 2$ . The leaf surface is rendered similarly, using only boundary particles, and coloured by height for clarity. A hair on the leaf surface is visible under the translucent droplet.

We initialised the simulation with a sphere of fluid particles just above a hair of the wheat leaf, with zero impact velocity, to study how the droplet settles onto the surface. Figure 15 shows the shape of the droplet at 50 ms, once it has lost momentum and settled on the leaf. We can visualise the contact line with a contour line on the density isosurface, where the approximate minimum distance to the leaf surface is $2\Delta x$ . Specifically, figure 16 shows the contact line visualised as the set of points

(4.1) \begin{align} \Bigg \{ \boldsymbol{x} \,:\, \sum _{j \,\in \,\mathcal{I}_{{f}}} m_j W(\boldsymbol{x} - \boldsymbol{x}_j; H) = \frac {\rho _0}{2} \Bigg \} \cap \left \{ \boldsymbol{x} \,:\, \min _{\boldsymbol{y} \in \mathcal{L}} \| \boldsymbol{x} - \boldsymbol{y} \| = 2 \Delta x \right \}, \end{align}

where $\mathcal{L}$ is the set of leaf surface points. Note the highly irregular contact line following the natural curvature of the wheat leaf. This droplet shape could be used to measure the wetted area of the leaf surface or serve as an initial condition for an evaporation model in future work.

Figure 16. The contact line of the droplet on a wheat leaf depicted in figure 15, visualised by computing the set of points on the density isosurface that are $2\Delta x$ from the surface.

4.4. Discussion on pairwise-force profiles

The exact design of the pairwise-force profile $f_{\textit{ij}}$ is in need of more investigation before the PF-SPH model can be applied as readily as, for example, the CSF formulation. Literature on the effect of the force profile is scarce. Tartakovsky & Panchenko (Reference Tartakovsky and Panchenko2016) report some analytical results predicting the surface tension from $f_{\textit{ij}}(\tau )$ and $s_{\textit{ij}}$ in multiphase simulations, but we could not verify these results in our single phase simulations. Our force profile causes most fluid–fluid particle interactions to be attractive, which could cause excess numerical stress in the fluid. In the simulations tested, the mostly attractive pairwise force has the somewhat desirable effect of keeping the particle distribution ordered – for example, in a high velocity impact event – which ensures that interpolation error is low. More investigation is needed to quantify the dissipative effect of the pairwise forces before this model is applied to viscosity-dependent scenarios.