3.2.1 Point-particle approaches

For particles smaller than the Eulerian grid cells, the particle surface is no longer geometrically resolved. Under these conditions one frequently assumes that it suffices to treat the particle as a mass point with a virtual spatial extent that is smaller than the grid cell size (Reference BalachandarBalachandar, 2009; Reference Balachandar and EatonBalachandar & Eaton, 2010). Since the particle surface is not resolved by the computational mesh, a simplified approach has to be employed to obtain the relevant contributions to the hydrodynamic force acting on the particles

(3.5)\begin{equation} \boldsymbol{F}_{h,p} = \boldsymbol{F}_{d,p} + \boldsymbol{F}_{l,p} + \boldsymbol{F}_{f,p} + \boldsymbol{F}_{a,p} + \boldsymbol{F}_{b,p} + \cdots, \end{equation}

viz. the drag $\boldsymbol {F}_{d,p}$ and lift forces $\boldsymbol {F}_{l,p}$, the pressure gradient force from the undisturbed flow $\boldsymbol {F}_{f,p}$, the added mass force $\boldsymbol {F}_{a,p}$ and the Basset history force $\boldsymbol {F}_{b,p}$, as derived by Reference Maxey and RileyMaxey and Riley (1983). Depending on the particular flow field under consideration, some of these contributions may be negligible in size, so that (3.5) can be further simplified, or additional terms may have to be accounted for, as indicated by the $(\cdots )$ notation. Since the flow around the particles is not spatially resolved, closure models are needed for the hydrodynamic force contributions. Semi-empirical expressions for the individual terms are typically derived from simplified situations, such as a sphere settling under its own weight in a quiescent, viscous fluid (Reference BassetBasset, 1890; Reference BoussinesqBoussinesq, 1903; Reference OseenOseen, 1927). For denser systems of volume fractions larger than $10^{-3}$, corrections exist that account for higher particle concentrations (e.g. Reference Bogner, Mohanty and RüdeBogner, Mohanty, & Rüde, 2015; Reference Tenneti, Garg and SubramaniamTenneti, Garg, & Subramaniam, 2011). In case of two-way or four-way coupling, different schemes have been proposed to compute the particle forcing in the Navier–Stokes equation (3.1) (e.g. Reference Ferrante and ElghobashiFerrante & Elghobashi, 2003). A robust method that converges under grid refinement has been put forward by Reference Capecelatro and DesjardinsCapecelatro and Desjardins (2013),

(3.6)\begin{equation} f_d(\boldsymbol{x}) = \sum_{p = 1}^M \boldsymbol{F}_{h,p} G(\boldsymbol{x}-\boldsymbol{x}_p), \end{equation}

where $G(\boldsymbol {x}-\boldsymbol {x}_p)$ is a kernel, often Gaussian, with units inverse volume and $\boldsymbol {x}_p$ is the centre of particle $p$.

3.2.2 Particle-resolved direct numerical simulations

As reviewed by Reference Biegert, Vowinckel and MeiburgBiegert, Vowinckel, and Meiburg (2017a), there are numerous frameworks to conduct pr-DNS. A particularly successful and popular one has been the immersed boundary method (IBM) with direct forcing (Reference Biegert, Vowinckel and MeiburgBiegert et al., 2017a; Reference Kempe and FröhlichKempe & Fröhlich, 2012a; Reference Lucci, Ferrante and ElghobashiLucci, Ferrante, & Elghobashi, 2010; Reference UhlmannUhlmann, 2005). For this method, the underlying Newton–Euler equations for spherical particles (3.3) and (3.4) are solved, where the hydrodynamic stresses on the particle surface are spatially resolved to yield the hydrodynamic stresses and torques as

(3.7)\begin{equation} \boldsymbol{F}_{h,p} = \oint_{\varGamma_p} \boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{n}\,{\text{d}\varGamma} \end{equation}

and

(3.8)\begin{equation} \boldsymbol{T}_{h,p} = \oint_{\varGamma_p} \boldsymbol{r}\times(\boldsymbol{\tau}\boldsymbol{\cdot}\boldsymbol{n})\,{\text{d}\varGamma}, \end{equation}

respectively. In contrast to point-particle methods, a key feature of the IBM is that the hydrodynamic forces and torques $\boldsymbol {F}_{h,p}$ and $\boldsymbol {T}_{h,p}$, respectively, become a direct result of the fluid–particle coupling. The fluid acts on the particles through the hydrodynamic stress tensor $\boldsymbol {\tau }$, where $\boldsymbol {r}$ represents the vector from the particle centre to a point on the surface $\varGamma ^p$ and $\boldsymbol {n}$ is the unit normal vector pointing outwards from that point. The net force and torque acting on the particle centre of mass due to collisions and contacts are again given by $\boldsymbol {F}_{c,p}$ and $\boldsymbol {T}_{c,p}$, respectively.

Vice versa, the fluid at an arbitrary location $\boldsymbol {x}$ experiences drag due to particle motion via $\boldsymbol {f}_d = \boldsymbol {f}_{ibm}$,

(3.9)\begin{equation} \boldsymbol{f}_{ibm}(\boldsymbol{x}) = \rho_f\sum_{p}^{N_{tot}} \sum_{l}^{N_{l,p}}\delta_h(\boldsymbol{X}_l-\boldsymbol{x}) \frac{\boldsymbol{U}_{l,p}^d-\boldsymbol{U}_{l,p}}{\Delta t}V_{l,p}, \end{equation}

where $N_{tot}$ is the number of particles in the domain and $N_{l,p}$ is the number of Lagrangian markers employed to discretize the surface of particle $p$. These Lagrangian markers are placed on the particle in equidistant spacing of the order of the grid cell size to discretize its surface. To couple the mesh of these Lagrangian markers to the fluid grid, a regularized Dirac delta function $\delta _h$ is introduced that spreads the forcing from the Lagrangian marker point on the particle surface to the Eulerian grid (Reference Roma, Peskin and BergerRoma, Peskin, & Berger, 1999; Reference Yao, Biegert, Vowinckel, Köllner, Meiburg, Balachandar and FringerYao et al., 2022). Furthermore, $\boldsymbol {U}^d_{l,p}$ is the desired velocity at the particle surface, $\boldsymbol {U}_{l,p}$ is the interpolated fluid velocity computed prior to the forcing correction and $V_{l,p}$ is the volume element of the thin shell around particle $p$ that is associated with Lagrangian marker $l$, and that is located inside a given fluid grid cell. Hence, $\boldsymbol {f}_{ibm}$ acts as a correction to the fluid velocity to enforce the no-slip condition on the particle surface.

3.4.1 The JKR-type models

