2.1 Fully coupled grain-resolving simulations

We solve the unsteady Navier–Stokes equations for an incompressible Newtonian fluid, given by

(2.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}\boldsymbol{u})=-\frac{1}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D708}_{f}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\boldsymbol{f}_{\mathit{IBM}},\end{eqnarray}$$

along with the continuity equation

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

on a uniform rectangular grid with grid cell size $\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=\unicode[STIX]{x0394}z=h$ . Here, $\boldsymbol{u}=(u,v,w)^{\text{T}}$ designates the fluid velocity vector in Cartesian components, $p$ denotes the pressure, $\unicode[STIX]{x1D708}_{f}$ is the kinematic viscosity, $t$ the time and $\boldsymbol{f}_{\mathit{IBM}}$ represents an artificial volume force introduced by the immersed boundary method (IBM; Uhlmann Reference Uhlmann2005; Kempe & Fröhlich Reference Kempe and Fröhlich2012b ). This volume force, which acts on the right-hand side of (2.1) in the vicinity of the interphase boundaries, connects the motion of the particles to the fluid phase. We integrate equations (2.1) and (2.2) by a third-order low-storage Runge–Kutta (RK) scheme and a finite differencing approach in time and space, respectively. The pressure is treated with a direct solver based on fast Fourier transforms.

Note that gravity has been omitted from the equation of motion for the fluid (2.1), because the contribution from hydrostatic pressure is not relevant to the problems presented in the following. Gravity is, however, explicitly accounted for in the equations of motion for the particles. Within the framework of the IBM, we calculate the motion of each individual spherical particle by solving an ordinary differential equation for its translational velocity $\boldsymbol{u}_{p}=(u_{p},v_{p},w_{p})^{\text{T}}$ ,

(2.3) $$\begin{eqnarray}m_{p}\frac{\text{d}\boldsymbol{u}_{p}}{\text{d}t}=\underbrace{\oint _{\unicode[STIX]{x1D6E4}_{p}}\boldsymbol{\unicode[STIX]{x1D749}}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}A}_{=\,\boldsymbol{F}_{h,p}}+\underbrace{V_{p}(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})\boldsymbol{g}}_{=\,\boldsymbol{F}_{g,p}}+\,\boldsymbol{F}_{c,p},\end{eqnarray}$$

and its angular velocity $\boldsymbol{\unicode[STIX]{x1D74E}}_{p}=(\unicode[STIX]{x1D714}_{p,x},\unicode[STIX]{x1D714}_{p,y},\unicode[STIX]{x1D714}_{p,z})^{\text{T}}$ ,

(2.4) $$\begin{eqnarray}I_{p}\frac{\text{d}\boldsymbol{\unicode[STIX]{x1D74E}}_{p}}{\text{d}t}=\underbrace{\oint _{\unicode[STIX]{x1D6E4}_{p}}\boldsymbol{r}\times (\boldsymbol{\unicode[STIX]{x1D749}}\boldsymbol{\cdot }\boldsymbol{n})\,\text{d}A}_{=\,\boldsymbol{T}_{h,p}}+\,\boldsymbol{T}_{c,p}.\end{eqnarray}$$

Here, $m_{p}$ denotes the particle mass, $\unicode[STIX]{x1D6E4}_{p}$ the fluid–particle interface, $\boldsymbol{\unicode[STIX]{x1D749}}$ the hydrodynamic stress tensor, $\unicode[STIX]{x1D70C}_{p}$ the particle density, $V_{p}$ the particle volume, $g$ the gravitational acceleration, $I_{p}=8\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{p}R_{p}^{5}/15$ the moment of inertia, and $R_{p}$ the particle radius. Furthermore, the vector $\boldsymbol{n}$ represents the outward-pointing normal on the interface $\unicode[STIX]{x1D6E4}_{p}$ , $\boldsymbol{r}=\boldsymbol{x}-\boldsymbol{x}_{p}$ is the position vector of the surface point with respect to the centre of mass $\boldsymbol{x}_{p}$ of a particle, and $\boldsymbol{F}_{c,p}$ and $\boldsymbol{T}_{c,p}$ indicate the force and torque due to particle collisions, respectively. For the sake of brevity, we denote the hydrodynamic force and torque as $\boldsymbol{F}_{h,p}$ and $\boldsymbol{T}_{h,p}$ , respectively, and the gravitational force as $\boldsymbol{F}_{g,p}$ . We use an RK scheme that subdivides the three-step procedure of the fluid into a total of 15 substeps per fluid time step to integrate the particles’ equations of motion (2.3) and (2.4) in time. It was shown by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) that this is a necessity to resolve short-range effects of lubrication forces in time.

The fluid–particle interaction was validated in Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) by comparing our simulation results to experimental data for a settling sphere in a large container (Mordant & Pinton Reference Mordant and Pinton2000), and for a particle settling above a wall (Ten Cate et al. Reference Ten Cate, Nieuwstad, Derksen and Van den Akker2002), yielding excellent agreement.

2.2 Cohesionless particle–particle interaction

The computational approach for modelling cohesionless particle–particle interactions is described in detail in Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ), and validation results are provided for normal and oblique binary collisions, as well as for the collective motion of a sediment bed sheared by a Poiseuille flow. In order to keep this paper self-contained, we provide a brief summary in the following.

The particle–particle interaction comprises short-range hydrodynamic effects due to lubrication forces $\boldsymbol{F}_{l}$ , as well as forces acting in the normal and tangential directions for direct particle contact, denoted as $\boldsymbol{F}_{n}$ and $\boldsymbol{F}_{t}$ , respectively. The resulting collision force on particle $p$ is the sum off all these effects,

(2.5) $$\begin{eqnarray}\boldsymbol{F}_{c,p}=\mathop{\sum }_{q,q\neq p}^{N_{p}}(\boldsymbol{F}_{l,pq}+\boldsymbol{F}_{n,pq}+\boldsymbol{F}_{t,pq})+\boldsymbol{F}_{l,pw}+\boldsymbol{F}_{n,pw}+\boldsymbol{F}_{t,pw},\end{eqnarray}$$

where the subscripts $pq$ and $pw$ indicate interactions with particle $q$ or a wall, respectively. In what follows, we present the algebraic expressions for particle–particle interaction only. Analogous formulations for particle–wall interactions can be found in Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ). Consistent with the findings of Cox & Brenner (Reference Cox and Brenner1967), we model the unresolved component of the lubrication forces in our simulations as

(2.6) $$\begin{eqnarray}\boldsymbol{F}_{l,pq}=\left\{\begin{array}{@{}ll@{}}-{\displaystyle \frac{6\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D708}_{f}R_{eff}^{2}}{\max (\unicode[STIX]{x1D701}_{n},\unicode[STIX]{x1D701}_{min})}}\boldsymbol{g}_{n}\quad & 0<\unicode[STIX]{x1D701}_{n}\leqslant 2h,\\ 0\quad & \text{otherwise},\end{array}\right.\end{eqnarray}$$

where $\boldsymbol{g}_{n}=\boldsymbol{u}_{p}-\boldsymbol{u}_{q}$ is the relative velocity of the two colliding particles. To prevent the unresolved lubrication forces from diverging to infinity with decreasing gap size, the force is limited by $\unicode[STIX]{x1D701}_{min}$ , which can be interpreted as a surface roughness of the particles. The value of $\unicode[STIX]{x1D701}_{min}=3\times 10^{-3}R_{m}$ was calibrated in Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) for particle–wall collisions to match the rebound trajectories of the experiments by Gondret, Lance & Petit (Reference Gondret, Lance and Petit2002). The mean radius becomes $R_{m}=R_{p}$ and $R_{m}=(R_{p}+R_{q})/2$ for particle–wall and particle–particle interactions, respectively. The effective radius $R_{eff}$ is defined as $R_{eff}=R_{p}R_{q}/(R_{p}+R_{q})$ , where we set $R_{q}=\infty$ for particle–wall collisions. We also note that (2.6) only accounts for the part that cannot be resolved with the IBM as $\unicode[STIX]{x1D701}_{n}$ becomes smaller than $2h$ . We have conducted detailed tests repeating the test case of Ten Cate et al. (Reference Ten Cate, Nieuwstad, Derksen and Van den Akker2002) for the particle sizes of interest here. Our results show that the particles rapidly decelerate as soon as they come as close as $\unicode[STIX]{x1D701}_{n}=2D_{p}$ , illustrating that nearly all of the work required to squeeze the fluid out of the gap is fully resolved in our simulations.

Direct particle contact is accounted for by a normal and a tangential component of the collision force. The repulsive normal component is represented by a nonlinear spring–dashpot model for the normal direction,

(2.7) $$\begin{eqnarray}\boldsymbol{F}_{n,pq}=-k_{n}|\unicode[STIX]{x1D701}_{n}|^{3/2}\boldsymbol{n}-d_{n}\boldsymbol{g}_{n,cp},\end{eqnarray}$$

where $\boldsymbol{g}_{n,cp}$ denotes the normal component of the relative velocity at the surface contact point (Kempe & Fröhlich Reference Kempe and Fröhlich2012a ). Furthermore, $k_{n}$ and $d_{n}$ represent stiffness and damping coefficients that are adaptively calibrated for every collision as described by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ), in order to yield a prescribed restitution coefficient $e_{dry}=-u_{out}/u_{in}$ . Here, $u_{out}$ and $u_{in}$ indicate the normal components of the relative particle speed immediately after and right before the particle impact, respectively. The forces in the tangential direction are modelled by a linear spring–dashpot model capped by the Coulomb friction law as

(2.8) $$\begin{eqnarray}\boldsymbol{F}_{t,pq}=\min (-k_{t}\boldsymbol{\unicode[STIX]{x1D73B}}_{t}-d_{t}\boldsymbol{g}_{t,cp},\Vert \unicode[STIX]{x1D707}\boldsymbol{F}_{n}\Vert \boldsymbol{t}),\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ represents the coefficient of friction between the two surfaces and $\boldsymbol{\unicode[STIX]{x1D73B}}_{t}$ is the tangential displacement integrated over the time interval for which the two particles are in contact. The tangential stiffness and damping coefficients $k_{t}$ and $d_{t}$ are adapted to account for zero-slip rolling or sliding according to the Coulomb friction law (Thornton, Cummins & Cleary Reference Thornton, Cummins and Cleary2013). For all of the simulations to be presented in the following, we have chosen $e_{dry}=0.97$ and $\unicode[STIX]{x1D707}=0.15$ , which is a common parametrization for silicate materials (e.g. Joseph et al. Reference Joseph, Zenit, Hunt and Rosenwinkel2001; Joseph & Hunt Reference Joseph and Hunt2004; Vowinckel, Kempe & Fröhlich Reference Vowinckel, Kempe and Fröhlich2014; Biegert et al. Reference Biegert, Vowinckel and Meiburg2017a ; Vowinckel et al. Reference Vowinckel, Nikora, Kempe and Fröhlich2017a ,Reference Vowinckel, Nikora, Kempe and Fröhlich b ).

3.1 Physical background

