2.1. Electrostatics

At the micro-particle scale, physicochemical forces become significant relative to hydrodynamic forces. Derjaguin–Landau–Verwey–Overbeek (DLVO) theory (Derjaguin & Landau Reference Derjaguin and Landau1941; Verwey & Overbeek Reference Verwey and Overbeek1948) is applied ubiquitously to describe these forces for solid objects interacting in a fluid medium, with its physical validity down to separation distances approximately $0.1\times 10^{-9}\ {\rm m}$ proven extensively (Israelachvili Reference Israelachvili2011b). DLVO theory states that the total electrostatic potential, $V_{e}$, is the superposition of the attractive London–van der Waals (LvdW) potential, $V_{LvdW}$, the repulsive electric double layer (EDL) potential, $V_{EDL}$, and the repulsive Born potential, $V_{BR}$:

(2.1)\begin{equation} V_{e} = V_{LvdW} + V_{EDL} + V_{BR}. \end{equation}

In the present work, these potentials are resolved fully according to their analytical expressions, as opposed to a simplified cohesion force model (Vowinckel et al. Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019), and are incorporated directly into the DEM as forces, as outlined in § 2.2. The primary advantage of this direct numerical simulation (DNS) approach over a simplified model is that it remains physically valid down to the physical limit of Brownian motion (i.e. particle sizes approximately ${1}\ {\mathrm {\mu }}{\rm m}$, depending on flow conditions). The remainder of this section presents the analytical DLVO forces and parameters employed in this work, and then analyses the physical attributes of these DLVO systems.

The LvdW potential arises when atomic dipole forces, which are weak individually, interact simultaneously between two macroscopic bodies. In a typical system where two different materials interact in a fluid medium, the interaction strength is quantified by the system Hamaker constant, $A_{132}$, which is calculated from the individual vacuum Hamaker constants of the interacting materials (Khilar & Fogler Reference Khilar and Fogler1998). In this work, subscript 1 denotes the particles, subscript 2 denotes the channel surface, and subscript 3 denotes the fluid. Accordingly, the LvdW interaction force between a sphere (particle) of diameter, $d$, and a flat wall, separated by a distance, $\zeta$, depends on $A_{132}$ (Hamaker Reference Hamaker1937),

(2.2)\begin{equation} F_{LvdW,p-w} ={-}\frac{A_{132}}{12}\,\frac{d^3}{(\zeta+d)^2\zeta^2}, \end{equation}

where the negative sign signifies attraction in keeping with convention. Similarly, the interaction force between two spheres (particles) is described by $A_{131}$,

(2.3)\begin{equation} F_{LvdW,p-p} ={-}\frac{A_{131}}{6}\,\frac{d^6}{\left((d+\zeta)^2-d^2\right)^2(d+\zeta)^3}. \end{equation}

It is noted that (2.3) is obtained by differentiating the original expression for potential $V_{LvdW}$ (Hamaker Reference Hamaker1937) with respect to $\zeta$, while (2.2) is obtained by additionally taking the limit of $V_{LvdW}$ as the diameter of one sphere approaches infinity. In this work, $A_{132}$ and $A_{131}$ are the only electrostatic parameters that are varied. Four distinct values are used in this work (presented in table 1), which represent a wide range of materials (Ackler, French & Chiang Reference Ackler, French and Chiang1996; Israelachvili Reference Israelachvili2011a).

Table 1. The electrostatic parameters and their values employed in this study, with resulting maximum adhesion forces $F_{a,p-w}$ and $F_{a,p-p}$. The physical values presented here are intended to provide an understanding of how $F_{a,p-w}$ and $F_{a,p-p}$ originate; however, all results are subsequently presented in terms of the non-dimensional parameters $Ad_w$ and $Ad_p$ (see (2.6) and (2.7)).

Two separate potentials, $V_{EDL}$ and $V_{BR}$, both of which arise due to the dissociation of ions within a fluid, act to oppose the attractive LvdW potential. Here, their origins and force equations are summarised. However, for an in-depth description of the electrostatic potentials that arise in liquids, the interested reader is referred to Israelachvili (Reference Israelachvili2011b). The surface of a solid object becomes negatively charged via a chemical process when submerged in a liquid. The surface then attracts a cloud of dissociated counter-ions, referred to as the EDL, some of which are bound to the surface in what is known as the Stern layer.

When the ion clouds of two surfaces overlap, they mutually repulse, which is the origin of the EDL potential. The decay of this repulsive energy with $\zeta$ is governed chiefly by the characteristic length of the EDL, known as the inverse Debye length, $\kappa$, and scales with the absolute fluid permittivity, $\varepsilon$, surface potential, $\psi$, and $d$. In this work, a simplification of the EDL force when $\psi \lesssim {25}\ {\rm mV}$ is employed for two identical spheres,

(2.4)\begin{equation} F_{EDL,p-p} = {\rm \pi}\,{d} \varepsilon \kappa \psi^2\,{\rm e}^{-\kappa \zeta}, \end{equation}

and a sphere and a flat wall,

(2.5)\begin{equation} F_{EDL,p-w} = 2{\rm \pi} \,{d} \varepsilon \kappa \psi^2\,{\rm e}^{-\kappa \zeta}, \end{equation}

noting that these lose accuracy for $\zeta < 1/\kappa$ but are useful for efficiently describing a wide variety of situations (Israelachvili Reference Israelachvili2011b). The Debye length depends on the type of electrolyte and the molar concentration of dissolved salt, $C_m$. For NaCl solutions, $\kappa = 0.73\times 10^{8}\ {\rm m}^{-1}\,\sqrt {1\ {\rm m}^{3}\ {\rm mol}^{-1}\times C_mz^2}$ and $z=1$. All EDL parameters used in this work are held constant at the values indicated in table 1.