As mentioned earlier, the JKR theory of Reference Johnson, Kendall and RobertsJohnson et al. (1971) is appropriate for soft materials with moderate to high surface energy. In fact, the JKR model was first developed for rubber materials (Reference Johnson and GreenwoodJohnson & Greenwood, 1997). Nevertheless, it has been found useful for DEM simulations to investigate, e.g. the mixing, and (de-)agglomeration of fine powders as reviewed by Reference Chen and ElliottChen and Elliott (2020). One model compatible with DEM was put forward by Reference Chokshi, Tielens and HollenbachChokshi, Tielens, and Hollenbach (1993) and implemented in the framework of CFD-DEM simulations by Reference MarshallMarshall (2009). Subsequently, this approach was successfully applied to simulate the break-up of particle agglomerates in linear shear flow (Reference Dizaji and MarshallDizaji & Marshall, 2017) and isotropic turbulence (Reference Chen and ElliottChen & Elliott, 2020).

According to Reference MarshallMarshall (2009), the normal contact forces between particles $p$ and $q$ are expressed as

(3.20)\begin{equation} F_{n,pq} = F_{n,pq}^{E} + F_{n,pq}^{D}, \end{equation}

where $F_{n,pq}^{D}$ is a damping force to set the restitution coefficient to zero (Reference Dizaji and MarshallDizaji & Marshall, 2017) and $F_{n,pq}^{E}$ represents the forces from van der Waals attraction and elastic deformation. These forces are nonlinearly coupled to the radius of the contact region as a function of time, $a(t)$, via

(3.21)\begin{equation} F_{n,pq}^{E} = {-}4F_c\biggl[ \left(\frac{a}{a_0}\right)^3-\left(\frac{a}{a_0}\right)^{3/2}\biggr], \end{equation}

where $F_c$ it the critical normal pull-off force and $a_0$ is the equilibrium contact region radius. The area of the contact region is computed by solving

(3.22)\begin{equation} \delta_n = 6 \delta_c\biggl[ 2\left(\frac{a}{a_0}\right)^2-\frac{4}{3}\left(\frac{a}{a_0}\right)^{1/2}\biggr], \end{equation}

where $\delta _n=-\zeta _n$ is the particle overlap of the two interacting spheres. The pull-off force is the force needed to pull particles apart if they are bonded by cohesive forces only, i.e. $\delta _n = \delta _c$. Owing to the nonlinear coupling of contact forces and adhesive forces through $a(t)$, this yields

(3.23)\begin{equation} F_c = 3{\rm \pi} \gamma R_r \end{equation}

and $\delta _c$ is the critical distance for which necking of the two particles remains possible in the contact region. Note that (3.21) replaces the terms of (3.10) due to the nonlinear coupling of cohesive force and surface deformation during the collision of the soft spheres. Similarly, the torque induced by cohesion upon contact is rewritten as

(3.24)\begin{equation} T_{c,p} = {-}k_R \zeta_t, \end{equation}

where $k_R=4F_c(a/a_0)^{3/2}$ is the tangential stiffness for cohesive materials and $\zeta _t$ is the same tangential displacement that enters (3.18). With this type of contact modelling, particles remain in contact until the tensile normal forces exceed the critical pull-off force $F_c$ and the particles move apart beyond $\delta _c$. Even though the JKR model has been derived for soft materials with $\lambda _T\gg 1$, Reference Johnson and GreenwoodJohnson and Greenwood (1997) claim that similarly good results can be obtained for materials with a smaller Tabor parameter.

A somewhat similar approach to the JKR model is the so-called hit-and-stick assumption. Under this assumption, particles will always stick together and form larger aggregates as they come into contact. It was pointed out by Reference Chen, Li and MarshallChen et al. (2019) that such an assumption is only valid if there are exclusively low energy collisions in the system. The hit-and-stick assumption has been a popular tool for dust aggregation (Reference Okuzumi, Tanaka and SakagamiOkuzumi, Tanaka, & Sakagami, 2009) and even for planet formation in protoplanetary disks (Reference Güttler, Blum, Zsom, Ormel and DullemondGüttler, Blum, Zsom, Ormel, & Dullemond, 2010), but these types of studies are usually done for DEM simulations in dry conditions only.

To model the potential energy suddenly released to the particles upon detachment in pr-DNS, Reference DerksenDerksen (2014) proposed a square-well potential for agglomerated particles. This potential applies to two particles bonded by cohesive forces until the tensile force pulls them apart and the potential energy is converted into kinetic energy. Reference Delenne, El Youssoufi, Cherblanc and BénetDelenne, El Youssoufi, Cherblanc, and Bénet (2004) proposed a contact force model that introduces cohesive bonds that establish a rigid bridging of the particles. These bonds fail as soon as the yield load is reached and the usual contact and friction laws are used.

3.4.2 The DMT-type models

Even though the JKR model has been applied with success to the deagglomeration of cohesive particles, many studies deal with materials that are quite stiff, such as silica grains and clay flocculi described in § 2.2. In addition, the nonlinear coupling of forces due to collision and cohesion poses challenges to the numerical stability of classical CFD-DEM solvers. For these reasons, many studies have adopted the additive DMT model that is more in line with the underlying idea of the Newton–Euler equations (3.3) and (3.4). For example, Reference van Wachem, Thalberg, Remmelgas and Niklasson-Björnvan Wachem, Thalberg, Remmelgas, and Niklasson-Björn (2017) model collisions by a spring-dashpot system and introduce an additional cohesive contact force $F_{vdW}$ that acts normal to the surface and is proportional to the particle overlap $\delta$ at a finite range $\delta _0$ representing the equilibrium distance for cohesive contact. For larger overlaps, the cohesive forces are equivalent to the critical pull-off force $F_c/2$.

For the DMT model, $F_c$ is computed via the Hertzian contact theory, which yields $F_c=4{\rm \pi} R_r \gamma$, where again $\gamma =A_h/(24{\rm \pi} \delta _0^2)$ is the surface energy due to van der Waals forces. Compared with the JKR model (3.23), this slightly larger value is the direct consequence of the additive and nonlinear coupling of the DMT- and JKR-framework approaches, respectively. A detailed derivation of $F_c$ for both frameworks can be found in Reference Marshall and LiMarshall and Li (2014). Owing to the additive nature of the DMT model, this cohesive force $F_{vdW}$ can subsequently be added to the resulting force from particle interaction and contact (3.10).

Similarly, Reference Yao and CapecelatroYao and Capecelatro (2021) followed Reference Gu, Ozel and SundaresanGu, Ozel, and Sundaresan (2016) and defined cohesive forces as $F_{vdW}\propto \zeta _2^{-2}$ that act in a finite gap in between particles and becomes constant for particle overlap. These forces are added to the classic spring-dashpot model for contact, so that they can directly enter (3.10). Full details are provided in Reference Gu, Ozel and SundaresanGu et al. (2016) and Reference Yao and CapecelatroYao and Capecelatro (2021). Reference Yao and CapecelatroYao and Capecelatro (2021) investigated the deagglomeration of cohesive aggregates by isotropic turbulence and determined $0.19\leq \lambda _T\leq 0.98$. They found that their results are not sensitive to the choice of using either a DMT or a JKR model for the cohesive forces.