To derive a cohesive force model suitable for the framework of the IBM, we start with the classical Derjaguin–Landau–Verwey–Overbeek (DLVO) theory (Derjaguin & Landau Reference Derjaguin and Landau1941; Verwey & Overbeek Reference Verwey and Overbeek1948). This theory was derived for colloids and is based on the assumption that there are two dominant short-range forces that can be interpreted as opposing potentials surrounding particles with grain sizes in the micro- to nanometre range. On one hand, there exists a repulsive force when equally charged surfaces are in close proximity. On the other hand, as one particle causes correlations in the fluctuating polarization of a nearby particle surface, an attractive force is generated. The former effect is usually called the repulsive ‘double-layer’ (DL) force, while the latter effect is commonly referred to as van der Waals (vdW) force. These forces become important for gap sizes $\unicode[STIX]{x1D701}_{0}<\unicode[STIX]{x1D701}_{n}<\unicode[STIX]{x1D701}_{\infty }$ , where $\unicode[STIX]{x1D701}_{0}$ defines the microscopic size of surface asperities and $\unicode[STIX]{x1D701}_{\infty }$ is the distance for which these forces decay to zero (Israelachvili Reference Israelachvili1992). The repulsive DL force and the attractive vdW force due to polarization scale as $F_{rep}\propto \text{e}^{-\unicode[STIX]{x1D701}_{n}}$ and $F_{att}\propto \unicode[STIX]{x1D701}_{n}^{-2}$ , respectively. Note that the quadratic scaling of $F_{att}$ applies to radii much larger than $\unicode[STIX]{x1D701}_{0}$ (Kosinski & Hoffmann Reference Kosinski and Hoffmann2010; Breuer & Almohammed Reference Breuer and Almohammed2015), while linear scaling of vdW forces has been reported for cylinders of smaller size (Israelachvili Reference Israelachvili1992). A qualitative sketch of the DLVO theory is given in figure 1(a). We elaborate further on the shape of this figure in appendix A. The superposition of the two potentials yields a net force as a function of the gap size $\unicode[STIX]{x1D701}_{n}$ . Depending on the properties of the particles’ surface charge and the salinity of the ambient fluid, this net force exhibits several distinct characteristics. Hence, figure 1(a) displays the characteristics of silica particles with small to medium surface charge (Wu, Ortiz & Jerolmack Reference Wu, Ortiz and Jerolmack2017) and rather low salinity: (i) For larger gap sizes, there exists an outer range where attractive forces dominate over repulsive forces. (ii) Within this outer range the net attractive force displays a maximum. (iii) In the inner range with a maximum, a force barrier dominated by the repulsive double-layer force can be found. (iv) Finally, the attractive forces diverge to infinity for gap sizes smaller than $\unicode[STIX]{x1D701}_{0}$ . The latter condition is equivalent to merging two separate objects into one, although evidence suggests that this condition is never reached for rough surfaces, since asperities on the particle surfaces prevent them from coming into such close contact (Parsons, Walsh & Craig Reference Parsons, Walsh and Craig2014).

Figure 1. Schematic of short-range forces versus gap size: (a) force profiles according to the DLVO theory; (b) the model ansatz given by (3.2). The red vertical lines in (a) indicate the range considered by the model ansatz shown in (b).

The DLVO theory holds only for gap sizes in the nanometre range. Since we cannot resolve this length scale in simulations involving hundreds of micrometre-size particles, we instead employ a simplified algebraic expression for the vdW forces that reproduces its integral properties. The most common expression for vdW forces scales linearly with the particle diameter, as reviewed by Visser (Reference Visser1989) and Israelachvili (Reference Israelachvili1992),

(3.1) $$\begin{eqnarray}\boldsymbol{F}_{vdW}=\frac{A_{H}R_{eff}}{6\unicode[STIX]{x1D701}_{0}^{2}}\boldsymbol{n}.\end{eqnarray}$$

Owing to its simplicity and ease of implementation and interpretation, this scaling assumption has been popular in discrete element methods (DEMs) and point-particle approaches, as well as for experimental analysis (e.g. Ye, van der Hoef & Kuipers Reference Ye, van der Hoef and Kuipers2004; Pandit, Wang & Rhodes Reference Pandit, Wang and Rhodes2005; Righetti & Lucarelli Reference Righetti and Lucarelli2007; Breuer & Almohammed Reference Breuer and Almohammed2015). However, such methods do not resolve the gap size, and cohesive forces effectively act only when particles come into contact. This approach is hence equivalent to lowering the restitution coefficient of the inelastic collision, in line with the underlying hard-sphere collision model. The hard-sphere model modifies the velocity right after the impact as $u_{out}=-eu_{in}$ , where $e$ is the restitution coefficient of the collision and $u_{in}$ denotes the normal component of the impact velocity. Hence, accounting for cohesive forces using the hard-sphere model involves manipulating $u_{out}$ and introducing thresholds for particle escape, but it prevents quantifying the intergranular stresses or work required for floc breakup.

Recently, various models for cohesive collisions were tested in the framework of a DEM by Thornton, Cummins & Cleary (Reference Thornton, Cummins and Cleary2017). The authors report that piecewise cohesive force models that distinguish between approach and rebound give rise to unphysical properties. In particular, this approach leads to an overestimation of the tensile forces during the rebound process. As a consequence, flocculation may be overestimated because of unphysical sticking conditions at high impact velocities. Hence, studies treating the particles as mass points that collide according to a hard-sphere model are difficult to interpret physically, because this approach does not capture and quantify the forces associated with cohesive particle–particle interaction. Consequently, it becomes impossible to distinguish between the different contributions to particle–particle interactions, such as repulsive collision, tangential friction and lubrication forces as outlined in § 2.2, as repulsive collision forces and cohesive forces are lumped into the inelastic restitution coefficient.

Derksen (Reference Derksen2014) carried out particle-resolving simulations using a finite-size square-well potential to account for flocculation. This square-well potential considers two particles attached as soon as $u_{in}$ falls below a critical threshold and, vice versa, to break up if the escape velocity lies above a critical threshold. For the case of floc breakup, this model converts the effect of cohesion from potential to kinetic energy. This treatment represents an improvement over point-particle approaches, although it does not allow for the space-resolved computation of the cohesive forces and stresses acting on flocculated particles. Similarly, Gu et al. (Reference Gu, Ozel and Sundaresan2016) proposed a cohesive force model that scales inversely with $\unicode[STIX]{x1D701}_{n}$ . We tested this approach and found it to be prone to numerical instabilities when dealing with rather stiff particles, as it introduces large attractive forces for small gap sizes that are discontinuously shut off at a minimal gap size.

In order to computationally simulate realistic cohesive sediment dynamics, we aim to resolve in space and time the following three phases of particle–particle interaction, cf. figure 2: (a) particle approach/flocculation, (b) capture/steady-state contact, and (c) separation in the presence of external forces. In doing so, we will obtain detailed information on the work performed by the interparticle forces. We furthermore aim for a computational model that recovers the original DEM scheme proposed by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) for cohesionless grains. These goals will be achieved by the approach to be described in the following.

Figure 2. Computational scenario for the binary interaction of two cohesive particles.

3.2 Cohesive force model

To exploit the advantages offered by the soft-sphere model of Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) for grain-resolving simulations, we develop an approach that is consistent with the DLVO theory as sketched in figure 1(a), with an attractive interparticle force within the interval $2~\text{nm}\leqslant \unicode[STIX]{x1D701}_{n}\leqslant 10~\text{nm}$ that has a local maximum at $\unicode[STIX]{x1D701}_{n}\approx 4~\text{nm}$ . Note that we do not wish to resolve the layer of $\unicode[STIX]{x1D701}_{n}<2~\text{nm}$ , but instead consider this to be part of the surface roughness. Ideally, the cohesive forces would decay to zero for $\unicode[STIX]{x1D701}_{n}=0$ as the repulsive forces are already accounted for through (2.7). This can be accomplished by the ansatz of a parabolic spring force with the following properties: (i) it decays to zero as the gap size goes to zero; (ii) it has a maximum at a gap width orders of magnitude smaller than the particle diameter; and (iii) it decays to zero for larger gap sizes, without any discontinuous jumps, cf. figure 1(b). These characteristics are incorporated by the mathematically simple model

(3.2) $$\begin{eqnarray}\boldsymbol{F}_{coh}=\left\{\begin{array}{@{}ll@{}}-k_{coh}(\unicode[STIX]{x1D701}_{n}^{2}-\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D706})\boldsymbol{n}\quad & 0<\unicode[STIX]{x1D701}_{n}\leqslant \unicode[STIX]{x1D706},\\ 0\quad & \text{otherwise},\end{array}\right.\end{eqnarray}$$

where $k_{coh}$ denotes the stiffness constant and $\unicode[STIX]{x1D706}$ represents the range over which the cohesive force is smeared. This length scale can be interpreted as a Debye length, which is typically of the order of several micrometres in DEM simulations (Mari et al. Reference Mari, Seto, Morris and Denn2014). As will be shown in § 4.3, the simulation results are insensitive to the exact value of $\unicode[STIX]{x1D706}$ . A reasonable choice that will allow us to resolve the cohesive force computationally, while limiting it to a range much smaller than the particle size, is $D_{50}/\unicode[STIX]{x1D706}=20$ , where $D_{50}$ is the median grain size of a polydisperse ensemble of particles. Particle-resolving simulations typically employ a resolution of 20 grid cells of size $h$ per diameter, so that choosing $\unicode[STIX]{x1D706}\approx h$ and utilizing the substepping routine as proposed by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) guarantees a proper resolution of cohesive effects in space and time. This modelling approach is also consistent with the experimental observations of Delenne et al. (Reference Delenne, El Youssoufi, Cherblanc and Bénet2004), who coated rods of $D=8~\text{mm}$ in diameter with epoxy resin to glue them together. These rods where then put under tension to determine the cohesive forces. It was found in this study that the cohesive force increases with gap size $\unicode[STIX]{x1D701}_{n}$ to a maximum at $D/\unicode[STIX]{x1D701}_{n}\approx 80$ . For larger gap sizes, the force decreases and eventually the rods detach at $D/\unicode[STIX]{x1D701}_{n}\approx 40$ . For the present study, we have chosen the latter value as a reference for the largest particles of the considered polydisperse particle mixtures.

We determine the stiffness $k_{coh}$ of the model by preserving the energy contained in the vdW forces,

(3.3a ) $$\begin{eqnarray}E_{coh}=E_{vdW}.\end{eqnarray}$$

According to the DLVO theory, the vdW forces (3.1) are defined within the interval $\unicode[STIX]{x1D701}_{n}=[\unicode[STIX]{x1D701}_{0},\infty ]$ , where the lower boundary $\unicode[STIX]{x1D701}_{0}=0.2~\text{nm}$ is taken from Israelachvili (Reference Israelachvili1992). On the other hand, (3.2) is defined over the interval $\unicode[STIX]{x1D701}_{n}=[0,\unicode[STIX]{x1D706}]$ . Taking the endpoints of the respective intervals as the integration limits yields