The short-range Born repulsion arises due to the atomic repulsion of the surface ions residing in two objects’ Stern layers. This short-range atomic repulsion prevents $F_{LvdW}$ from reaching infinity, instead creating the primary attractive force minimum that is characteristic of electrostatic systems under most conditions. The hard-sphere approximation assumes that the electrostatic force is immediately cut off at some minimum atomic separation and is a common choice for attenuating $F_{LvdW}$ in numerical simulations (Mitchell & Leonardi Reference Mitchell and Leonardi2016; Lohaus, Perez & Wessling Reference Lohaus, Perez and Wessling2018; Vowinckel et al. Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019). A more realistic approximation applies the assumption of pairwise additivity to integrate the repulsive potential of two macroscopic bodies, similar to Hamaker's derivation of the LvdW potential (Feke et al. Reference Feke, Prabhu, Adin Mann and Adin Mann1984). This results in a gradual increase of the Born repulsion with decreasing separation distance. The complete expression for the Born potential obtained from pairwise integration can be found in Feke et al. (Reference Feke, Prabhu, Adin Mann and Adin Mann1984), but is not included here for brevity. In this work, the approach of Feke et al. (Reference Feke, Prabhu, Adin Mann and Adin Mann1984) is applied with a moderate hardness scaling $F_{BR} \propto \zeta ^{-6}$ and Stern layer thickness $\sigma =0.5\times 10^{-9}\ {\rm m}$.

A difficulty arises when quantifying DLVO forces for physical systems, namely, that the origin of the EDL force is at the surfaces of the counter-ions in the Stern layer, whereas the LvdW force is measured from the centre of the co-ions within the surface of the object (Israelachvili Reference Israelachvili2011b). The resulting force profiles are highly sensitive to these origins, and shifts of just $0.2\times 10^{-9}\ {\rm m}$ in the EDL force can completely eliminate the primary minimum. As discussed in § 2.2, the macroscopic frictional response of the system is highly sensitive to the normal and tangential particle contact forces. Consequently, clogging results depend strongly on the way that electrostatics is treated at contact. In this work, electrostatics at contact is treated by selecting a separation distance $\zeta /\sigma =0.6$ from which the electrostatic force is decreased linearly to $F_e=0$ at $\zeta =0$, and maintained at $F_e=0$ for $\zeta \leqslant 0$. This distance is within the primary minimum and is commensurate with minimum contact distances reported in the literature (Visser Reference Visser1995; Israelachvili Reference Israelachvili2011a). In this way, there is no electrostatic force acting when two computational surfaces overlap, which helps to resolve soft-sphere collision accurately in the DEM outlined in § 2.2. Quantifying the uncertainty of the macroscopic friction to this treatment of electrostatics, as well as other contributing factors, is discussed in detail in § 2.2.

The variation of the net electrostatic force, $F_e = F_{LvdW} + F_{EDL} + F_{BR}$, with respect to separation distance determines whether conditions are favourable or unfavourable for particle–particle and particle–wall attachment. As outlined in § 1, particle–wall attachment in particular has a key influence on clogging (Sendekie & Bacchin Reference Sendekie and Bacchin2016; Kang et al. Reference Kang, Kim and Kim2018) and the behaviour of electrostatic systems in general (Elimelech et al. Reference Elimelech, Gregory, Jia and Williams1995; Khilar & Fogler Reference Khilar and Fogler1998; Israelachvili Reference Israelachvili2011b; Chequer et al. Reference Chequer, Carageorgos, Naby, Hussaini, Lee and Bedrikovetsky2021). One of the primary aims of this work is to therefore determine how particle–wall attachment, and to a lesser degree particle–particle attachment, affect planar channel clogging. The values that influence these primarily are, respectively, the maximum particle–wall adhesion force, $F_{a,{p-w}}$, and the maximum particle–particle adhesion force, $F_{a,{p-p}}$. To quantify the strength of these adhesion forces in relation to the hydrodynamic drag force, two non-dimensional adhesion numbers are defined, one for particle–wall adhesion,

(2.6)\begin{equation} Ad_w=\frac{F_{a,{p-w}}}{3{\rm \pi}\rho_f\nu d \langle v \rangle}, \end{equation}

and the other for particle–particle adhesion,

(2.7)\begin{equation} Ad_p=\frac{F_{a,{p-p}}}{3{\rm \pi}\rho_f\nu d \langle v \rangle}, \end{equation}

where $\rho _f$ and $\nu$ are the fluid density and kinematic viscosity, respectively, and $\langle v \rangle = Gw^2/(12\rho _f\nu )$ is the superficial mean channel velocity (based on momentum conservation of a pure fluid in a planar channel, where $G$ is the applied fluid pressure gradient).

Figure 2 plots the force profiles for the two particle–particle Hamaker constants and four particle–wall Hamaker constants employed in this work, for a representative physical particle size $d={5}\ {\mathrm {\mu }}{\rm m}$. Noting that $F_e$ scales linearly with $d$, the range of physical particle sizes and Hamaker constants in this work results in a range of $F_{a,{p-w}}$ and $F_{a,{p-p}}$ values, which are presented in table 1. The physical values presented in this section are intended to provide the reader with an understanding of where the electrostatic forces originate; however, it is emphasised that the results in § 4 are presented in terms of only the non-dimensional quantities $Ad_w$ and $Ad_p$ for complete generalisation. For all parameter combinations, however, the primary minimum of the DLVO potential is attractive, corresponding to points of zero force on the force profiles, and also resulting in an attractive force minimum for all. This means that particles will remain attached to a surface when trapped within the potential minimum, unless acted on by a sufficiently large external force. There also exists a repulsive barrier in the potential past the primary minimum, which prevents particles from attaching in the primary minimum unless acted on by a sufficiently large external force. Instead, they can become trapped in the significantly weaker secondary minimum, leading to rolling or sliding along the surface rather than attachment at asperities (Chequer et al. Reference Chequer, Carageorgos, Naby, Hussaini, Lee and Bedrikovetsky2021). The magnitude of the energy barrier is influenced strongly by the fluid salinity; however, here $C_m$ is held constant for the aforementioned reasons.