A recent study by Reference Yu, Yu and BalachandarYu et al. (2022) followed the example of Reference Yao and CapecelatroYao and Capecelatro (2021) but, in addition, these authors claim that adhesive forces according to the JKR model need to be included. They follow the reasoning by Reference Parteli, Schmidt, Blümel, Wirth, Peukert and PöschelParteli et al. (2014) and propose to add an additional force during contact that is proportional to the radius of the contact region $a$ (Reference KendallKendall, 1971). Therefore, this approach represents a mixture of DMT and JKR modelling and it contradicts the rationale put forward by Reference Marshall and LiMarshall and Li (2014) and Reference Yao and CapecelatroYao and Capecelatro (2021), who argued that the JKR model implies that collision and cohesion do not represent additive effects. Nevertheless, Reference Yu, Yu and BalachandarYu et al. (2022) have successfully used this model to generate data for particles flocculating in isotropic turbulence to quantify rates of aggregation and disaggregation, respectively.

3.4.3 Cohesive force model compatible with point-particle and particle-resolved simulations

For geometrically resolved CFD-DEM applications, i.e. pr-DNS, the modelling approach poses some additional difficulties that go beyond the DEM approaches summarized in §§ 3.4.1 and 3.4.2. Owing to the high resolution of pr-DNS, the computational model needs to properly reflect the following three phases of particle–particle interaction: (a) particle approach/flocculation, (b) capture/steady state contact and (c) separation/deagglomeration in the presence of external forces. Such a model provides detailed information on the work performed by the inter-particle forces upon approach and separation. To this end, Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) propose a model that follows the DMT theory, where collision and cohesion are treated as additive effects. Hence, the computational model recovers the original DEM scheme by Reference Biegert, Vowinckel and MeiburgBiegert et al. (2017a) for cohesionless grains in the framework of IBM simulations of particle-laden flows.

To guarantee the compatibility with the soft-sphere DEM of Reference Biegert, Vowinckel and MeiburgBiegert et al. (2017a), Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) developed a model that does not alter the original algorithm for contact and collision and is consistent with the DLVO theory as sketched in figure 2. In particular, Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) aim to provide an algebraic expression that mimics the secondary maximum of the DLVO curve, where flocculation is possible. According to the DLVO theory, the secondary maximum yields an attractive inter-particle force within the interval $2\,\text {nm} \leq \zeta _n \leq 10\,\text {nm}$ that has a local maximum at $\zeta _n \approx 4\,\text {nm}$. Length scales of $\zeta _n < 2\,\text {nm}$ are not considered as they are part of the surface roughness. To this end, the model is designed such that cohesive forces decay to zero for $\zeta _n = 0$, while the repulsive forces act during contact through (3.17). Such a model is well represented by a parabolic spring force

(3.25) \begin{equation} \boldsymbol{F}_{coh} = \begin{cases} -k_{coh}(\zeta_n^2 - \zeta_n \lambda) \boldsymbol{n} & 0, < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}, \end{cases} \end{equation}

where $k_{coh}$ denotes the stiffness constant. As desired, (3.25) has 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. As an important parameter, $\lambda$ represents the range over which the cohesive force is smeared. This measure becomes necessary because length scales of a few nanometres are typically not resolved in CFD.

As was shown by Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b), the simulation results are insensitive to the exact value of $\lambda$. For pr-DNS, $\lambda$ has to be large enough so that cohesive forces can be numerically resolved, but much smaller than the particle size, e.g. $D_{50} / \lambda = 20$, where $D_{50}$ is the median grain size of a polydisperse grain size distribution. Owing to the substepping routine proposed by Reference Biegert, Vowinckel and MeiburgBiegert et al. (2017a), $\lambda$ can be as small as the grid cell size $h$, where a typical resolution of pr-DNS is 20 grid cells per diameter. Such a choice still guarantees a proper resolution of cohesive effects in space and time. This modelling approach is also consistent with the experimental observations of Reference Delenne, El Youssoufi, Cherblanc and BénetDelenne et al. (2004), who used the same data to derive their cohesion model outlined in § 3.4.1. In these experiments, rods of $D = 8$ mm in diameter were coated 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 $\zeta _n$ to a maximum at $D / \zeta _n \approx 80$. For larger gap sizes, the force decreases and eventually, the rods detach.

Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) determine the stiffness $k_{coh}$ of the model by preserving the energy contained in the van der Waals forces (Reference IsraelachviliIsraelachvili, 1992)

(3.26)\begin{equation} F_{vdW} = \frac{A_h R_{eff}}{6 \zeta_0^2}, \end{equation}

which yields

(3.27)\begin{equation} k_{coh} = \frac{A_h R_{eff}}{\zeta_0 \lambda^3 } \end{equation}

to obtain the dimensional cohesive force model

(3.28) \begin{equation} \boldsymbol{F}_{coh} = \begin{cases} -\dfrac{A_h R_{eff}}{\zeta_0 \lambda^3 }(\zeta_n^2 - \zeta_n \lambda) \boldsymbol{n}, & 0 < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}. \end{cases} \end{equation}

This form provides an expression for cohesive forces in pr-DNS via an additional force term in the collision model (3.10)

(3.29)\begin{equation} \boldsymbol{F}_{c,p} = \sum_{q,\ q \neq p}^{N_p} \left( \boldsymbol{F}_{l,pq} + \boldsymbol{F}_{n,pq} + \boldsymbol{F}_{t,pq} + \boldsymbol{F}_{{coh},pq}\right) + \boldsymbol{F}_{l,pw} + \boldsymbol{F}_{n,pw} + \boldsymbol{F}_{t,pw} + \boldsymbol{F}_{{coh},pw}, \end{equation}

where $pq$ and $pw$ indicate collisions of particle $p$ with another particle of index $q$ or the wall, respectively. Such a framework can be used in conjunction with the DEM of Reference Biegert, Vowinckel and MeiburgBiegert et al. (2017a) without any further modification and the hydrodynamic forces and the particle weight, $\boldsymbol {F}_{h,p}$ and $\boldsymbol {F}_{g,p}$, respectively, remain a result of the original IBM. Such a model resolves the particle bonding process in space and time and, at the same time, it uses the DMT rationale to retain the distinction between the individual inter-particle force components via (3.29).

