3.1. Normal reaction

When a particle comes into contact with a flat solid surface, we define the overlap as (see figure 1 a)

(3.1) \begin{equation} \delta = R_c -(z_p-z_s), \end{equation}

Figure 1.

Definition of the overlap between two solids in contact. (a) Definition of the overlap between a spherical particle and a flat substrate. (b) Definition of the overlap between two spherical particles.

where $R_c$ is the mechanical radius of the particle, $z_p$ is the $z$ -coordinate of the centre of the particle and $z_s$ is the $z$ -coordinate corresponding to the top of the solid substrate. The contact force is modelled with a spring–dashpot model. For a perfectly elastic contact without any dissipation, the normal contact force is given by $F^0_n=k_n\delta$ , where $k_n$ is the spring stiffness for overlap in the normal direction. But since energy is generally dissipated during realistic collisions, a damping term is added to the expression such that the magnitude of the normal reaction force becomes

(3.2) \begin{equation} F_n=k_n\delta + \gamma _nv_n. \end{equation}

Here, the second term is a damping term, where $\gamma _n$ is the damping constant for contact in the normal direction and $v_n$ is the magnitude of the relative velocity between the particle and the surface in the direction perpendicular to the surface. The normal reaction force on the particle acts perpendicular to the flat surface. In our model, we only account for the force on the particle due to the flat substrate and ignore the consequence of the equal and opposite force on the substrate due to the particle. The rationale for ignoring the force on the substrate is that we only consider solid substrates that are much larger and heavier than the particle. Thus, we keep the position of the substrate fixed and do not evolve it when a net force acts on it.

When two spherical particles $i$ and $j$ come into contact, we define the overlap between them as (figure 1 b)

(3.3) \begin{equation} \delta _{ij} = R_i + R_{\!j} - |\boldsymbol{X}_{\!i}-\boldsymbol{X}_{\!j}|, \end{equation}

where $\boldsymbol{X}_{\!i}$ and $\boldsymbol{X}_{\!j}$ are the position vector of the centre of mass of particle $i$ and $j$ , respectively. When the overlap is positive ( $\delta \gt 0$ ), the magnitude of this force acting each particle is given by

(3.4) \begin{equation} F_n=k_n\delta _{ij}+ \gamma _nv_{ij}. \end{equation}

Here, $v_{ij}=|(\boldsymbol{v}_{\!i}-\boldsymbol{v}_{\!j})\boldsymbol{\cdot }\boldsymbol{n}_{\!ij}|$ is the component of the relative velocity that points along the direction joining their centre of masses, where $\boldsymbol{n}_{\!ij}=(\boldsymbol{X}_{\!i}-\boldsymbol{X}_{\!j})/|\boldsymbol{X}_{\!i}-\boldsymbol{X}_{\!j}|$ . Both the particles experience an equal normal reaction force, but the force acts in the opposite direction. The direction of the force on the particles points towards their respective centre of mass. For particle $i$ , the normal force reaction force vector is $F_n\boldsymbol{n}_{\!ij}$ , whereas for particle $j$ , it is $-F_n\boldsymbol{n}_{\!ij}$ , where $F_n$ is given by (3.4).

3.2. Sliding friction

Sliding friction becomes active when two solid surfaces slide relative to each other. Here, we use the approach proposed by Luding (Reference Luding2008) to model sliding friction as arising from a tangential spring that resists relative sliding motion between the two solid surfaces at the contact point. For an initially stationary particle, the spring restores the particle to its initial position when the applied force is below the Coulomb threshold, $F_C=\mu F_n$ . When the applied force exceeds this limit, the particle can move continuously while experiencing friction. In the following, we briefly outline the algorithm that we used to implement these principles. A complete description is provided in the paper by Luding (Reference Luding2008).

When a contact is active (i.e. positive overlap), we use a linear spring–dashpot model to compute the sliding friction force acting tangential to the contact,

(3.5) \begin{equation} \boldsymbol{F}^0_t(t)=-k_t\boldsymbol{\xi }(t)-\gamma _t\boldsymbol{v}_{\!t}(t), \end{equation}

where $k_t$ is the tangential (or sliding) spring stiffness, $\boldsymbol{\xi}$ is the extension of the tangential spring, $\gamma _t$ is the tangential (or sliding) damping constant and $\boldsymbol{v}_{\!t}$ is the relative velocity between the contacting surfaces at the contact point.

When the restoring force due to the tangential spring is below the Coulomb threshold $(|\boldsymbol{F}_{\!t}^0|\leq F_C)$ , the particle is in the static regime and the restoring force balances the applied force. In the opposite case when $|\boldsymbol{F}_{\!t}^0|\gt F_C$ , the sliding regime becomes active and the magnitude of the friction force is set to $F_C$ . By accounting for these two cases, the tangential force applied to the particle can be written as

(3.6) \begin{equation} \boldsymbol{F}_{\!t}=\begin{cases} \boldsymbol{F}_{\!t}^0, & |\boldsymbol{F}_{\!t}^0|\leq F_C,\\[5pt] F_C\frac {\boldsymbol{F}_{\!t}^0}{|\boldsymbol{F}_{\!t}^0|}, & |\boldsymbol{F}_{\!t}^0|\gt F_C. \end{cases} \end{equation}

The extension of the tangential spring is incremented for the next time step, as follows:

(3.7) \begin{equation} \boldsymbol{\xi }(t+\Delta t)=\begin{cases} \boldsymbol{\xi }(t)+\boldsymbol{v}_{\!t} \Delta t, & |\boldsymbol{F}_{\!t}^0|\leq F_C,\\\\ -\frac {1}{k_t}\left (F_C\frac {\boldsymbol{F}_{\!t}^0}{|\boldsymbol{F}_{\!t}^0|}+\gamma _t\boldsymbol{v}_{\!t}\right )\!, & |\boldsymbol{F}_{\!t}^0|\gt F_C.\\ \end{cases} \end{equation}

In this expression, when the first condition is met (static regime), the extension of the tangential spring is incremented for the next time step. When the second condition is met (sliding regime), the extension of the spring is defined such that in the next time step, the particle experiences a friction force equal to the Coulomb threshold, $\boldsymbol{F}_{\!t}^0(t+\Delta t)\approx F_C\boldsymbol{F}_{\!t}^0(t)/|\boldsymbol{F}_{\!t}^0(t)|$ .

The torque that arises due to the tangential friction forces is

(3.8) \begin{equation} \boldsymbol{T}_{\!t}=(\boldsymbol{X}_{\!C}-\boldsymbol{X})\times \boldsymbol{F}_{\!t}, \end{equation}

where $\boldsymbol{X}_{\!C}$ is the position vector of the contact point and $\boldsymbol{X}$ is the position vector of the centre of the particle.

3.3. Rolling friction

In addition to sliding friction, objects also experience rolling friction (or rolling resistance) that restrains their rolling motion. In practice, rolling friction can arise due to roughness on the surface of the particle, surface deformation and asymmetries in the shape of the contact when particles roll.

The model for rolling friction follows a similar line of reasoning as the model for the sliding friction, as explained in detail in Luding (Reference Luding2008). We introduce three new parameters for the rolling friction model: the rolling spring stiffness $k_r$ , the rolling damping constant $\gamma _r$ and the coefficient of rolling friction $\mu _r$ . To model rolling friction, we follow the same procedure as described above for sliding friction, except for two important differences. First, for the relative velocity at the contact, we use the relative rolling velocity rather than the relative tangential velocity.

Second, the tangential force that is computed is merely a fictive force used to compute the torque. Only the torque is used in Newton’s second law for angular motion. The fictive tangential force does not contribute to the equation of motion for the translational motion because rolling resistance only affects the angular motion.

3.4. Net contact force

The net contact force acting on a particle when it makes contact with another solid body is obtained by adding the normal reaction, sliding friction and rolling friction:

(3.9) \begin{equation} \boldsymbol{F}_c=\boldsymbol{F}_{\!n}+\boldsymbol{F}_{\!t}. \end{equation}

Rolling friction is not included in the above expression because it only contributes to the torque and not to the net force. The total torque on the particle is given by

(3.10) \begin{equation} \boldsymbol{T}_{\!c}=\boldsymbol{T}_{\!t}+\boldsymbol{T}_{\!r}, \end{equation}

where $\boldsymbol{T}_{\!r}$ is the torque contribution due to rolling friction. Torques due to normal reaction do not enter this expression because normal forces always point towards the centre of the particle for spherical particles and thus do not generate torques.

3.5. Benchmarking contact forces

We tested the normal reaction force model by analysing the centre of mass position of a particle that was released from rest from an initial height, $h_0$ (figure 2 a). In the absence of damping ( $\gamma _n=0$ ), the collision was perfectly elastic, as expected. The sum of the kinetic and gravitational potential energy was conserved. The speed $v_0$ just before impact agrees with the prediction obtained by equating the loss in gravitational potential energy to the gain in kinetic energy, $v_{0}=\sqrt {2gh_0}$ , where $g$ is the acceleration due to gravity. Since energy is conserved during the impact, this speed is equal to the speed just after impact. The centre of mass bounced to the same height after every impact (figure 2 b, solid blue line). Between successive impacts, the centre of mass height varied with time as expected according to the equation of motion, $h(t)=v_0 t-gt^2/2=\sqrt {2gh_0}t-gt^2/2$ .