Figure 2. Profiles of the total electrostatic force, $F_e$, with respect to the separation distance $\zeta$ for the two $A_{131}$ and four $A_{132}$ used in this work. An example physical particle size $d={5}\ {\mathrm {\mu }}{\rm m}$ is used; however, the primary force minima remain attractive ($Ad_w$ and $Ad_p$ remain positive) for all $d$ tested in this work since $F_e$ scales linearly with $d$. (a,b) Close-range interactions (small $\zeta$) to illustrate the maximum adhesion force for particle–particle, $F_{a,{p-p}}$, and particle–wall, $F_{a,{p-p}}$, interactions, respectively. (c,d) Long-range interactions (large $\zeta$) to illustrate the secondary force minimum for particle–particle and particle–wall interactions, respectively. Note the significant difference in magnitude between the $F_e$ axis scales, indicating the significantly lower magnitude of the secondary minimum.

2.2. Fluid–solid coupling

To explicitly capture the mechanical and physicochemical forces acting on each particle, the DEM is employed in this work. Through the open-source DEM granular media solver LIGGGHTS (Kloss et al. Reference Kloss, Goniva, Hager, Amberger and Pirker2012), the movement of each solid particle is updated at each discrete time step according to a velocity Verlet integration of Newton's second law of motion. The equations of motion of particle $A$ are

(2.8)\begin{equation} \left.\begin{gathered} m_A\,\frac{{\rm d}\boldsymbol{u}_A}{{\rm d}t} = \sum_{B, \zeta<0}\left(\boldsymbol{F}_{c,AB}\right) + \sum_{C, 0<\zeta<\zeta_c}\left(\boldsymbol{F}_{e,AC}\right) + \boldsymbol{F}_{h,A},\\ I_A\,\frac{{\rm d}\boldsymbol{\omega}_A}{{\rm d}t} = \sum_{B, \zeta<0}\left(\frac{d}{2}\,\boldsymbol{n} \times \boldsymbol{F}_{c,AB}\right), \end{gathered}\right\}\end{equation}

where $B$ therefore represents the set of all particles and walls in contact with particle $A$ ($\zeta <0$), and $C$ represents the set of all particles and walls within the electrostatic cutoff distance $\zeta _c$, but not in contact ($0<\zeta <\zeta _c$). Hence $\boldsymbol {F}_{c,AB}$ accounts for each collision force, $\boldsymbol {F}_{e,AC}$ accounts for electrostatic interactions and $\boldsymbol {F}_{h,A}$ accounts for hydrodynamic force.

The LvdW, EDL and Born forces are implemented as separate forces in the DEM, but are written in (2.8) as a total electrostatic force. The cutoff distance, which is incorporated to reduce computational load, is set to $\zeta _c=100\times 10^{-9} \textrm {m}$, at which point $\boldsymbol {F}_{e,AC}$ reduces to approximately $10^{-5}$ of its maximum value in all cases, as shown in figure 2. Here, this ‘outer’ cutoff distance is distinguished from the ‘inner’ treatment of electrostatics at contact. As described in § 2.1, the latter linearly decreases $\boldsymbol {F}_{e,AC}$ to 0 at contact from an ‘inner’ distance $\zeta /\sigma =0.6$, so as to not interfere with contact force calculations.

The collision force can be further decomposed into the normal and tangential forces, $\boldsymbol {F}_{c,AB} = \boldsymbol {F}_n + \boldsymbol {F}_t$, which are resolved via a nonlinear soft-sphere Hertz contact model. Each comprises a stiffness term subtracted from a damping term, which are nonlinearly dependent on the overlap and relative velocity of two bodies, respectively, for which full expressions can be found in Kloss et al. (Reference Kloss, Goniva, Hager, Amberger and Pirker2012). The terms are also dependent on the Young's modulus, Poisson ratio and coefficient of restitution, which are set to values $5\times 10^{6}$, 0.5 and 0.8, respectively. In § 2.3, it is verified that the collision model and parameter values reproduce the correct contact behaviour.

As outlined in § 1, friction has a strong influence on arch stability and consequently on the clogging probability. Its numerical implementation is therefore important and requires close attention. In this work, friction is implemented for both particle–particle and particle–wall collisions via the Coulomb friction law,

(2.9)\begin{equation} \boldsymbol{F}_t \leqslant \mu \boldsymbol{F}_n, \end{equation}

where $\boldsymbol {F}_t$ is truncated according to the microscopic Coulomb friction coefficient, $\mu$. This microscopic friction coefficient is not equivalent to the macroscopic friction coefficient of the bulk medium, but is a numerical tool used to reproduce it. The macroscopic frictional response of the medium is a property of the system as a whole, and is achieved through resolution of both the normal and tangential forces in all particle contacts in the DEM. The coefficient $\mu$ impacts the latter directly, and therefore the bulk friction.

It is noted that particles are able to roll freely along surfaces, and that rolling friction is not included in this implementation. This is reflected in the absence of an additional torque component in (2.8), as opposed to more sophisticated models that distinguish between rolling and sliding. More complex friction models, however, introduce additional numerical parameters that require calibration to experimental data (Syed, Tekeste & White Reference Syed, Tekeste and White2017). This would be in addition to other numerical parameters including the electrostatic inner cutoff distance (introduced in § 2.1) and the lubrication inner cutoff distance (discussed in detail in § 2.3). In general, this type of rigorous calibration would require comparison to experiments of bulk properties, such as rheology or pressure drop in a confined channel, which consider electrostatics and surface roughness as parameters. However, to the best of the authors’ knowledge, such data have not yet appeared in the literature. To circumvent this problem, this work instead opts for a friction model with only one parameter ($\mu$) and performs a sensitivity analysis with respect to this parameter. Throughout this work, $\mu =1$, which represents a high roughness and particle irregularity, but has the advantage of producing clogging in the majority of cases. In the sensitivity analyses, it is reduced to $\mu =0.2$. For simplification, $\mu$ is set equal for particle–particle and particle–wall collisions.