(3.3b ) $$\begin{eqnarray}\int _{0}^{\unicode[STIX]{x1D706}}-k_{coh}(\unicode[STIX]{x1D701}_{n}^{2}-\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D706})\,\text{d}\unicode[STIX]{x1D701}_{n}=\int _{\unicode[STIX]{x1D701}_{0}}^{\infty }\frac{A_{H}R_{eff}}{6\unicode[STIX]{x1D701}_{n}^{2}}\,\text{d}\unicode[STIX]{x1D701}_{n},\end{eqnarray}$$

so that

(3.3c ) $$\begin{eqnarray}k_{coh}=\frac{A_{H}R_{eff}}{\unicode[STIX]{x1D701}_{0}\unicode[STIX]{x1D706}^{3}}.\end{eqnarray}$$

Note that we have replaced $R_{p}$ by $R_{eff}$ on the right-hand side of (3.3b ) to account for polydisperse particle sizes. Substituting this expression for $k_{coh}$ into (3.2) provides the final expression for the dimensional cohesive force model

(3.4) $$\begin{eqnarray}\boldsymbol{F}_{coh}=\left\{\begin{array}{@{}ll@{}}-{\displaystyle \frac{A_{H}R_{eff}}{\unicode[STIX]{x1D701}_{0}\unicode[STIX]{x1D706}^{3}}}(\unicode[STIX]{x1D701}_{n}^{2}-\unicode[STIX]{x1D701}_{n}\unicode[STIX]{x1D706})\boldsymbol{n}\quad & 0<\unicode[STIX]{x1D701}_{n}\leqslant \unicode[STIX]{x1D706},\\ 0\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

This form enables us to account for cohesive forces in our current IBM–DEM framework via an additional force term in the collision model (2.5),

(3.5) $$\begin{eqnarray}\boldsymbol{F}_{c,p}=\mathop{\sum }_{q,q\neq p}^{N_{p}}(\boldsymbol{F}_{l,pq}+\boldsymbol{F}_{n,pq}+\boldsymbol{F}_{t,pq}+\boldsymbol{F}_{coh,pq})+\boldsymbol{F}_{l,pw}+\boldsymbol{F}_{n,pw}+\boldsymbol{F}_{t,pw}+\boldsymbol{F}_{coh,pw}.\end{eqnarray}$$

Equation (3.5) will be the basis of the simulations to be discussed in §§ 4 and 5. Note that this approach treats particles individually and is not modified for the interaction of clusters. Furthermore, no changes to the computation of $\boldsymbol{F}_{h,p}$ and $\boldsymbol{F}_{h,g}$ are necessary, as these are a direct result of the IBM. The main advantage of this approach lies in its ability to resolve the particle bonding process in space and time as cohesive forces act over a finite size shell surrounding each particle. At the same time, it retains the distinction between the individual interparticle force components via (3.5), which allows for an in-depth analysis of the different effects governing the particle–particle interaction. Note that (2.3) and (3.5) together take a form that resembles the minimal flocculation model proposed by Vicsek et al. (Reference Vicsek, Czirók, Ben-Jacob, Cohen and Shochet1995), which was developed for self-propelled organisms. While our particles are passive, their weight $\boldsymbol{F}_{g}$ provides a preferred direction of motion (Toner, Tu & Ramaswamy Reference Toner, Tu and Ramaswamy2005), the collision force $\boldsymbol{F}_{c}$ aligns their motion through the competition of repulsion and cohesion, and the hydrodynamic force $\boldsymbol{F}_{h}$ introduces a forcing that can cause particles to flocculate or to break up.

The dimensional form (3.4) still requires the proper parametrization of the Hamaker constant $A_{H}$ , which implies all of the difficulties mentioned in the introduction. This issue will be addressed via the rescaling to be discussed in § 3.3.

3.3 Non-dimensionalization of the cohesive force model

For the purpose of conducting numerical simulations, we wish to capture the effects of cohesive forces by means of a non-dimensional similarity parameter. To this end, we render the Navier–Stokes equation (2.1) and the particle equation of motion (2.3) dimensionless in appendix B. For the Navier–Stokes equation, the only dimensionless parameter to appear is the Reynolds number (Biegert et al. Reference Biegert, Vowinckel, Ouillon and Meiburg2017b ). For the particle equation of motion (2.3), the cohesive and lubrication forces combined can be viewed as a spring–dashpot system for two interacting particles.

As derived in appendix B, by choosing the buoyancy velocity $u_{s}=\sqrt{g^{\prime }D_{50}}$ (appendix C), the characteristic time scale $\unicode[STIX]{x1D70F}_{s}=D_{50}/u_{s}$ and the characteristic mass $m_{50}=\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x03C0}D_{50}^{3}/6$ , the characteristic force scale for particles settling under gravity in an otherwise quiescent fluid becomes the specific weight $m_{50}g^{\prime }$ . Here, $D_{50}$ is the median diameter of an ensemble of polydisperse particles, $g^{\prime }=(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})g/\unicode[STIX]{x1D70C}_{f}$ denotes the reduced gravity, and $g$ represents the gravitational acceleration. To write the algebraic expression for cohesive forces (3.4) in dimensionless form, we define a cohesive number as

(3.6) $$\begin{eqnarray}\mathit{Co}=\frac{\text{max}(\Vert \boldsymbol{F}_{coh,50}\Vert )}{m_{50}g^{\prime }}.\end{eqnarray}$$

It represents the ratio of the cohesive force maximum for particles of diameter $D_{50}$ to the characteristic gravitational force scale of the problem. A similar characteristic number was used by Sun et al. (Reference Sun, Xiao and Sun2018) in the framework of a point-particle approach. A complete derivation for the origin of $\mathit{Co}$ within the present numerical framework is provided in appendix B.

By design, (3.4) has its maximum at $\unicode[STIX]{x1D701}_{n}=\unicode[STIX]{x1D706}/2$ , so that we immediately obtain, for $R_{eff}=D_{50}/2$ ,

(3.7) $$\begin{eqnarray}\text{max}(\Vert \boldsymbol{F}_{coh,50}\Vert )=-\frac{A_{H}}{\unicode[STIX]{x1D701}_{0}}\frac{D_{50}}{2\unicode[STIX]{x1D706}^{3}}\left(\frac{\unicode[STIX]{x1D706}^{2}}{4}-\frac{\unicode[STIX]{x1D706}^{2}}{2}\right)=\frac{A_{H}}{\unicode[STIX]{x1D701}_{0}}\frac{D_{50}}{8\unicode[STIX]{x1D706}}.\end{eqnarray}$$

After specifying $\text{max}(\Vert \boldsymbol{F}_{coh,50}\Vert )$ for a given problem, we combine (3.7) with (3.4) to obtain the cohesive force as

(3.8) $$\begin{eqnarray}\boldsymbol{F}_{coh}=-\frac{8\,\text{max}(\Vert \boldsymbol{F}_{coh,50}\Vert )}{D_{50}}\frac{R_{eff}}{\unicode[STIX]{x1D706}^{2}}(\unicode[STIX]{x1D701}_{n}^{2}-\unicode[STIX]{x1D706}\unicode[STIX]{x1D701}_{n})\boldsymbol{n}.\end{eqnarray}$$

Owing to the smearing of the cohesive effects over the range $\unicode[STIX]{x1D706}$ , the model now scales with $\boldsymbol{F}_{coh}\propto \unicode[STIX]{x1D706}^{-2}$ rather than $\boldsymbol{F}_{coh}\propto \unicode[STIX]{x1D701}_{0}^{-2}$ . It is important to note, however, that this quantity does not reflect a tunable parameter but a constant property chosen in accordance with the length scales of the physical problem.

It is then convenient to define dimensionless quantities (denoted by tilde) as

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}_{coh}=m_{50}g^{\prime }\,\tilde{\boldsymbol{F}}_{coh}, & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle R_{eff}=D_{50}\tilde{R}_{eff}, & \displaystyle\end{eqnarray}$$

(3.9c ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D706}=D_{50}\tilde{\unicode[STIX]{x1D706}}, & \displaystyle\end{eqnarray}$$

(3.9d ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}_{n}=D_{50}\tilde{\unicode[STIX]{x1D701}}_{n}, & \displaystyle\end{eqnarray}$$

$m_{50}g^{\prime }$

(3.10) $$\begin{eqnarray}\tilde{\boldsymbol{F}}_{coh}=\left\{\begin{array}{@{}ll@{}}-\mathit{Co}\,{\displaystyle \frac{8\,\tilde{R}_{eff}}{\tilde{\unicode[STIX]{x1D706}}^{2}}}(\tilde{\unicode[STIX]{x1D701}}_{n}^{2}-\tilde{\unicode[STIX]{x1D701}}_{n}\tilde{\unicode[STIX]{x1D706}})\boldsymbol{n}\quad & 0<\unicode[STIX]{x1D701}_{n}\leqslant \unicode[STIX]{x1D706},\\ 0\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

The stiffness of our cohesive force model thus becomes $\tilde{k}_{coh}=8\mathit{Co}\tilde{R}_{eff}/\tilde{\unicode[STIX]{x1D706}}^{2}$ , so that the cohesive forces for a given physical system scale linearly with the cohesive number and the effective radius of the two colliding particles, which is consistent with the considerations of Visser (Reference Visser1989), Lick et al. (Reference Lick, Jin and Gailani2004) and Righetti & Lucarelli (Reference Righetti and Lucarelli2007). The characteristics are meant to represent rough macroscopic particles, i.e. $D_{p}>2~\unicode[STIX]{x03BC}\text{m}$ , in saline water. A comprehensive translation of our modelling approach to various physical systems is given in appendix A.

To rewrite the unresolved lubrication forces in dimensionless form, we substitute dimensional quantities in (2.6) by

(3.11a ) $$\begin{eqnarray}\boldsymbol{F}_{lub}=m_{50}g^{\prime }\,\tilde{\boldsymbol{F}}_{lub}=\frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x03C0}D_{50}^{3}}{6}\,g^{\prime }\,\tilde{\boldsymbol{F}}_{lub}\end{eqnarray}$$

and

(3.11b ) $$\begin{eqnarray}\boldsymbol{g}_{n}=u_{s}\tilde{\boldsymbol{g}}_{n}=\sqrt{g^{\prime }D_{50}}\,\tilde{\boldsymbol{g}}_{n}.\end{eqnarray}$$

Combining (2.6) with (3.11) then yields

(3.12) $$\begin{eqnarray}\tilde{\boldsymbol{F}}_{lub}=-36\frac{\unicode[STIX]{x1D708}_{f}}{\sqrt{g^{\prime }D_{50}}\,D_{50}}\frac{\tilde{R}_{eff}^{2}\tilde{\boldsymbol{g}}_{n}}{\max (\tilde{\unicode[STIX]{x1D701}}_{n},\tilde{\unicode[STIX]{x1D701}}_{min})}\end{eqnarray}$$

so that we obtain the dimensionless form