Figure 2.

Particle of radius $R_p$ bouncing off a solid surface. (a) When the spherical solid particle makes contact with the solid surface, the contact force is modelled by a spring–dashpot system. Contact forces are only active when the particle and the surface overlap. (b) Centre of mass height of the particle (relative to the top of the solid surface) as a function of time with and without damping. In the absence of damping ( $k_n=2.5, \gamma _n=0$ ), the particle bounces to the same height each time (solid blue line), indicating that energy is conserved. When there is damping ( $k_n=2.5, \gamma _n=1.0$ ), energy is lost during the contact and the maximum height decreases with each contact (dashed black line). The damped rebound height agrees well with the predicted height ((3.12), red points). The dotted green horizontal line indicates the solid surface.

When the collision was damped ( $\gamma _n\gt 0$ ), the centre of mass height decayed with time (figure 2 b, dashed line), as expected for situations where energy is lost during a collision. The theoretical rebound height can be predicted from the coefficient of restitution. In the DEM model, the coefficient of restitution is given by (Luding Reference Luding2008)

(3.11) \begin{equation} r=v_{n+1}/v_n=\exp (\eta _n/t_c), \end{equation}

where $v_n$ and $v_{n+1}$ are the impact speed after impact number $n$ and $(n+1)$ , respectively. The term $\eta _n$ is related to the damping constant according to $\eta _n = \gamma _n/(2m_p)$ and $t_c=\pi /\omega_{c}$ , where $\omega_{c} =\sqrt {(k_n/m_p - \eta _n^2)}$ . By balancing the kinetic energy just after impact with the potential energy when the particle reaches its maximum height and using (3.11) to relate the impact speed of the nth impact to the impact speed of the previous impact, we obtain a relation for the rebound height reached after the nth impact:

(3.12) \begin{equation} h_n = r^{2n} h_0. \end{equation}

The rebound height obtained in our numerical simulations (figure 2 b, dashed line) shows excellent agreement with the height predicted by (3.12) (figure 2 b, red points).

Figure 3.

Motion of a particle down an inclined plane for different friction scenarios: (a) in the absence of sliding friction, the particle slides without any rolling (the green dot on the circumference of the particle is a marker to track the angular motion of the particle); (b) when sliding friction is non-zero, the particle rolls with negligible slip at the contact between the particle and the flat surface; (c) when both sliding and rolling friction are non-zero, the particle rolls slower; (d) trajectory of the marker point (scatter points) when the particle rolls with negligible slip, for example in (b) and (c); The trajectory is described by the cycloid equation (black line), as expected for pure rolling.

Figure 3. Long description

Panel A: A series of images showing a particle sliding down an inclined plane without any rolling friction. The green dot on the circumference of the particle is a marker to track the angular motion of the particle. Panel B: A series of images showing a particle rolling down an inclined plane with non-zero sliding friction but zero rolling friction. The green dot on the circumference of the particle is a marker to track the angular motion of the particle. Panel C: A series of images showing a particle rolling down an inclined plane with both sliding and rolling friction. The green dot on the circumference of the particle is a marker to track the angular motion of the particle. Panel D: A scatter plot showing the trajectory of the marker point when the particle rolls with negligible slip. The scatter points represent the marker position, and the black line represents the cycloid equation.

To benchmark the sliding and rolling friction, we analysed the motion of a sphere down an inclined plane with an inclination angle of $45^\circ$ under gravity (figure 3). In the absence of sliding friction ( $\mu =0$ ) and tangential damping ( $\gamma _t=0$ ), the particle slid without any rolling (figure 3 a). The particle moved with a constant acceleration and no energy was lost. When the coefficient of sliding friction was increased to a finite value ( $\mu \gt 0$ ), friction force acted at the contact and opposed sliding. The particle began to roll (figure 3 b). In steady state, the particle displayed an almost pure rolling motion, with a small slip ratio (He et al. Reference He, Sandu, Mousavi, Shenvi, Braun, Kruger and Els2020), $s=\omega_{p} R_p/v_{CM}-1\approx -0.001$ , where $\omega_{p}$ is the angular velocity, $R_p$ is the particle radius, and $v_{CM}$ is the centre of mass velocity of the particle. This expression for the slip ratio quantifies the amount of slip between the particle and the surface at the contact point, with $s\to 0$ corresponding to pure rolling (zero slip) and $s \to -1$ corresponding to pure sliding (i.e. when $\omega_{p} =0)$ .

When both rolling friction and sliding friction were non-zero ( $\mu \gt 0$ , $\mu _r\gt 0$ , $\mu _r\lt \mu$ ), the particle still rolled, but with a smaller acceleration since energy is now lost due to both sliding and rolling friction (figure 3 c).

When a particle rolls with negligible slip at the contact point, the trajectory traced by a marker point on the circumference of the particle is expected to follow the cycloid equations given by the following parametric equations:

(3.13) \begin{align} x_{m}(t) &= R_p[\alpha_{p} (t)-\sin \alpha_{p} (t)], \quad z_{m}(t) =R_p[1-\cos \alpha_{p}(t)]+z_0. \end{align}

Here, $x_m$ and $z_m$ are the $x$ and $z$ coordinates (labelled in figure 3 c) traced by the marker point. $z_0$ is the vertical position of the top of the solid surface on which the particle rolls and $\alpha_{p}$ is the angular position of the marker point with respect to the centre of the particle. The cycloid equation provides an excellent fit to the trajectory traced by the marker point in figure 3(d), as expected for a spherical particle for which the sliding friction exceeds the rolling friction.

5.1. No-slip and no-penetration at solid boundaries

To achieve the no-slip and no-penetration boundary conditions at the solid wall, we modify the distribution functions of fluids $a$ and $b$ using the halfway bounce-back scheme as follows:

(5.1) \begin{equation} f^a_i(\boldsymbol{x}_{\!f},t)=f^a_{\bar {i}}(\boldsymbol{x}_b,t)+6w_i\rho ^a(\boldsymbol{x}_{\!f},t)\frac {\boldsymbol{u}_{\!w}(\boldsymbol{x}_{\!w},t)\boldsymbol{\cdot }\boldsymbol{c}_{\!i}}{c_s^2}, \end{equation}

(5.2) \begin{equation} f^b_i(\boldsymbol{x}_{\!f},t)=f^b_{\bar {i}}(\boldsymbol{x}_b,t)+6w_i\rho ^b(\boldsymbol{x}_{\!f},t)\frac {\boldsymbol{u}_{\!w}(\boldsymbol{x}_{\!w},t)\boldsymbol{\cdot }\boldsymbol{c}_{\!i}}{c_s^2}, \end{equation}

where $\boldsymbol{u}_{\!w}$ is the velocity of the particle wall and the subscript $\bar {i}$ denotes the discrete velocity direction opposite to $i$ . This update is performed after the streaming step. After streaming, the distribution function will have propagated into the solid. $\boldsymbol{x}_b$ are locations of solid nodes that neighbour at least one fluid node. When computing the wall velocity, both the translational and rotational motion are considered as follows:

(5.3) \begin{equation} \boldsymbol{u}_{\!w}=\boldsymbol{u}_{\!p}+\boldsymbol{\omega_{p} }\times (\boldsymbol{x}_{\!w}-\boldsymbol{x}_{\!p}), \end{equation}

where $\boldsymbol{x_p}$ and $\boldsymbol{u_p}$ is the position and velocity of the centre of mass of the particle, respectively. On fixed solid substrate, $\boldsymbol{u}_{\!w}=0$ since the solid is stationary.

5.2. Transferring momentum from fluid to particle

When fictive fluid particles bounce back on the surface of solid particles, they transfer momentum to the particle. We use the momentum exchange method to obtain the momentum transferred to the solid particle (Wen et al. Reference Wen, Zhang, Tu, Wang and Fang2014). Other methods that can used for the force coupling between the fluid and the particle include the immersed boundary method (Peskin Reference Peskin1977) or the stress integration method. Here, we choose the momentum exchange method because of its simplicity and its intuitive physical interpretation. With this method, momentum is transferred along each of the 19 discrete velocity directions. The momentum transferred to the particle along the $i$ direction per time step can be obtained from the distribution functions according to

(5.4) \begin{equation} \Delta \boldsymbol{p}(\boldsymbol{x}_{\!w})= \sum _{i=1}^{19}\left [ (\boldsymbol{c}_{\bar {i}}-\boldsymbol{u}_{\!w})f_{\bar {i}}^{*}-(\boldsymbol{c}_{\!i}-\boldsymbol{u}_{\!w})f_i \right ]\!, \end{equation}