Finally, the hydrodynamic force contribution from the suspending fluid, $\boldsymbol {F}_{h,A}$, is resolved using a previously described numerical framework (Di Vaira et al. Reference Di Vaira, Łaniewski Wołłk, Johnson, Aminossadati and Leonardi2022) that couples a single-relaxation-time LBM with the DEM via the partially saturated method (PSM) (Noble & Torczynski Reference Noble and Torczynski1998), using the open-source LBM code-base TCLB (Łaniewski-Wołłk & Rokicki Reference Łaniewski-Wołłk and Rokicki2016). Hydrodynamic torque is not applied to particles in order to reduce computational complexity. The primary effect of hydrodynamic torque would be to decrease arch stability, just as the inclusion of rolling torque would increase arch stability. Therefore, also like the decision to not include rolling friction in favour of computational simplicity, hydrodynamic torque is treated as an additional calibration parameter that is accounted for in the sensitivity analysis to $\mu$. While it would also influence surface attachment, § 3.4 demonstrates that the numerical simplifications and geometry reproduce valid critical attachment velocities.

The LBM is employed in this work for two primary reasons. First, the PSM can resolve exactly the momentum exchange between fluids and particles (within the numerical accuracy of the LBM), including when the gap between solid boundaries is less than one lattice spacing (Strack & Cook Reference Strack and Cook2007). While there also exist other methods of particle-resolved DNS, past studies of hydrodynamic particle clogging have applied only classical finite element CFD methods that approximate the solid boundary momentum coupling with drag force correlations (Sun et al. Reference Sun, Xu, Pan, Shi, Geng and Cai2019), the more complex fictitious domain method (Mondal et al. Reference Mondal, Wu and Sharma2016), or a simplified two-fluid Eulerian–Eulerian approach (Li et al. Reference Li, Qiu, Zhong, Zhao and Huang2020). Second, the LBM uses a local and explicit time-stepping scheme that retains the properties of implicit integration (Łaniewski-Wołłk Reference Łaniewski-Wołłk2017), increasing the ease and efficiency of parallel computation on CPU and GPU architectures compared to finite element methods.

Throughout this work, a kinematic LBM lattice viscosity $\nu ^*=1/6$ is used, and the LBM grid spacing, $\delta _x$, is set such that the particle diameters are represented by ten lattice nodes, $d/\delta _x=10$. A number of DEM sub-iterations, $n_{sub}$, may be computed for every LBM iteration, depending on the time resolution required to capture electrostatic forces during a particle rebound. This issue is explored in detail in § 2.3, along with other important considerations of contact modelling. For a general discussion of the numerical error and stability associated with the spatial and temporal discretisation of this LBM–DEM implementation, the reader is referred to the earlier works of Di Vaira et al. (Reference Di Vaira, Łaniewski Wołłk, Johnson, Aminossadati and Leonardi2022) and Han, Feng & Owen (Reference Han, Feng and Owen2007).

2.3. Verification

The accuracy of fluid–solid momentum exchange for the PSM is well established (Noble & Torczynski Reference Noble and Torczynski1998; Strack & Cook Reference Strack and Cook2007; Owen, Leonardi & Feng Reference Owen, Leonardi and Feng2011), while validation of the present implementation, including second-order convergence and bulk suspension behaviour in channel flows, has been demonstrated previously (Di Vaira et al. Reference Di Vaira, Łaniewski Wołłk, Johnson, Aminossadati and Leonardi2022). As described in § 1, however, hydrodynamic clogging is fundamentally dependent on the formation of arches and electrostatic attachment, both of which depend on individual particle contacts and the interplay of hydrodynamic and electrostatic forces. Specifically, the coupled LBM–DEM model must correctly resolve: particle collisions according to the rebound regime; short-range hydrodynamics; the dynamics of electrostatic–hydrodynamic forces near contact; and electrostatic forces as two surfaces separate.

First, the degree of rebound of a particle is dictated by the Stokes number, $St$, which may be interpreted as the ratio of the particle inertia to fluid viscous effects. Below a critical Stokes number, $St_{cr}$, approximately 10–15, viscous dissipation is large enough such that no rebound occurs (Legendre et al. Reference Legendre, Zenit, Daniel and Guiraud2006). With increasing $St$ above $St_{cr}$, however, a particle will rebound to an increasing degree, and the DEM contact time of a numerical model must be calibrated to ensure an adequate fluid response and the correct rebound trajectories (Rettinger & Rüde Reference Rettinger and Rüde2022). The task here, then, is to first determine the rebound regime(s) of the simulations in this work, and then ensure that the correct behaviour is reproduced. To do so, the Stokes number is defined as

(2.10)\begin{equation} St=\frac{1}{9}\,\frac{\rho_p}{\rho_f}\,\frac{\langle v \rangle d}{\nu}, \end{equation}

where $\rho _p$ is the particle density. Based on this definition, $St=0.026$–0.84 for the range of simulations conducted in this work (note that these values would be reduced further if the reductions in $\langle v \rangle$ due to the presence of solids were taken into account). These low $St$ are indicative of the micron-scale particles studied here, and suggest that particles should not physically rebound. To verify this, a simplified version of the pendulum experiment of Yang & Hunt (Reference Yang and Hunt2006) is conducted. Two spheres with $d/\delta _x=10$ are positioned at an initial $x$ separation distance $3d$ at the centre of a periodic fluid domain of size $10d \times 5d \times 5d$ ($x \times y \times z$). A constant acceleration is applied in the $x$ direction to sphere 1, which is switched off just prior to contact with sphere 2. The contact is characterised by the velocities of the two particles at contact, $u_{c,1}$ and $u_{c,2}$, via the binary Stokes number, $St_B=(1/9)(\rho _p/\rho _f)(u_{c,1} - u_{c,2})d/\nu$. For $St_B=0.4$, figure 3(a) illustrates that the two particles remain touching after contact, verifying that the suspending fluid provides the necessary energy dissipation to prevent particle rebound at lattice resolution $d/\delta _x=10$.

Figure 3. Verification of particle contacts for two separate test cases, including (a) the simplified pendulum test, with the trajectories of each particle relative to the point of contact, $x_c$, and time of contact, $t_c$, showing that rebounding does not occur, and (b) the hydrodynamic force, $F$, normalised by the Stokes drag force, $F_{St}$, for two approaching spheres at small separation distances, demonstrating that the PSM with smooth spheres at resolution $d/\delta _x=10$ is equivalent to the analytical solution (sum of Stokes drag and lubrication (2.11)) for particles of effective roughness $\eta _e/d=0.034$.

