Flocculation Dynamics of Cohesive Sediment in Turbulent Flows Using CFD-DEM Approach

Abstract

Two-phase computational fluid dynamics - discrete element method (CFD-DEM) framework has gained attention in cohesive sediment transport due to its capability of resolving particle-particle interactions and capturing the time evolution of individual flocs and hence the flocculation dynamics of cohesive sediment in turbulent flows. For cohesive sediments of size smaller than the Kolmogorov length scale, the point-particle approach is commonly used, in which the flow around particles is not fully resolved, and the hydrodynamic force on particles is parameterized by the drag law. The accuracy of floc dynamics, aggregation, breakup, and reshaping therefore strongly depends on force parameterization of individual point-particles that make up the floc. In this chapter, we review recent advances in the state-of-art two-phase CFD-DEM model approach on cohesive sediment transport and make recommendation for future research.

1. Introduction

Cohesive sediment transport is important for many geoscience and engineering applications. Cohesive sediments can adsorb pollutants (e.g., heavy metals, microplastics, and nutrients) and are a concern for water quality [1, 2]. The transport of estuarine and nearshore cohesive sediment plays an important role in coastal processes and the functioning of healthy ecosystems [3]. Cohesive sediment particles can aggregate with organic particles, such as bacteria, phytoplankton, and algae to form large flocs, marine snow [4], or algal-sediment aggregates [5]. The transport of flocs is therefore one of the key processes in the biogeochemical cycling of the ocean. Modeling the dynamics and transport of suspended sediment is essential to calculate sediment budgets and to provide relevant knowledge for the investigation of biogeochemical cycles. During transport, cohesive sediments undergo a complex flocculation process of aggregation, breakup, and reshaping, which determines the floc size and shape spectrum that in turn affects their sedimentation rate. Our ability to confidently predict the suspended sediment transport critically depends on how accurately we are able to parameterize the shape and size distribution of the cohesive sediments and hence their density and settling velocity.

Population balance equation (PBE) approach has been widely used to model flocculation dynamics in suspended cohesive sediment, in which floc properties change with ambient flow conditions due to either aggregation or breakup processes [6, 7]. By assuming fractal entities of flocs, the aggregation and breakup processes can be modeled as power-law relations. With recent advancements in imaging techniques, Spencer et al. [8] proposed a five-level hierarchy of floc aggregation, in which the floc structures are scale dependent. For microflocs ($O\left(10\mathrm{\mu m}\right)$), electrochemical forces and bridging mechanism by biological activities are of equal importance. For macroflocs ($O(100\mathrm{\mu m}$–$1\mathrm{mm})$), the backbone mechanism due to protruding filamentous bacteria influences the size and shape of macroflocs. The filaments can provide structural integrity by increasing the tensile strength of the flocs. Through the anchoring mechanism, the abundance and morphotype of bacteria affect the shapes of megaflocs ($O\left(1\mathrm{mm}\right)$). Due to the scale dependence of flocs structures, fractal theory cannot represent or predict the behavior of flocs in natural environments. A model framework that can resolve interactions among particles to better characterize the floc structure needs to be developed for more accurate prediction of floc properties in engineering or geoscience applications.

The advantage of the PBE approach is that it can scale to large systems of practical interest. In addition, PBE approach can capture the complex behavior of floc dynamics and predict the bimodal or multimodal floc size distributions that have been observed in the field [9, 10], compared to the single-class model using averaged floc size. However, the accuracy of the PBE approach entirely relies upon the fidelity of the closure models of aggregation and breakup processes, which have generally been motivated by experimental measurements of floc size and settling velocity. The video camera systems used in previous studies also have difficult in tracking individual flocs and directly observe aggregation and breakup processes in 3D [11, 12]. Numerical simulations using CFD-DEM approach can fill this gap to develop accurately models floc-floc and flow-floc interactions and hence improve PBE model accuracy on the prediction of cohesive sediment transport.