where $f^*$ denotes the distribution function after bounce-back and $f$ denotes the distribution function before bounce-back. The summation runs over all discrete velocity directions. The change of momentum per time step gives the force on the particle. The total hydrodynamic force and torque on the particle is obtained by summing the local force at all solid wall nodes,

(5.5) \begin{equation} \boldsymbol{F}_{\!\eta} =\sum _{W\textit{all}} \Delta \boldsymbol{p}(\boldsymbol{x}_{\!w}), \end{equation}

(5.6) \begin{equation} \boldsymbol{T}_{\!\eta} =\sum _{W\textit{all}} (\boldsymbol{x}_{\!w}-\boldsymbol{x}_{\!p})\times \Delta \boldsymbol{p}(\boldsymbol{x}_{\!w}). \end{equation}

Here, the summation runs over all solid wall nodes.

5.3. Refilling of new fluid nodes that emerge on the receding side of a particle

When particles move, nodes that were previously in the solid region may be converted to fluid nodes. The distribution functions at these fresh fluid nodes must be assigned since we do not evolve the distribution function inside the solid. To assign distribution functions at these fresh fluid nodes, we calculate the corresponding equilibrium distribution function from density of the neighbouring fluid nodes and the macroscopic velocity of the solid wall as follows (Caiazzo Reference Caiazzo2008):

(5.7) \begin{equation} f^{a,eq}_i(\boldsymbol{x}_{\!n})=\rho ^a(\boldsymbol{x}_{\!n})w_i\left (1+\frac {\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{c}_{\!i}}{c_s^2}+\frac {(\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{c}_{\!i})^2}{2c_s^4}-\frac {\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{u}(\boldsymbol{x}_{\!n})}{2c_s^2}\right )\!, \end{equation}

(5.8) \begin{equation} f^{b,eq}_i(\boldsymbol{x}_{\!n})=\rho ^b(\boldsymbol{x}_{\!n})w_i\left (1+\frac {\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{c}_{\!i}}{c_s^2}+\frac {(\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{c}_{\!i})^2}{2c_s^4}-\frac {\boldsymbol{u}(\boldsymbol{x}_{\!n})\boldsymbol{\cdot }\boldsymbol{u}(\boldsymbol{x}_{\!n})}{2c_s^2}\right )\!, \end{equation}

where $\boldsymbol{x}_{\!n}$ denotes the position of a fluid node that neighbours a solid particle node. The local density at the fresh fluid node is computed by taking an average of the local density at all neighbouring fluid nodes that lie in any of the 19 directions of the D3Q19 velocity set,

(5.9) \begin{equation} \rho ^{a,b}(\boldsymbol{x}_{\!n})=\frac {1}{N_{\!f}}\sum ^{N_{\!f}}_{i \ \in \ \textit{ Fluid}} \rho ^{a,b}(\boldsymbol{x}_{\!n}+\boldsymbol{c}_{\!i}\Delta t,t), \end{equation}

where $N_{\!f}$ is the number of fluid nodes surrounding the fresh fluid node that is being considered. The summation excludes neighbouring solid nodes and freshly created fluid nodes. This averaging procedure does not lead to local instability in density because the two fluids have equal density in our model. The macroscopic velocity is obtained from the particle velocity as

(5.10) \begin{equation} \boldsymbol{u}(\boldsymbol{x}_{\!n})=\boldsymbol{u}_{\!p}+\boldsymbol{\omega_{p} }\times (\boldsymbol{x}_{\!n}-\boldsymbol{x}_{\!p}). \end{equation}

Here $\boldsymbol{v}_p$ is the particle centre of mass velocity. Note that while the fluid is constrained to a lattice grid, the position of the particle is continuous and is not set by the fluid grid.

5.4. Benchmarking hydrodynamic force

To validate the accuracy of the hydrodynamic force calculation, we simulated the migration of a neutrally buoyant rigid sphere in a cylindrical Poiseuille flow, as schematically illustrated in figure 5(a). In this simulation, the sphere is expected to eventually reach an equilibrium position within the cylindrical tube, a phenomenon known as the Segrè–Silberberg effect (Karnis et al. Reference Karnis, Goldsmith and Mason1966).

For this benchmark, we used the following parameter values in simulation units: tube radius $R_t=40$ , particle radius $R_p=12.2$ , fluid density $\rho _{\!f}=1.05$ , fluid kinematic viscosity $\eta /\rho _{\!f}=0.457$ and average velocity at inlet $v_{{in}}=0.0113$ . Initially, the particle was positioned at $l/R_t = 0.21$ , where $l$ is the distance from the central axis of the tube as shown in figure 5 a). These values were chosen to match the experimental parameters used by Karnis et al. (Reference Karnis, Goldsmith and Mason1966) to allow a direct comparison with their experiments. The corresponding experimental parameters are $R_t = 2\,\mathrm{mm}$ , $R_p = 0.61\,\mathrm{mm}$ , $\rho _{\!f}=1.05 \times 10^3\,\mathrm{kg\,m^{-3}}$ , $\eta /\rho _{\!f}=1.14 \times 10^{-5}\,\mathrm{m^2\,s^{-1}}$ and $v_{{in}}=5.6\,\mathrm{mm\,s^{-1}}$ .

At both the inlet and outlet of the tube, we applied pressure boundary conditions applied using the anti-bounce-back scheme to modify the distribution functions as follows (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2016):

(5.11) \begin{equation} f_{\bar {i}}(\boldsymbol{x}_w, t + \Delta t) = -f_i^*(\boldsymbol{x}_w, t) + \frac{2w_i}{c_s^2} \left [ p_w + {\rho _w}{c_s^2} \left (\frac {(\boldsymbol{c}_{\!i} \boldsymbol{\cdot }\boldsymbol{u}_{\!w})^2}{2c_s^4} - \frac {u_w^2}{2c_s^2} \right )\right ]\!. \end{equation}

Here, $p_w$ , $\rho_w$ and $\boldsymbol{u}_w$ correspond to the fluid pressure, density, and velocity at the inlet/outlet boundaries, respectively. As shown in figure 5(b), our numerical results show good agreement with previous experimental (Karnis et al. Reference Karnis, Goldsmith and Mason1966) and numerical (Wen et al. Reference Wen, Zhang, Tu, Wang and Fang2014) results. Therefore, errors that propagate into the hydrodynamic force calculation due to the refilling procedure and the staircase approximation of the particle shape in the bounce-back procedure do not affect the particle dynamics significantly.

Figure 5.

Benchmarking the hydrodynamic force by investigating the migration of a spherical particle in a Poiseuille flow in a cylindrical tube (Segrè–Silberberg effect). (a) Schematic of the simulation set-up. The Poiseuille flow was generated by applying a pressure inlet/outlet boundary condition. (b) Motion of the particle perpendicular to the cylindrical axis as a function of time. Our numerical results (black line) agree well with previous experimental results by Karnis et al. (Reference Karnis, Goldsmith and Mason1966) and previous numerical results by Wen et al. (Reference Wen, Zhang, Tu, Wang and Fang2014).

5.5. Wetting boundary condition

When the interface between the two fluids contacts a solid surface (either the particle or the flat solid substrate), it forms a contact angle that depends on the surface energy of the solid. We use a geometric wetting boundary condition that allows us to independently tune the contact angle. The order parameter gradient at the fluid nodes next to a solid boundary within the diffuse interface region is modified such that it results in the prescribed contact angle (Xu et al. Reference Xu, Liu and Valocchi2017; Akai et al. Reference Akai, Bijeljic and Blunt2018). Two possible vectors, $\boldsymbol{n}_+$ and $\boldsymbol{n}_-$ are computed for the modified gradient direction,

(5.12) \begin{equation} \boldsymbol{n}_{\pm }=\left (\cos \pm \theta _p-\frac {\sin \pm \theta _p \cos \alpha ^*}{\sin \alpha ^*}\right )\boldsymbol{n}_s+\frac {\sin \pm \theta _p}{\sin \alpha^* }\boldsymbol{n}^*. \end{equation}

Here, $\boldsymbol{n}_s$ is the unit normal vector that points outwards from the surface of the solid, $\alpha ^*=\cos ^{-1}(\boldsymbol{n}_s\boldsymbol{\cdot }\boldsymbol{n}^*)$ , and $\boldsymbol{n}^*$ is the initial direction of the gradient of the order parameter calculated before applying the correction. Here $\theta _p$ is the prescribed contact angle. The vector that is closest to $\boldsymbol{n}^*$ is chosen as the gradient direction. The direction of the gradient is then updated (while keeping the magnitude constant) to obtain the prescribed contact angle.

Figure 6.