For pr-DNS, it is convenient to write the dimensional form (3.28) in terms of non-dimensional quantities. This especially applies to the proper parameterization of the Hamaker constant $A_h$ with respect to particle inertia. To this end, suitable non-dimensional numbers are needed to define cohesive effects for arbitrary systems. As usual, one can render the Navier–Stokes equation (3.1) dimensionless by choosing a proper reference length scale, velocity and density ($L_0$, $U_0$ and $\rho _0$, respectively), which yields the Reynolds number $Re=U_0 L_0 / \nu _f$ as a dimensionless parameter (Reference Biegert, Vowinckel, Ouillon and MeiburgBiegert et al., 2017b) where again $\nu _f$ represents the kinematic viscosity. Applying the same logic to the particle equation of motion (3.3), non-dimensional values emerge for the lubrication force and the cohesive force scales. By writing the algebraic expression for cohesive forces (3.28) in dimensionless form, Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) obtained the cohesive number

(3.30)\begin{equation} Co = \frac{\text{max}(\| \boldsymbol{F}_{{coh},50} \|)}{U_0^2L_0^2\rho_0} \end{equation}

as the ratio of the cohesive force maximum to the characteristic inertial force scale of the problem under consideration. Similar characteristic numbers have been introduced for CFD-DEM applications in the form of an adhesion parameter $Ad=\gamma / (\langle u_p'\rangle ^2 R_p^2 \rho _p)$ as the ratio of surface energy due to the mean kinetic energy based on the granular temperature of the particles (e.g. Reference Dizaji and MarshallDizaji & Marshall, 2017; Reference Li and MarshallLi & Marshall, 2007; Reference Yao and CapecelatroYao & Capecelatro, 2021) or the bond number $Bo$ as the ratio of cohesive forces to the particle weight (e.g. Reference Sun, Xiao and SunSun, Xiao, & Sun, 2018). In this context, $Bo$ should not be confused with the Bond number, which represents the ratio of gravitational forces to surface tension forces and is unrelated to cohesive bonding. This dimensionless number is named after Wilfrid Noel Bond and is commonly used to characterize the shape of bubbles in a viscous flow (Reference Clift, Grace and WeberClift, Grace, & Weber, 2005). For the applications presented in the present review, the choices of the relevant velocity and length scales, $U_0$ and $L_0$ that enter (3.30), very much depend on the problem under considerations as will be detailed in § 3.5 below

To determine the numerator of (3.30), Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) compute the maximum cohesive force at $\zeta _n = \frac {1}{2} \lambda$ to obtain, for $R_{eff} = D_{50} / 2$,

(3.31)\begin{equation} \text{max}(\| \boldsymbol{F}_{{coh},50} \|) = {-}\frac{A_h}{\zeta_0} \frac{D_{50}}{2 \lambda^3} \left(\frac{\lambda^2}{4} - \frac{\lambda^2}{2}\right) = \frac{A_h}{\zeta_0} \frac{D_{50}}{8 \lambda}, \end{equation}

which yields the dimensionless cohesive force as

(3.32)\begin{equation} \tilde{\boldsymbol{F}}_{coh} = \begin{cases} - Co \dfrac{8 R_{eff}}{\lambda}^2(\tilde{\zeta}_n^2 - \zeta_n \lambda) \boldsymbol{n}, & 0 < \zeta_n \leq \lambda, \\ 0 & \text{otherwise}. \end{cases} \end{equation}

The stiffness of the cohesive force model is therefore determined as $k_{coh} = 8 {Co} R_{eff} / \lambda ^2$, where the tilde indicates dimensional values for the effective radius and the cohesive force range. As desired, the cohesive forces for a given physical system scale linearly with the cohesive number and the effective radius of the two interacting particles, a scaling that is consistent with the considerations of Reference VisserVisser (1989), Reference Lick, Jin and GailaniLick, Jin, and Gailani (2004) and Reference Righetti and LucarelliRighetti and Lucarelli (2007). According to § 2.2, such a model is meant to represent rough macroscopic (non-Brownian) particles. This model is based on principles of the DLVO theory, and at the same time it is compliant with the DMT theory. Given the advantages outlined above, such a model is especially useful for flocculation processes in pr-DNS (Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al., 2019a, Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019b; Reference Zhu, He, Zhao, Vowinckel and MeiburgZhu et al., 2022), but it can also be applied to point-particle approaches undergoing continuous agglomeration and break-up in a quasisteady state (Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al., 2021, Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020).

4.1.1 Flocculation in cellular model flows

An instructive setting for investigating the three-way coupled dynamics of flocculation and break-up is given by the Taylor-Green cellular flow shown in figure 5(a), with fluid velocity field ${\boldsymbol u}_f = (u, v)^{\rm T}$,

(4.1a,b)\begin{equation} u = \frac {U_0}{\rm \pi} \sin\left(\frac{{\rm \pi} x}{L}\right) \cos\left( \frac{{\rm \pi} y}{L}\right),\quad v_f = {-}\frac {U_0}{\rm \pi} \cos\left( \frac{{\rm \pi} x}{L}\right) \sin\left( \frac{{\rm \pi} y}{L}\right). \end{equation}

Here $L$ and $U_0$ represent the characteristic length and velocity scales of the vortex flow, which are used to render the problem dimensionless. This steady, doubly periodic flow field can be viewed as a simple analytical model of turbulence at the Kolmogorov scale, and it has been successfully employed for investigating several fundamental aspects of non-cohesive particle–vortex interactions (Reference Bergougnoux, Bouchet, Lopez and GuazzelliBergougnoux, Bouchet, Lopez, & Guazzelli, 2014; Reference MaxeyMaxey, 1987).

Figure 5. (a) Streamlines of the doubly periodic background flow. (b) Typical floc configuration made up of spherical primary particles, with individual flocs distinguished by colour.

A typical floc configuration for $N_p = 50$ particles in a flow domain of size $L_x \times\, L_y = 2 \,\times\, 2$ is shown in figure 5(b). All primary particles have identical diameter $D_p = 0.1$ and density $\rho _s = 1$, with $St = 0.1$ and $Co = 5 \,\times\, 10^{-4}$. Details regarding the non-dimensionalization are provided in Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al. (2020). Initially the particles are at rest and randomly distributed throughout the domain. When the distance between two particles is less than $\lambda /2$, we consider them as part of the same floc. We then track the number of flocs $N_f$ as a function of time, with an individual particle representing the smallest possible floc. Figure 6(a) demonstrates that the number of flocs decreases rapidly from its initial value $N_p$ due to flocculation, before levelling off around a constant value $N_{f,min}$ that reflects a stable balance between aggregation and breakage. We can fit the transient variation of $N_f(t)$ by an exponential function of the form

(4.2)\begin{equation} N_f(t) = (N_{p} - N_{f,min}) \,{\rm e}^{b t} + N_{f,min}, \end{equation}