The two-phase CFD-DEM coupling framework has gained popularity in studies of sediment transport [13, 14], where sediment particles are modeled as dispersed phase by the discrete element method (DEM) [15, 16]. Turbulence-sediment interactions are modeled by momentum transfer between the continuous fluid phase and the dispersed sediment phase. Two approaches, namely the particle-resolved (PR) and the point-particle (PP) approaches, have been widely used. The PR simulations resolve the flow around individual particles and have been used to investigate important mechanisms in flocculation dynamics, such as aggregation and differential settling [17, 18]. The PR approach is computationally expensive and limited to simulations with a small number of particles (a few thousand), which may not be sufficient to capture mesoscale processes and obtain converged statistics. On the other hand, the PP approach can efficiently implement millions of particles [19]. However, important hydrodynamic interactions at the microscale, such as flow blockage and sheltering by neighboring particles, are not resolved [20].

In this chapter, we will review the recent advancements in the two-phase CFD-DEM approach to cohesive sediment transport studies. The review is structured as follows: Section 2 outlines the fundamental methodology of the CFD-DEM approach, followed by a discussion of its applications to flow modeling and future research (Section 3). We provide a summary of this approach in Section 4.

2. Two-phase CFD-DEM coupling approach for cohesive sediment transport in turbulent flows

2.1 Two-phase modeling approach

Two-phase modeling approach has been widely used in sediment transport. Depending on the volumetric concentration of sediment ($\varphi$) and Stokes number ($\mathit{St}$) defined as the ratio of the particle response time to the characteristic timescale of the turbulent flow, different computational approaches has been developed to turbulent multiphase flow [21], including the dusty gas, equilibrium Eulerian, Eulerian, Lagrangian point-particle, and fully resolve approaches. For dilute flow, one-way coupled or two-way coupled approaches are often adopted. In one-way coupled approach, only the effects of the continuous fluid phase on particles are considered, while the reverse effect is neglected. On the other hand, in two-way coupled approach, the fluid-particle interaction term is also introduced in the fluid momentum equation. When interactions among particles (collisions) become important, four-way coupled approach should be used. For sufficiently dilute flow, the fluid does not feel the presence of particles and particles move in response to the fluid; one-way coupling can be used. For large particles of $\mathit{St}>1$, point-particle Lagrangian approach is uniquely suited. For small particles, the Eulerian approaches are more computationally efficient. In the limit of $\mathit{St}<0.2$, the equilibrium Eulerian approximation (or mixture theory) can be adopted to study fine sediment suspension [22, 23]. For $0.2\le \mathit{St}\le O\left(1\right)$, both the Eulerian two-fluid and Lagrangian point-particle approaches are comparable. For cohesive sediment transport, the volumetric concentration is typically dilute and particle-particle interactions are important; the point-particle Lagrangian approach is better suited.

Two-phase CFD-DEM approach has gained popularity in flocculation dynamics of cohesive sediment in turbulent flows [17, 18, 24, 25, 26]. For primary particles of size smaller than the Kolmogorov length scale ($\eta$), the local average method of Anderson and Jackson [27] can be used to directly derive the fluid momentum equation as follows (equation set I in [28]):

$$\frac{\partial \alpha }{\partial t}+\nabla \cdot \left(\alpha \mathbf{u}\right)=0,$$E1

and

$$\frac{\partial \alpha \mathbf{u}}{\partial t}+\nabla \cdot \left(\alpha \mathbf{uu}\right)=-\frac{1}{\rho }\nabla p+\nabla \cdot \mathit{\tau }+{\mathbf{F}}_{\mathit{\text{inter}}}+{\mathbf{F}}_t,$$E2

where $\mathbf{u}$ is the fluid velocity, $\alpha$ is the fluid volume fraction, $\mathit{\tau }$ is the viscous stress tensor, and $p$ is the fluid-phase pressure. ${\mathbf{F}}_{\mathit{\text{inter}}}$ is the interphase force exchange term, and ${\mathbf{F}}_t$ is the external force to drive the turbulent flow. For homogeneous isotropic turbulence, a linear forcing term is commonly implemented [29]. The particle-fluid interaction force is the sum of all particle-fluid interaction forces acting on individual particles, including the drag force, pressure gradient force, viscous force due to fluid shear stress, virtual mass force, Basset force, and lift force [28]. In one-way coupled approach, ${\mathbf{F}}_{\mathit{\text{inter}}}$ term is neglected.