(3.13) $$\begin{eqnarray}\tilde{\boldsymbol{F}}_{lub}=\left\{\begin{array}{@{}ll@{}}-{\displaystyle \frac{36}{Re}}{\displaystyle \frac{\tilde{R}_{eff}^{2}\tilde{\boldsymbol{g}}_{n}}{\max (\tilde{\unicode[STIX]{x1D701}}_{n},\tilde{\unicode[STIX]{x1D701}}_{min})}}\quad & 0<\unicode[STIX]{x1D701}_{n}\leqslant 2h,\\ 0\quad & \text{otherwise}.\end{array}\right.\end{eqnarray}$$

Hence lubrication forces scale with the inverse of the particle Reynolds number $Re=D_{50}u_{s}/\unicode[STIX]{x1D708}_{f}$ , while cohesive forces scale with the cohesive number $\mathit{Co}=\text{max}(\Vert \boldsymbol{F}_{coh,50}\Vert )/(m_{50}g^{\prime }).$ By quantifying these two dimensionless similarity parameters, the physical system is thus fully specified.

3.4 Cohesive force model validation

The key idea behind the cohesive number introduced in (3.6) in § 3.3 is to define a critical threshold of $\mathit{Co}=1$ , for which the cohesive forces balance the specific weight of a particle of size $D_{50}$ . To test and verify this behaviour, we consider a simple test case of two particles in quiescent fluid. Particle $p$ is held fixed, and particle $q$ is placed right below it (cf. figure 2). Particle $q$ wants to settle as a result of gravity, whereas the cohesive force acts to keep it attached to particle $p$ . The governing Reynolds number is based on the quantities introduced in § 3.3, which yields $Re=u_{s}D/\unicode[STIX]{x1D708}_{f}$ . We chose the Reynolds number for the test case from the experimental observations of Lick et al. (Reference Lick, Jin and Gailani2004), who reported that cohesive forces start to alter the erosion behaviour of silicate particles with a density of $\unicode[STIX]{x1D70C}_{p}=2650~\text{kg}~\text{m}^{-3}$ and a diameter in the range of $140~\unicode[STIX]{x03BC}\text{m}\leqslant D_{p}\leqslant 390~\unicode[STIX]{x03BC}\text{m}$ , which are submerged in water with a density and kinematic viscosity of $\unicode[STIX]{x1D70C}_{f}=1000~\text{kg}~\text{m}^{-3}$ and $\unicode[STIX]{x1D708}_{f}=10^{-6}~\text{m}^{2}~\text{s}^{-1}$ , respectively. Choosing a representative value of $D_{p}=242~\unicode[STIX]{x03BC}\text{m}$ yields $Re=u_{s}D/\unicode[STIX]{x1D708}_{f}=15.1$ . Since only two particles are involved, we define $D_{50}=(D_{p}+D_{q})/2$ . Initially, particle $q$ is at rest and the size of the vertical gap with respect to the fixed particle $p$ is set to $\unicode[STIX]{x1D701}_{n}=\unicode[STIX]{x1D706}/2$ , where we chose $D_{p}/\unicode[STIX]{x1D706}=20$ . Particle $q$ is then free to move. Note that, for this case of particle $q$ interacting with the fixed particle $p$ , the effective radius becomes $R_{eff}=R_{q}$ . The median particle size $D_{50}$ is discretized by 20 grid cells per diameter. In order to investigate whether or not particle $q$ will detach from particle $p$ , a relatively small computational domain of $L_{x}\times L_{y}\times L_{z}=2.5D_{50}\times 5D_{50}\times 2.5D_{50}$ suffices, with gravity acting in the negative $y$ -direction. No-slip walls are imposed in the $y$ -direction, whereas the $x$ - and $z$ -directions are being treated as periodic.

For $\mathit{Co}<1$ , the weight of particle $q$ is larger than the maximum of the cohesive force at $\unicode[STIX]{x1D701}_{n}=\unicode[STIX]{x1D706}/2$ , so that we expect particle $q$ to detach from particle $p$ . This is confirmed by figure 3(a), which shows the gap size versus non-dimensional time $t/\unicode[STIX]{x1D70F}_{s}$ . In addition, the detachment happens more slowly as the cohesive number approaches its critical value. For the critical cohesive number $\mathit{Co}=1$ , the particles remain stationary at a constant gap size $\unicode[STIX]{x1D701}_{n}=\unicode[STIX]{x1D706}/2$ . Increasing the cohesive number even further causes particle $q$ to move closer towards particle $p$ , as the maximum of the attractive forces increases. Since the cohesive force approaches zero as the gap decreases, particle $q$ will find an equilibrium position within the interval $0\leqslant \unicode[STIX]{x1D701}_{n}\leqslant \unicode[STIX]{x1D706}/2$ at which the cohesive force is balanced by the weight of the particle. Also note that, even if we do not explicitly resolve the nanoscale over which cohesive forces typically act, the smearing of the cohesive potential over $\unicode[STIX]{x1D706}$ bonds particles at a constant gap size of $\unicode[STIX]{x1D701}_{n}\leqslant D_{p}/40$ , which is sufficiently small to reproduce physically realistic behaviour. Particles can remain a finite distance apart during flocculation under the influence of divergent forces (Israelachvili Reference Israelachvili1992; Thornton et al. Reference Thornton, Cummins and Cleary2017), and they can experience friction while in direct contact during collisions.

Figure 3. Gap size versus time: (a) different runs with varying cohesive number; (b) different runs with varying ratio of the two particle radii with $\mathit{Co}=1$ .

Corresponding tests can be carried out for polydisperse particles. Here we set the cohesive number to $\mathit{Co}=1$ while varying the ratio of the two particle radii, $R_{p}/R_{q}$ . The median grain size $D_{50}=(D_{p}+D_{q})/2$ serves a basis for calculating the similarity parameters $Re$ and $\mathit{Co}$ . The results are shown in figure 3(b). If the radius ratio is smaller than unity, particle $q$ is larger than particle $p$ and its weight causes particle $q$ to detach. On the other hand, if particle $q$ is smaller than particle $p$ they stay in contact at a constant gap size $\unicode[STIX]{x1D701}_{n}<\unicode[STIX]{x1D706}/2$ . This test case validates the arguments underlying (3.10) for polydisperse particles.

The above analysis demonstrates that our computational model allows for the precise control of cohesive forces. Since Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017a ) showed that the present computational approach also yields excellent agreement with experimental data for cohesionless particles, we expect it to reproduce the settling dynamics of cohesive particles with high fidelity.

4.1 Drafting, kissing, tumbling

Figure 4. Cohesionless particles undergoing DKT: (a) kissing at $t/\unicode[STIX]{x1D70F}_{s}=10$ , and (b) tumbling at $t/\unicode[STIX]{x1D70F}_{s}=50$ . Contours show the downward fluid velocity component.

In order to assess the influence of cohesive forces under simplified conditions, we focus on the classical DKT experiment of two particles settling under gravity, which has been explored in depth for cohesionless grains (e.g. Fortes et al. Reference Fortes, Joseph and Lundgren1987; Glowinski et al. Reference Glowinski, Pan, Hesla, Joseph and Priaux2001). The initial configuration is shown in figure 2. Two particles with a density greater than that of the ambient fluid are placed above each other in a tank of quiescent fluid, with an initial gap size that is substantially larger than the distance of the short-range lubrication and cohesive forces. As the particles are released and settle under the influence of gravity, the trailing particle is drafted by the wake of the leading particle, so that it experiences reduced drag. The two particles touch (or kiss, figure 4 a) and subsequently rearrange themselves (tumble) into a side-by-side configuration. In the absence of cohesive effects, hydrodynamic forces eventually push them apart, and they start to separate laterally (figure 4 b).

Figure 5. DKT results for equal-sized spheres with different cohesive numbers showing separation width $\unicode[STIX]{x1D701}_{n}$ for (a) drafting phase, (b) tumbling phase, and (c) steady state.

To investigate the DKT scenario for cohesive particles, we employ the same physical set-up and the same numerical parameters as described in § 3.4, i.e. $Re=15.1$ , $D_{p}/h=20$ and $D_{50}/\unicode[STIX]{x1D706}=20$ . Since we are dealing with grains that should resemble the characteristic of natural silt, the Reynolds number of our system is substantially lower compared to the studies of Fortes et al. (Reference Fortes, Joseph and Lundgren1987) and Glowinski et al. (Reference Glowinski, Pan, Hesla, Joseph and Priaux2001). Decreasing the Reynolds number results in slower dynamics of the interacting particles and increases the time scales to observe the DKT motion. Hence, we choose a rather long computational domain of $L_{x}\times L_{y}\times L_{z}=120D_{50}\times 5D_{50}\times 5D_{50}$ , with gravity acting in the negative $x$ -direction. The boundary conditions assume periodicity in $x$ , and free slip in $y$ and $z$ , and everything is at rest initially. Owing to the low Reynolds number, the domain length in the $x$ -direction is sufficiently large for particles to establish a steady-state configuration. We begin by considering two equal-sized spheres that are initially placed at $\boldsymbol{x}_{p}=(x_{q}+D_{50}+2h,y_{q}+0.5h,z_{q})^{\text{T}}$ , which is sufficiently close to trigger the DKT behaviour, but far enough apart for lubrication and cohesive forces to be unimportant initially. We remark that the simulation results do not depend on the exact initial conditions, as long as $\unicode[STIX]{x1D701}_{n}\leqslant D_{p}$ , so that DKT is initiated. The initial horizontal offset of $0.5h$ in the $y$ -coordinate triggers the physical instability, leading to the particle rearrangement during the kissing phase.

The results for the monodisperse case are presented in figure 5(a,b), which show the gap width $\unicode[STIX]{x1D701}_{n}$ of the two particles as a function of time. Consistent with the classical observations by Fortes et al. (Reference Fortes, Joseph and Lundgren1987), for $\mathit{Co}=0$ the cohesionless particles approach each other (figure 5 a), touch for approximately 20 time units and then separate in a tumbling behaviour (figure 5 b). Increasing the cohesive forces speeds up the drafting phase, and slows down or prevents the subsequent separation of the particles. The fact that the particles remain in steady-state contact for $\mathit{Co}<1$ indicates that hydrodynamic forces are not as effective in pulling them apart as gravity was in § 3.4. Figure 5(c) demonstrates that already a cohesive number value of $\mathit{Co}=0.35$ suffices to maintain the steady-state bond between the particles.

Figure 6. DKT results for spheres of different size with $\mathit{Co}=1$ and different ratios of $R_{p}/R_{q}$ : (a) $R_{p}/R_{q}=0.6$ , (b) $R_{p}/R_{q}=0.7$ , (c) $R_{p}/R_{q}=1$ and (d)  $R_{p}/R_{q}=4$ . Circles indicate the initial and final particle configuration, respectively.

4.2 Settling of cohesive particles of different size