Benchmark of the wetting boundary condition between the particle and the fluids by varying the prescribed contact angle and measuring the corresponding contact angle obtained in the simulations. The snapshots show the final positions of the particle relative to the interface between the two fluids for each contact angle. The line denotes the trend that should be obtained if there is perfect agreement between the prescribed and measured angles. The data points lie close to the line, showing that the wetting boundary condition accurately models contact angles between $30^\circ$ and $150^\circ$ (inclusive). For these simulations, the domain size is $200\times 200\times 200$ and $R_p=25$ . The viscosities of the two liquids are equal (relaxation times, $\tau _a=\tau _b=1.0$ ) and density the fluid density is $\rho_{f} =0.05$ . The interfacial tension between the two fluids is $\sigma =0.02$ . The filled triangles, crosses and empty triangles correspond to when the particle was initialised above, at and below its predicted equilibrium position relative to the interface to demonstrate contact angle hysteresis. Here CM stands for centre of mass. Contact angle hysteresis arises because the fluid sees the particle as being rough due to the lattice Boltzmann grid. The image processing procedure to obtain the snapshots shown here is described in Supplementary figure 1 are available at https://doi.org/10.1017/jfm.2026.11865.

We benchmarked the wetting boundary condition by investigating the equilibrium position of a particle at a liquid–fluid interface for a range of prescribed contact angles, as shown in figure 6. For this benchmark, the viscosity of the two fluids was set to be equal. The fluid parameters were chosen such that the Ohnesorge number, which characterises the ratio between the interfacial forces and the viscous forces, was ${Oh}=\eta /\sqrt {\rho_{f} R_p\sigma }=0.075$ . Initially, the particle was positioned exactly halfway across the liquid–fluid interface. Gravity was neglected. As the simulation evolved, the capillary force shifted the particle upward (when $\theta _p\gt 90^\circ$ ) or downward (when $\theta _p\lt 90^\circ$ ) until the interface relaxed back to its equilibrium shape and there was no net capillary force on the particle. By measuring the final position reached by the centre of the particle, we can obtain the actual measured contact angle between the particle and the lower fluid and compare it with the prescribed contact angle to test the implementation of the wetting boundary condition. Based on circular geometry, the measured angle (in radians) is related to the centre of the particle according to

(5.13) \begin{equation} \theta _{\textit{mes}}=\frac {\pi }{2}-\arcsin \left (\frac {z_p-z_{\textit{int}}}{R_p}\right )\!, \end{equation}

where $z_p$ is the final vertical position of the particle’s centre of mass and $z_{\textit{int}}$ is the position of the liquid–fluid interface.

We obtain a good agreement between the prescribed contact angle and the measured contact angle (figure 6), especially for contact angles $\gt 30^\circ$ when the particle is initialised close to its predicted equilibrium position (in which case the difference between the predicted and measured angles $\approx 6\,\%$ ). These results do not depend on the viscosity ratio between the two fluids as long as the particle is allowed to reach its final steady-state position.

We also observed contact angle hysteresis (Butt et al. Reference Butt2022) when the particle was initialised away from its equilibrium position. When the particle was initialised above the interface such that it had to move downward to reach its final position, the measured angle corresponds to the advancing contact angle, as shown in figure 6 (blue filled triangles). In contrast, when the particle was initialised below the interface such that it had to move upward to reach its final position, the measured angle was lower and corresponds to the receding contact angle, as shown in figure 6 (empty red triangles). Contact angle hysteresis is expected because the particle does not appear perfectly smooth from the LBM perspective. Instead, it appears rough because the solid nodes are confined to a lattice grid. Roughness leads to contact angle hysteresis and prevents the particle from reaching the equilibrium position that it would have reached if it had been perfectly smooth.

5.6. Capillary force

When particles are at the liquid–fluid interface, they experience capillary force due to the action of surface tension at the three-phase contact line. The capillary force depends on several factors, including the particle shape and size, the interfacial tension between the two fluids, and the contact angle between the particle and the two fluids.

In principle, the capillary force can be obtained by integrating the surface tension vector along the three-phase contact line,

(5.14) \begin{equation} \boldsymbol{F}_{\!\sigma} =\oint _{\textit{CL}}\boldsymbol{\sigma }\boldsymbol{\cdot }\mathrm{d}\boldsymbol{l}, \end{equation}

where $\boldsymbol{l}$ is the line element and the integral is performed around the contact line. However, performing this integral is not straightforward in LBM. Since we use a diffuse interface model, this integral has to be modified to account for the fact that the contact line is not a sharp line but a ring with a finite thickness. Furthermore, the integral must be discretized and converted into a summation over a set of Lagrangian points spanning the surface of the particle. In this work, we use the capillary force model by Zhang et al. (Reference Zhang, Zhang, Liu and Jia2022). By taking into account the above considerations, Zhang et al. (Reference Zhang, Zhang, Liu and Jia2022) proposed that the capillary force on a particle can be computed as

(5.15) \begin{equation} \boldsymbol{F}_{\!\sigma} =\sum _i\frac {9}{2}k \beta [1-\phi (\boldsymbol{x_i})^2]\Delta A_i \sin \theta \ \boldsymbol{m}(\boldsymbol{x}_{\!i}). \end{equation}

Here, the summation runs across all Lagrangian points on the surface of the particle that lie within the diffuse interface region where $-1\lt \phi \lt 1$ , $k$ is a lattice-dependent geometric constant that links the segregation parameter $\beta$ to the equilibrium interface thickness. For the D3Q19 discrete velocity set used in our method, $k\approx 0.134$ (Riaud et al. Reference Riaud, Zhao, Wang, Cheng and Luo2014; Zhang et al. Reference Zhang, Zhang, Liu and Jia2022). The term $\Delta A_i=R_p^2\sin \alpha _i\Delta \alpha \Delta \varphi$ is the area element on the spherical particle surface at the Lagrangian point $i$ , where $\alpha _i$ is the polar angle of the Lagrangian point and $\varphi$ is the azimuthal angle. Throughout this paper, we set $\Delta \alpha =\Delta \varphi =2^\circ$ . The vector $\boldsymbol{m}(\boldsymbol{x}_{\!i})=(\boldsymbol{n}(\boldsymbol{x}_{\!i})\times \boldsymbol{n}_s(\boldsymbol{x}_{\!i}))\times \boldsymbol{n}(\boldsymbol{x}_{\!i})$ is the unit vector tangential to the liquid–fluid interface at the Lagrange point $\boldsymbol{x}_{\!i}$ .

When the capillary force is distributed unequally around the particle, for example, due to a difference in contact angle on different sides of the particle, this leads to a capillary torque (Naga et al. Reference Naga, Vollmer and Butt2021b ). This consideration is implicitly embedded in the model.

5.7. Benchmarking capillary force

Next, we focus on benchmarking the capillary force when a particle is at the liquid interface. The capillary force plays a key role in understanding how drops interact with particles on surfaces, which we will discuss in the next section.

We benchmarked the capillary force by measuring the force required to detach a particle from a liquid–fluid interface, as shown in figure 7. For this benchmark, we initialised the particle halfway across the interface and allowed it to reach its equilibrium position, which is a function of the contact angle. We then moved the particle upward at constant velocity and measured the capillary force due to the action of surface tension at the three-phase contact line (figure 7 c). When the particle was hydrophobic ( $\theta _p\gt 90^\circ$ ), it detached cleanly from the interface (figure 7 a), and no liquid residue was left on the particle. In contrast, with hydrophilic particles ( $\theta _p\lt 90^\circ$ ), the interface became unstable and ruptured, leaving a liquid residue on the particle (figure 7 b).

During the detachment, the magnitude of the capillary force increased (figure 7 c) and reached a maximum. The measured maximum force can be compared with the following analytical prediction (Scheludko & Nikolov Reference Scheludko and Nikolov1975):

(5.16) \begin{equation} F_{\textit{max}}=2\pi \sigma R_p\cos ^2\left (\frac {\theta _p}{2}\right )\!. \end{equation}

Here, this expression assumes that the contact angle maintains a constant value $\theta _p$ as the three-phase contact line slides on the particle. It also neglects the effects arising from the rupture of the capillary bridge. To test our capillary force model, we compared the maximum measured forces with the predictions given by (5.16) for a range of contact angles. Our simulation results agree well with the analytical prediction for $30^\circ \leq \theta _p \leq 150^\circ$ , as shown in figure 7(d). For small contact angles ( $\theta _p=30^\circ$ ), the error in the capillary force can be as large as 30 % for small particles with radius $10$ lattice units. The error can be reduced to within 5 % by increasing the particle radius to $25$ lattice units. Two factors contribute to this error. First, the prescribed contact angle that we use in the analytical prediction (5.16) is not exactly equal to the actual contact angle on the particle, as shown in figure 6 where discrepancies can be as large as $\approx 16\,\%$ for the worst cases. Second, our capillary force model computes the capillary force assuming that the particle is a perfect sphere, but from the LBM perspective, the particle is jagged because it must conform to the discrete lattice grid. This explains why the agreement is worse for smaller particles in figure 7(d). For small particles the assumption of a perfect sphere becomes less realistic due to the jagged edges becoming more significant. Therefore, when choosing the particle size, we must make a compromise between the computational cost required and the accuracy of the capillary force model.