Second, a large repulsive hydrodynamic force is created as fluid is forced out of the small closing gap of two approaching surfaces. This force, commonly termed the lubrication force, must be captured correctly by numerical models. An approximation of the normal lubrication force for two spheres of equal diameter is given by (Brenner Reference Brenner1961; Izard, Bonometti & Lacaze Reference Izard, Bonometti and Lacaze2014)

(2.11)\begin{equation} \boldsymbol{F}_{lub} ={-}\frac{6{\rm \pi}\rho_f\nu\bar{\boldsymbol{u}}}{\zeta + \eta_e}\left(\frac{d}{4}\right)^2, \end{equation}

where $\bar {\boldsymbol {u}}$ is the relative velocity, and $\eta _e$ is the effective roughness of the spheres. Incorporating an effective roughness height in this manner is consistent with experimental observations (Izard et al. Reference Izard, Bonometti and Lacaze2014), and prevents the force from diverging as the separation distance $\zeta \rightarrow 0$. However, while $\eta _e$ is of the same order of magnitude as the physical roughness, it is not exactly equivalent to any standard measure, such as the mean or root mean square of asperity heights (Mongruel et al. Reference Mongruel, Chastel, Asmolov and Vinogradova2013).

In general, the LBM correctly predicts the total hydrodynamic force (drag plus lubrication) for $\zeta > (2/3)\delta _x$, but then resolves only part of the lubrication force when $\zeta < (2/3)\delta _x$ (Ladd & Verberg Reference Ladd and Verberg2001; Rettinger & Rüde Reference Rettinger and Rüde2022). The same problem faces all fluid numerical solvers, therefore a lubrication correction force can be added within distances where the solver underestimates the total hydrodynamic force. To account for roughness, or to prevent divergence, the correction force is typically cut off at a separation distance that is calibrated against experimental results (e.g. particle rebound trajectories for the case $St>St_{cr}$; Rettinger & Rüde Reference Rettinger and Rüde2022). In this work, the lubrication correction is not added because at the resolution used, the element size and roughness of the particles are in the same range. In fact, the hydrodynamic force predicted by the PSM matches the analytical hydrodynamic force of rough particles. To demonstrate this, two smooth particles are approached at constant velocity, and the total hydrodynamic force acting on each particle is measured. Figure 3(b) demonstrates that for particles with $d/\delta _x=10$, the force obtained from the LBM closely matches the analytical force (Stokes drag plus lubrication) for $\eta _e/d=0.034$ over all separation distances. In other words, the computational resolution used throughout this work is representative of particles with effective roughness $\eta _e/d=0.034$. For a particle size $d=5\times 10^{-6}\ \textrm {m}$, this equates to $\eta _e=0.17\times 10^{-6}\ \textrm {m}$, which is commensurate with the mean or root mean square physical roughness of many materials, such as steel, glass and polymers, which commonly fall into the range $0.1{-}0.5\times 10^{-6}\ \textrm {m}$ (Gondret, Lance & Petit Reference Gondret, Lance and Petit2002; Yang & Hunt Reference Yang and Hunt2006). It is reiterated, however, that $\eta _e$ and the physical roughness are not exactly equal. Any calibration of or sensitivity analysis in relation to $\eta _e$ (via changing the grid resolution or adding a lubrication correction) is considered outside the scope of this work. Instead, a sensitivity analysis is performed with respect to $\mu$ as discussed in § 2.2.

Third, the implementation of DLVO forces and resulting dynamic behaviour with hydrodynamic forces near contact is verified. As illustrated in figure 4(a), two spheres with $d/\delta _x=10$ are positioned at an initial separation distance $\zeta _0=0.62\sigma$, at the centre of a periodic fluid domain of size $5d \times 5d \times 5d$ in a fully coupled LBM–DEM simulation. The full inter-particle DLVO force, as well as a constant gravitational body force, $F_g$, are applied to one of the spheres only, such that $F_a/F_g=1.22$. The other sphere is fixed. Based on a force balance, the moving sphere should reach an equilibrium position $\zeta =0.66\sigma$ when $F_e=F_g$. Figure 4(b) shows that the spheres initially move apart and then oscillate to $\zeta =0.66\sigma$ due to damping from the fluid. In other words, the spheres behave as a nonlinear spring–damper system. For reference, the non-damped DEM-only simulation is also plotted. Overall, this verifies first that the DEM correctly calculates the analytical DLVO forces presented in § 2.1, and second that the LBM coupling introduces no spurious forcing. It is recognised, however, that at such small separation distances, which are much less than the lattice spacing ($\zeta \approx \delta _x 10^{-4}$), the system dynamics is highly sensitive to the lattice spacing. Analysis of the fluid–particle dynamic response in relation to the computational resolution is outside the scope of this work. However, this issue has hitherto been given insufficient attention in CFD modelling that incorporates electrostatics, and should be a focus of future research.

Figure 4. Verification test for two spheres in close proximity, showing (a) the initial test configuration (sphere sizes and separation distance are not to scale), where the top sphere is fixed and the total DLVO force, $F_e$, and a body force, $F_g$, are applied to the bottom sphere, and (b) that an equilibrium position of $\zeta =0.66\sigma$ is attained, verifying correct implementation of hydrodynamic and DLVO forces.