where we evaluate $N_{f,min}$ as the average number of flocs during the equilibrium stage $50 \leqslant t \leqslant 200$. A characteristic flocculation time scale $t_{char}$ can then be defined as the time it takes for the number of flocs to decrease from its initial value $N_p$ to $N_{f,char} = N_p /2$. Hence, the corresponding characteristic time can be calculated as $t_{char} = {\rm ln} [(N_p /2 - N_{f,min}) / (N_p - N_{f,min})] / b$.

Figure 6. (a) Typical evolution of the number of flocs $N_f$ as a function of time. (b) Probability density function (PDF) of the floc size distribution during the equilibrium stage $50 \leqslant t \leqslant 200$. Simulation parameters are $N_p = 50$, $D_p = 0.1$, $\rho _s = 1$, $St = 0.1$ and $Co = 5 \times 10^{-4}$.

Figure 6(b) shows the statistical floc size distribution during the equilibrium stage $50 \leqslant t \leqslant 200$, where the floc size $N_{p,local}$ denotes the number of particles in a floc. Here $N_{f,local}$ refers to the number of the flocs of the same size. We find that the floc size distribution closely follows a log-normal distribution, consistent with previous experimental observations (Reference Bouyer, Line and Do-QuangBouyer et al., 2004; Reference Hill, Boss, Newgard, Law and MilliganHill, Boss, Newgard, Law, & Milligan, 2011; Reference Verney, Lafite, Burn-Cottan and le HirVerney, Lafite, Burn-Cottan, & le Hir, 2011b).

Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al. (2020) present a detailed discussion of the flocculation dynamics as a function of the governing dimensionless parameters. In particular, they observe that the number of flocs during the equilibrium stage $N_{f,min}$ decreases for increasing $Co$, along with the characteristic flocculation time. They find that the flocculation time $t_{char}$ has a minimum around $St \sim O(1)$, which reflects the well-known optimal coupling between particle and fluid motion when the particle response time is of the same order as the characteristic time scale of the flow (Reference Wang and MaxeyWang & Maxey, 1993). This value is larger than the critical value $St_c=1/(8{\rm \pi} )=0.04$ derived by Reference MassotMassot (2007) that marks the transition for which particle inertia allows cohesionless grains to cross between different cellular flow fields (Reference Yao and CapecelatroYao & Capecelatro, 2018). Under these conditions, particles rapidly accumulate near the edges of the vortices, which facilitates the formation of flocs. The authors furthermore discuss the average number of primary particles per floc, the mean floc size and the associated fractal dimension of the floc. These quantities are useful for parameterizing the terms in population balance models, cf. the pioneering work by Reference LevichLevich (1962) and current extensions (Reference Maggi, Mietta and WinterwerpMaggi et al., 2007; Reference Shin, Son and LeeShin et al., 2015; Reference Verney, Lafite, Burn-Cottan and le HirVerney et al., 2011b) and simplified versions (Reference Lee, Toorman, Molz and WangLee, Toorman, Molz, & Wang, 2011; Reference Shen, Lee, Fettweis and ToormanShen, Lee, Fettweis, & Toorman, 2018; Reference Son and HsuSon & Hsu, 2008, Reference Son and Hsu2009; Reference WinterwerpWinterwerp, 1998).

4.1.2 Flocculation in homogeneous isotropic turbulence

As a next step, we consider three-way coupled simulations with 10 000 primary particles in homogeneous isotropic turbulence (Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al., 2021). As discussed by those authors, the problem is fully characterized by a turbulent Reynolds number $Re$, a characteristic parameter of the random turbulent forcing process $D_s$, the dimensionless particle diameter $D_p$, the density ratio $\rho _s$ and the cohesive number $Co$. Here $Re$ and $D_s$ can equivalently be expressed by the shear rate $G$ of the turbulence, and a suitably defined particle Stokes number $St$. The ratio of the Kolmogorov length scale to the primary particle size takes values up to 3.3 in the simulations of Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021). Based on the observations by Reference Bosse, Kleiser and MeiburgBosse, Kleiser, and Meiburg (2015) that particle loading can modify the turbulence statistics even for volume fractions as low as $10^{-5}$, we expect two-way coupling effects to have an impact on the flocculation dynamics even in moderately dilute flows. Furthermore, even for globally dilute flows the local volume fraction inside a floc will be $O(1)$, resulting in significant shielding effects, so that the one-way coupled assumption generally will not hold within a floc. Nevertheless, as we will discuss in the following, three-way coupled simulations such as those by Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) are able to reproduce several experimentally observed statistical features of the flocculation dynamics.

Figure 7(a) shows the evolution of the number of flocs $N_f(t)$ with time for a representative case with $Co = 1.2 \times 10^{-7}$, $St=0.06$, $G = 0.62$, $\rho _s = 2.65$ and $\eta /D_p = 2.25$. Similar to the above cellular flow simulations of Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al. (2020), $N_f$ decreases rapidly with time from its initial value, before levelling off around a constant value $N_{f,eq}$ that reflects a stable equilibrium between aggregation and breakage. Figure 7(b) shows separately the number of flocs with $N_{p} = 1, 2, 3$ and more than three primary particles. While the number of two- and three-particle flocs initially grows quickly, it soon reaches a peak and subsequently declines, as more flocs of larger sizes form.

Figure 7. (a) Temporal evolution of the number of flocs $N_f$. The vertical dashed line divides the simulation into the flocculation and equilibrium stages. (b) Number of flocs containing $N_p$ primary particles. The number of flocs with a single particle rapidly decreases from its initial value of $N_f = 10\,000$. The numbers of flocs with two or three particles intially grow and subsequently decay, as increasingly many flocs with three or more particles form. All results are for $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $G = 0.62$, $\rho _s = 2.65$ and $\eta /D_p = 2.25$.

To obtain insight into the dynamics of floc growth and breakage, it is useful to keep track of those flocs that maintain their identity over a suitably specified time interval $\Delta T$, those that add additional primary particles while keeping all of their original ones, and all those that have undergone a breakage event during the time interval. Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) present detailed results in this regard, which show that all three fractions reach statistically steady states.

The simulations furthermore enable us to track the floc size $D_f$ as the Feret diameter and its fractal dimension $n_f$ with time. Figure 8 shows the evolution of $D_f$ and $n_f$ over time for a representative floc. During the time interval $200 \leqslant t \leqslant 210$, hydrodynamic forces deform the floc so that it becomes more compact, which reduces $D_f$ and increases $n_f$. Later on, near $t \approx 240$, the floc is being stretched, with corresponding changes in $D_f$ and $n_f$.

Figure 8. Temporal evolution of the characteristic diameter $D_f$ and the fractal dimension $n_f$ of a typical floc that maintains its identity over the time interval considered. Three instants are marked by vertical dashed lines, and the corresponding floc shapes are shown. In response to the fluid forces acting on it, the floc first changes from a slightly elongated to a more compact shape, and subsequently to a more strongly elongated one. The floc with seven primary particles is taken from the case with governing parameters $Co = 1.2 \times 10^{-7}$, $St = 0.1$, $G = 0.91$.