Figure 7.

Benchmarking the capillary force by detaching a particle from a liquid–fluid interface. (a) Simulation snapshots show the detachment process of a hydrophobic particle ( $R_p=25$ lattice units, $\theta _p=110^\circ$ ) as the particle moves upward at a constant velocity ( $v_{p}=0.01$ ). (b) Detachment of a hydrophilic particle ( $R_p=25$ , $\theta _p=50^\circ$ ) from an interface. A small amount of liquid residue remains on the particle after the detachment. (c) Vertical component of the normalised capillary force on the particle against normalised time during the detachment. This force curve corresponds to the snapshots shown in (a) for $\theta _p=110^\circ$ . On the $x$ -axis, time is normalised by the time taken for the particle to move a distance of one particle radius, $R_p/v_{p}$ . The minimum of the force curve is given by $F_{\textit{max}}$ . (d) Here $F_{\textit{max}}$ against $\theta _p$ for contact angles between $30^\circ$ and $150^\circ$ . The results are in good agreement with the analytical prediction (dashed red line), in particular when using large particles ( $R_p=25$ ).

Figure 7. Long description

Panel A: Simulation snapshots show the detachment process of a hydrophobic particle as it moves upward at a constant velocity. The snapshots are labeled with different normalized times. Panel B: Simulation snapshots show the detachment process of a hydrophilic particle from an interface, with a small amount of liquid residue remaining on the particle after detachment. The snapshots are labeled with different normalized times. Panel C: A line graph shows the vertical component of the normalized capillary force on the particle against normalized time during the detachment process. The x-axis is labeled ‘t/(R_p/v_p)’ and the y-axis is labeled ‘F_σ/(σR_p)’. The graph shows a curve with a minimum point labeled ‘−F_max/(σR_p)’. Panel D: A scatter plot shows the relationship between the maximum force and contact angle for different particle sizes and surface tensions. The x-axis is labeled ‘Contact angle (deg.)’ and the y-axis is labeled ‘F_max/(σR_p)’. The plot includes data points for different particle sizes and surface tensions, with a dashed red line representing the analytical prediction.

5.8. Other technical considerations

In LBM, there must always be at least one fluid node separating two particles or a particle and a solid substrate for the hydrodynamic force calculated by the momentum exchange method to be accurate (Jiang et al. Reference Jiang, Liu, Chen and Tsuji2022; Zhang et al. Reference Zhang, Zhang, Liu and Jia2022). However, DEM requires a finite overlap between solids for contact forces to become active. To harmonise the differing requirements of LBM and DEM, we use a DEM particle radius (which we call the mechanical contact radius, $R_c$ ) that is larger than the LBM radius (which we call the effective hydrodynamic radius, $R_p$ ), following the approach proposed by Boutt et al. (Reference Boutt, Cook and Williams2011), Jiang et al. (Reference Jiang, Liu, Chen and Tsuji2022) and Fei et al. (Reference Fei, Qin, Wang, Huang, Wen, Zhao, Luo, Derome and Carmeliet2023). With this approach, there can still be contact from the perspective of DEM while ensuring that there is always at least one fluid node separating the solid nodes between two particles or between a particle and the solid substrate. The difference between the mechanical and hydrodynamic radii, $\delta R=R_c-R_p$ , is set to one lattice unit in our simulations. As suggested by Jiang et al. (Reference Jiang, Liu, Chen and Tsuji2022), this approach qualitatively mimics a similar effect to lubrication forces, since the contact forces start to become active when there is still a thin liquid film separating the hydrodynamic radii of the two solid surfaces. Importantly, for the particle removal scenarios that we will focus on in § 6, a quantitative treatment of lubrication forces is not necessary because there is direct contact between the particle and the surface, as shown in previous experiments with hydrophobic particles on hydrophobic surfaces (Naga Reference Naga2021).

To resolve the contact force dynamics, the integration time step $\Delta t$ must be smaller than the duration $t_c$ of a collision between a particle and a solid surface or between two particles. The duration of a collision is given by

(5.17) \begin{equation} t_c=\frac {\pi }{\omega_{c} }, \end{equation}

where $\omega_{c} =\sqrt {({k_n}/{m_p})- ( {\gamma _n}/{2m_p} )^2}$ is the angular frequency of the damped oscillations that arise when a particle collides with a solid substrate. For the collision between two particles, $m_p$ must be replaced by the reduced mass of the particles, $m_{ij}=m_im_{\!j}/(m_i+m_{\!j})$ . The term within the square root of the denominator must be positive for $t_c$ to be a real number. This imposes a constraint on the maximum damping constant that can be used for a given spring constant and particle mass, $\gamma _n\lt 2\sqrt {k_nm_p}$ . In LBM, the integration time step is typically set to unity. The integration time step for the DEM model can also be set to unity as long as $\Delta t\ll t_c$ . Throughout this paper, we use $\Delta t \approx t_c/45$ .

6.1. Simulation set-up and parameters

To investigate how a drop removes a particle from a surface, we initialised the drop as a hemisphere on a horizontal solid surface and positioned the particle on the surface in front of the drop. For all the simulations presented in this section, the domain size was 4 mm $\times$ 1.81 mm $\times$ 7.1 mm ( $400\times 181\times 71$ lattice units) and the particle radius was $R_p=100$ μm (10 lattice units). The drop radius was $R_d=5R_p$ , except in § 6.4 where we varied the radius from $2R_p$ to $6R_p$ .

To move the drop along the horizontal surface, we applied a horizontal force $f_x$ per unit volume to the drop along $x$ . On the particle, we applied a downward force $g_p$ (henceforth called the gravitational force on the particle) to account for the combined effects of gravity and intermolecular adhesion forces when the particle is in contact with the surface. Our model also accounts for the reduction in the effective gravitational force on the particle due to buoyancy when the particle enters the liquid drop. Periodic boundary conditions were applied along all three Cartesian directions.

In this paper, we systematically investigate the effects of varying the coefficients of friction between the particle and the surface, the radius of the drop, and the force applied to move the drop, expressed in terms of the dimensionless Bond number, ${Bo}=f_xR_d/\sigma$ . We focus on hydrophobic particles ( $90^\circ \leq \theta _p \leq 95^\circ$ ) on hydrophobic surfaces ( $90^\circ \leq \theta _s \leq 110^\circ$ ). The full list of parameters is provided in table 1.

Table 1.

Description and values of the parameters used when simulating the collision between a drop and a particle on a flat solid surface. The conversion factor for the length scale is chosen such that the particle radius is 100 μm, and the conversion factors for the mass and time are chosen such that the surface tension matches the surface tension between air and water (72 mN m⁻¹) and the dynamic viscosity of the drop matches that of water (1 mPa s) for some of the simulations. All other quantities can be converted from the conversion factors for length, mass and time because the SI units of these quantities are combinations of these three quantities. The viscosity ratio between the two fluids is 10. The contact angles are chosen to match the experimental parameters in Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). When comparing with experiments, the gravitational force on the particle was chosen such that the friction force on the particle has the same order of magnitude as the capillary force between the drop and the particle, as is the case in Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). Varying the damping constants $\gamma _n$ , $\gamma _t$ and $\gamma _r$ in the range indicated above did not noticeably affect the results. Therefore, when conducting the systematic investigations in figures 12, 13 and 14, they were set to zero.

Table 1. Long description

A table with 18 rows and 3 columns comparing simulation parameters for fluid dynamics experiments. The columns are labeled Property, Symbol Value in simulation units, and Value in SI units. The table lists various properties such as time unit, length unit, fluid density, surface tension, dynamic viscosity, contact angles, particle radius, gravitational force, total horizontal force, particle density, spring stiffness, rolling spring stiffness, friction coefficients, damping constants, and their corresponding values in simulation units and SI units. Each row provides specific values for these properties.

Our simulations correspond to Bond numbers in the range of $0.16\leq {Bo}\leq 1.44$ . Within this range, the capillary number ( ${Ca}=\eta _bU/\sigma$ , where $\eta _b$ is the dynamic drop viscosity, $\sigma$ is the surface tension and $U$ is the drop centre of mass velocity) is of the order of $10^{-2}$ or lower, indicating that the capillary force on the particle dominates the viscous force. Furthermore, the particle Reynolds number ( ${Re}=\rho _{\!f} R_p U/\eta _b$ , where $\rho _{\!f}$ is the liquid density) is of order $10^{-1}$ or lower, indicating that our simulations are in the laminar flow regime. For particle removal by water drops on hydrophobic surfaces, the interaction is a competition between the capillary force and the friction force (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). To characterise the interplay between the friction force and the capillary force, we introduce the dimensionless number $\mu g_p/(\sigma R_p)$ , where the numerator is the friction force that the particle experiences when sliding on a (dry) solid surface and the denominator is a measure of the capillary force on the particle. In the following, we varied this parameter between $0\leq \mu g_p/(\sigma R_p)\leq 26$ .