Finally, the small separation distances over which DLVO forces act mean that the DEM time resolution, determined by $n_{sub}$, must be sufficiently fine that the electrostatic interaction of two bodies is captured as they separate after contact. At the same time, however, a finite limit on $n_{sub}$ is required, since the computational expense scales with $n_{sub}$. The required time resolution depends on the relative velocity of the interaction: for a particle contacting and rebounding from a wall, for example, it will either escape the electrostatic force minimum or remain trapped and attached to the wall. If the DEM time step is too large, however, then the distance rebounded in a single time step will exceed the distance at which the electrostatic minimum occurs ($\zeta \approx 0.72\sigma$) and the electrostatic force will not be seen by the particle. Depending on the contact velocity, this may mean the difference between a particle remaining trapped or escaping. Therefore, to determine the required $n_{sub}$, the same simulation set-up as for the previous verification test is used; however, in this case, the stationary particle is replaced with a wall (figure 5a). An initial velocity, $V_0$, is applied to the particle, and it then overlaps with the wall and rebounds due to the DEM contact force (figure 5b). This set-up reflects a ‘worst case’ scenario, where a particle is already attached to a wall and is impacted by another particle travelling in the flowing fluid, and the hydrodynamic force is imparted completely to the wall particle, perpendicular to the wall, via $V_0$. Four different adhesion numbers ($Ad_w=4.5,9,18,70$) are tested, reflecting the full range of simulations in this work, where $V_0$ replaces the characteristic mean channel velocity in (2.6). The initial separation distance is $\zeta _0=0.6\sigma$. The resulting particle trajectories are shown in figures 5(c)–5(f), noting that the large negative $\zeta$ values are due to the non-physical overlap of the DEM. First, for $Ad_w=70$, it can be concluded that no sub-stepping is required, since the electrostatic adhesion force is sufficiently large for $n_{sub}=1$ to adequately match $n_{sub}=100$. For the lower $Ad_w=18$, however, upon rebounding, the particle should remain captured by the wall, according to the high-resolution case $n_{sub}=100$. For $n_{sub}=1$, the DEM time resolution fails to reproduce this, while for $n_{sub}=2$ and $n_{sub}=3$, the trajectory is significantly different from $n_{sub}=100$. The case $n_{sub}=4$, however, may be regarded as an adequate resolution of the contact. As $Ad_w$ is decreased further to 9, the particle escapes in all cases; however, it is clear that $n_{sub}=8$ is required to approximate adequately the high-resolution case (the larger scale of the $y$-axis should be noted). Finally, for the lowest $Ad_w=4.5$, $n_{sub}=8$ again gives an adequate approximation. According to these results, for the following simulations of this work, a maximum $n_{sub}=8$ is employed for $Ad_w \lesssim 9$, a minimum $n_{sub}=1$ is employed for $Ad_w \gtrsim 70$, and for $9 < Ad_w < 70$, $n_{sub}$ is varied in a linear manner from 8 to 1.

Figure 5. Verification test for a particle contacting a wall to determine the required number of DEM sub-iterations, $n_{sub}$, as a function of $Ad_w$. (a) The initial test configuration (particle size and separation distance are not to scale), where the total DLVO force, $F_e$, and an initial velocity, $V_0$, are applied to the particle. (b) The particle contacting the wall, and the resulting DEM contact force, $F_c$, where the overlap is accentuated to convey the non-physical overlap of the DEM. The resulting particle trajectories in (c–f) demonstrate that increasing DEM time resolution is required to capture electrostatic forces as $Ad_w$ decreases, up to maximum resolution $n_{sub}=8$, where (c) $Ad_w=70$, (d) $Ad_w=18$, (e) $Ad_w=9$, and (f) $Ad_w=4.5$.

3.1. Numerical test cell

The numerical test cell depicted in figure 6 is employed in this work to replicate clogging in planar channels and obtain the clogging data presented in § 4. The key feature of this test cell is the ridges of the channel walls, which mimic a planar channel with long-wavelength roughness. These ridges have two functions. First, they provide surface asperities to which particles can attach with electrostatics since, even if it was included, the torque generated by the inclusion of rolling friction would be insufficient to cause attachment at the desired pressure gradients. The modelling of attachment pressure gradients is discussed at length in § 3.4. Second, the asperities induce the transverse particle motion required to cause arching, since clogging does not occur in perfectly smooth channels. This has been shown experimentally (Barree & Conway Reference Barree and Conway2001; Ray et al. Reference Ray, Lewis, Martysevich, Shetty, Walters, Bai and Ma2017) and is explained by the fact that walls push particles away since pressure is higher on the side of a particle closest to a wall (Feng, Hu & Joseph Reference Feng, Hu and Joseph1994). Consequently, even if particles could attach to a perfectly smooth numerical wall via electrostatics, they would never contact the wall in the first place. Overall, in numerical implementations, the absence of realistic surface roughness may be regarded as a numerical artefact that negates attachment and arching, even in suspensions approaching the maximum particle packing limit.

Figure 6. Planar channel test cell with ridged surfaces for clogging simulations. Fluid is driven by a constant pressure gradient in the $z$ direction, and fluid boundaries are periodic in the $y$ and $z$ directions. Particles are injected continuously in the region $z \in [l/10,l/5]$, and particle boundaries are periodic in the $y$ direction only. The shaded surfaces represent the walls that impart an electrostatic force on the particles.

The suspending fluid is driven by a constant body force in the $z$ direction, akin to a constant pressure gradient $G=\textrm {d}P/\textrm {d}z$. The fluid is Newtonian, with particle to fluid density ratio $\rho _p/\rho _f=3.5$; however, no gravitational force is applied such that particles do not settle transverse to the flow.

The inlet and outlet boundaries of the fluid are periodic; however, particles exit the simulations once they reach the end of the computational domain. In this way, the inlet of the fluid experiences the flow effects of particles that are close to the outlet; however, the particles do not re-enter at the inlet. This set-up is required since the pressure gradient is modelled by a constant body force in the LBM simulation, which provides a quick flow profile development in a small computational domain, in comparison to prescribed pressure boundaries at the inlet and outlet. The particles are not periodic since a random distribution is desired at insertion, in which the incoming particles are not influenced by or coupled to the rest of the simulation. Particles are injected at the inlet of the channel in the region $z \in [l/10,l/5]$ to allow development of a smooth fluid profile. They are added constantly to the injection region such that the bulk solid volume fraction, $\bar {\phi }$, is maintained in this region. The vertical direction, $y$, is periodic for both particles and fluid. A study of the domain dependence in the $z$ and $y$ directions is presented in § 3.3. The clogging probability is shown to be approximately independent of the domain dimensions for $l/d=20$ and $h/d=15$, which are the values used subsequently throughout this work.