Figure 9 shows representative results for different Stokes numbers $St$ of the primary particles. Figure 9(a) indicates that the equilibrium value of $\bar N_{p}$ increases for smaller $St$. The evolution of the average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles, shown in figure 9(b), demonstrates that smaller Stokes numbers result in more compact flocs. By varying the governing parameters, Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) find that, as a general trend, weaker turbulence, lower Stokes numbers and higher cohesive numbers result in larger and more compact flocs.

Figure 9. Temporal evolution of (a) the average number of primary particles per floc $\bar N_p$, and (b) the average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles, for different Stokes number values $St$, with $Co = 1.2 \times 10^{-7}$, $G = 0.91$ and $\eta /D_p = 1.85$. During the equilibrium stage, the number of primary particles per floc and the fractal dimension increase for smaller Stokes numbers.

In order to discuss the floc size distribution during the equilibrium stage, we sort the flocs into bins of width $\Delta (D_f/D_p) = 0.7$. Figure 10(a) shows that, for all values of the turbulent shear $G$, the size distribution peaks at the smallest flocs and then decreases exponentially with the floc size. The decay rate is largest for the strongest turbulence, confirming our earlier observation that strong turbulence breaks up large flocs and reduces the average floc size. This finding is consistent with the experimental observations by Reference Braithwaite, Bowers, Nimmo Smith and GrahamBraithwaite et al. (2012) in an energetic tidal channel. Corresponding results for different $St$ values display a similar trend.

Figure 10. Floc size distribution during the equilibrium stage, obtained by sorting all flocs into bins of constant width $\Delta (D_f/D_p) = 0.7$. (a) Results for different shear rates $G$, with $Co = 1.2 \times 10^{-7}$ and $St = 0.06$, during the time interval $1000 \leqslant t \leqslant 4000$. (b) Results for different cohesive numbers $Co$, with $St = 0.02$ and $G = 0.29$, for the time interval $15\,000 \leqslant t \leqslant 19\,000$.

We refer the reader to Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) for further detailed information regarding the preferential alignment of the flocs with the principal strain directions of the turbulent velocity field, as well as with the local vorticity vector. The authors observe that elongated flocs and the vorticity vector are strongly aligned with the direction corresponding to the largest eigenvalue of the Lagrangian stretching tensor. This alignment is consistent with corresponding previous findings by Reference Parsa, Guasto, Kishore, Ouellette, Gollub and VothParsa et al. (2011) and Reference Ni, Ouellette and VothNi, Ouellette, and Voth (2014).

Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) furthermore discuss how information obtained from the simulations can be employed towards the development of flocculation models. Based on their observation that the diameter and the average fractal dimension of the flocs are typically related by a power law expression, the authors proceed to fit the exponent of this relationship from the simulation results, which avoids the assumption of a constant fractal dimension in earlier investigations (Reference Kuprenas, Tran and StromKuprenas et al., 2018; Reference WinterwerpWinterwerp, 1998; Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al., 2020). Consequently, the authors are able to formulate a new flocculation model, with variable fractal dimension.

In order to relate the cellular flow results to those for homogeneous isotropic turbulence, we note that one can view the cellular fields described in § 4.1.1 (Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al., 2020) such that each cell represents an idealized eddy with the size of the Kolmogorov length scale $L_0$ and a maximum velocity $U_0$. For the case of homogeneous isotropic turbulence, the turbulent eddy size is set by the characteristic time scale $T_0=L_0/U_0$ that defines the magnitude of the force term $\mathbf {f}_t$ entering (3.1). Using this framework, Reference Zhao, Vowinckel, Hsu, Köllner, Bai and MeiburgZhao et al. (2020) and Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) consider primary particles with diameter $D_p = 5 \,\mathrm {\mu }{\rm m}$ and demonstrated good agreement with the experiments of Reference Kuprenas, Tran and StromKuprenas et al. (2018) and Reference Maggi, Mietta and WinterwerpMaggi et al. (2007), respectively. For the case of homogeneous isotropic turbulence, Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) selected $T_0 = 7.81 \times 10^{-5}\,{\rm s}$ for the random turbulent forcing process. By choosing $L_0$, $T_0$ and $\rho _f = 1000\,{\rm kg}\,{\rm m}^{-3}$ as the characteristic length, time and density scales, one can obtain the characteristic velocity scale $U_0 = L_0/T_0 = 8\,{\rm m}\,{\rm s}^{-1}$, which is similar to values employed in previous investigations (Reference Chen and LiChen & Li, 2020; Reference Chen, Li and MarshallChen et al., 2019). The definition of (3.36) yields much smaller values of the cohesive number. For example, Reference Zhao, Pomes, Vowinckel, Hsu, Bai and MeiburgZhao et al. (2021) reported values in the range $10^{-10}\leq Co \leq 10^{-7}$. A conversion between the two definitions of the cohesive number (3.35) and (3.36) requires a choice for the value of the gravitational acceleration. For example, the standard value of $g=9.81$ m s$^{-2}$ yields a conversion factor of $3.7\times 10^8$.

4.2.1 Hindered settling