Figure 1a shows the CFD-DEM simulation results of flocculation dynamics of cohesive sediment in homogeneous isotropic turbulence. Initially, 200,000 mono-dispersed spherical particles were introduced into the flow field. Due to collision, flocs formed and grew in size till an equilibrium floc size distribution developed when the aggregation process balanced with the breakup process. The isosurfaces represent the turbulent coherent structures. By turning off the linear forcing term, the decaying turbulent flow relaminarized to quiescent flow. Figure 1b shows the gravitational settling of flocs in still water with gravitational acceleration pointing downward. Preferential orientation occurred with major axis of flocs aligning with the direction of gravity for large flocs. In CFD-DEM simulations, sediment is modeled as dispersed phase, and the evolution of floc properties can be directly obtained from model results (Figure 1b). Flocs can be identified by the connectivity of spheres, and the floc size can be characterized by the number of primary particles (spheres) consisting the floc. Floc structure can then be derived from the relative position of particles to the center of the floc, defined as

Figure 1.

(a) Flocculation dynamics of cohesive sediment in homogeneous isotropic turbulence. (b) Preferential orientation of flocs with respect to gravity in gravitational sedimentation with effects of long-range hydrodynamic interactions.

$${\mathbf{x}}_f=\frac{{\sum }_{i=1}^N{\mathbf{x}}_{p,i}}{N},$$E3

where ${\mathbf{x}}_f$ is the center of the floc, ${\mathbf{x}}_p$ is the center of the particles consisting the floc, $i$ is the label of the particle, and $N$ is the total number of the particles in the floc. To characterize the structure of flocs, fractal dimension and gyration radius have been commonly used. The fractal dimension $N_f$ is defined as

$$N={\left(\frac{L_f}{L_p}\right)}^{n_f},$$E4

where $L_f$ is the floc size, defined as the floc size along the longest principle axis, and $L_p$ is the size of the primary particle.

The gyration radius is defined as the root-mean-square distance of particles from the floc center as

$$R_g=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{\left|{\mathbf{x}}_{p,i}-{\mathbf{x}}_f\right|}^2}.$$E5

By turning off the linear forcing term, turbulence was allowed to decay over time, and the flow relaminarized to a quiescent state. Figure 1b shows the gravitational settling of flocs in the final still water with gravitational acceleration pointing downward. Preferential orientation occurred with the major axis of flocs aligning with the direction of gravity for large flocs.

2.2 Discrete element method for cohesive sediment

Interactions among primary particles play an important role in aggregation and breakup processes of cohesive sediment. To resolve particle-particle interactions, the sediment phase is modeled by the discrete element method (DEM), in which motions of individual particles are tracked following

$$\frac{d{\mathbf{x}}_i}{\mathit{dt}}={\mathbf{v}}_i,$$E6

$$m_i\frac{d{\mathbf{v}}_i}{\mathit{dt}}=\sum {\mathbf{F}}_i,$$E7

and

$$I_i\frac{d{\mathit{\omega }}_{\mathbf{i}}}{\mathit{dt}}=\sum {\mathbf{T}}_i.$$E8

where $\mathbf{x}$ is the position vector, $\mathbf{v}$ is the particle velocity vector, $\mathbf{F}$ is the force vector, and $m$ is the mass of the particle. The subscript “$i$” is the particle label. The force on particle “$i$” is the sum of the collision force (${\mathbf{F}}_c$) between particle $i$ and all other particles $j$, the hydrodynamic force (${\mathbf{F}}_f$), and the gravitational force (${\mathbf{F}}_g$) as ${\mathbf{F}}_i={\sum }_j{\mathbf{F}}_{c,\mathit{ij}}+{\mathbf{F}}_{f,i}+{\mathbf{F}}_{g,i}$. $I$ is the moment of inertia, $\omega$ is the angular velocity of the particle, and $\mathbf{T}$ is the hydrodynamic torque on the particle.

2.2.1 Cohesive force: DLVO and XDLVO theories

In CFD-DEM approach, the stickiness of cohesive sediment is modeled by pairwise interaction forces, which are given as functions of distance between particles. Parteli et al. [30] showed both types of attractive forces, namely, adhesion and nonbonded van der Waals forces, need to be taken into account to accurately model the particle interactions [30]. DLVO theory [31, 32] has been widely used in cohesive sediment studies to model the electrochemical interactions among particles, in which nonbonded attractive van der Waals and repulsive electrostatic double layer interactions are considered. In CFD-DEM simulations, the repulsive electrostatic interactions are often neglected and the van der Waals force is given by