To investigate the impact of polydispersity on the settling behaviour of cohesive particles, we repeat the above simulations for a fixed cohesive number value of $\mathit{Co}=1$ , while varying the ratio of the particle radii $R_{p}/R_{q}$ in the range $0.25\leqslant R_{p}/R_{q}\leqslant 4$ . For $R_{p}/R_{q}<1$ , the lower particle is larger and tends to settle faster than the upper, trailing particle. However, already for a ratio of $R_{p}/R_{q}=0.6$ , the wake of the leading particle is sufficiently strong to draft the trailing particle into the DKT motion. Corresponding behaviour was also reported for two-dimensional simulations of cohesionless settling circular disks by Wang, Guo & Mi (Reference Wang, Guo and Mi2014), who found that there exists a critical value for $R_{p}/R_{q}$ to initiate DKT for a smaller particle trailing a larger one. To analyse the particle positions relative to each other, we define the vector connecting the particle centres as $\boldsymbol{r}_{pq}=\boldsymbol{x}_{p}-\boldsymbol{x}_{q}$ . We then plot the distance between these two centres, as well as the horizontal and vertical separation components,

(4.1a ) $$\begin{eqnarray}\displaystyle & \displaystyle \Vert \boldsymbol{r}_{pq}\Vert =\unicode[STIX]{x1D701}_{n}+D_{50}=\sqrt{(x_{p}-x_{q})^{2}+(y_{p}-y_{q})^{2}+(z_{p}-z_{q})^{2}}, & \displaystyle\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D701}_{h}=\sqrt{(y_{p}-y_{q})^{2}+(z_{p}-z_{q})^{2}}, & \displaystyle\end{eqnarray}$$

and

(4.1c ) $$\begin{eqnarray}\unicode[STIX]{x1D701}_{v}=\sqrt{(x_{p}-x_{q})^{2}},\end{eqnarray}$$

respectively. Note that particles touch when $\Vert \boldsymbol{r}_{pq}\Vert /D_{50}=1$ . The temporal evolution of these quantities is shown in figure 6. For the ratio $R_{p}/R_{q}=0.6$ , which was found to be the approximate threshold for drafting, the particles initially separate but then quickly approach each other and touch. The horizontal distance between the particle centres is seen to increase slowly throughout the simulation (figure 6 a), which suggests that the pair rotates into an oblique configuration, although the simulation time is too short for a quasi-steady state to be reached. For $R_{p}/R_{q}=0.7$ this rotation occurs more rapidly, and a quasi-steady oblique configuration emerges (figure 6 b). When the particle radii ratio is close to, but not equal to, unity, we are in the regime for which cohesionless particles repeatedly undergo DKT interactions as reported by Shao, Liu & Yu (Reference Shao, Liu and Yu2005).

For equal-sized particles the vertical distance between the centres decays to zero, while their horizontal distance approaches the particle diameter (figure 6 c), indicating that the particles align horizontally while touching. This observation is consistent with the classical kissing behaviour of the monodisperse case described in the literature (Fortes et al. Reference Fortes, Joseph and Lundgren1987; Glowinski et al. Reference Glowinski, Pan, Hesla, Joseph and Priaux2001). Increasing the ratio even further, so that the trailing particle $p$ becomes larger than the leading particle $q$ , causes the two particles to swap positions by rotating around each other, since particle $p$ has a bigger settling velocity than particle $q$ (figure 6 d). They subsequently align approximately vertically with $\unicode[STIX]{x1D701}_{v}\approx \Vert \boldsymbol{r}_{pq}\Vert$ , and only a small horizontal separation width. Corresponding observations of larger trailing particles swapping positions with smaller leading ones were reported for cohesionless grains in both two- and three-dimensional simulations (Wang et al. Reference Wang, Guo and Mi2014; Liao et al. Reference Liao, Hsiao, Lin and Lin2015). It is interesting to note that, even though particles do not bond for $R_{p}/R_{q}\leqslant 0.5$ , they do form a lasting bond for the more disparate size ratio of $R_{p}/R_{q}=4$ via this swapping mechanism. This demonstrates that the ability of the particles to form flocs strongly depends on their initial configuration. The quasi-steady settling configuration of cohesive particle pairs undergoing DKT is sketched in figure 7. Particle pairs of very different sizes tend to align vertically (figure 7 d), while the alignment becomes increasingly horizontal as the particle sizes approach each other.

Figure 7. Sketch of the asymptotically steady settling configuration of polydisperse, cohesive particles, as illustrated in figure 6: (a) $R_{p}/R_{q}=0.6$ , (b) $R_{p}/R_{q}=0.7$ , (c)  $R_{p}/R_{q}=1$ and (d) $R_{p}/R_{q}=4$ . The colour coding of the connecting lines corresponds to figure 6.

Figure 8. Settling velocity in the DKT scenario as a function of time for monodisperse cohesionless and polydisperse cohesive particles as shown in figures 6 and 7.

The mechanism governing the particle orientation has important implications for the settling speeds of the interacting particles. Figure 8 compares the situations addressed in figures 6 and 7 by showing the settling velocity $u_{pq}=0.5(u_{p}+u_{q})$ of the two interacting particles $p$ and $q$ as a function of time. Note that the settling velocity is normalized by $u_{s}$ , which is identical for all of the cases shown here since we only change the ratio of $R_{p}/R_{q}$ but not the sum $R_{p}+R_{q}$ . This allows for a direct comparison of the settling velocities for the different radii ratios. Figure 8 shows only the first passage of the two particles through the periodic domain, in order to exclude any effects from perturbations that might be caused by the particle wake flows. Owing to the rather long domain employed for this study, this is equivalent to more than $140\unicode[STIX]{x1D70F}_{s}$ even for the fastest-settling particles. As soon as the monodisperse, cohesionless grains tumble apart (at $t/\unicode[STIX]{x1D70F}_{s}\approx 40$ ), their settling speed decreases. As expected, introducing cohesive forces increases the settling velocity, as the particles remain in contact, which reduces their total drag. This is true for both monodisperse and polydisperse particle configurations. For polydisperse particles, the settling speed is ultimately governed by the larger, leading particle. The kink for $R_{p}/R_{q}=4$ at $t/\unicode[STIX]{x1D70F}_{s}\approx 11$ reflects the situation of the two particles swapping their leading/trailing configuration. Subsequently, these two particles acquire the largest settling velocity among all of the cases presented here. As $R_{p}/R_{q}$ approaches unity, the settling velocity decreases. These results clearly illustrate the effect of polydispersity on particle settling speeds.

Figure 9. Quasi-steady configuration of cohesive particle pairs undergoing DKT for different particle radii ratios $R_{p}/R_{q}$ : (a) orientation angle $\unicode[STIX]{x1D703}$ ( $\unicode[STIX]{x1D703}=90^{\circ }$ , horizontal alignment; $\unicode[STIX]{x1D703}=0^{\circ }$ , vertical alignment), and (b) settling speed. The dotted horizontal line in (b) indicates the undisturbed settling velocity.

The relationship between the quasi-steady particle alignment and the settling speed is displayed in figure 9 for the parameter range $1\leqslant R_{p}/R_{q}\leqslant 4$ . When this ratio exceeds 2, the particles are aligned approximately vertically, while oblique configurations are observed for smaller ratios, as indicated by the orientation angle $\cos \unicode[STIX]{x1D703}=\unicode[STIX]{x1D701}_{v}/\Vert \boldsymbol{r}_{pq}\Vert$ . The orientation angle is seen to decrease approximately linearly between $1\leqslant R_{p}/R_{q}\leqslant 2$ .

To investigate whether the particle settling is accelerated, we estimate the undisturbed settling velocity $u_{ra}$ by using Rayleigh’s drag equation (appendix C). Figure 9(b) displays the settling velocities $u_{p}$ and $u_{q}$ normalized by their respective undisturbed settling velocity $u_{ra}$ . Equal-sized particles, i.e. $R_{p}/R_{q}=1$ , are seen to settle with a velocity that is slightly higher than their undisturbed settling velocity. For increasing ratios $R_{p}/R_{q}$ , we find that the settling of the smaller particle $q$ is substantially accelerated by the stable bond, whereas the larger particle $p$ still settles approximately with its undisturbed settling velocity $u_{ra}$ . Hence the centre of mass of a cohesive particle pair with a stable bond settles more rapidly than that of two cohesionless particles.

4.3 Sensitivity of cohesive force range

Figure 10. Influence of the cohesive force range $\unicode[STIX]{x1D706}$ on the settling velocity of two interacting particles: (a) $R_{p}/R_{q}=1$ and (b) $R_{p}/R_{q}=4$ .

While the model proposed in § 3 replaces the empirical constants $A_{H}$ and $\unicode[STIX]{x1D701}_{0}$ by the physically meaningful dimensionless cohesive number $\mathit{Co}$ , it still contains the dependence on the cohesive force range $\unicode[STIX]{x1D706}$ . To test the sensitivity of the simulation results on this parameter, we conducted simulations with $R_{p}/R_{q}=1$ and $R_{p}/R_{q}=4$ and different values of $\unicode[STIX]{x1D706}$ . The results are shown in figure 10. For the equal-sized particles, the settling velocities for all values of $\unicode[STIX]{x1D706}$ collapse onto a single curve. The same is true for $R_{p}/R_{q}=4$ and $\unicode[STIX]{x1D706}\leqslant 3h$ . However, for $\unicode[STIX]{x1D706}=4h$ the particles start to separate at $t/\unicode[STIX]{x1D70F}_{s}\approx 20$ . The reason for this detachment is that $\unicode[STIX]{x1D706}$ is now equal to the radius $R_{q}$ of the smaller particle, so that the range of the cohesive force is no longer much smaller than the particle radius, which violates the assumptions underlying the model of § 3. In summary, we find that, as long as $\unicode[STIX]{x1D706}$ is significantly smaller than the smallest particle radius, the simulation results are independent of $\unicode[STIX]{x1D706}$ . The analysis presented in this section also provides evidence that decreasing the steady-state separation distance by increasing the cohesive number does not affect the hydrodynamics of the two interacting particles.

5.1 Computational set-up

In order to explore the influence of cohesive forces on the sedimentation process of a large, polydisperse ensemble of particles, we reproduce the experiments by te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015), albeit on a smaller spatial scale. These authors investigated the hindered settling of silt particles, with diameters in the range $2~\unicode[STIX]{x03BC}\text{m}\leqslant D_{p}\leqslant 63~\unicode[STIX]{x03BC}\text{m}$ . To this end, we place a polydisperse mixture with a homogenous particle volume fraction of $15\,\%$ in a tank of quiescent fluid (figure 11 a). As before, it is convenient to define the reference velocity based on the buoyancy velocity $u_{s}=\sqrt{g^{\prime }D_{50}}$ , where $D_{50}$ denotes the median grain size of the entire particle size distribution. The characteristic time scale based on the buoyancy velocity and the median diameter then becomes $\unicode[STIX]{x1D70F}_{s}=D_{50}/u_{s}$ .