6.2. Push–pull versus enter–exit scenarios

To test whether our numerical model captures the key interactions necessary to study particle removal by drops, we simulated the collision between a drop and a particle on a surface and compared our results with previous experiments (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). The experiments were performed with water drops and spherical hydrophobic glass particles ( $R_p\approx 115\,{\unicode{x03BC}} \mathrm{m}$ ) on a hydrophobic polydimethylsiloxane surface, as described in Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). Briefly, the collision between the drop and the particle was imaged using an inverted laser scanning confocal microscope. The drop was fixed within the field of view of the microscope using a metallic cantilever while the surface moved at controlled speeds. With this set-up, it is possible to image the collision between the drop and the particle over distances of several centimetres, which is much larger than the millimetric field of view of the microscope. The experimental imaging provided horizontal slices of the collision, as shown in figure 8(a). Furthermore, the horizontal force $F_x$ acting on the drop and the particle along $x$ could also be measured directly from the deflection $\Delta x$ of the cantilever and using Hooke’s law, $F_x=k_{s}\Delta x$ , where $k_{s}$ is the spring constant of the cantilever. The contact angle between the drop and the surface was between $80^\circ$ and $120^\circ$ , and the contact angle between the drop and the particle was between $75^\circ$ and $100^\circ$ . The range spanned by these angles is due to contact angle hysteresis caused by roughness, chemical heterogeneities and surface adaptation effects. For simplicity, our numerical model does not include static contact angle hysteresis. However, it still captures dynamic contact angle hysteresis.

Figure 8.

Push–pull scenario. (a) Experiments showing the push–pull scenario when a water drop collides with a hydrophobic particle on a hydrophobic surface at 50 ${\unicode{x03BC}} \mathrm{m\,s^{-1}}$ . The snapshots show the bottom view of the collision in a horizontal plane through the centre of the particle. Reproduced from Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ) under a CC BY 3.0 licence. (b) Simulations of a drop colliding with a particle with low friction, as seen in a horizontal plane through the centre of the particle. The velocity field is overlaid on top, with the arrows corresponding to velocity in the drop’s centre-of-mass frame. Green line in $t_3$ indicates where two images have been stitched together when the drop is at the periodic boundary. (c) Experimental force (along $x$ ) curve for the force on the particle. Data from Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). (d) Force that drop exerts on particle (i.e. capillary + hydrodynamic forces), obtained with our LBM–DEM method. Different colours show force components in different directions. When the particle is frictionless, the forces along the $x$ and $y$ directions are zero (dashed lines) because there is no resistance to lateral motion. The drop still exerts a vertical force on the particle due to the asymmetric shape of the liquid meniscus around the particle. In Supplementary figure 2, we split the forces into contributions due to the hydrodynamic force and the capillary force.

Experimentally, two scenarios were reported for the collision between a water drop and a hydrophobic particle. In the first scenario, which we call the push–pull scenario, the particle moved around the drop’s footprint and remained attached to the rear drop–air interface (figure 8 a). In the second scenario (enter–exit scenario), the particle entered and exited the drop (figure 9 a). Both the push–pull and enter–exit scenarios were imaged directly using laser scanning confocal microscopy (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). However, due to the limited temporal resolution of this technique, the imaging was restricted to a single 2-D horizontal plane. By tuning the friction coefficients, we could successfully reproduce these two scenarios using our LBM–DEM method (figures 8 b and 9 b). Our numerical method also allows us to visualise the interaction in three dimensions for the first time (figure 10, Supplementary movies 1 and 2).

Figure 9.

Enter–exit scenario. (a) Experiments showing the enter–exit scenario when a water drop collides with a hydrophobic particle on a hydrophobic surface at 500 ${\unicode{x03BC}} \mathrm{m\,s^{-1}}$ . The snapshots show the bottom view of the collision in a horizontal plane through the centre of the particle. Reproduced from Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ) under a CC BY 3.0 licence. (b) Simulations of a drop colliding with a particle with high friction, as seen in a horizontal plane through the centre of the particle. The velocity field is overlaid on top, with the arrows corresponding to velocity in the drop’s centre of mass frame. In both experiments and simulations, the particle enters and exits the drop. (c) Experimental force curve for the force (along $x$ ) acting on the particle during the collision. Data from Naga et al. (Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). (d) Force that drop exerts on particle (i.e. capillary + hydrodynamic forces), obtained with our numerical method. The force components in all three Cartesian directions can be obtained, as shown by the different colours. In both experiments and simulations, the force is maximum when the drop crosses the drop–air interface, highlighting that capillary forces dominate. Here $F_z$ remains negative after the particle exits because a small amount of liquid is left behind with the particle, leading to an attractive capillary force between the particle and the surface.

Figure 9. Long description

Panel A: A series of four images showing the bottom view of a water drop colliding with a hydrophobic particle on a hydrophobic surface at 500. The snapshots are taken in a horizontal plane through the center of the particle. The images are labeled t0, t1, t2, and t3, showing the progression of the collision. Panel B: Simulations of a drop colliding with a particle with high friction, viewed in a horizontal plane through the center of the particle. The velocity field is overlaid on top, with arrows corresponding to velocity in the drop’s center of mass frame. The images are labeled t0, t1, t2, and t3, showing the progression of the collision. Panel C: A line graph showing the experimental force curve for the sum of the force acting on the drop and particle during the collision. The x-axis is labeled Normalized time, t/(Rd/U) and the y-axis is labeled F/(σRp). The graph shows two peaks, indicating the force is maximum when the drop crosses the drop-air interface. Panel D: A line graph showing the force that the drop exerts on the particle, obtained with a numerical method. The x-axis is labeled Normalized time, t/(Rd/U) and the y-axis is labeled F/(σRp). The graph shows three lines representing the force components in all three Cartesian directions, labeled Fx, Fy, and Fz. The force is maximum when the drop crosses the drop-air interface, highlighting that capillary forces dominate.

Figure 10.

Two types of collision scenarios between a drop and a particle on a horizontal surface. (a) Push–pull scenario. When the friction between the particle and the surface is small compared with the maximum capillary force that the drop can exert on the particle, the particle moves along the circumference of the drop and remains attached to the rear. In the example shown here, the particle moves clockwise around the drop because it is initially positioned slightly to the right-hand side of the drop. When there is no initial offset, the drop pushes the particle, and the latter remains in front of the drop for much longer. The insets show that the particle rolls when it is pulled by the drop. Frictionless particles roll anticlockwise, following the flow in the drop (arrows show fluid flow as viewed in the drop’s centre of mass frame). Frictional particles roll in the opposite direction (clockwise) due to friction between the particle and the surface. (b) Enter–exit scenario. When the friction force exceeds the maximum capillary force, the particle penetrates the drop, travels across its base and detaches from the drop at the rear side.

To simulate the two scenarios, we positioned a spherical particle on the substrate in front of the drop. The contact angle between the drop and the particle and between the drop and the surface were both set to $90^\circ$ . We then applied a constant force to the drop to move it towards the particle. We chose the applied force, drop viscosity and surface tension such that the capillary force between the drop and the particle dominates the hydrodynamic force, as is typically the case in experiments using water drops (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). This regime corresponds to when the capillary number of the drop is significantly less than unity. In these simulations, the capillary number was of the order of ${Ca}\approx 10^{-2}$ . Our numerical method successfully reproduced the push–pull scenario in the limit where the friction coefficients $\mu$ and $\mu _r$ were both small, such that the capillary force between the drop and the particle exceeds the friction force between the particle and the surface. The enter–exit scenario was reproduced in the opposite limit when the friction force exceeds the capillary force.

In the push–pull scenario, the particle goes around the drop because any slight offset in the initial alignment between the particle and the drop relative to the direction of motion causes the particle to move sideways due to the convex shape of the front of the drop. Once at the rear of the drop, the particle remains attached to the drop–air interface. The particle always ends up in the same final configuration at the rear, regardless of the initial offset between the drop and the particle.

To compare the experiments and simulations quantitatively, we analysed the forces acting on the particle. In the experiments, the force acting on the particle alone was obtained by subtracting the total measured force (on drop $+$ particle) from the force acting on the drop alone before it collided with the particle. In the simulations, the force that the drop exerts on the particle can be obtained directly. While current experimental methods only provide the force acting along the direction of motion ( $x$ ), our numerical method provides the force in all three Cartesian directions ( $x$ , $y$ , $z$ , as defined in figure 8 a time $t_0$ ).

The experimental and numerical force curves show good agreement (figures 8 c and 8 d, respectively), indicating that the numerical model captures all the key physical interactions. The force curve along $x$ is characterised by two broad plateaus, corresponding to when the particle is in front and behind the drop. The first plateau is due to the force exerted by the drop as it pushes the particle, overcoming friction forces. The width of the first plateau depends on the amount of time that the particle spends in front of the drop before moving to the side. When there is an initial offset in the alignment of the drop and the particle (as is the case in figures 8 b and 8 d), the force in the $y$ direction is non-zero (orange curve in figure 8 d). This non-zero force causes the particle to move around the drop perimeter. When the particle is at the side (time $t_2$ ), the $x$ component of the force decreases to zero since at that instant, the drop does not push the particle in the forward direction. The second plateau in the horizontal force corresponds to when the particle is pulled at the rear of the drop. Interestingly, the drop also exerts a vertically downward force when pulling the particle (green curve in figure 8 d).

