Computing sessile droplet shapes on arbitrary surfaces with a new pairwise force smoothed particle hydrodynamics model

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

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

which, as we shall see, lends itself naturally to discretisation using smoothed particle hydrodynamics.

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 (1)-(2) is

where $\bm{\mathrm{x}}_{ij}=\bm{\mathrm{x}}_{i}-\bm{\mathrm{x}}_{j}$ for notational convenience, $W_{ij}=W(\bm{\mathrm{x}}_{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 (5). 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 $\mathrm{D}/\mathrm{D}t$ to the total derivative $\mathrm{d}/\mathrm{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, 1995) as the SPH kernel, which we parameterise by $H$, its radius of support, also called the kernel cutoff radius:

Wendland functions were developed independently of SPH for their smoothness properties, but have been found to give accurate and stable SPH results (Dehnen and Aly, 2012). 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 $\bm{\mathrm{F}}^{\mathrm{(pf)}}$, which will be the main focus of this work, the discretisation summarised above is well established in the literature: see, for example, Monaghan (2005). 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 $\bm{\mathrm{F}}^{\mathrm{(pf)}}$ that yield stable and accurate simulations of droplets.

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. (2011) 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 (3) and (4) has a time complexity of $\mathcal{O}(N^{2})$. This is made more efficient by pre-computing a neighbour list

for each fluid particle $i$. Any sum over $j=1,\dots,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. (2011) for CPU implementations, or for GPU implementations: Domínguez et al. (2011); Crespo et al. (2011, 2015); Domínguez et al. (2022). Simulations were run using up to 64 processors in parallel.

Surface tension and wetting are both effects of intermolecular forces (Bormashenko, 2017). A molecule at an interface misses approximately half of its interactions with neighbours when compared to 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, 1924). The basis for these forces in SPH is empirical; nevertheless, in Section 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 (4) via the term $\bm{\mathrm{F}}_{i}^{\mathrm{(pf)}}$:

where $s_{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_{ij}(\|\bm{\mathrm{x}}_{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 and Panchenko, 2016), SPH kernels (Shigorina et al., 2017; Kordilla et al., 2013), or part of a cosine curve (Tartakovsky and Meakin, 2005; Nair and Pöschel, 2018). For simplicity, we will instead use polynomials for the force profile, with some intuitive constraints at key points. Those constraints are:

The last constraint, controlling the zero crossing of the force profile, is intentionally different for the fluid-fluid ($\mathrm{ff}$) and the fluid-solid ($\mathrm{fs}$) interactions. 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. Fourth-degree polynomials are simple to construct from these five constraints, and so we have:

Figure 2 shows these force profiles over the range of possible particle distances.

The parameter $s_{ij}$ in equation (7) controls the strength of the pairwise force and, like $f_{ij}(\tau)$, depends on the labels of the particles $i$ and $j$: fluid or solid. We denote the cohesive force strength $s_{\mathrm{ff}}$ (fluid-fluid) and the adhesive force strength $s_{\mathrm{fs}}$ (fluid-solid). Note the inclusion of the length scale $H$ in the form of $\bm{\mathrm{F}}_{i}^{\mathrm{(pf)}}$ in equation (7): this is a departure from established pairwise-force methods (Kordilla et al., 2013; Nair and Pöschel, 2018; Shigorina et al., 2017; Tartakovsky and Panchenko, 2016; Tartakovsky and Meakin, 2005; Arai et al., 2020; Liu et al., 2006), which have a scale dependency on $H$. Our approach ensures that the units of $s$ are that of an interfacial tension. Taking $f_{ij}(\tau)$ to be dimensionless, the units of equation (7) are

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

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, 2017):

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

where $d$ is the spatial dimension (taken here to be $d=3$), $V(r)$ is the volume of a sphere of radius $r$, and $\bm{\mathrm{F}}_{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<R$ is used in place of the actual droplet radius $R$ to exclude the region near the interface (Tartakovsky and Meakin, 2005). 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\leq 4H$, we take $r=R-2H$.

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 $\bm{\mathrm{F}}^{\mathrm{(pf)}}$. We then allow the droplet to reach equilibrium over $200\,\mathrm{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_{\mathrm{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_{\mathrm{ff}}$ reveals a simple linear relationship in Figure 5, namely

This is consistent with the dimensional analysis in Section 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\cdot 10^{-5}\,\mathrm{m}$ and as large as $8\cdot 10^{-5}\,\mathrm{m}$ and found it to be independent of the resolution, as expected.

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 (1879) (with a more succinct derivation given by Landau and Lifshitz (1987)) to be

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 1ms.

Then we ‘stretch’ the spherical droplet into an ellipsoid with the coordinate transform

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.

The simulation then proceeds, with the diameter of the droplet oscillating over 3ms as shown in Figure 6 for $s_{\mathrm{ff}}=0.00156$ (corresponding to $\sigma=47\,\mathrm{mN/m}$). Despite the fluid having zero viscosity in the simulation, the oscillations are clearly damped – the amplitude decreases with each oscillation. Nair and Pöschel (2018) 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 138Hz. With equation (11), we can calculate the surface tension as $\sigma=47\,\mathrm{mN/m}$ when $s_{\mathrm{ff}}=0.00156$. Figure 5 shows the surface tensions calculated this way, plotted against the cohesive force $s_{\mathrm{ff}}\in[0,0.002]$, corroborating the linear relationship from the Laplace pressure test in Section 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_{\mathrm{ff}}$.

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_{\mathrm{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, 2017). 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\,\mathrm{kg/m^{3}}$, viscosity $0.89\,\mathrm{mPa\cdot s}$, and artificial speed of sound $100\,\mathrm{m/s}$. 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\,\mathrm{ms}$ under zero gravity, without interacting with the surface. Then, we switch on gravity ($g=-9.81\,\mathrm{m/s^{2}}$), allowing the droplet to fall and spread across the solid surface. After $50\,\mathrm{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(\bm{\mathrm{x}})\rangle$:

where $\mathcal{I}_{\mathrm{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., 2016) 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

where $\theta_{\mathrm{CA}}$ is the contact angle, $Z_{\text{shift}}$ is a vertical shift of the data position, and $\mathrm{mindist}(R,Z;\theta_{\mathrm{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>0$ is a small constant to avoid division by zero.

To construct the distance function $d(R,Z;\theta_{\mathrm{CA}})$ we must first compute the Y-L droplet shape for a given volume $V_{0}$, surface tension $\sigma$, and contact angle $\theta_{\mathrm{CA}}$. Following Danov et al. (2016) (using equations first reported by Hartland and Hartley (1976, p. 40)), we solve the system of ordinary differential equations

from the “initial condition” $r=z=\theta=V=0$ until $\theta=\theta_{\mathrm{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}$.

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

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,

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}(\Psi)$ 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,

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_{\text{low}}$ and $\theta_{\text{high}}$, such that $\theta_{\text{CA}}\in[\theta_{\text{low}},\theta_{\text{high}}]$ and $L(\theta_{\text{low}},Z_{\text{shift}})=L(\theta_{\text{high}},Z_{\text{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_{\mathrm{fs}}/s_{\mathrm{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_{\text{CA}}$ result in only small changes in the droplet profile, the feasible range $\theta_{\text{high}}-\theta_{\text{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_{\text{CA}}$ via $s_{\mathrm{fs}}$.

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_{\text{CA}}$ in the PF-SPH scheme. We can calculate the required cohesive strength $s_{\mathrm{ff}}=\sigma/30.96$ from equation (10). Then we can interpolate the data from Figure 10 to find the required ratio $s_{\mathrm{fs}}/s_{\mathrm{ff}}$, and, knowing $s_{\mathrm{ff}}$, calculate $s_{\mathrm{fs}}$.

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 et al., 2021)). 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 setup of Dupuis and Yeomans (2005) in simulating a $3\,\mathrm{\mu L}$ droplet on a surface featuring square pillars, as shown in Figure 12. The parameters used in these simulations are given in Table 3.

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\,\mathrm{m/s}$, allowing it to overcome the energy barrier discussed by Dupuis and Yeomans (2005) 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.

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. (2014) 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_{\mathrm{fs}}$ across the surface, depending on whether a boundary particle is hydrophobic (low wettability, high $\theta_{\text{CA}}$) or hydrophilic (high wettability, low $\theta_{\text{CA}}$). The fluid properties for this simulation are given in Table 4. For the hydrophobic surface type, we used $s_{\mathrm{fs}}/s_{\mathrm{ff}}=1.4$ for an equilibrium contact angle of $\theta_{\mathrm{hydrophobic}}\approx 155^{\circ}$, and for the hydrophilic surface type, $s_{\mathrm{fs}}/s_{\mathrm{ff}}=3.5$ for $\theta_{\mathrm{hydrophilic}}\approx 45^{\circ}$.

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_{\mathrm{fs}}/s_{\mathrm{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 (12) 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_{\mathrm{hydrophobic}}$ and $\theta_{\mathrm{hydrophilic}}$. The side views reveal that the contact angle the droplet makes is vastly different on the hydrophilic substrate sections as compared to 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.

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., 2015). 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 microCT scan using a radial basis function partition of unity method as described in related work: Whebell et al. (2023). 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}(\bm{\mathrm{x}})<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.

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 50ms, 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

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.

The exact design of the pairwise force profile $f_{ij}$ is in need of more investigation before the PF-SPH model can be applied as readily as, for example, the continuum surface force formulation. Literature on the effect of the force profile is scarce. Tartakovsky and Panchenko (2016) report some analytical results predicting the surface tension from $f_{ij}(\tau)$ and $s_{ij}$ in multiphase simulations, but we could not verify these results in our single phase simulations.

We have designed the fluid-fluid pairwise force curve to have its zero crossing at approximately the ‘rest distance’ of the particles. To do this, we imagined the particles as spheres, packed optimally 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 of diameter $\Delta x$ 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 was the geometric motivation behind the location of the zero crossing of the force profile in (7).

Clearly, however, this will lead to most fluid-fluid pairwise force interactions being 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.