$${\mathbf{F}}_{\mathit{vdW}}=\left\{\begin{array}{ll}\frac{A_H{R}_{\mathit{eff}}}{6D_{\mathit{min}}^2}{\mathbf{e}}_n& \text{if}{\delta }_n>0,\\ \frac{A_H{R}_{\mathit{eff}}}{6{\left({\delta }_n-D_{\mathit{min}}\right)}^2}{\mathbf{e}}_n& \text{if}-D_{\mathit{max}}\le {\delta }_n\le 0,\\ 0,& \text{if}{\delta }_n<-D_{\mathit{max}},\end{array}\right.$$E9

in which ${\delta }_n$ is the overlap distance (Figure 2), and ${\mathbf{e}}_n$ is the normal unit vector connecting centers of two spheres. $A_H$ is the Hamaker constant. $D_{\mathit{min}}$ is introduced to avoid the singularity of the equation at ${\delta }_n=0$. As the surface of the particle is not smooth, $D_{\mathit{min}}$ can also be interpreted as the surface roughness as there is always a minimal distance between two particles at contact. $D_{\mathit{max}}$ is the cutoff distance of the van der Waals interaction. The effective radius is defined as $R_{\mathit{eff}}=\frac{R_i{R}_j}{R_i+R_j}$ for two contacting particles of radius $R_i$ and $R_j$, respectively.

Figure 2.

(a) Schematic of a particle-particle collision handled by the soft-sphere discrete element method. The overlap distance ${\delta }_n$ represents the deformation in the contact region. (b) The spring-dashpot soft-sphere contact model in discrete element method.

The classic DLVO theory models the electrochemical interactions among particles and assumes the particle surfaces to be chemically inert, which is not valid for organic particles, such as bacteria. As the organic particles play important role in floc properties [8] and hence the flocculation dynamics, effects from microbial communities need to be considered in cohesive sediment transport. To model biological interactions, DLVO theory has been extended as XDLVO [33, 34] to include additional mechanisms, such as bridging effects and acid-base energy from hydrogen bonding [35, 36, 37, 38]. DLVO and XDLVO theories provide the ideal framework for modeling the complex cohesive particle interactions. In natural environments, several different mechanisms simultaneously contribute to the attractive force among cohesive sediment particles (i.e., stickiness) and affect the floc properties. To predict flocculation dynamics more accurately, measurement on particle interaction forces becomes important. With Atomic Force Microscopy (AFM), particle level force as a function of distance between two particles can be measured and directly implemented in DEM simulations [39, 40, 41]. However, detailed measurements on particle interaction force are still limited; systematic studies to quantify the interaction forces between particles under different ionic strengths, surface coasting, pH, biomaterials, and so forth need to be conducted.

2.2.2 Adhesive contact mechanics

In DEM framework, sediment particles are modeled as soft spheres, allowing small overlaps between two particles in contact. Contact mechanics models are often used to account for the collisional behavior of particles, which provide a variety of options for the normal, tangential, rolling, and twisting forces resulting from the contact. These forces are hence implemented to model the complex interactions among the cohesive sediment particles. To model adhesion of particles, two conflicting approximations have been developed, namely the Johnson-Kendall-Roberts (JKR) approximation [42] and the Derjaguin-Muller-Toporov (DMT) approximation [43]. In JKR model, adhesion forces outside the area of contact are neglected and elastic stresses at the edge of the contact are infinite. In DMT model, the geometry is assumed to be Hertzian and hence the adhesion force could not deform the surface even the adhesion is taken into account [44]. The dispute was about the magnitude of the pull-off force. Tabor [45] introduced a parameter $\mu$ [45], such that for small $\mu$, DMT should be used and for large $\mu$, JKR should be used. Recently, Greenwood [46] pointed out the errors in DMT thermodynamic method and argued that Hertzian geometry does not occur [46]. Therefore, we will focus on the JKR model in the following discussion.

For soft clay particles, the Johnson-Kendall-Roberts (JKR) model gives

$${\mathbf{F}}_{\mathit{ne},\mathit{jkr}}=\left(\frac{4{\mathit{Ea}}^3}{3R}-2\pi a^2\sqrt{\frac{4\mathit{\gamma E}}{\mathit{\pi a}}}\right){\mathbf{e}}_n,$$E10

where $a$ is the radius of the contact zone and is related to the overlap ${\delta }_n$ according to

$${\delta }_n=\frac{a^2}{R}-2\sqrt{\frac{\mathit{\pi \gamma a}}{E}},$$E11

in which $E$ is the Young’s Modulus, $R$ is the radius of the particle, and $\gamma$ is the surface energy density. The overlap between particle “$i$” and particle “$j$” is given as ${\delta }_n=R_i+R_j-\mid {\mathbf{x}}_i-{\mathbf{x}}_j\mid$. JKR model allows for a tensile force beyond contact (${\delta }_n<0$), up to a maximum of $3\mathit{\pi \gamma R}$. When two particles are not in contact initially, they will not experience this force until they come into contact (${\delta }_n>0$); then, as they move apart, they experience a tensile force up to $3\mathit{\pi \gamma R}$ till they lose contact. Therefore, this force can be used to define the yield strength of the floc.

In addition to the adhesive elastic force, a viscoelastic damping force model is used

$${\mathbf{F}}_{n,d}=-{\eta }_D{\mathbf{V}}_n,$$E12

where ${\eta }_D$ is the viscoelastic damping coefficient and ${\mathbf{V}}_n$ is the relative velocity along the direction of the unit vector ${\mathbf{e}}_n$. The total normal force is the sum of the adhesive elastic JKR and viscoelastic damping terms

$${\mathbf{F}}_n={\mathbf{F}}_{\mathit{ne},\mathit{kjr}}+{\mathbf{F}}_{n,d}.$$E13

The tangential force (${\mathbf{F}}_t$) is given by the Mindlin no-slip model [47] as

$${\mathbf{F}}_t=-min\left({\mu }_t{F}_{n0},|-k_{t}a\mathit{\xi }+{\mathbf{F}}_{t,d}|\right){\mathbf{e}}_t,$$E14

in which ${\mu }_t$ is the friction coefficient, $k_t$ is the elastic constant for tangential contact, and $\mathit{\xi }$ is the tangential displacement accumulated during the entire duration of the contact. The unit vector ${\mathbf{e}}_t$ is in the relative tangential velocity direction. ${\mathbf{F}}_{t,d}$ is the damping term for the tangential force, which follows the same general form as the normal damping force (Eq. (13)) but uses the relative velocity along the direction of the tangential vector ${\mathbf{e}}_t$. The normal force value $F_{n0}$ used to compute the critical force is given as

$$F_{n0}=\mid {\mathbf{F}}_n+2{\mathbf{F}}_{\mathit{\text{pulloff}}}\mid =\mid {\mathbf{F}}_{\mathit{ne},\mathit{jkr}}+6\mathit{\pi \gamma R}\mid .$$E15

Cohesive sediment flocs experience restructuring in turbulent flow, in which the relative position of primary particles consisting the flocs changes without breakup [48]. To account for floc restructuring, a rolling friction model of a pseudo-force formulation [49] is often implemented. The rolling pseudo-force model does not contribute to the total force on either particle but acts only to induce an equal and opposite torque on the particle, which allows adjustment of rolling displacement of the contacting pair and hence floc restructuring. The rolling friction model allows the adjustment of rolling displacement of the contacting pair. The rolling pseudo-force is computed analogously to the tangential force as

$${\mathbf{F}}_{\mathit{\text{roll}},0}=k_{\mathit{\text{roll}}}{\mathit{\xi }}_{\mathit{\text{roll}}}-{\gamma }_{\mathit{\text{roll}}}{\mathbf{v}}_{\mathit{\text{roll}}},$$E16

in which $k_{\mathit{\text{roll}}}$ is the elastic constant for rolling, ${\gamma }_{\mathit{\text{roll}}}$ is the damping constant for rolling, ${\mathit{\xi }}_{\mathit{\text{roll}}}$ is the rolling displacement, and ${\mathbf{v}}_{\mathit{\text{roll}}}$ is the relative rolling velocity [50].

A Coulomb friction criterion truncates the rolling pseudo-force if it exceeds a critical value:

$${\mathbf{F}}_{\mathit{\text{roll}}}=min\left({\mu }_{\mathit{\text{roll}}}{F}_{n0},|{\mathbf{F}}_{\mathit{\text{roll}},0}|\right){\mathbf{e}}_k,$$E17

where ${\mathbf{e}}_k$ is the direction of the pseudo-force. The rolling pseudo-force does not contribute to the total force on either particle, but acts only to induce an equal and opposite torque on each particle according to

$${\mathbf{T}}_{\mathit{\text{roll}},i}=\frac{R_i{R}_j}{R_i+R_j}\mathbf{n}\times {\mathbf{F}}_{\mathit{\text{roll}}},$$E18

$${\mathbf{T}}_{\mathit{\text{roll}},j}=-{\mathbf{T}}_{\mathit{\text{roll}},i}.$$E19

In DEM simulations, the time step is constrained by Young’s modulus of particles. To accelerate the simulation, particle stiffness is often reduced by several orders of magnitude. Numerous studies show stiffness can be reduced by orders of magnitude without altering collisional behavior of particles [51] and bulk transport rate of sediment [52]. With adhesive JKR model, the adhesive forces need to be scaled appropriately with a reduced stiffness [53, 54]. In addition, in point-particle Lagrangian approach, the flow around particles is not resolved; the lubrication force can affect the collisional behavior of the particles by decelerating two approaching particles. The enhanced dissipation from the lubrication forces can be modeled by a reduced effective restitution coefficient [55], which is defined as the ratio of the magnitude of relative velocity of particles after collision to the magnitude of the velocity before collision.

2.2.3 Hydrodynamic force: long-range hydrodynamic interactions

Maxey-Riley-Gatignol equation describes the motion of a small rigid sphere in a nonuniform flow [56] and is widely used in point-particle Lagrangian approach. The hydrodynamic force on a particle moving in the flow is the sum of the quasi-steady drag force, undisturbed flow force (sum of pressure gradient force and viscous force due to undisturbed fluid shear stress), added mass force, Basset history force, and lift force. For cohesive sediment transport studies, the Stokes number of the flocs is generally small. In most previous studies, only the drag force and pressure gradient force are considered for simplicity. The total hydrodynamic force on particle “$i$” is therefore given as

$${\mathbf{F}}_{\mathit{hd},i}={\mathbf{F}}_{d,i}+{\mathbf{F}}_{p,i},$$E20

in which ${\mathbf{F}}_d$ and ${\mathbf{F}}_p$ are the quasi-steady drag force and undisturbed flow force, respectively.

The undisturbed flow force experienced by the particle is calculated as

$${\mathbf{F}}_{p,i}=\left(-\nabla p+\nabla \cdot \tau \right)V_{p,i},$$E21

where the pressure gradient and stress divergence are interpolated to the particle center. In the current formulation, only the hydrostatic pressure is used to calculate the force $F_p$ for simplicity. The buoyancy force due to the hydrostatic pressure is ${\mathbf{F}}_{p,i}=-{\rho }_f\mathbf{g}{V}_{p,i}$, where $\mathbf{g}$ is the gravitational acceleration vector and $V_p=\pi D_p^3/6$ is the volume of the particle.

The accuracy of model prediction on particle motion strongly depends on the drag model. In the context of cohesive sediments, the Reynolds number of the primary particle (2–6 μm) and flocculi (6–50 μm) is smaller than 1. In a turbulent flow, the Reynolds number of aggregates that are smaller than the Kolmogorov scale is expected to be smaller than unity [57, 58]. Since most flocs tend to be smaller than the Kolmogorov scale, the hydrodynamic drag will be evaluated based on Stokes flow approximation for cohesive sediment transport.

The Free Draining Approximation (FDA) has been widely used to study deformation and breakup of aggregates in shear flow due to its simplicity and computational efficiency [59, 60]. In FDA, Stokes drag law is used to calculate the hydrodynamic force of individual particles that make up the floc, neglecting the hydrodynamic interactions between particles within the floc. Previous studies found FDA can yield qualitatively the same results as studies that fully resolve the flow around particles [61, 62]. In FDA, the shielding effect of neighbors is neglected, and hence, the drag force is overpredicted. However, it can be expected that the drag force on a particle will be influenced by its neighbors—those that are attached to it as part of the floc and those that are part of a nearby floc in case of floc-floc interaction. These effects are currently ignored, and their importance must be carefully examined.

The drag force ${\mathbf{F}}_d$ on particle “$i$” is given as

$${\mathbf{F}}_{d,i}=\frac{1}{2}\rho C_{D}A\mid \mathbf{u}-{\mathbf{v}}_i\mid \left(\mathbf{u}-{\mathbf{v}}_i\right).$$E22

$\rho$ is the fluid density, and $A=\pi D_p^2/4$ is the projected area of the spherical particle with $D_p$ as the diameter of the spherical particle. For very dilute flow with sediment concentration $\varphi \ll 0.1\%$, the standard drag coefficient $C_D$ for an individual particle is used, which is given as

$$C_D=\left\{\begin{array}{ll}\frac{24}{{\mathit{Re}}_p}\left(1+0.15{\mathit{Re}}_p^{0.687}\right),& \text{if}{\mathit{Re}}_p<1000\\ 0.44,& \text{otherwise},\end{array}\right.$$E23

where ${\mathit{Re}}_p=\mid \mathbf{u}-\mathbf{v}\mid D_p/\nu$ is the particle Reynolds number, which is mostly likely to be quite small.

In FDA approach, effects of neighbor particles are oversimplified by a correction to the drag coefficient, which is a function of volumetric concentration of sediment ($\varphi$). It has been claimed that this simplification works reasonably for sand transport with moderate particle Reynolds number as hindered settling. However, in the Stokes limit, the drag force on two nearby particles is smaller than the Stokes drag of a single particle due to the shielding effect [63]. The shielding effect is represented by ${\lambda }_s=1-\lambda$, with $\lambda =F_i/6\pi \mu {\mathit{aU}}_{\infty }$, and $F_i$ is the drag force on the particle.

While a particle-resolved simulation is computationally out of reach, traditional point-particle approach with free draining assumption completely ignores the shielding effect. Different approaches have been developed to model the shielding effects by neighbor particles. At the level of two approaching particles in an unbounded fluid, their hydrodynamic interaction can be completely described in terms resistance functions [64]. For a system of $N$ interacting particles, Brady [65, 66] introduced Stokesian Dynamics (SD) and accelerated Stokesian Dynamics (ASD), which have been applied to study breakage of colloidal aggregates in shear flow [67, 68]. In SD approach, the particles are modeled as rigid spheres, and long-range hydrodynamic interactions are modeled accurately except for close-contact of particles. Stokesian Dynamics simulations are computationally expansive as the computer time scales with total number of particles as $N^3$. SD approach has only been applied to study aggregation, breakup, and restructuring of aggregates (flocs) under simple shear in laminar flow.

There have been many theoretical studies on the sheltering effects by neighbor particles at Stokes limit. Filippov developed an exact approach using spherical harmonics to calculate the hydrodynamic force and torque on aggregates of spherical particles in the Stokes limit [69]. However, the exact method of Filippov is also computationally expensive and scales as $N^3$ for a floc containing $N$ particles. On the other hand, data-driven approach has been quite successful in predicting the shielding effects by neighbor particles. Kim and Lee [70] developed an empirical formula for the shielding effect of neighbor particles, which unlike Filippov’s exact approach does not require solving a large linear system. Ma et al. [71] developed binary (pairwise) and trinary interaction models enabled by graph neural network to accurately account for the effects of neighboring particles in a computationally efficient manner.

Figure 3 shows numerical simulation results using two different drag models. By including the shielding effect from neighbor particles, the Kim-Lee (K-L) model is able to predict the dependence of floc settling velocity as a function of floc size (represented as the number of primary particles $n_f$ consisting the floc), while the floc settling velocity is constant (single particle settling velocity) in the FDA approach. In addition, the shielding effect from neighbor particles can also affect the breakup and restructuring by turbulence, as it shields the particle from external flow, allowing the growth of larger flocs. The equilibrium floc size distributions show a significant difference. The K-L model predicts a bimodal distribution with two peaks, while the FDA model predicts a single peak. The first peak appears at around $n_f=10$ in both cases, while a second peak shows at $n_f$ around 50 with K-L model.

Figure 3.

(a) Floc settling velocity. (b) Floc size distribution.

3. Applications and future research

The advantage of the point-particle Lagrangian approach is its capability to resolve interactions among particles, which makes it the ideal candidate to study the flocculation dynamics of cohesive sediment. In many previous studies, the particle-unresolved DNS-DEM approach has been implemented, in which all scales up to the Kolmogorov length scale are fully resolved; however, flows around individual particles are not resolved. Zhao et al. applied the unresolved DNS-DEM model to study flocculation dynamics under homogeneous isotropic turbulence and identified the equilibrium floc size distribution depends on the cohesive number, which is defined as the ratio of floc strength to the turbulence intensity [72]. They also proposed a new flocculation model with a variable fractal dimension that can predict the temporal evolution of the floc size and shape. Zhao et al. [73] applied the one-way coupled CFD-DEM model framework and found that the intermediate shear gives rise to the largest flocs and proposed a simple model for the equilibrium floc size [73]. Yu et al. [24] applied the one-way coupled CFD-DEM model to investigate the reshaping of flocs by turbulence, which leads to more compact flocs at equilibrium [48]. They identified the cohesion number defined as the ratio of floc strength to the turbulent shear rate controls the reshaping mechanisms. With low cohesion number, the breakup-regrowth dominates the reshaping of turbulence dominates as turbulent shear can effectively break the flocs. With high cohesion number, flocs are resilient to small-scale turbulent shear; floc restructuring by hydrodynamic drag plays an important role, in which flocs reshape without breakup. The one-way coupled simulation results have also been used to develop aggregation and breakup kernels [24, 48] by tracking the time evolution of the flocs (Figure 4). Yao and Capecelatro [26] performed two-way coupled model simulations to study the breakup of aggregates in homogeneous isotropic turbulence [26]. They also showed the cohesion number plays a key role in the breakup process and conducted a scaling analysis to develop a phenomenological model of the breakup process.

Figure 4.

(a) The probability density function (PDF) that a floc of size $n_f$ becomes a floc of size $n_f=8$ in a nondimensional time period of 2 (blue circles) and the PDF that a floc of size $n_f=8$ evolves into a floc of size $n_f$ over the same period (red pluses). (b, c) the same as frame (a) for $n_f=16$ and 24.

CFD-DEM approach tracks motions of individual particles; therefore, we can study how turbulence affects floc structures. Different geometric measures have been used to characterize the floc structure, including fractal dimension, gyration radius, and sphericity. Previous studies [24, 48] have identified the compaction of flocs in turbulent flow, in which the gyration radius decreases and fractal dimension increases for flocs of given size (defined by the number of particles in the floc). For instance, the most stable structure for trimers (flocs consisting of three particles) is the equilateral triangle (Figure 5) with the interior angle of the triangle formed by three particles of $\pi /3$. At early stage of the simulation, there were many linear trimers with interior angle of 0 or $\pi$; however, we observed the dominant peak at $\pi /3$ at the equilibrium stage, suggesting the trimers evolved into the most stable and compact structure by turbulence.

Figure 5.

(a) Flocs at early stage of simulation, (b) flocs at the equilibrium stage, and (c) the corresponding distribution of the interior angles of trimers (a floc consisting of three primary particles).

Looking at the future, there is great potential in point-particle Lagrangian modeling framework that can resolve the complex particle-particle and turbulence-particle interactions. More observational data sets on force versus distance between particles are needed to accurately model interactions among particles, including effects of temperature, ionic strength, bacterial strength, and so on. Previous studies did not consider the long-range hydrodynamic interactions among particles mainly due to computational cost. Efficient algorithms to account for the sheltering effects from neighbor particles need to be developed. With both numerical simulation results and observational data sets, detailed scaling analyses need to be conducted to develop parameterization of important flocculation processes, including aggregation, breakup, and restructuring processes. However, we must recognize the role of high-resolution laboratory experiments as a ground truth for model validation, as well as the role of theoretical frameworks to develop more accurate phenomenological models of flocculation dynamics with more accurate representations of physical processes.

4. Conclusion

Two-phase CFD-DEM framework has gained traction in cohesive sediment transport due to its capability of resolving particle-particle interactions. In this chapter, we have presented an overview on recent advances in the state-of-art two-phase CFD-DEM approach on cohesive sediment transport in turbulent environments. We outlined the fundamental methodology of CFD-DEM approach, including models for the CFD-DEM coupling and interparticle interactions, which is followed by a quick discussions of the applications to cohesive sediment transport. The two-phase CFD-DEM approach tracks motions of individual particles and hence can provide valuable information of floc evolution under turbulent flow. These information can be used to develop more accurate parameterizations of aggregation and breakup models for population balance approach for more precise prediction of cohesive sediment transport in natural environments. At last, we identified research gaps on the drag force parameterization, which could guide future research. To overcome the flaws in the widely used free-draining approximation, shielding effects from neighbor particles need to be considered.

Acknowledgments

XY gratefully acknowledges supports by the U.S. Army Engineer Research and Development Center (W912HZ-21-2-0035).