Macroscopically, the push–pull scenario could be reproduced even without accounting for friction. When both $\mu$ and $\mu _r$ were set to zero, the scenario looked similar to the one shown in figure 8(a,b). However, despite the macroscopic similarity, there are significant differences in the force curves (dashed lines in figure 8 d show force curves for a frictionless particle). Since no force is required to push or pull a frictionless particle, the forces along the $x$ and $y$ directions were always zero (dashed blue and orange lines). The absence of a horizontal force also meant that the drop did not deform as much when pushing or pulling a frictionless particle. However, a small vertical force was still present (dashed green line) due to the asymmetric shape of the water meniscus around the particle. Therefore, although ignoring friction leads to a similar macroscopic outcome, it is crucial to include friction to reproduce the forces on the particle accurately.

An important difference between the simulations and experiments is the presence of random roughness on the surface and on the particle in experiments. Roughness causes the particle to roll unevenly on the surface, which leads to fluctuations in the plateaus in the force curve when the drop pushes and pulls the particle, as seen in figure 8(c). Roughness also causes the particle to move around the drop perimeter even when the particle is initially positioned in perfect alignment with the drop. Random fluctuations in the particle motion due to roughness cause the particle to break the left–right symmetry and go around the drop. In contrast, in the simulations, the particle can be pushed in front of the drop until the left–right symmetry is broken by numerical noise (Connington & Lee Reference Connington and Lee2012) in the particle position. For the comparison with experiments in figure 8(d), we introduced a small misalignment by positioning the particle slightly to the right of the drop, to favour clockwise motion as it moves around the drop perimeter.

For the enter–exit scenario, the force curve (for force along $x$ ) is characterised by two peaks, shown in figure 9(c,d). The maxima of these peaks correspond to when the particle crosses the front and rear drop–air interface. The time interval between the maximum of the first peak and the beginning of the second peak corresponds to the time taken for the particle to travel across the drop footprint, $\Delta t\approx 2L/U$ , where $L$ is the length of the drop along the direction of motion and $U$ is the relative velocity between the drop and the particle.

When the particle moved across the footprint of the drop, the force reduced significantly compared with when the particle crossed the interface. This reduction in force is because when the particle is fully submerged, the capillary force is zero, and the particle only experiences hydrodynamic force due to viscous drag. For slow speeds, as was the case in the experiments where the speed was $O(10^{-4}\,\mathrm{m\,s^{-1}})$ , the viscous drag force is small (estimated to be $O(10^{-9}\,\mathrm{N})$ (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a )) compared with the capillary force when the particle is at the liquid–air interface. Consequently, the force reduced to $\approx 0$ between the two peaks in the experimental force curve. In the simulations, the reduction in force was less significant because we used a higher drop speed $O(1\,\mathrm{m\,s^{-1}})$ to keep computational costs tractable. Due to the higher drop speed in the simulations, the force along $x$ does not decrease completely to zero when the particle moves inside the drop. Note that, despite the four orders of magnitude difference in the absolute drop speed, the Reynolds number of the particle only differs by one order of magnitude in the experiments ( ${Re}\approx 10^{-2}$ ) and simulations ( ${Re}\approx 10^{-1}$ ) because the difference in speed is compensated by having a kinematic drop viscosity $\eta _b/\rho _{\!f}$ that is around three orders of magnitude higher in simulations. Crucially, in both experiments and simulations, the same hierarchy of forces applies, with the capillary force being typically larger than the viscous force ( ${Ca}\ll 1$ ), which is itself larger than the inertial force ( ${Re}\lt 1$ ).

Our simulations further indicate that the drop exerted a vertical force on the particle when the liquid–air interface moved across the particle. Figure 9(d) shows that the vertical force was negative throughout the interaction, except when the particle first reached the rear drop–air interface due to dewetting of the liquid on the particle. The vertical force also remained negative even after the particle exited the drop (time $t_3$ in figure 9 d) due to the presence of a small amount of liquid being left under the particle after the detachment. This remnant liquid leads to capillary adhesion between the particle and the surface.

In addition to enabling 3-D visualisation and providing force measurement in all directions, another benefit of the numerical method is that it allows us to obtain flow profiles. In both the push–pull and enter–exit scenarios, the flow in the drop followed a clockwise rolling motion when viewed in the centre of mass frame of the drop, as shown in figure 11(a,b). Supplementary movies 3, 4 and 5, 6 show flow profiles in vertical and horizontal slices across the centre of the particle, respectively. Interestingly, the flow direction in the drop is opposite to the direction in which friction causes the particle to rotate. When the drop pulls the particle, flow promotes anticlockwise rotation of the particle, whereas friction between the particle and the surface promotes clockwise rotation, as shown in the inset at the bottom of figure 10(a). Consequently, frictional particles roll clockwise, driven by the friction force between the particle and the surface, whereas frictionless particles roll anticlockwise, driven by the fluid flow in the drop. This effect reinforces the conclusion that even when friction is small, it cannot be ignored since it plays a crucial role in determining the microscopic motion of the particle.

Figure 11.

Velocity profile in the drop and surrounding phase for the push–pull (a) and enter–exit (b) scenarios. These velocity profiles are taken in a vertical slice close to the centre of mass of the drop. The velocities shown are in the centre of mass frame of the drop. In (a), there is a small but finite friction force. In (b), there is a large friction force. A satellite drop becomes entrained by the particle due to the breakup of the capillary bridge when it exits the drop.

Figure 12.

(a) Outcomes of the collision between a drop and a particle as $\mu$ and $\mu _r$ are varied independently. Here, ${Bo}=1.25$ and $0\leq \mu g_p/(\sigma R_p)\leq 5.2$ . The green circles correspond to the case where the particle remains attached to the drop, as shown in figure 10(a). The red squares correspond to when the particle detaches from the drop. (b) Variation of the rotational velocity of the particle with the friction coefficient. The rotational velocity is normalised by the centre of mass velocity of the particle. A ratio of 1 indicates pure rolling, whereas a ratio of 0 indicates pure sliding. Even when the rolling friction is set to zero ( $\mu _r=0$ ), pure rolling is not the preferred mode of motion, and there remains some relative sliding between the particle and the surface up to $\mu =0.8$ . When $\mu _r=\mu$ , the particle slides rather than rolls. The rotational velocity can be negative when $\mu =0$ , indicating that the particle rotates anticlockwise. This anticlockwise rotation is driven by the anticlockwise flow of liquid inside the drop.

6.3. Transition from push–pull to enter–exit scenario

We showed that the friction force between the particle and the surface determines whether the collision follows the push–pull or enter–exit scenarios. Yet, there are several other factors that may influence these scenarios. First, there are two types of friction that enter this problem: sliding friction and rolling friction. How do these two types of friction influence the collision scenario? Furthermore, even when the rolling and sliding friction coefficients are fixed, the collision scenario may depend on the drop size or the tilt angle of the surface. In the following, we investigate the effects of varying the coefficients of sliding and rolling friction, the drop size and the force applied to the drop.

To investigate the role of sliding and rolling friction, we varied the coefficients of sliding and rolling friction ( $\mu$ and $\mu _r$ , respectively) independently between 0 and 1 and mapped out whether the particle remained attached to the drop (figure 12 a). In this set of simulations, we kept everything else constant, including the contact angles, particle weight and force applied to the drop. It turns out that the particle could always be captured when at least one of the friction coefficients was zero (figure 12 a). However, this was no longer the case when the friction coefficients were finite.

As shown in figure 12(a), the outcome depends on both the sliding and rolling friction coefficients. However, for successful particle capture, it is more important for the sliding friction to be small because the distribution of outcomes is not symmetric about the line $\mu _r=\mu$ . A higher number of successful capture cases lies above the $\mu _r=\mu$ line, indicating that having a high $\mu$ is more detrimental to capture than having a high $\mu _r$ . For example, when $\mu =1.0$ , the particle is never captured for any $\mu _r\gt 0$ . However, when $\mu _r=1.0$ , the particle can still be captured as long as $\mu \lt 0.6$ . Therefore, although both sliding and rolling friction are important to predict whether a particle can be captured, the sliding friction plays a larger role in determining the outcome.