In order to explore the influence of cohesive forces on the sedimentation process, Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019a, Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019b) simulated the settling of an ensemble of 1261 polydisperse particles for the experimental conditions reported by Reference te Slaa, van Maren, He and Winterwerpte Slaa et al. (2015). The geometry of the particles was fully resolved, so that the fluid–particle coupling does not require any closures (§ 3.2.2). The particle–particle interactions were modelled the same way as described in § 4.1. A polydisperse mixture with a homogeneous particle volume fraction of $15\,\%$ was placed in a tank of quiescent fluid. To characterize the physical problem under consideration, it is convenient to define a reference velocity based on the buoyancy of a representative particle, which yields $u_{b}=\sqrt {g'D_{50}}$, where $g'=(\rho _p-\rho _f)g/\rho _f$ and $D_{50}$ denotes the median grain size of the particle size distribution. The characteristic time scale based on the buoyancy velocity and the median diameter then becomes $\tau _s = D_{50} / u_{b}$. Using a density ratio of $\rho _p/\rho _f=2.6$, the particle Reynolds number becomes $Re=u_b D_{50}/\nu _f = 1.35$. The grain sizes obey a cumulative log-normal distribution, with a maximum size ratio of ${\max \{D\} / \min \{D\}=4}$. The computational domain is of size $L_x \times L_y \times L_z = 13.1 D_{50} \times 40.0 D_{50}\times 13.1 D_{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. Two different values of the Cohesion number ${Co}$ are considered: (i) cohesionless grains with ${Co}=\max (\| \boldsymbol {F}_{{coh},50} \|) / m_{50} g'=0$, and (ii) strongly cohesive sediment with ${Co} = 5$. For both simulations, the particles are released from rest in quiescent fluid, and subsequently settle under the influence of gravity. Collisions are inelastic with $e_{p} = 0.97 < 1$, and subject to friction according to (3.18). Full details regarding the simulation set-up can be found in Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019a).

The impact of cohesive forces on the settling behaviour is illustrated by figure 11. During the early stages the particle configurations are very similar for both 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 horizontally averaged concentration profiles shown as contours in figure 12, as a function of time.

Figure 11. Particle configurations during the settling process for $Co=0$ (a) and $Co=5$ (b). Left column: $t=17.6 \tau _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 \tau _s$. The particle colouring reflects the vertical particle velocity. The cohesive sediment is seen to settle more rapidly than its non-cohesive counterpart.

Figure 12. Contours of the horizontally averaged particle volume fraction $\phi _v$ over time: (a) cohesionless sediment and (b) cohesive sediment (image taken from Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al., 2019a).

In agreement with the experimental observations of Reference Been and SillsBeen and Sills (1981), the volume fraction contours for the two simulations are nearly identical until $t = 120\tau _s$, as cohesive particles have not yet had sufficient time to flocculate. Subsequently, two distinct fronts emerge for both simulations. The upper front marks the transition between the clear fluid and the suspended sediment, whereas the second front shows the transition between the suspended sediment and the sediment bed. At $t \approx 350\tau _s$, the two fronts merge into one for the simulation data of the cohesive sediment. This point in time is called the point of contraction and marks the transition between hindered settling and consolidation (Reference WinterwerpWinterwerp, 2002). The cohesionless sediment, on the other hand, has not yet reached the point of contraction by the end of the simulation.

As the particles settle, they displace fluid at the bottom of the tank and generate an upward counterflow that is sufficiently strong in the current simulations to sweep smaller particles upward. This represents one reason for the reduced settling velocity of the smaller grains, and for the separation of the different grain sizes in very tall water columns (Reference te Slaa, van Maren, He and Winterwerpte Slaa et al., 2015). Figure 13 illustrates this effect by showing the final volume fraction profiles for the smallest, intermediate and largest third of the particles. The figure demonstrates that small and intermediate cohesive particles settle much more rapidly than their non-cohesive counterparts, consistent with the observation of Reference Lick, Jin and GailaniLick et al. (2004), who found intermediate size particles to be most strongly affected by cohesive forces. These types of grains have their peak concentration in the interval $3 D_{50} < y < 10 D_{50}$. In contrast, large cohesive grains have a lower volume fraction near the bottom of the tank than large cohesionless grains (figure 13c). 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 (Reference Mehta, Hayter, Parker, Krone and TeeterMehta et al., 1989).

Figure 13. Particle volume fraction profile for different particle radii at $t=480 \tau _s$: (a) small particles with $D_p \leq D_{33}$, (b) medium sized particles in the range $D_{33} < D_p \leq D_{66}$, and (c) large particles with $D > D_{66}$. Note the different horizontal axis scalings for the individual frames. The results in (a,b) were smoothed by a moving average with a filter width of $1.5 D_{50}$ for clarity.

Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) compared their simulation results for the effect of cohesion on the settling velocity of polydisperse particles to established empirical relations for silt. They computed the settling behaviour of the particle phase by means of an averaging operator applied to the Eulerian fluid grid that evaluates instantaneous snapshots of the particle velocity distribution $(\phi _v, v_p)$ to obtain horizontally averaged values of $\overline {\langle v_p \rangle }$ and the settling velocity $v_\infty$ (Reference Clift, Grace and WeberClift et al., 2005) of a single particle in still fluid at the same particle Reynolds number, where the angular brackets denote horizontal averaging and $\phi _v$ 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 this quantity still depends on the particle diameter. Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) computed its value in an iterative fashion using the force balance of a particle and empirical correlations for the drag coefficients. The ratio $\overline {\langle v_p \rangle }/v_\infty$ 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 Reference Richardson and ZakiRichardson and Zaki (1954),

(4.3)\begin{equation} \frac{\overline{\langle v_p \rangle}}{v_{\infty} } = \left(1-\frac{\phi}{\phi_s}\right)^n, \end{equation}

who experimentally investigated the settling behaviour of non-cohesive grains. Here, $\phi _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 Reference Dankers and WinterwerpDankers and Winterwerp (2007), 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. Hence, these authors introduced another critical volume fraction $\phi _{max}$ with the property $\phi _s<\phi _{max}$, where the transition from $\phi _s$ to $\phi _{max}$ is governed by gel collapse and self-weight consolidation rather than settling. Due to its simplicity, (4.3) has been very popular in hydraulic engineering. However, as described by Reference te Slaa, van Maren, He and Winterwerpte Slaa et al. (2015), this function is known to underestimate the settling velocity for higher concentrations. Instead, these authors have used the hindered settling function of Reference WinterwerpWinterwerp (2002),

(4.4)\begin{equation} \frac{\overline{\langle v_p \rangle}}{v_\infty} = \frac{\left(1-\dfrac{\phi}{\phi_s}\right)^m (1-\phi)}{\left(1-\dfrac{\phi}{\phi_{max}}\right)^{-({5}/{2})\phi_{max}}}, \end{equation}

where the numerator represents the effects of the counterflow and the buoyancy, respectively, while the denominator reflects the increased viscosity of dense suspensions according to Reference Krieger and DoughertyKrieger and Dougherty (1959). Here, $m$ is an empirical exponent, which plays a similar role to the parameter $n$ in (4.3). It was shown by Reference te Slaa, van Maren, He and Winterwerpte Slaa et al. (2015) that (4.3) and (4.4) both yield good agreement with experimental results for settling coarse silt.

Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019a) compared their simulation results to the hindered settling functions as parameterized by Reference te Slaa, van Maren, He and Winterwerpte Slaa et al. (2015) in figure 14(a) to validate their simulations. The results agreed well with the two empirical hindered settling functions, and they demonstrated the enhanced settling velocity due to the cohesive forces for all concentration values. It was shown by Reference Vowinckel, Biegert, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019a) that choosing $\phi _s = \phi _{max}$ yields good results. This is in contrast to the observations of Reference Dankers and WinterwerpDankers and Winterwerp (2007), who found for mud that $\phi _s<\phi _{max}$, but since the simulations of Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) deal with coarse silt, consolidation does not play a large role and it is justified to choose $\phi _s=\phi _{max}$. 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. Reference Houssais, Ortiz, Durian and JerolmackHoussais, Ortiz, Durian, & Jerolmack, 2016). Based on these observations, Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) propose to not calibrate $\phi _s$ but to parameterize critical volumetric concentrations by the maximum value of their concentration data. Using this data, they can immediately parameterize $\phi _s = \phi _{max}=0.7$, which is in line with experimental and computational studies of polydisperse particle packings (Reference Desmond and WeeksDesmond & Weeks, 2014; Reference Sohn and MorelandSohn & Moreland, 1968). 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 employed grid fully resolves the volume fraction over intervals smaller than the particle diameter. Using this physically based parameterization, one obtains a much better fit of (4.4) to the data (dotted line in figure 14b), which illustrates the importance of including the effects of the counterflow, of buoyancy, and of the increased viscosity in the formulation of the hindered settling function, even during the consolidation phase. On the other hand, (4.3) underestimates the settling velocity of cohesive sediment, but is very close to the settling behaviour of the simulated cohesionless sediment. This was expected, as (4.3) was derived for cohesionless sediments in the first place.