The ridges on each channel wall are symmetrically opposed, and the channel width, $w$, is defined as the perpendicular distance from an electrostatic surface to the surface of an opposing ridge. All particles in a single simulation have identical sizes (i.e. all suspensions are monodisperse), yielding a single value of $w/d$ for each simulation. While this work is focused on monodisperse suspensions, it has been shown that for polydisperse suspensions (comprising particles of multiple sizes), it is the largest particles that define whether clogging will occur (STIM-LAB, Inc. 1995; Xu et al. Reference Xu, Sun, Cai and Geng2020). Consequently, the study of monodisperse suspensions is a meaningful simplification.

For ease of computational implementation, wall electrostatic forces are implemented on each side of the channel (in the $z$–$y$ plane) using a smooth electrostatic plane without ridges. Its position therefore coincides with the recessed part of the physical wall, as indicated by the shaded surfaces in figure 6, and the force is in the $x$ direction only. The choice of ridge geometry is relatively arbitrary; however, ridge heights $a/d=0.1$ provide sufficient moment arms for the walls’ normal electrostatic adhesion force to balance with the hydrodynamics of the suspending fluid (figure 7). Overall, these modelling simplifications adequately reproduce surface attachment, which is investigated in detail in § 3.4. Nevertheless, dependence of surface attachment and clogging on the ridge geometry and electrostatic wall implementation warrants further investigation, particularly when translating the computational domain to nonlinear realistic textures.

Figure 7. Torque balance of hydrodynamic and electrostatic forces at a surface asperity.

3.2. Determining clogging probability

To determine the clogging probability, $P$, this work effectively employs the commonly used method (Guariguata et al. Reference Guariguata, Pascall, Gilmer, Sum, Sloan, Koh and Wu2012; Mondal et al. Reference Mondal, Wu and Sharma2016; Chen et al. Reference Chen, Ho-Minh and Tan2018; Sun et al. Reference Sun, Xu, Pan, Shi, Geng and Cai2019) of conducting multiple Bernoulli trials at a fixed combination of parameters, and dividing the number of trials that clog, $N_j$, by the total number of trials, $N_t$, i.e. $P = N_j/N_t$. Each trial uses a unique and random particle injection seed to achieve numerical randomness. In § 4, however, a methodology is detailed for predicting the variation of the discrete $P$ with a free parameter using a continuous logarithmic regression model, and results are interpreted in terms of this predictive model.

The approach of obtaining a probability using a finite number of trials, obtained over a finite time, comprises two aspects of uncertainty. First, for what total simulation time, $t_s$, should each simulation be run; and second, how many trials are required to give a sufficiently accurate $P$? This second aspect is addressed in an uncertainty quantification in the Appendix, while the first aspect is addressed in the following.

Apart from the fact that minimising $t_s$ minimises computational expense, constraining $t_s$ is necessary since, in theory, $P \rightarrow 1$ as $t_s \rightarrow \infty$ (Koivisto & Durian Reference Koivisto and Durian2017). However, as noted in § 1, there are also critical values of parameters for which the probability of clogging is observable. Using the example of particle concentration, there exists a critical concentration above which clogging is observed ($P>0$), but below which it is not. Predicting $P$ in the limit as $t_s \rightarrow \infty$ is also crucial if numerical models are to be applied at physically realistic time and length scales. In this respect, as discussed in § 1, $P$ follows an exponential probability distribution function with time (Lafond et al. Reference Lafond, Gilmer, Koh, Sloan, Wu and Sum2013; Marin et al. Reference Marin, Lhuissier, Rossi and Kähler2018; López et al. Reference López, Hernández-Delfin, Hidalgo, Maza and Zuriguel2020). This modelling aspect, however, is considered outside the scope of the present work.

Instead, this study determines $t_s$ by defining a dimensionless distance,

(3.1)\begin{equation} \bar{D}=\frac{\displaystyle\sum\nolimits_{i} {\langle v \rangle}_i \delta_t}{d}, \end{equation}

where ${\langle v \rangle }_i$ is the mean velocity of all particles with non-zero velocity, at each time step $i$. Particles that have arched locally while other regions of the domain are still flowing (as shown in § 3.3) are therefore not considered in the calculation of ${\langle v \rangle }_i$, nor are particles that have attached to the channel wall when electrostatics is present. Consequently, $\bar {D}$ is effectively how many particle diameters the suspension would travel if not clogged. Letting each simulation run for the same $\bar {D}$ gives consistency to the measurement of $N_j$.

A value $\bar {D}=200$ gives sufficient time for full flow development. As illustrated in figure 8, $\langle v \rangle$ initially decreases as more particles enter the channel and the porosity effectively decreases. Complete flattening of the transient velocity curves indicates the point where the channel has become fully occupied. This behaviour occurs independently of any arching.

Figure 8. Variation in the mean velocity of non-arched particles, $\langle v \rangle$, with dimensionless suspension distance $\bar {D}$ (see (3.1)). For these five selected simulations at these parameters ($w/d=2.4$, $St=0.19$, $\bar {\phi }=0.24$), full suspension clogging does not occur within $\bar {D}=200$; however, dips in the velocity indicate periods of increased localised arching. This demonstrates how the flow has a transient development period as the channel gradually becomes fully occupied by particles.

3.3. Domain dependence

The use of a periodic particle boundary in the $y$ direction imposes a dependence on the domain height, while a dependence on the domain length also exists since the number of ridges increases with $l$. Here, these dependencies are quantified by measuring $P$ over a range of normalised domain dimensions, $h/d$ and $l/d$. Ten trials are conducted at each parameter combination. Only the non-electrostatic scenario is analysed.

First, considering dependence on $h$, it must be recognised that for a planar channel, localised arching occurs first prior to complete clogging. This is depicted in both figure 9(a) and supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.214. It is emphasised, however, that complete channel clogging does not necessarily follow from localised arching, since non-stable arches can disintegrate. If the domain height is insufficiently large, then an occurrence of a local arch could result in arching immediately propagating over the entire domain height (complete channel clogging) due to interactions of particles over the periodic $y$ boundary. This immediate propagation would be more indicative of a confined channel (such as a cylinder), where the flow is bounded on all sides. Therefore, for a planar channel simulation to be $h$-independent, local arches must be able to form independently of one another. The difference between immediate and gradual arch propagation can be characterised by plots of the number of particles exiting the simulation ( figures 9b,c). For $h/d=4$ (figure 9b), the cumulative particle exit count becomes constant instantaneously, while for $h/d=15$ (figure 9c), the rate of particles exiting decreases gradually before becoming constant.