To understand when the particle rolls or slides on the surface, we analysed the angular velocity of the particle for the cases when it remained attached at the rear of the drop (figure 12 b). When both the sliding and rolling friction were zero, the angular velocity of the particle was slightly negative, indicating an anticlockwise rotation in the direction of the flow in the drop, as discussed earlier in figure 10. However, as soon as the sliding friction was non-zero, the angular velocity switched sign, indicating a clockwise rotation driven by sliding friction between the particle and the surface. When the sliding and rolling friction coefficients were equal, sliding was the preferred mode of motion. In this case, the preference for sliding can be understood to be because there are additional sources of rolling resistance that enhance the rolling friction beyond the sliding friction, including the rolling resistance experienced by the particle when it rotates against the opposing liquid flow inside the drop and due to the resistive capillary torque (Naga et al. Reference Naga, Vollmer and Butt2021b ) associated with the particle rotating at a liquid–air interface. Increasing $\mu$ beyond $\mu _r$ reduced the relative amount of sliding between the particle and the surface, until the particle started to display a pure rolling motion. Pure rolling corresponds to when the rotational velocity and translational velocity become equal, $\omega_{p} R_p=v_p$ , such that there is no relative sliding between the particle and the surface at the point of contact. For a constant sliding friction coefficient, the amount of rolling relative to sliding decreases with an increase in the rolling friction coefficient (the orange compared with the blue curve in figure 12 b).

Experimentally, it was observed that a spherical hydrophobic particle rolls on a hydrophobic surface without any noticeable sliding when it is pulled by a water drop (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ). In experiments, rolling was deduced by tracking the position of defects on the particle. Based on the findings shown in figure 12(b), our numerical results reveal that, in general, the amount of rolling relative to sliding depends on the relative size of the coefficients of rolling and sliding friction. However, since in experiments, the rolling friction coefficient for spherical particles is typically much smaller than the sliding friction coefficient, spherical particles are much more likely to roll.

6.4. Effects of drop size and speed

In practical scenarios and cleaning applications, we are likely to encounter drops of different sizes moving at different speeds. Thus, it is important to understand how the drop size and speed affect its ability to capture a particle. To understand how the removal mechanism is influenced by the drop size, we varied the drop size while keeping the particle size constant (figure 13). In all cases, we applied the same force per unit volume to move the drop. Therefore, the total force applied to the drop was proportional to its volume. This procedure mimics drops of different volumes sliding down a vertical surface under the action of gravity.

The drop size influences its ability to capture as well as transport the particle (figure 13). When the drop was too small ( $R_d/R_p=2$ , figure 13 a), it was unable to displace the particle because the total force exerted by the drop was insufficient to overcome the friction force on the particle. As the drop volume increased ( $R_d/R_p=3,4$ , figure 13 b,c), the drop successfully displaced the particle by pushing it forward. Here, within the simulated time scale, the particle remained in front of the drop rather than going around (as shown previously in figure 8) because it was positioned in perfect alignment with the drop. Therefore, a larger drop improves removal efficiency because it moves faster and can transport the particle over a larger distance within the same time interval.

Figure 13.

(a–d) Time-lapse showing the effect of increasing the drop radius $R_d$ relative to the particle radius, $R_p=100\,{\unicode{x03BC}} \mathrm{m}$ . In all cases, $\mu g_p/(\sigma R_p)=5.23$ , $\mu =\mu _r$ , and the drop Bond number varies between $0.16 \leq {Bo}\leq 1.44$ . The force applied per unit drop volume is the same in all cases, therefore changes in ${Bo}$ are due to changes in drop radius. Translucent images indicate earlier times. The particle was positioned in perfect alignment with the drop. Therefore, within the simulated time scale, it remained in front of the drop and did not go around. In (d), the particle enters and exits the drop. In this case, satellite droplets are left behind when the particle detaches from the main drop. (e) Position of the particle as a function of time for the cases shown in (a–d). Time is normalised using the initial drop radius and the characteristic drop speed corresponding to the enter–exit case when $R_d/R_p=6$ . When $R_d/R_p=2$ , the force exerted by the drop is insufficient to push the particle. For the largest drop ( $R_d/R_p=6$ ), the particle only moves when it crosses the front and rear drop–air interface. There are four steps in the curve because the drop collides with the particle a second time after crossing the periodic boundary of the simulation domain.

However, drops that are too large are no longer efficient at transporting particles. As shown in figure 13(d), when $R_d/R_p=6$ , the particle entered and exited the drop. When the particle enters and exits the drop, it can no longer be displaced efficiently because it only moves during the brief periods when it interacts with the liquid–air interface (figure 13 e). This transition to the enter–exit scenario can be understood because larger drops move faster, causing the dynamic advancing and receding contact angles at the front and rear of the drop to change, as seen by the greater drop deformation in figure 13(d) compared with figure 13(c). Changes in the contact angles influence the capillary force between the drop and the particle (Naga Reference Naga2021) and can therefore influence the outcome of the collision. Therefore, particles can only be efficiently transported by drops that are neither too large nor too small, highlighting that there is an optimal drop size. In general, the optimal drop size will depend on the friction coefficients, contact angles and the surface tension of the liquid. All these quantities must be taken into account when optimising the drop size to enable efficient particle capture and removal on surfaces.

Figure 14.

Varying the force per unit volume applied to the drop. For all cases here, the drop volume is constant, $\mu g_p/(\sigma R_p)=5.23$ , $\mu =\mu _r$ , and the particle is perfectly aligned with the drop. (a) Force acting on the particle along $x$ for different drop speeds. For the lowest three ${Bo}$ , the drop pushes the particle, and therefore the force curves have a single plateau, corresponding to the first plateau in the push–pull scenario. For the highest two ${Bo}$ , the particle enters and exits the drop. The two peaks in the force curves correspond to when the particle crosses the front and rear liquid–air interface, respectively. (b) Vertical force (along $z$ ) acting on the particle corresponding to the five cases in (a). The drop exerts a downward force on the particle when the particle crosses the front and rear liquid–air interface. In both (a) and (b), the time is normalised by $R_d/U$ , where $U$ is the approximate drop speed for the two enter–exit scenarios shown here.

Figure 14. Long description

Panel A: A line graph shows the force acting on the particle along the x-axis for different drop speeds. The x-axis represents normalized time, t/(R_d/U), and the y-axis represents the force, F_x/(σR_p). The graph includes multiple lines corresponding to different Bond numbers (Bo), with labels indicating Bo values of 1.25, 1.00, 0.75, 0.50, and 0.25. The lines show different force behaviors, including push and enter-exit scenarios. Panel B: Another line graph shows the vertical force (along the z-axis) acting on the particle for the same cases as in Panel A. The x-axis represents normalized time, t/(R_d/U), and the y-axis represents the force, F_z/(σR_p). The graph includes labels for push and enter-exit scenarios, with lines corresponding to different Bond numbers. The vertical force shows peaks when the particle crosses the front and rear liquid-air interface.

Besides optimising the drop size, it is also important to consider how the tilt angle of the surface affects the collision scenario for a fixed drop volume. Varying the tilt angle changes the force acting on the drop parallel to the surface. Therefore, to understand the effect of varying the tilt angle, we varied the force applied to the drop parallel to the surface. In figure 14, we compare five cases corresponding to different forces applied to the drop, where the applied force is indicated by the Bond number. For these simulations, the drop and particle were positioned in perfect alignment. Therefore, when the particle did not enter the drop, it remained in front of the drop within the simulated time scale. As the applied force was increased, the collision scenario transitioned from the particle being pushed by the drop to it entering and exiting the drop. For the lowest three applied forces, the particle was pushed by the drop, whereas for the highest two forces, it entered and exited the drop. Comparing the force curves for the first three and last two cases separately allows us to understand how changes in the force applied to the drop influence the force on the particle for the two different scenarios.

When the particle remained in front of the drop, a higher force was required to push the particle as the force applied to the drop increased (figure 14 a). The increase in the horizontal force required to push the particle was correlated with an increase in the downward vertical force acting on the particle (figure 14 b). Since the friction force experienced by the particle is proportional to the normal reaction force between the particle and the surface, the increase in the downward force causes the friction to increase. As a consequence, when the drop exerts a downward force on the particle (e.g. due to a change in the contact angle or in the flow profile), it also needs to exert a larger horizontal force to push the particle.

As the applied force increased, the horizontal force required to push the particle increased until the drop could no longer push the particle. For the highest two applied forces ( ${Bo}=1.00$ and $1.25$ ), the particle entered and exited the drop. The horizontal force on the particle as it entered the drop was almost the same for these two cases (figure 14 a, first peak of the top two curves). Once the particle had fully crossed the front drop–air interface, it experienced hydrodynamic drag before reaching the rear liquid–air interface. When the particle reached the rear drop–air interface, the horizontal force initially decreased as it snapped into contact with the interface before increasing again as it detached from the interface. When the particle exited the drop, the horizontal force on the particle was almost the same for ${Bo}=1.00$ and $1.25$ . Overall, these results indicate that the tilt angle of the surface can affect particle capture because the capillary force exerted by the drop–air interface is limited by a threshold value determined by the surface tension and contact angles (Naga et al. Reference Naga, Kaltbeitzel, Wong, Hauer, Butt and Vollmer2021a ), and cannot increase indefinitely to overcome friction.