Consistent with the experiments, we choose a Reynolds number of $Re=1.35$ . As in the experiments of te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015), the computational grain sizes obey a cumulative log-normal distribution $(1/2)+(1/2)\operatorname{erf}[(\ln (D_{p}-\unicode[STIX]{x1D707}))/\sqrt{2}\,\unicode[STIX]{x1D70E}]$ around the median diameter $D_{50}$ , with the arithmetic moments $\unicode[STIX]{x1D707}=-1.33$ and $\unicode[STIX]{x1D70E}=0.34$ (figure 11 c). This yields a total of 1261 particles with maximum size ratio of $\max \{D\}/\text{min}\{D\}=4$ , which is smaller than it was in the experiments of te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015). Our computational approach of particle-resolved simulations, however, requires us to resolve even the smallest grain size by at least eight grid cells per diameter (Uhlmann Reference Uhlmann2005). Hence, it was concluded from § 4 that a size ratio of 4 is plenty to account for polydispersity but remain computationally feasible. The small deviations from the analytical log-normal distribution stem from the fact that we slightly rearranged the particle distribution due to the initial random particle placement. This became necessary as the rather small computational domain yielded large variations in volume fraction over the vertical extent of the domain. To smooth out the horizontally averaged volume fraction profile, we applied a two-step procedure: first we removed larger particles from $y$ -locations with higher concentrations, and subsequently we replaced them with exactly the same volume of a few smaller particles in $y$ -locations with lower concentrations at random $x$ - and $z$ -positions. This procedure yields an almost uniform particle volume fraction $\unicode[STIX]{x1D719}_{v}=V_{p}/V_{0}\approx 0.155$ (figure 11 b), where $V_{p}$ and $V_{0}$ denote the volume occupied by the particles and the computational domain size, respectively. Note that this rearrangement of particles would not have been required for a much larger tank and many more particles, but this would have been prohibitively expensive computationally. The computational domain is of size $L_{x}\times L_{y}\times L_{z}=13.1D_{50}\times 40.0D_{50}\times 13.1D_{50}$ , with gravity pointing in the negative $y$ -direction. We assume periodic boundary conditions in the $x$ - and $z$ -directions, respectively, along with a no-slip condition at the bottom wall and a free-slip condition at the top wall. The median particle size is discretized by $D_{50}/h=18.25$ grid cells.

Figure 11. (a) Initial particle distribution, (b) initial particle volume fraction distribution, and (c) cumulative distribution function of the particle diameter, along with the log-normal distribution.

Table 1. Parameters for simulations of a large ensemble, where $L_{x}$ , $L_{y}$ , $L_{z}$ indicate the domain size, and $s$ represents the standard deviation of the grain size.

Three simulations were performed for different values of the cohesion number $\mathit{Co}$ : (i) cohesionless grains with a cohesive number $\mathit{Co}=\max (\Vert \boldsymbol{F}_{coh,50}\Vert )/m_{50}g^{\prime }=0$ , (ii) mildly cohesive sediment with $\mathit{Co}=1$ , and (iii) strongly cohesive sediment with $\mathit{Co}=5$ . For all simulations, the particles are released from rest in quiescent fluid, and subsequently settle under the influence of gravity. The key simulation parameters are listed in table 1. The particle collisions are inelastic with $e_{dry}=0.97<1$ , and they experience friction through (2.8).

5.2 Hindered settling behaviour

The impact of cohesive forces on the settling behaviour is illustrated by figure 12. During the early stages, the particle distributions are very similar for all three simulations. Over the course of the simulations, however, the cohesive sediment is seen to settle faster than its non-cohesive counterpart. This qualitative observation is confirmed by the concentration profiles of figure 13. At $t=17.6\unicode[STIX]{x1D70F}_{s}$ , when the particle phase has its maximum kinetic energy (cf. § 5.3), the profiles for all three simulations remain nearly identical, as cohesive forces have not yet had sufficient time to cause a noticeable change. As time progresses, the concentration profiles remain very similar in the dilute region near the top of the tank, where the volume fraction remains below $5\,\%$ so that particle–particle interactions are negligible (Capart & Fraccarollo Reference Capart and Fraccarollo2011), cf. figure 13(b). In the lower part of the tank, differences begin to emerge, as cohesive forces result in the formation of flocs with larger settling speeds, so that particles accumulate at the bottom of the tank more quickly. Also, note that the undulations in the profiles are milder for larger cohesive forces. At the final simulation time (figure 13 c), cohesive sediment has a lower volume fraction at the very bottom of the tank at ( $0\leqslant y\leqslant 2D_{50}$ ), as compared to cohesionless grains. This reflects the impact of cohesive forces on the consolidation process, as larger cohesive forces yield stable flocs, whereas cohesionless sediment rearranges itself into a denser configuration under the weight of the overlying sediment (Been & Sills Reference Been and Sills1981). Above the dense layer of sediment at the very bottom of the tank ( $0\leqslant y\leqslant 5D_{50}$ ), there exists a layer of loosely flocculated particles with a lower volume fraction, which increases in depth over time. This observation holds for both $\mathit{Co}=1$ and $\mathit{Co}=5$ , and is consistent with experimental observations of freshly sedimented flocs with larger pore spaces above older sediments (e.g. Winterwerp Reference Winterwerp2001). At the final simulation time, we observe a sharp decrease in particle volume fraction near $y\approx 10D_{50}$ for the cohesive sediment, which is more pronounced for the simulation with larger cohesive forces. The cohesionless sediment, on the other hand, shows higher volume fractions at $y>10D_{50}$ , as a result of the unfinished settling process.

As the particles settle, they replace fluid at the bottom of the tank and generate an upward counterflow. For the current simulations, this counterflow is sufficiently strong to sweep smaller particles upwards. It represents one reason for hindered settling and for the separation of the grain sizes for very large water columns (te Slaa et al. Reference te Slaa, van Maren, He and Winterwerp2015). This effect is illustrated in figure 14, which shows the final volume fraction profiles for the smallest, intermediate and largest third of the particles. Figure 14 demonstrates that small and intermediate cohesive particles settle much more rapidly than their non-cohesive counterparts, consistent with the observation of Lick et al. (Reference Lick, Jin and Gailani2004), who found intermediate-sized particles to be most strongly affected by cohesive forces. These types of grains have their peak concentration in the interval $3D_{50}<y<10D_{50}$ . In contrast, large cohesive grains have a lower volume fraction near the bottom of the tank than large cohesionless grains (figure 14 c). As a consequence, the effect of size segregation is less pronounced for cohesive grains, as flocculation results in particles of different sizes settling with the same velocity (Mehta et al. Reference Mehta, Hayter, Parker, Krone and Teeter1989).

Figure 12. Particle configurations during the settling process. (a–f)  $Co=0$ , (g–l) $Co=1$ , (m–r)  $Co=5$ . Left column: $t=17.6\unicode[STIX]{x1D70F}_{s}$ , which corresponds to the time at which the particle phase has its maximum kinetic energy. From left to right, the columns are separated by time intervals of $72.5\unicode[STIX]{x1D70F}_{s}$ . The grey shading reflects the vertical particle velocity. The cohesive sediment is seen to settle more rapidly than its non-cohesive counterpart (see also supplementary movie available at https://doi.org/10.1017/jfm.2018.757).

Figure 13. Horizontally averaged particle volume fractions during the settling process: (a) at $t=17.6\unicode[STIX]{x1D70F}_{s}$ , (b) at $t=162.6\unicode[STIX]{x1D70F}_{s}$ and (c) at $t=380.1\unicode[STIX]{x1D70F}_{s}$ .

Figure 14. Particle volume fraction profile for different particle radii at $t=380.1\unicode[STIX]{x1D70F}_{s}$ : (a) small particles with $D_{p}\leqslant D_{33}$ , (b) medium-sized particles in the range $D_{33}<D_{p}\leqslant D_{66}$ , and (c) large particles with $D>D_{66}$ . Note the different horizontal axis scalings for the individual panels. The results in (a) and (b) were smoothed by a moving average with filter width of $1.5D_{50}$ for clarity.

To validate our simulation results, we compare the effect of cohesive forces on the settling of our polydisperse particles to established empirical relations describing hindered settling of silt. We can compute the settling behaviour of the particle phase in a double-averaged sense (Vowinckel et al. Reference Vowinckel, Nikora, Kempe and Fröhlich2017a ). To this end, we apply an averaging operator to our Eulerian fluid grid that evaluates instantaneous snapshots of the particle velocity distribution $(\unicode[STIX]{x1D719}v_{p})$ , where $\unicode[STIX]{x1D719}$ is the part of a cell occupied by solids and $v_{p}$ is the settling velocity of the particle taking up this space. Note that we assume rigid-body motion and zero rotation for this analysis. This yields

(5.1) $$\begin{eqnarray}\langle v_{p}\rangle (y,t)=\frac{\displaystyle \int _{0}^{L_{z}}\int _{0}^{L_{x}}\unicode[STIX]{x1D719}(x,y,z,t)\,v_{p}(x,y,z,t)\,\text{d}x\,\text{d}z\,\text{d}t}{\displaystyle \int _{0}^{L_{z}}\int _{0}^{L_{x}}\unicode[STIX]{x1D719}(x,y,z,t)\,\text{d}x\,\text{d}z\,\text{d}t},\end{eqnarray}$$

where the angular brackets denote horizontal averaging. These data are then evaluated for binned values of $\unicode[STIX]{x1D719}$ using

(5.2) $$\begin{eqnarray}\overline{\langle v_{p}\rangle }(\unicode[STIX]{x1D719})=\frac{1}{N_{t}N(\unicode[STIX]{x1D719})}\mathop{\sum }_{N_{t}}\mathop{\sum }_{N(\unicode[STIX]{x1D719})}\langle v_{p}\rangle (y,t),\end{eqnarray}$$

where $N(\unicode[STIX]{x1D719})$ is the number of samples recorded for a given $\unicode[STIX]{x1D719}$ , and $N_{t}$ is the number of the evaluated datasets in time outputted with an interval of $\unicode[STIX]{x1D6E5}_{t}=0.25\unicode[STIX]{x1D70F}_{s}$ . Time averaging indicated by the overbar is performed from $t_{s}=240\unicode[STIX]{x1D70F}_{s}$ to $t_{e}=380\unicode[STIX]{x1D70F}_{s}$ , which is the end of the simulation as displayed in figure 13(c). The averaging time interval was chosen to start well after the initial stir-up phase, as will be discussed in detail below. The sample size hence comprises $N_{t}=560$ datasets over a total time of $140\unicode[STIX]{x1D70F}_{s}$ , which is long enough to obtain statistically meaningful data of the well-developed settling behaviour. Similarly, we evaluate the undisturbed settling velocity $v_{ra}$ (cf. appendix C) for the instantaneous particle distribution $\unicode[STIX]{x1D719}$ to compute $\overline{\langle v_{ra}\rangle }$ . As a result, the ratio $\overline{\langle v_{p}\rangle }/\overline{\langle v_{ra}\rangle }$ is still a function of the volume fraction and can therefore be used to compare with the hindered settling functions available in the literature.

One of the first hindered settling functions was proposed by Richardson & Zaki (Reference Richardson and Zaki1954),

(5.3) $$\begin{eqnarray}\frac{\overline{\langle v_{p}\rangle }}{\,\overline{\langle v_{ra}\rangle }\,}=\left(1-\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{s}}\right)^{n},\end{eqnarray}$$