Figure 14. Settling velocity normalized with the undisturbed settling velocity. Comparison to empirical relationships (4.3) of Reference Richardson and ZakiRichardson and Zaki (1954) (RZ) and (4.4) of Reference WinterwerpWinterwerp (2002) (W).

4.2.2 Submerged cohesive granular collapse

Submerged non-cohesive granular collapses can be classified into three different categories (free fall, inertial and viscous), depending on the particle and fluid properties (Reference Bougouin and LacazeBougouin & Lacaze, 2018; Reference Jing, Yang, Kwok and SobralJing, Yang, Kwok, & Sobral, 2019). The relevant regime is determined by two dimensionless numbers, the Stokes number $St$ and the density ratio $r$,

(4.5)$$\begin{gather} St = \frac{1}{18\sqrt{2}}\frac{\sqrt{\rho_p(\rho_p-\rho_f)gD_p^3}}{\nu_f\rho_f}, \end{gather}$$

(4.6)$$\begin{gather}r = \sqrt{\frac{\rho_p}{\rho_f}}. \end{gather}$$

For $St\gg 10$ and $r\gg 4$, the collapse is in the free-fall regime, for $St\gg 2.5r$ and $r\ll 4$, it is in the inertial regime, and otherwise, it is in the viscous regime. For cohesive collapses, an additional parameter $Co$ arises in the form of (3.35).

Using the method by Reference Vowinckel, Withers, Luzzatto-Fegiz and MeiburgVowinckel et al. (2019b) summarized in § 4.2.1 above, Reference Zhu, He, Zhao, Vowinckel and MeiburgZhu et al. (2022) discuss particle-resolved cohesive simulation results for both shallow and tall columns, with representative results shown in figure 15. Consistent with earlier observations for non-cohesive collapses (Reference Lee, Huang and YuLee, Huang, & Yu, 2018; Reference Meruane, Tamburrino and RocheMeruane, Tamburrino, & Roche, 2010), cohesive collapses proceed through an acceleration or collapse stage, a constant front velocity stage and a deceleration stage. More cohesive collapses travel over shorter distances than their less cohesive counterparts, and with a smaller front velocity. For $Co=50$ with shallow columns and $Co=250$ with tall columns, the cohesive force is sufficiently large so that the column no longer collapses. This finding is qualitatively similar to previous experimental and numerical observations for dry cohesive granular collapse (Reference Artoni, Santomaso, Gabrieli, Tono and ColaArtoni, Santomaso, Gabrieli, Tono, & Cola, 2013; Reference Langlois, Quiquerez and AllemandLanglois, Quiquerez, & Allemand, 2015; Reference Santomaso, Volpato and GabrieliSantomaso, Volpato, & Gabrieli, 2018).

Figure 15. Evolution of granular collapse for $a=1$ and $Co=0$ and 30. From left to right: $t=5$, 20 and 30. (a–c) Magnitude of the angular velocity $\|\boldsymbol \omega _p\|$ for $Co=0$; (d–f) magnitude of the translational velocity $\|\boldsymbol {u_p}\|$ for $Co=0$; (g–i) $\|\boldsymbol \omega _p\|$ for $Co=30$; (j–l) $\|\boldsymbol {u_p}\|$ for $Co=30$. The red lines indicate the location of the failure surface. The black arrows represent vectors of the average particle velocity.

The angular and translational particle velocity magnitudes are displayed in figure 15. Figure 15(d,j) shows that the velocity magnitude in the $x$ and $y$ directions remains very small in the lower left corner of the columns, as the particles remain approximately at rest there. For both cohesive and non-cohesive (Reference Sun, Zhang, Wang and LiuSun, Zhang, Wang, & Liu, 2020) columns, during the acceleration stage particles near the upper right corner slide down along an inclined failure surface (indicated by a red line in frames 15d,j). The failure surface is defined as the contour where $\|\boldsymbol {u_p}\| = 0.05\|\boldsymbol {u_p}\|_{max}$ (Reference Lacaze and KerswellLacaze & Kerswell, 2009), with $\|\boldsymbol {u_p}\|_{max}$ denoting the maximum translational velocity at the same time. By comparing figure 15(d,j), we note that cohesive forces corresponding to $Co=30$ elevate the location of the failure surface, resulting in the growth of the region of stationary particles in the lower left corner. The angular velocity of the particles remains quite small near the failure surface (cf. frames 15a,b,g,h), which indicates that the particles slide, rather than roll, past each other (Reference Xu, Dong and DingXu, Dong, & Ding, 2019).

Based on the above results, Reference Zhu, He, Zhao, Vowinckel and MeiburgZhu et al. (2022) are able to establish scaling laws for the quasisteady front velocity and runout length as functions of $a$ and $Co$. They proceed to analyse the deposit morphology as well as the energy budget of the granular collapse events, in terms of the potential and kinetic energy components of the fluid and particulate phases, respectively.

The particle-resolving simulations provide complete information on the evolution of the individual cohesive force bonds, so that these can be analysed in detail. For the initial configurations with $Co=10$ as well as $a=1$ and 8.6, figure 16(a,c) indicates all cohesive bonds between individual particles by straight line segments that connect the particle centres. Figure 16(b,d) shows those of the initial cohesive bonds that have survived until the end of the collapse. We observe a few clear differences between the two aspect ratios. As discussed earlier, for $a=1$, the particles in the lower left corner hardly move at all, so that many of the initial cohesive bonds between them survive the entire collapse process. Near the left wall, and in the entire upper (pink) layer, quite a few of the bonds also survive, whereas this is not true for the lower sections of the deposit near the front. For the particles in that section, the collapse process destroys most of their initial cohesive bonds. For $a=8.6$, on the other hand, cohesive bonds primarily survive in a very small section in the lower left corner, and along the entire top layer of the deposit, including at the very front. For the entire interior section of the final deposit, almost all initial cohesive bonds are destroyed during the course of the collapse.

Figure 16. Cohesive bonds at the initial time for $a=1$, $Co=10$ (a) and $a=8.6$, $Co=25$ (c), respectively. For the same two flows, (b,d) shows those cohesive bonds that have stayed intact during the entire collapse process until the final time $t=60$.