Figure 9. Illustrating the concept of localised arching and gradual arch propagation for planar channels, showing (a) the existence of localised arching, which does not necessarily lead to complete channel clogging. (b,c) Plots of the particle exit count versus time for five trials, each at $w/d=1.8$, $St=0.1$, $\bar {\phi }=0.07$, and (b) $h/d=4$ or (c) $h/d=15$. A gradual decrease in the rate of particles exiting for $h/d=15$ indicates gradual arch propagation leading to channel clogging. For $h/d=4$, on the other hand, the number of particles exiting the domain becomes constant instantaneously, indicating immediate arch propagation and clogging due to particle interactions over the periodic $y$ boundary.

Figure 10(a) shows that the clogging probability becomes independent of the domain height at approximately $h/d=20$, for both of the $w/d$ values tested. As $h$ decreases, $P$ increases due to increased periodic particle interaction, and therefore more occurrences of immediate arch propagation. In the subsequent results in this paper, a domain height $h/d=15$ is employed to minimise computational expense while eliminating most of the periodic height effects.

Figure 10. Dependence of $P$ on (a) the domain height, $h/d$, and (b) the domain length, $l/d$, where comparison is made between the original dimensionless distance, $\bar {D}$, and the alternative length-based dimensionless distance $\bar {D}_l$. Ten trials are conducted to obtain each data point.

The clogging probability increases as the channel length increases because the number of ridges, and hence opportunities for arching, becomes larger. If the number of ridges is too low, then there will be insufficient opportunities for arching. An intuitive way of accounting for this dependence on length is to define an alternative dimensionless distance,

(3.2)\begin{equation} \bar{D}_l=\frac{\displaystyle\sum\nolimits_{i} {\langle v \rangle}_i \delta_t}{d}\,\frac{l}{d}, \end{equation}

which effectively decreases the simulation time as the channel length increases. Figure 10(b) plots the dependence of $P$ on the channel length, with the simulation time determined by both the original dimensionless distance, $D$, and the new length-based dimensionless distance, $\bar {D}_l$. Evidently, when using the original dimensionless distance, $P$ becomes independent of the domain length at $l/d \gtrsim 16$. The length-based dimensionless distance, on the other hand, over-compensates for the additional channel length/ridges. This is evidenced by the decreasing $P$ for increasing $l/d$ at $l/d \gtrsim 20$, which is a result of the decreasing simulation time. These results suggest that for $l/d \gtrsim 16$, the addition of more ridges does not influence the clogging probability. Further, $l/d=12$ clearly provides insufficient opportunities for arching. Consequently, $l/d=20$ is used throughout this work, which corresponds to five ridges along the channel length, as depicted in figure 6.

3.4. Surface attachment

To validate that the asperity geometry defined in § 3.1 gives a physical representation of surface attachment, the critical mean fluid velocities at which attachment occurs, $\langle v \rangle ^*$, are here determined for a range of parameters. To do so, a single particle is initialised at a surface asperity, as in figure 7, and $\langle v \rangle$ is decreased in subsequent simulations until the particle remains attached to the wall.

First, figure 11(a) shows that $\langle v \rangle ^*$ increases approximately linearly with the magnitude of the particle–wall adhesion force, $F_{a,{p-w}}$. The $\langle v \rangle ^*$ values are greater for $w/d=2.4$ compared to $w/d=1.8$ since the larger channel results in a lower velocity at the wall, and therefore lower hydrodynamic drag, for an equivalent mean velocity.

Figure 11. The critical mean fluid velocity for electrostatic surface attachment, $\langle v \rangle ^*$. Dependence on (a) the particle–wall adhesion force, $F_{a,{p-w}}$, and non-dimensional channel width, $w/d$, at $d={5}\ {\mathrm {\mu }}\textrm {m}$, and (b) the particle diameter, $d$, at $A_{132}=8.4\times 10^{-21}\ \textrm {J}$ and $w/d=1.8$.

Next, the physical particle size is varied at a constant Hamaker constant to give an indication of how surface attachment changes as hydrodynamic forces dominate electrostatic forces. Both the lattice resolution and $w/d$ are also held constant as $d$ is increased. Figure 11(b) demonstrates that $\langle v \rangle ^*$ decreases with increasing $d$ in a nonlinear fashion. This trend has been predicted previously theoretically (Yang et al. Reference Yang, Siqueira, Vaz, You and Bedrikovetsky2016; You et al. Reference You, Yang, Badalyan, Bedrikovetsky and Hand2016) and observed experimentally (Chequer et al. Reference Chequer, Carageorgos, Naby, Hussaini, Lee and Bedrikovetsky2021). The nonlinear relationship can be attributed to the scaling of the particle forces, where the electrostatic attractive force increases approximately linearly with $d$, while the hydrodynamic drag increases quadratically with both $d$ and the channel velocity.

The range of critical attachment velocities obtained here ($0.017 \leqslant \langle v \rangle ^* \leqslant 0.202$) is commensurate with experimental values obtained previously at similar electrostatic parameter values (Chequer et al. Reference Chequer, Carageorgos, Naby, Hussaini, Lee and Bedrikovetsky2021). It is emphasised, however, that the ridge height employed in the current study is non-physical compared to the size of the experimental surface asperities (${\approx }1\times 10^{-8}\ \textrm {m}$), and instead may be viewed as an additional numerical parameter on which $\langle v \rangle ^*$ is dependent. Since physically valid $\langle v \rangle ^*$ and correct trends in relation to other parameters have been achieved, however, the choice of ridge height in this study, along with the decisions to implement a simplified friction model and neglect hydrodynamic torque, can be considered valid.