where $\unicode[STIX]{x1D719}_{s}$ and $n$ are empirical parameters describing the volume fraction of a freshly deposited sediment bed and the particle size and shape, respectively. As argued by Dankers & Winterwerp (Reference Dankers and Winterwerp2007), cohesive sediment such as mud with a significant amount of clay deposits at the bottom in a gel-like structure with a volume fraction that is lower than the maximum possible volume fraction $\unicode[STIX]{x1D719}_{max}$ . Hence, these authors define $\unicode[STIX]{x1D719}_{s}<\unicode[STIX]{x1D719}_{max}$ based on measurements of the settling velocity of mud to separate the effects of hindered settling from the consolidation of the deposited sediment. Owing to its simplicity, equation (5.3) has been very popular in hydraulic engineering. However, as described by te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015), this function is known to underestimate the settling velocity for higher concentrations. Instead, these authors have used the hindered settling function of Winterwerp (Reference Winterwerp2002),

(5.4) $$\begin{eqnarray}\frac{\overline{\langle v_{p}\rangle }}{\,\overline{\langle v_{ra}\rangle }\,}=\frac{\left(1-{\displaystyle \frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{s}}}\right)^{m}(1-\unicode[STIX]{x1D719})}{\left(1-{\displaystyle \frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{max}}}\right)^{-5\unicode[STIX]{x1D719}_{max}/2}},\end{eqnarray}$$

where the numerator represents the effects of the counterflow and the increased buoyancy, respectively, while the denominator reflects the increased viscosity of dense suspensions according to Krieger & Dougherty (Reference Krieger and Dougherty1959). Here, $m$ is an empirical exponent, which plays a similar role to the parameter $n$ in (5.3). It was shown by te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015) that equations (5.3) and (5.4) both yield good agreement with experimental results for settling coarse silt. The best agreement for (5.3) was shown for the parameter values $n=4.91$ and $\unicode[STIX]{x1D719}_{s}=1$ , while (5.4) performs best with $m=1$ , $\unicode[STIX]{x1D719}_{s}=0.5$ and $\unicode[STIX]{x1D719}_{max}=0.65$ .

Figure 15. Double-averaged settling velocity normalized with the undisturbed settling velocity. Comparison to empirical relationships (5.3) of Richardson & Zaki (Reference Richardson and Zaki1954) (RZ) and (5.4) of Winterwerp (Reference Winterwerp2002) (W). (a) Parametrization according to te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015). (b) Same parametrization except for choosing $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x1D719}_{max}=0.7$ according to the simulation data of figure 13(c).

We compare our double-averaged simulation results obtained from (5.2) to the hindered settling functions as parametrized by te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015) in figure 15(a) to validate our simulations. Our results agree well with the two empirical hindered settling functions, and they demonstrate the enhanced settling velocity due to the cohesive forces for all concentration values. Unlike the hindered settling functions, the simulated settling velocities do not approach the undisturbed value for very low volume fractions. This can be attributed to the finite size of the computational domain and the limited number of particles. As can be seen in figure 13, these low volume fractions are typically found in the top part of the tank, where the smallest particles are still accelerating with a very low Reynolds number. Nevertheless, the agreement is remarkable, which demonstrates the ability of the current simulation approach to produce physically realistic results.

Surprisingly, the hindered settling function of Richardson & Zaki (Reference Richardson and Zaki1954) as parametrized by te Slaa et al. (Reference te Slaa, van Maren, He and Winterwerp2015) does not underestimate the settling velocities of higher volume fractions, and the agreement of (5.3) with our data seems to be even better than for (5.4). This is because the parameter $\unicode[STIX]{x1D719}_{s}=1$ was calibrated for best fit to match experimental results. Hence, figure 15(a) serves as validation of our simulation approach. However, this parametrization is not in line with the definition of $\unicode[STIX]{x1D719}_{s}$ as given by Dankers & Winterwerp (Reference Dankers and Winterwerp2007). Moreover, it is immediately obvious that the solid content of a freshly deposited sediment bed cannot reach this value. Hence, to improve the parametrization of (5.3) and (5.4), we propose to choose $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x1D719}_{max}$ since we do not deal with mud but with coarse silt. The consolidation of the sediment will continue to squeeze out water from the bed until the particle packing jams. This process will maintain a counterflow over very long time scales (e.g. Houssais et al. Reference Houssais, Ortiz, Durian and Jerolmack2016). Based on these observations, we propose to not calibrate  $\unicode[STIX]{x1D719}_{s}$ but to parametrize critical volumetric concentrations by the maximum value of our concentration data as shown in figure 13(c). Using these data, we can immediately parametrize $\unicode[STIX]{x1D719}_{s}=\unicode[STIX]{x1D719}_{max}=0.7$ , which is in line with experimental and computational studies of polydisperse particle packings (Sohn & Moreland Reference Sohn and Moreland1968; Desmond & Weeks Reference Desmond and Weeks2014). Note that this value is larger than the volume fraction of a randomly closed packing of monodisperse spheres, since we are dealing with polydisperse particles and the operator (5.1) fully resolves the volume fraction over intervals smaller than the particle diameter. Using this physically based parametrization, we obtain a much better fit of (5.4) to our data (dotted line in figure 15 b), which illustrates the importance of including the effects of the counterflow, the buoyancy and the increased viscosity in the formulation of the hindered settling function even during the consolidation phase. On the other hand, (5.3) underestimates the settling of cohesive sediment, but is very close to the settling behaviour of the simulated cohesionless sediment. This was expected, as (5.3) was derived for cohesionless sediments in the first place.

5.3 Energetics of enhanced settling

We now analyse the energetics of the sedimentation process, with a focus on the conversion of the initial potential energy of the particles into kinetic energy, and on the modulation of this process by hydrodynamic and collision forces. Corresponding studies of the energetics of gravity and turbidity currents have provided useful information for the fluid phase that can facilitate the development of simplified models (e.g. Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005; Konopliv & Meiburg Reference Konopliv and Meiburg2016). By integrating the particle equation of motion (2.3) along its trajectory, we obtain for the energy budget of a single particle

(5.5) $$\begin{eqnarray}\displaystyle & & \displaystyle \underbrace{{\textstyle \frac{1}{2}}m_{p}(u_{p}^{0})^{2}}_{=\,E_{k,p}^{0}}+\underbrace{V_{p}(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})gy_{p}^{0}}_{=\,E_{y,p}^{0}}\nonumber\\ \displaystyle & & \displaystyle \quad =\underbrace{{\textstyle \frac{1}{2}}m_{p}\boldsymbol{u}_{p}^{2}}_{=\,E_{k,p}(t)}+\underbrace{V_{p}(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})gy_{p}}_{=\,E_{y,p}(t)}-\underbrace{\int _{0}^{t}\boldsymbol{F}_{h,p}\boldsymbol{u}_{p}\,\text{d}t^{\ast }}_{=\,W_{h,p}(t)}-\underbrace{\int _{0}^{t}\boldsymbol{F}_{c,p}\boldsymbol{u}_{p}\,\text{d}t^{\ast }}_{=\,W_{c,p}(t)}.\end{eqnarray}$$

The left-hand side represents the energies of the initial (conservative) budget $E_{p}^{0}=E_{k,p}^{0}+E_{y,p}^{0}$ , where $E_{k,p}^{0}$ and $E_{y,p}^{0}$ denote the initial kinetic and potential energies, respectively. The right-hand side represents the budget $E_{p}(t)=E_{k,p}(t)+E_{y,p}(t)-W_{h,p}(t)-W_{c,p}(t)$ as a function of time, where $W_{h,p}(t)$ and $W_{c,p}(t)$ indicate the work performed on the particle up to time $t$ by the non-conservative hydrodynamic and collision forces $\boldsymbol{F}_{h,p}$ and $\boldsymbol{F}_{c,p}$ , respectively. To compute this work, these non-conservative forces are integrated over the entire particle path. The hydrodynamic forces and collision forces modify the total conservative energy $E_{k,p}+E_{y,p}$ .

To assess the energy budget, we store particle information such as position, velocity, hydrodynamic and collision forces every 1000 time steps. We can then compute the potential and kinetic energy, along with work performed by the hydrodynamic and collision forces, based on the integration scheme:

(5.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle E_{k}^{n}=\mathop{\sum }_{N_{p}}\frac{1}{2}m_{p}(\boldsymbol{u}_{p}^{n})^{2}, & \displaystyle\end{eqnarray}$$

(5.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle E_{y}^{n}=\mathop{\sum }_{N_{p}}V_{p}(\unicode[STIX]{x1D70C}_{p}-\unicode[STIX]{x1D70C}_{f})gy_{p}^{n}, & \displaystyle\end{eqnarray}$$

(5.6c ) $$\begin{eqnarray}\displaystyle & \displaystyle W_{h}^{n}=W_{h}^{n-1}+\frac{\unicode[STIX]{x0394}t}{4}\mathop{\sum }_{N_{p}}[\boldsymbol{F}_{h,p}^{n}+\boldsymbol{F}_{h,p}^{n-1}][\boldsymbol{u}_{p}^{n}+\boldsymbol{u}_{p}^{n-1}], & \displaystyle\end{eqnarray}$$

(5.6d ) $$\begin{eqnarray}\displaystyle & \displaystyle W_{c}^{n}=W_{c}^{n-1}+\frac{\unicode[STIX]{x0394}t}{4}\mathop{\sum }_{N_{p}}[\boldsymbol{F}_{c,p}^{n}+\boldsymbol{F}_{c,p}^{n-1}][\boldsymbol{u}_{p}^{n}+\boldsymbol{u}_{p}^{n-1}], & \displaystyle\end{eqnarray}$$

$n$

$p$

$n-1$

$n$

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\frac{|E^{0}-E(t)|}{E^{0}}\end{eqnarray}$$

over time, where $E^{0}=E_{y}^{0}$ and $E(t)=E_{k}(t)+E_{y}(t)-W_{h}(t)-W_{c}(t)$ . This error remained below 0.0045 for all times, which is sufficiently small to establish confidence in the results. We furthermore note that the error saturates at an almost constant level for larger times and is not expected to grow any further as the particles gradually come to rest at the bottom of the tank towards the end of the simulation.

Figure 16. Time evolution of the mechanical energy budget of all particles in the flow, normalized by the initial energy $E^{0}=E_{k}^{0}+E_{y}^{0}$ : (a) potential energy, (b) work performed by hydrodynamic forces, (c) kinetic energy, and (d) work due to collision forces.

The results of the energy analysis are shown in figure 16. Integrated over all particles, the contributions of the kinetic energy and the work performed by collision forces are orders of magnitude smaller than the contributions of potential energy and the work of the hydrodynamic forces, which indicates that the initially available potential energy $E_{p}$ is primarily used to overcome the viscous drag force as the particles settle (figure 16 a,b). During the initial stage $t<80\unicode[STIX]{x1D70F}_{s}$ , the curves for cohesive and non-cohesive sediment are nearly indistinguishable. Subsequently, however, as the cohesive sediment forms flocs and settles out more rapidly, its potential energy decays faster and it performs more work against the viscous drag forces.

During the initial stage, the kinetic energy data in figure 16(c) collapse for all simulations, with a distinct peak at $t=17.6\unicode[STIX]{x1D70F}_{s}$ and a subsequent exponential decay. This behaviour can be attributed to the fact that particles initially are distributed throughout the entire domain, so that particles close to the bottom wall immediately begin to feel the presence of the confinement. These particles will never accelerate towards their undisturbed settling velocity, and as soon as there are more particles decelerating than accelerating, $E_{k}$ starts to decay. Hence, the evolution of $E_{k}$ reflects the behaviour of a dissipative dynamical system released from rest, with an initial supply of potential energy, so that its dynamics resemble the temporal evolution of the fluid kinetic energy for a lock-release turbidity current propagating in a channel (Necker et al. Reference Necker, Härtel, Kleiser and Meiburg2005) or spreading radially in a basin (Francisco et al. Reference Francisco, Espath, Laizet and Silvestrini2018). The difference in settling velocities of the larger and smaller particles during the initial stage results in the strong stirring and mixing of both the fluid and the particles. During the interval $100\unicode[STIX]{x1D70F}_{s}<t<270\unicode[STIX]{x1D70F}_{s}$ , the kinetic energy is larger for the cohesive sediment, reflecting its higher settling velocity. As particles begin to deposit, this effect becomes less and less prominent since fewer particles remain in suspension and $E_{k}$ eventually approaches zero for all simulations. During the final simulation stages, the kinetic energy is slightly larger for the cohesionless sediment, since almost all of the cohesive sediment has already settled out.

This behaviour is also reflected in the work performed by the particles against the collision forces (figure 16 d). Even though the forces acting on particles $p$ and $q$ through equation (3.5) must be opposite and equal, the total work they perform against the collision forces is non-zero. The cohesive forces modify the amount of work performed by the particles during collisions, as they tend to align the particles and reduce their velocity difference. This observation is more pronounced for $\mathit{Co}=5$ . Since cohesive sediment settles out more quickly than cohesionless sediment, the late stages of the cohesive simulations see more collisions of flocculated particles with deposited particles, so that the work performed against the collision forces is larger for cohesive sediments.

Figure 17. Average velocity of the particle centre of mass $\langle v_{p}\rangle$ as a function of time: (a) all particles, (b) $D_{p}>D_{66}$ , (c) $D_{33}<D_{p}\leqslant D_{66}$ , and (d) $D_{p}\leqslant D_{33}$ .

The above demonstrates that cohesive forces modify the processes by which potential energy $E_{y}$ is converted into kinetic energy $E_{k}$ . As a result, the effective settling rate of the particle ensemble is altered, as reflected by the vertical velocity component of the centre of mass,

(5.8) $$\begin{eqnarray}\langle v_{p}\rangle =\frac{1}{\displaystyle \mathop{\sum }_{p=1}^{N_{p}}M_{p}}\displaystyle \mathop{\sum }_{p=1}^{N_{p}}M_{p}v_{p},\end{eqnarray}$$

cf. figure 17(a). Figure 17(b–d) display the ensemble-averaged velocity $\langle v_{p}\rangle$ , conditioned by particle size in the same way as in figure 14. These data confirm that the enhanced kinetic energy of the cohesive sediment during the time interval $100\unicode[STIX]{x1D70F}_{s}\leqslant t\leqslant 270\unicode[STIX]{x1D70F}_{s}$ can mainly be attributed to faster settling velocities, rather than enhanced horizontal velocity fluctuations. After peaking at $t=17.6\unicode[STIX]{x1D70F}_{s}$ , the settling process slows down most noticeably for cohesionless sediment. Consistent with our earlier observations, the settling of medium and small cohesive grains is accelerated most strongly by cohesive forces (figure 17 c,d). We note that for small cohesionless grains the ensemble-averaged settling velocity decays almost to zero at $t=110\unicode[STIX]{x1D70F}_{s}$ , as a significant fraction of them are swept upwards by the counterflow. Subsequently these smaller particles settle towards the bottom, reaching a constant settling velocity over time. By contrast, small cohesive grains attach to larger ones and consequently settle more rapidly. In addition, the settling velocity of the smaller particles follows the decelerating behaviour of the larger ones.

The above observations confirm the enhanced settling of cohesive sediment. We can quantify this effect by computing the relative increase in the settling velocity,

(5.9) $$\begin{eqnarray}\unicode[STIX]{x0394}\langle v_{p}\rangle =\frac{\langle v_{p}\rangle (\mathit{Co})-\langle v_{p}\rangle (\mathit{Co}=0)}{\langle v_{p}\rangle (\mathit{Co}=0)},\end{eqnarray}$$

cf. figure 18. After the acceleration phase up to $t=75\unicode[STIX]{x1D70F}_{s}$ , the cohesive particles with $\mathit{Co}=1$ and $\mathit{Co}=5$ settle up to 24 % and 29 % faster than cohesionless sediment. Beyond $t>250\unicode[STIX]{x1D70F}_{s}$ , more and more of the cohesive particles have reached the sediment bed, so that the relative settling velocity increase $\unicode[STIX]{x0394}\langle v_{p}\rangle$ loses its meaning.

Figure 18. (a) Relative settling velocity increase as a function of time. (b) Relative enhancement of the decay of potential energy for different values of $\unicode[STIX]{x1D6FC}=E_{y}/E^{0}$ .

An alternative way of quantifying the enhanced settling behaviour of cohesive sediment relies on the time $T_{\unicode[STIX]{x1D6FC}}$ it takes for the initial potential energy to decay to a relative value of $\unicode[STIX]{x1D6FC}=E_{y}/E^{0}$ . We define the relative settling time reduction as

(5.10) $$\begin{eqnarray}\unicode[STIX]{x0394}T_{s}=\frac{T_{\unicode[STIX]{x1D6FC}}(\mathit{Co}=0)-T_{\unicode[STIX]{x1D6FC}}(\mathit{Co})}{T_{\unicode[STIX]{x1D6FC}}(\mathit{Co}=0)},\end{eqnarray}$$

and show corresponding computational results in figure 18(b). The relative settling time reduction is larger for lower values of $\unicode[STIX]{x1D6FC}$ , as cohesive sediment continues to settle out faster than cohesionless particles for later times. The speedup is most pronounced as we compare cohesionless sediment with the cohesive case $\mathit{Co}=1$ . A further increase of the cohesive forces to $\mathit{Co}=5$ results in only slightly more rapid settling.

We conclude that the energy budget analysis represents a suitable tool for clarifying the mechanisms by which cohesive forces accelerate the settling of particles, as observed in § 5.2 as well as in experiments (Mehta et al. Reference Mehta, Hayter, Parker, Krone and Teeter1989).

5.4 Preferred particle interaction configurations

We proceed to explore if the observations of § 4 for interacting particle pairs can help explain the origins of enhanced settling for large particle ensembles. There we had found that, while cohesionless grains tend to undergo the DKT process, cohesive particle pairs bond to each other in certain preferred geometrical configurations that depend on the ratio $R_{p}/R_{q}$ of the particle radii. With a ratio of $R_{p}/R_{q}\neq 1$ , i.e. a polydisperse size distribution, particles were seen to align obliquely or even vertically.

In the following, we define particle $p$ as having a larger $y$ -coordinate than particle $q$ , so that $y_{p}>y_{q}$ . To test how much of the behaviour observed for isolated cohesive particle pairs can still be found in the context of many settling particles, we analyse the probability density function (PDF) of the angles of interacting particles. Owing to the symmetry of the problem, we can consider the angle of the contact point with respect to the vertical coordinate of particle $p$ as

(5.11) $$\begin{eqnarray}\text{cos}\,\unicode[STIX]{x1D703}=\frac{y_{p}-y_{q}}{\Vert \boldsymbol{r}_{pq}\Vert }~,\end{eqnarray}$$

and we can define vertical, oblique and horizontal contacts as having angles $0^{\circ }<\unicode[STIX]{x1D703}\leqslant 22.5^{\circ }$ , $22.5^{\circ }<\unicode[STIX]{x1D703}\leqslant 67.5^{\circ }$ , and $67.5^{\circ }<\unicode[STIX]{x1D703}\leqslant 90^{\circ }$ , respectively. Furthermore, we distinguish between direct contact ( $\unicode[STIX]{x1D701}_{n}<0$ ) and cohesive bonding ( $0\leqslant \unicode[STIX]{x1D701}_{n}\leqslant \unicode[STIX]{x1D706}$ ).

Figure 19. Probability density function of the radii ratio with respect to their contact angles: (a) horizontal direct contact, (b) horizontal cohesive bonding, (c) vertical direct contact, (d) vertical cohesive bonding, (e) oblique direct contact, and (f) oblique cohesive bonding. Particles bonding horizontally are preferentially of similar size, while particles in vertical or oblique configurations tend to have different sizes. Interacting cohesive grains most often have the smaller particle trailing the larger one, whereas for cohesionless grains the larger particle tends to catch up with the smaller one from above.

Table 2. Median values of particle radii ratios $R_{p}/R_{q}$ at different contact angles and direct contact and cohesive bonding.

Figure 19 shows PDFs of the particle radii ratio $R_{p}/R_{q}$ for horizontal, vertical and oblique particle interactions, respectively. These PDFs were calculated from all contact points throughout the simulation between particle pairs with a vertical velocity greater than $5\,\%$ of the characteristic velocity $u_{s}=\sqrt{g^{\prime }D_{50}}$ . This conditioning by velocity is necessary to avoid counting particles that have already settled out. The median values of the PDFs displayed in figure 19 are summarized in table 2. Figure 19(a,b) confirm the observation of § 4.1 that particles interacting horizontally preferentially are of similar size, as the median value of the PDF is very close to unity. This holds for direct contact as well as for cohesive bonding, and for all three simulations. For vertical particle interactions, figure 19(c,d) show that cohesive particles interact very differently from cohesionless grains. While cohesive particle pairs tend to arrange themselves with the smaller particle trailing the larger one, cohesionless grains behave oppositely. This reflects the fact that for cohesionless particle pairs of different size the larger particle tends to settle more rapidly, so that it catches up with the smaller one from behind. For cohesive particle pairs, on the other hand, DKT causes the smaller particle to arrange itself in the wake of the larger one, so that the PDF has a distinct maximum for particle radii ratios below unity. For oblique contact angles (figure 19 e,f), the direct contact PDFs have maxima slightly below (above) unity for cohesive (cohesionless) sediment. For oblique cohesive bonding, on the other hand, the PDFs have median values clearly below unity, although somewhat larger than for vertical contact angles. In summary, we find that many of the features observed for isolated particle pairs in § 4.1 remain relevant within a larger ensemble of sedimenting polydisperse particles.