2.1. Particle motion in homogeneous isotropic turbulence

We consider the one-way coupled motion of suspended cohesive particles in three-dimensional, incompressible homogeneous isotropic turbulence. The motion of the single-phase fluid with constant density $\rho _f$ and kinematic viscosity $\nu$ is governed by

(2.1)$$\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_{\boldsymbol{f}} = 0, \end{gather}$$

(2.2)$$\begin{gather}\frac{\partial \boldsymbol{u}_{\boldsymbol{f}}}{\partial t} + (\boldsymbol{u}_{\boldsymbol{f}} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{u}_{\boldsymbol{f}} ={-} \frac{1}{\rho_f} \boldsymbol{\nabla} p + \nu \nabla^2 \boldsymbol{u}_{\boldsymbol{f}} + \boldsymbol F_{tur} , \end{gather}$$

where $\boldsymbol {u}_{\boldsymbol {f}} = (u_f, v_f, w_f)^\textrm {T}$ denotes the fluid velocity vector and $p$ indicates the hydrodynamic pressure. We employ the spectral approach of Eswaran & Pope (Reference Eswaran and Pope1988) to obtain the forcing term $\boldsymbol F_{tur}$, which generates and maintains statistically stationary turbulence, as implemented in Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015). Here, $\boldsymbol F_{tur}$ is non-zero only in the low-wavenumber band where the wavenumber vector $|\boldsymbol \kappa | < \kappa _f$, with $\kappa _f = 2.3\kappa _0$ and $\kappa _0 = 2 {\rm \pi}/ L_0$, with $L_0$ denoting the length of the physical domain. The origin $\boldsymbol \kappa = 0$ is not forced. In addition to the cutoff wavenumber $\kappa _f$, the random forcing process is governed by the dimensionless parameter $D_s = \sigma ^2 T_0 L_0^4 / \nu ^3$, where $\sigma ^2$ and $T_0$ indicate the variance and the time scale of the random process, respectively. Regarding the details of evaluating $\boldsymbol F_{tur}$ from $\kappa _f$ and $D_s$, we refer the reader to the original work by Eswaran & Pope (Reference Eswaran and Pope1988).

We approximate each primary suspended particle $i$ as a sphere moving with translational velocity $\boldsymbol {u}_{p,i} = (u_{p,i}, v_{p,i}, w_{p,i})^\textrm {T}$ and angular velocity ${\boldsymbol \omega }_{p,i}$. These are obtained from the linear and angular momentum equations

(2.3)$$\begin{gather} m_p \frac{\mathrm{d}\boldsymbol{u}_{p,i}}{\mathrm{d}t} = {\boldsymbol F}_{d,i} + \underbrace{\sum_{j=1,j \ne i}^{N}(\boldsymbol F_{con,ij} + \boldsymbol F_{lub,ij} + \boldsymbol F_{coh,ij})}_{{\boldsymbol F}_{c,i}} , \end{gather}$$

(2.4)$$\begin{gather}I_p \frac{\mathrm{d}{\boldsymbol \omega_{p,i}}}{\mathrm{d}t} = \underbrace{\sum_{j=1,j \ne i}^{N}(\boldsymbol T_{con,ij} + \boldsymbol T_{lub,ij})}_{{\boldsymbol T}_{c,i}} , \end{gather}$$

where the primary particle $i$ moves in response to the Stokes drag force $\boldsymbol F_{d,i} = -3 {\rm \pi}D_p \mu (\boldsymbol {u}_{p,i} - \boldsymbol {u}_{f,i})$, and the particle–particle interaction force $\boldsymbol F_{c,i}$. Buoyancy is not considered here, so that we can investigate the effects of particle inertia in isolation. We only consider primary particles that are larger than 2 $\mathrm {\mu }$m (cf. table 1), so that a suitably defined Péclet number measuring the relative importance of hydrodynamic and Brownian forces is sufficiently large for their Brownian motion to be negligible (Partheniades Reference Partheniades2009; Biegert et al. Reference Biegert, Vowinckel and Meiburg2017; Chen et al. Reference Chen, Li and Marshall2019; Vowinckel et al. Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019). Here, $\boldsymbol {u}_{p,i}$ indicates the particle velocity evaluated at the particle centre, $\boldsymbol {u}_{f,i} = \sum _{1}^{N_i}(\phi _{i,k} \boldsymbol {u}_{f,k})$ represents the fluid velocity averaged over the volume of particle $i$, where $N_{i}$ denotes the number of Eulerian grid cells covered by particle $i$, $\boldsymbol {u}_{f,k}$ is the fluid velocity at the centre of the grid cell $k$ and $\phi _{i,k}$ is the volume fraction of the particle $i$ in the grid cell $k$. We remark that the above implies that the diameter $D_p$ of the primary particle should be larger than the grid spacing $h$. This avoids the need for interpolating the fluid velocity within one grid cell, which would be required if $D_p < h$ (Chen et al. Reference Chen, Li and Marshall2019). Also, $m_p$ denotes the particle mass, $\mu$ the dynamic viscosity of the fluid and $N$ the total number of particles in the flow. We assume all particles to have the same diameter $D_p$ and density $\rho _p$. The parameter $\boldsymbol F_{c,i}$ accounts for the direct contact force $\boldsymbol F_{con,ij}$ in both the normal and tangential directions, as well as for short-range normal and tangential forces due to lubrication $\boldsymbol F_{lub,ij}$ and cohesion $\boldsymbol F_{coh,ij}$, where the subscript $ij$ indicates the interaction between particles $i$ and $j$. Also, $I_p = {\rm \pi}\rho _p D_p^5 / 60$ denotes the moment of inertia of the particle and $\boldsymbol T_{c,i}$ represents the torque due to particle–particle interactions, where we distinguish between direct contact torque $\boldsymbol T_{con,ij}$ and lubrication torque $\boldsymbol T_{lub,ij}$. Within a large floc, we account for all of the individual binary particle interactions.

Table 1. Non-dimensionalization employed in the present work: the characteristic values for length, velocity and density are $L_0 = 125D_p = 6.25 \times 10^{-4} \ \textrm {m}$, $U_0 = 8 \ \textrm {m}\ \textrm{s}^{-1}$ and $\rho _f = 1000 \ \textrm {kg}\ \textrm {m}^{-3}$, respectively.

The lubrication force $\boldsymbol F_{lub,ij}$ is accounted for based on Cox & Brenner (Reference Cox and Brenner1967) as implemented in Zhao et al. (Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020). We note that, although the present study is limited to monodisperse particles, polydisperse particle–particle interactions can be taken into account by an effective radius $R_{eff}=R_p R_q / (R_p + R_q)$, where $R_p$ and $R_q$ are the radii of two interacting spheres. Following Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017), the collision force $\boldsymbol F_{con,ij}$ is represented by a nonlinear spring–dashpot model in the normal direction, while the tangential component is modelled by a linear spring–dashpot model capped by the Coulomb friction law to account for zero-slip rolling or sliding of particles. We note that the tangential component of the contact force depends on the surface roughness, a prescribed restitution coefficient $e_{dry} = 0.97$ and a friction coefficient $e_{fri} = 0.15$ are implemented to yield adaptively calibration for every collision as described by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017). The cohesive force $\boldsymbol F_{coh,ij}$, which reflects the combined influence of the attractive van der Waals force and the repulsive electrostatic force, is based on the work of Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019), where additional details and validation results are provided. The model assumes a parabolic force profile, distributed over a thin shell surrounding each primary particle. Hence, the cohesive force between primary particles extends over a finite range, so that it is felt by the particles even before they come into direct contact. We consider two primary particles to be part of the same floc when their surface distance is smaller than half the range of the cohesive force, as implemented in Zhao et al. (Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020). We remark that, based on (2.3) and (2.4), the configuration of the primary particles within a floc can change with time in response to fluid forces, since the cohesive bonds are not rigid. Specifically, the contact points on the surface of the primary particles are not fixed, so that the primary particles can rotate individually within a floc.

2.2. Non-dimensionalization

In order to render the above governing equations dimensionless, we consider primary particles with diameter $D_p = 5 \ \mathrm {\mu } \textrm {m}$, which represents a typical value for clay or fine silt. The cubic computational domain has an edge length $L_0 = 125D_p = 6.25 \times 10^{-4} \ \textrm {m}$. As time scale of the random turbulent forcing process we select $T_0 = 7.81 \times 10^{-5} \ \textrm {s}$. By choosing $L_0$, $T_0$ and $\rho _f = 1000 \ \textrm {kg}\ \textrm {m}^{-3}$ as the characteristic length, time and density scales, we obtain the characteristic velocity scale $U_0 = L_0/T_0 = 8 \ \textrm {m}\ \textrm {s}^{-1}$, which is similar to values employed in previous investigations (Chen et al. Reference Chen, Li and Marshall2019; Chen & Li Reference Chen and Li2020). We employ $L_0$ and $U_0$ to define the turbulence Reynolds number $Re = L_0 U_0 \rho _f/ \mu$.

The dimensionless continuity and momentum conservation equations can then be expressed as

(2.5)$$\begin{gather} \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}_f =0, \end{gather}$$

(2.6)$$\begin{gather}\frac{\partial \tilde{\boldsymbol{u}}_f }{\partial \tilde t} + (\tilde {\boldsymbol u}_f \cdot \tilde{\boldsymbol{\nabla}}) \tilde {\boldsymbol u}_f ={-} \tilde{\boldsymbol{\nabla}} \tilde p + \frac{1}{Re} \tilde {\nabla}^2 \tilde {\boldsymbol u}_f + \tilde {\boldsymbol F}_{tur} \ , \end{gather}$$

while the dimensionless equations of motion for the primary cohesive particles take the form

(2.7)$$\begin{gather} \tilde m_p \frac{\mathrm{d} \tilde{\boldsymbol{u}}_{p,i}}{\mathrm{d} \tilde t} = \underbrace{- \frac{3 {\rm \pi}\tilde D_p (\tilde{\boldsymbol{u}}_{p,i} - \tilde{\boldsymbol{u}}_{f,i})}{Re}}_{\tilde{\boldsymbol{F}}_{d,i}} + \sum_{j=1,j \ne i}^{N}(\tilde{\boldsymbol{F}}_{con,ij} + \tilde{\boldsymbol{F}}_{lub,ij} + \tilde{\boldsymbol{F}}_{coh,ij} ) , \end{gather}$$

(2.8)$$\begin{gather}\tilde I_p \frac{\mathrm{d} \tilde {\boldsymbol \omega}_{p,i}}{\mathrm{d} \tilde t} = \sum_{j=1,j \ne i}^{N}(\tilde {\boldsymbol T}_{con,ij} + \tilde {\boldsymbol T}_{lub,ij}). \end{gather}$$

Here, dimensionless quantities are denoted by a tilde. The dimensionless particle mass is defined as $\tilde m_p = {\rm \pi}\tilde D_p^3 \tilde \rho _s / 6$, the moment of inertia $\tilde I_p = {\rm \pi}\tilde \rho _s \tilde D_p^5 / 60$ and the density ratio $\tilde \rho _s = \rho _p / \rho _f$. The dimensionless direct contact and lubrication forces, $\tilde {\boldsymbol {F}}_{con,ij}$ and $\tilde {\boldsymbol {F}}_{lub,ij}$, are accounted for based on Zhao et al. (Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020), while the dimensionless cohesive force $\tilde {\boldsymbol {F}}_{coh,ij}$ is defined as

(2.9)\begin{equation} \tilde{\boldsymbol{F}}_{coh,ij} = \left\{\begin{array}{@{}ll} - 4 Co \dfrac {\tilde \zeta_{n,ij}^2 - \tilde h_{co} \tilde \zeta_{n,ij}}{\tilde h_{co}^2} \boldsymbol n, & \tilde \zeta_{min} < \tilde \zeta_{n,ij} \leqslant \tilde h_{co} ,\\ 0, & \textrm{otherwise}. \end{array}\right. \end{equation}

Here, $\tilde \zeta _{min} = 0.0015 \tilde D_p$ and $\tilde h_{co} = 0.05 \tilde D_p$ represent the surface roughness of the particles and the range of the cohesive force, respectively. Also, $\boldsymbol n$ represents the outward-pointing normal on the particle surface, while $\tilde \zeta _{n, ij}$ is the normal surface distance between particles $i$ and $j$. The cohesive number $Co$ indicates the ratio of the maximum cohesive force $\|{\boldsymbol F_{coh,ij}}\|$ at $\tilde \zeta _{n,ij} = \tilde h_{co}/2$ to the characteristic inertial force

(2.10)\begin{equation} Co = \frac{\max(\|{\boldsymbol F_{coh,ij}}\|)}{U_0^2 L_0^2 \rho_f} = \frac{A_H D_p}{16 h_{co} \zeta_{0} } \frac{1}{U_0^2 L_0^2 \rho_f } , \end{equation}

where the Hamaker constant $A_H$ is a function of the particle and fluid properties, and the characteristic distance $\zeta _0 = 0.00025D_p$. Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019) provide representative values of various physicochemical parameters such as $A_H$, salt concentration and grain size of the primary particles for common natural systems. The present numerical approach for simulating the dynamics of cohesive sediment has been employed to predict the flocculation in simple vortical flow fields, and it was successfully validated with experimental data in our earlier work (Zhao et al. Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020).

To summarize, the simulations require as direct input parameters the turbulence Reynolds number $Re$, the characteristic parameter of the random turbulent forcing process $D_s$, the dimensionless particle diameter $\tilde D_p$, the total number of particles $N$, the density ratio $\tilde \rho _s$ and the cohesive number $Co$. As we will discuss below, $Re$ and $D_s$ can equivalently be expressed by the shear rate $G$ of the turbulence, cf. (3.1), and the Stokes number $St$ defined by (4.1). A list of the relevant dimensionless parameters is provided in table 1. We remark that due to computational limitations the simulations consider Kolmogorov scales that are somewhat smaller than typical field values, and turbulent shear rates that are larger than field values. Hence, the ratio of the Kolmogorov length scale to the primary particle size takes values up to 3.3 in the simulations, as compared with values up to $O(10)$ under typical field conditions. For convenience, the tilde symbol will be omitted henceforth.

3.1. Computational set-up

The triply periodic computational domain $\varOmega$ has a dimensionless size of $L_x \times L_y \times L_z = 1 \times 1 \times 1$, with the number of grid cells $N_x \times N_y \times N_z = 128 \times 128 \times 128$. This relatively modest number of grid points enables us to conduct the simulations over sufficiently long times for the flocculation and break-up processes to reach an equilibrium state (Tran, Kuprenas & Strom Reference Tran, Kuprenas and Strom2018), and it is in line with the earlier study of Chen et al. (Reference Chen, Li and Marshall2019). As mentioned above, we set the diameter $D_p$ of the primary particles moderately larger than the grid size $h = L_x/N_x$, at a constant value $D_p/h=1.024$.

Before introducing the particles into the flow, we simulate the single-phase turbulence until it reaches a statistically stationary state. Table 2 gives an overview of the physical parameters for the simulations conducted within the present investigation. Here, the Kolmogorov length scale and the root-mean-square velocity are defined as $\eta = 1/ (Re^3 \epsilon )^{1/4}$ and $u_{rms} = (2 k /3)^{1/2}$, respectively, where $\epsilon$ and $k$ denote the domain-averaged dissipation rate and kinetic energy of the fluctuations. The Taylor Reynolds number $Re_{\lambda } = \lambda u_{rms} Re$ of the turbulence is based on the Taylor microscale $\lambda = \sqrt {15} \, u_{rms}/(Re \, \epsilon )^{1/2}$. To provide a more complete quantitative description of the fluid shear, we define the vorticity fluctuation amplitude

(3.1)\begin{equation} G = \frac{1}{Re \, \eta^2} , \end{equation}

which can also be regarded as the turbulent shear rate. For additional details with regard to these quantities, we refer the reader to Pope (Reference Pope2001).

Table 2. Physical parameters of the single-phase turbulence simulations. As input parameters we specify the fluid Reynolds number $Re = L_0 U_0 / \nu$ and the characteristic parameter of the random turbulent forcing process $D_s = \sigma ^2 T_0 L_0^4 Re^3$. The simulation then yields the Taylor Reynolds number $Re_{\lambda } = \lambda u_{rms} Re$, the Kolmogorov scale $\eta$, the average root-mean-square velocity $u_{rms}$ and the shear rate $G = 1/ (Re \, \eta ^2)$. All of these output quantities are obtained by averaging over space and time, after a statistically stationary state has evolved.

3.2. Turbulence properties for different $Re_{\lambda }$

One key goal of the present investigation is to study the flocculation of primary particles whose diameter $D_p$ is smaller than the Kolmogorov length scale $\eta$. Since the particle diameter needs to be larger than the grid spacing, and the number of grid points is limited, suitable values of $\eta$ require a relatively low $Re_{\lambda }$. On the other hand, it is known that for $Re_{\lambda } \leq O(50)$ the turbulence may not be fully developed and isotropic (Mansour & Wray Reference Mansour and Wray1994). Hence this section presents a more detailed discussion of the turbulence properties for $Re_{\lambda } \leq O(50)$.

Figure 1 shows the time-dependent evolution of the box-averaged Kolmogorov length $\eta$, the root-mean-square velocity $u_{rms}$, the Taylor Reynolds number $Re_{\lambda }$ and the shear rate $G$ for cases Tur1 and Tur8, which have time-averaged Taylor Reynolds numbers of 9.72 and 50.34, respectively. Both cases are seen to reach statistically stationary states. We note that while case Tur8 results in $\eta /h = 0.6656$, Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015) demonstrated the validity of the current turbulent forcing approach even when the Kolmogorov length is smaller than the grid spacing. Snapshots of the vorticity modulus in a slice of the computational domain are shown in figure 2. They exhibit the intermittent multiscale patterns featuring eddies of different size along with thin filaments that are typical for turbulence.

Figure 1. Temporal evolution of box-averaged turbulence properties for cases Tur1 and Tur8 in table 2: (a) Kolmogorov length scale $\eta$; (b) root-mean-square velocity $u_{rms}$; (c) Taylor Reynolds number $Re_{\lambda }$; and (d) shear rate $G$. A statistically stationary state is seen to evolve for all quantities.

Figure 2. Representative snapshots of the vorticity modulus normalized by the vorticity fluctuation amplitude $G$, shown in the plane $z = 0.5$. (a) Case Tur1 and (b) case Tur8.

Figure 3 shows the temporal evolution of the domain-averaged magnitude of the velocity components $\langle |u_f| \rangle _{\varOmega }$, $\langle |v_f| \rangle _{\varOmega }$ and $\langle |w_f| \rangle _{\varOmega }$. During the statistically stationary state the three components are seen to oscillate around similar average values for both Tur1 and Tur8, which indicates that the flow is isotropic to a good approximation.

Figure 3. Temporal evolution of box-averaged magnitude of the fluid velocity components: (a) case Tur1 and (b) case Tur8. The flow is seen to be isotropic to a good approximation.

We define the instantaneous kinetic energy components in Fourier space, $E_{11}(\kappa )$, $E_{22}(\kappa )$ and $E_{33}(\kappa )$, as

(3.2a)$$\begin{gather} \int_0^\infty E_{11}(\kappa)\, \textrm{d} \kappa = \left\langle \frac{u_f \cdot u_f}{2} \right\rangle_{\varOmega} , \end{gather}$$

(3.2b)$$\begin{gather}\int_0^\infty E_{22}(\kappa)\, \textrm{d} \kappa = \left\langle \frac{v_f \cdot v_f}{2} \right\rangle_{\varOmega} , \end{gather}$$

(3.2c)$$\begin{gather}\int_0^\infty E_{33}(\kappa)\, \textrm{d} \kappa = \left\langle \frac{w_f \cdot w_f}{2} \right\rangle_{\varOmega} , \end{gather}$$

where $\kappa = |\boldsymbol \kappa |$ denotes the wavenumber. Figure 4 shows the time-averaged one-dimensional energy spectra. Only the wavenumbers below the cutoff wavenumber ($\kappa _f$, shown as vertical dashed lines in figure 4) are forced. The shapes of the energy spectra are in qualitative agreement with those obtained by Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015, p. 10) for higher values of $Re_{\lambda } \approx 60$. We conclude that the present forcing scheme yields statistically steady flow fields that are approximately isotropic for the current range of $Re_{\lambda }$-values.

Figure 4. Time-averaged one-dimensional energy spectra. The vertical dashed lines indicate the respective cutoff wavenumber of the turbulence forcing scheme, $\kappa _f \eta = 2.3 (2{\rm \pi} /L_x) \eta$. (a) Case Tur1 and (b) case Tur8. The spectra confirm that the statistically stationary flow fields are approximately isotropic.

4.1. One-way coupling

Once the single-phase turbulence reaches the statistically stationary regime, $N = 10\,000$ identical cohesive particles with diameter $D_p=0.008$ are randomly distributed throughout the domain, resulting in a particle volume fraction $\phi _p = 0.268\,\%$. Initially all particles are at rest and separated by a distance larger than the cohesive range $h_{co}$. To improve the statistics, we carry out repeated simulations for different random initial conditions, as the simulation results are statistically independent of the initial particle placement. The simulations to be discussed in the following are one-way coupled, so that the particles do not modify the background turbulence. Bosse, Kleiser & Meiburg (Reference Bosse, Kleiser and Meiburg2006) find that particle loading can modify the turbulence statistics even for volume fractions as low as $10^{-5}$, so that we expect two-way coupling effects to have an impact on the flocculation process even in moderately dilute flows. In addition, even for globally dilute flows the local volume fraction inside a floc will be $O(1)$, so that the one-way coupled assumption generally will not hold inside a floc. However, fully two-way coupled simulations for sufficiently many particles to obtain reliable statistical information, and for sufficiently long times to explore the balance between aggregation and breakup during the equilibrium stage, are not feasible on currently available supercomputers. Our assumption of one-way coupling hence limits the volume and mass fractions that we can reasonably consider. On the other hand, the current simulations and their comparisons to experimental observations are useful in that they help address the question as to which aspects of flocculation are governed by a one-way coupled dynamics, and which other aspects require a fully two-way coupled dynamics. As we will see below, for the range of physical parameters listed in table 1, even one-way coupled simulations are able to reproduce several experimentally observed statistical features of flocculation dynamics.

We adopt a multiscale time-stepping approach in which the fluid motion is calculated with a time step $\Delta t$ based on the criterion that the Courant–Friedrichs–Lewy number $\textrm {CFL} \leq 0.5$. The particle motion, on the other hand, is evaluated with a much smaller time step $\Delta t_p = \Delta t / 15$. Since the computational approach maintains a contact duration of $T_c = 10 \Delta t = 150 \Delta t_p$ (Biegert et al. Reference Biegert, Vowinckel and Meiburg2017), each particle collision is effectively resolved by 150 substeps, at the price of a marginal increase in the computational cost. The dynamics of the primary particles is characterized by the Kolmogorov-scale Stokes number

(4.1)\begin{equation} St = \frac{\rho_s}{18} \frac{D_p^2}{\eta^2} Re_{\eta} , \end{equation}

where the Kolmogorov Reynolds number $Re_{\eta } = \eta u_{rms} \, Re$. Since the particle diameter $D_p$ is constant throughout the present investigation, $St$ depends on the density ratio $\rho _s$ and the fluid properties. A particle with a small Stokes number tends to follow the fluid motion, while the dynamics of a particle with a large Stokes number is dominated by its inertia, so that it tends to continue along its initial direction of motion.

Table 3 summarizes the physical parameters of the simulations that we conducted. Following our analysis from § 2.2 and the examples given in Appendix A of Vowinckel et al. (Reference Vowinckel, Withers, Luzzatto-Fegiz and Meiburg2019), these values correspond to primary silica particles with a grain size of fine to medium silt in ocean water. In the following, we will investigate how the flocculation dynamics is influenced by the cohesive number $Co$, the Stokes number $St$ and the shear rate $G$. We remark that the density ratio $\rho _s$ and the size ratio $\eta /D_p$ are implicitly accounted for by $St$ and $G$.

Table 3. Physical parameters of the flocculation simulations. We separately investigate (bold) the influence of the cohesive number $Co$ (based on Flo1–5), the shear rate $G$ (Flo6–9) and the Stokes number $St$ (Flo10–13). The effects of $\rho _s$ and $\eta /D_p$ are implicitly accounted for by $St$ and $G$.

4.2. Flocculation and equilibrium stages

When the surface distance between two particles is smaller than half the range of the cohesive force, $h_{co}/2$, we consider these particles to be part of the same floc. Hence, in terms of a physical force balance breakage occurs when the net force pulling the particles apart is sufficiently strong to overcome the maximum of the cohesive force holding the particles together. An individual particle is considered to be the smallest possible floc. Figure 5(a) shows the evolution of the number of flocs $N_f(t)$ with time for the representative case Flo9, with $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $G = 0.62$, $\rho _s = 2.65$ and $\eta /D_p = 2.25$. As a result of flocculation, $N_f$ decreases rapidly with time from its initial value of 10 000, before levelling off around a constant value $N_{f,eq}$ that reflects a stable equilibrium between aggregation and breakage. This tendency of $N_f$ is consistent with our previous observation of flocculation in steady cellular flow fields (Zhao et al. Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020). Consequently, we can identify two pronounced stages of the flow, viz. an initial flocculation stage and a subsequent equilibrium stage. We define the end of the flocculation stage, i.e. the onset of the equilibrium stage, as the time $t_{eq}$ when $N_f$ first equals $N_{f,eq}$. Figure 5(b) shows separately the number of flocs with $N_{p} = 1, 2, 3$ and more than three primary particles. While the number of flocs with two or three particles initially grows quickly, they soon reach a peak and subsequently decline, as more flocs of larger sizes form. Toward the end of the flocculation stage, a stable equilibrium of the different floc sizes begins to emerge, although the distribution of flocs with different numbers of primary particles is still changing slowly.

Figure 5. (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 initially grow and subsequently decay, as increasingly many flocs with three or more particles form. (c) Temporal evolution of the fraction of flocs that maintain their identity ($\theta _{id}$), add primary particles ($\theta _{ad}$) or undergo breakage ($\theta _{br}$) over the time interval $\Delta T = 3$. All results are for case Flo9 with $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $G = 0.62$, $\rho _s = 2.65$, $\eta /D_p = 2.25$.

In order to gain insight into the dynamics of floc growth and breakage, we keep track of the evolution of three different types of flocs over a suitably specified time interval $\Delta T$: (a) those flocs that maintain their identity, i.e. they consist of the same primary particles at the start and the end of the time interval; (b) those that add additional primary particles while keeping all of their original ones; and (c) all others, i.e. all those who have undergone a breakage event during the time interval. We denote the fractions of these respective groups as $\theta _{id} = N_{f,id} / N_f$, $\theta _{ad} = N_{f,ad} / N_f$, and $\theta _{br} = N_{f,br} / N_f$. It follows that

(4.2)\begin{equation} \theta_{id} + \theta_{ad} + \theta_{br} = 1. \end{equation}

We found that a value of $\Delta T = 3$ is suitable for obtaining insight into the dynamics of the flocculation process, as it allows most of the flocs to maintain their identity during the time interval, while smaller but still significant numbers undergo primary particle addition or breakage. Figure 5(c) shows the evolution of $\theta _{id}$, $\theta _{ad}$ and $\theta _{br}$ for case Flo9. After an initial transient stage, all three fractions reach statistically steady states. Interestingly, even during the equilibrium stage when $N_f \approx \textrm {const.}$, we observe that $\theta _{ad} \ne \theta _{br}$. This reflects events such as when one floc breaks into three smaller parts, two of which then merge with other flocs. Here, the total number of flocs remains unchanged at three, in spite of only one break-up but two particle addition events.

4.3. Evolution of floc size and shape

While the number of primary particles in a floc, $N_{p}$, provides a rough measure of its size, flocs with identical values of $N_p$ can have very different shapes. In order to capture this effect, we define the characteristic diameter $D_f$ of the floc, also known as the Feret diameter, as

(4.3)\begin{equation} D_f = 2 \textrm{max}(\| \boldsymbol x_{p,i} - \boldsymbol x_c \|) + D_p,\quad 1 \leqslant i \leqslant N_{p}, \end{equation}

as well as its gyration diameter $D_g$ (Chen et al. Reference Chen, Li and Marshall2019),

(4.4)\begin{equation} D_g = \left \{ \begin{array}{ll} \displaystyle2\sqrt{\dfrac{1}{N_{p}} \sum_{i=1}^{N_{p}} \| \boldsymbol x_{p,i} - \boldsymbol x_c \|^2}, & N_{p} > 2 , \\ \sqrt{1.6}D_p, & N_{p} = 2 , \\ D_p, & N_{p} = 1 . \end{array}\right. \end{equation}

Here, $\boldsymbol x_{p,i}$ denotes the position of the centre of primary particle $i$, and the floc's centre of mass is evaluated as $\boldsymbol x_c = \sum _{i=1}^{N_{p}} \boldsymbol x_{p,i} / N_{p}$. While the characteristic diameter $D_f$ more closely represents the true spatial extent of the floc, the gyration diameter $D_g$ also accounts for the irregularity of the floc shape.

Following Khelifa & Hill (Reference Khelifa and Hill2006a,Reference Khelifa and Hillb), we then calculate the fractal dimension $n_f$ of the floc

(4.5)\begin{equation} n_f = \frac{\log N_{p}}{\log \dfrac{D_f}{D_p}} , \end{equation}

as a measure of its compactness. A dense, nearly spherical floc has $n_f \approx 3$, while for a linear floc $n_f \approx 1$. When $N_{p} = 1$ and $D_f/D_p = 1$, the above definition of the fractal dimension does not yield a finite value, and we set $n_f = 1$. In this way, the definition of the fractal dimension is continuous between $N_p=1$ and $N_p=2$. It is important to note that this differs from previous studies, which usually set $n_f = 3$ for this case (Khelifa & Hill Reference Khelifa and Hill2006a,Reference Khelifa and Hillb; Maggi et al. Reference Maggi, Mietta and Winterwerp2007; Son & Hsu Reference Son and Hsu2009).

For a typical floc with $N_p = 7$ that maintains its identity, figure 6 shows the evolution of $D_f$ and $n_f$ over time. 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, which modifies $D_f$ and $n_f$ in the opposite directions.

Figure 6. 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 case Flo10 with governing parameters $Co = 1.2 \times 10^{-7}$, $St = 0.1$, $G = 0.91$.

Figure 7 shows the evolution with time of the various floc size measures, for cases Flo6–9 in table 3 with different turbulent shear rates $G$. The other parameters are kept approximately constant at $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $\rho _s = 2.65$, and $2.24 \leqslant \eta /D_p \leqslant 2.28$. As can be seen from figure 7(a), a smaller shear rate results in a longer transient phase before the average number of primary particles per floc $\bar N_{p} = N / N_f$ reaches an equilibrium. A smaller value of $G$ furthermore gives rise to an equilibrium stage characterized by fewer flocs with more primary particles, since the weaker hydrodynamic stresses cannot break up the flocs as easily. Figures 7(a) and 7(b) indicate that both the average characteristic diameter $\bar D_f$ and the average gyration diameter $\bar D_g$ increase for smaller $G$. This is consistent with previous observations by other authors in both laboratory experiments (He et al. Reference He, Nan, Li and Li2012; Guérin et al. Reference Guérin, Coufort-Saudejaud, Liné and Frances2017) and river estuaries (Manning & Dyer Reference Manning and Dyer2002; Manning, Langston & Jonas Reference Manning, Langston and Jonas2010). Both $\bar D_g$ and $\bar D_f$ remain smaller than the Kolmogorov length scale $0.0179 \leqslant \eta \leqslant 0.0183$ for all cases. Since flocs with one or two primary particles have a constant fractal dimension $n_f = 1$, we evaluate the average fractal dimension $\bar n_{f,lar}$ from only those flocs with three or more particles. Figure 7(d) shows that $\bar n_{f,lar}$ increases for smaller shear rates, which demonstrates that for weaker turbulence the floc shape tends to be more compact. This finding is consistent with experimental observations by He et al. (Reference He, Nan, Li and Li2012), whereas previous numerical work by Chen et al. (Reference Chen, Li and Marshall2019) reports a constant value $\bar n_{f,lar} = 1.64$.

Figure 7. Temporal evolution of various floc size measures for different turbulent shear rates $G$, with $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $\rho _s = 2.65$ and $2.24 \leqslant \eta /D_p \leqslant 2.28$ (cases Flo6–9). (a) Average number of primary particles per floc $\bar N_p$. (b) Average characteristic floc diameter $\bar D_f$. (c) Average floc gyration diameter $\bar D_g$. (d) Average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles. Larger turbulent shear results in smaller flocs, with fewer primary particles and more elongated shapes.

Figure 8 discusses the floc growth during the very early flow stages, as a function of the turbulent shear rate $G$. As seen in figure 8(a), the evolution of $\bar D_f(t)$ can be closely approximated by an exponential function of the form

(4.6)\begin{equation} \bar D_f = b_1 (e^{b_2 t} - 1) + D_p , \end{equation}

where $b_1$ and $b_2$ represent fitting coefficients. Based on corresponding fits for different values of $G$, figure 8(b) displays the time-dependent floc growth rate $\textrm {d}\bar D_f /\textrm {d}t$ for different $G$. Consistent with the experimental observations by He et al. (Reference He, Nan, Li and Li2012), we find stronger shear to cause more rapid flocculation for $t < 5$. After this early stage the trends reverse, which reflects the fact that the equilibrium stage is reached faster for stronger turbulence. This agrees with the experimental findings by Braithwaite et al. (Reference Braithwaite, Bowers, Nimmo Smith and Graham2012), who also reported the equilibrium stage to emerge more quickly for stronger turbulence, due to more frequent floc collisions. We remark that the evolution of $\bar N_{p}$ (not shown) exhibits corresponding trends.

Figure 8. Early-stage flocculation rate for different turbulent shear rates $G$, with $Co = 1.2 \times 10^{-7}$, $St = 0.06$, $\rho _s = 2.65$ and $\eta /D_p \approx 2.26$ (cases Flo6–9). (a) The early-stage simulation results for $\bar D_f(t)$ can be accurately fitted by an exponential relation, as shown for the representative case Flo6 with $G = 1.49$. (b) The flocculation rate $\textrm {d}(\bar D_f)/\textrm {d}t$ obtained from the exponential fits of $\bar D_f(t)$. Initially flocs grow fastest in strong turbulence. Subsequently their growth rate decays, as the equilibrium stage is reached more rapidly for strong turbulence.

Figure 9 presents corresponding floc size results for different Stokes numbers, obtained from cases Flo10–13 in table 3. These simulations all employ the same turbulent flow Tur6, so that they have constant parameter values $Co = 1.2 \times 10^{-7}$, $G = 0.91$ and $\eta /D_p = 1.85$. The value of $St$ is varied by changing the density ratio $\rho _s$. Figures 9(a) and 9(b) indicate that the equilibrium values of both $\bar N_{p}$ and $\bar D_f$ increase for smaller $St$. This reflects the fact that cohesive forces become more dominant for smaller $St$, due to the lower drag force and the shorter particle response time. By again employing exponential fits for the early stages, we obtain the floc growth rate $\textrm {d}\bar D_f / \textrm {d}t$ for different $St$-values, as shown in figure 9(c). Initially flocs with intermediate Stokes numbers of $O(1)$ are seen to grow most rapidly, consistent with our earlier findings for two-dimensional cellular flows (Zhao et al. Reference Zhao, Vowinckel, Hsu, Köllner, Bai and Meiburg2020). This trend changes for $t > 20$, due to the later onset of the equilibrium stage for small Stokes numbers. The time evolution of the average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles is shown in figure 9(d). It demonstrates that smaller Stokes numbers result in more compact flocs.

Figure 9. Temporal evolution of various floc size measures, for different Stokes number values $St$, with $Co = 1.2 \times 10^{-7}$, $G = 0.91$ and $\eta /D_p = 1.85$ (cases Flo10–13). (a) Average number of primary particles per floc $\bar N_p$; (b) average characteristic floc diameter $\bar D_f$; (c) early-stage flocculation rate $\textrm {d}(\bar D_f)/\textrm {d}t$ obtained from exponential fits of $\bar D_f(t)$; and (d) average fractal dimension $\bar n_{f,lar}$ of flocs with three or more primary particles. During the equilibrium stage, the number of primary particles per floc, the characteristic floc diameter and the fractal dimension all increase for smaller Stokes numbers. Initially, flocs with $St \approx O(1)$ exhibit the fastest growth.

Figure 10 analyses the influence of the cohesive number $Co$ by comparing cases Flo1–5 in table 3. The other parameters are held constant at $St = 0.02$, $G = 0.29$, $\rho _s = 2.65$ and $\eta /D_p = 3.30$. We note that due to the small values of $St$ and $G$, the emergence of an equilibrium stage takes longer in these simulations. In fact, for case Flo5 with $Co = 1.2 \times 10^{-7}$, an equilibrium had not yet formed by $t=20\,000$, when the simulation terminated. Nevertheless, the simulations demonstrate the tendency of higher $Co$ to result in larger values of $\bar N_{p}$ during all phases of the flow, cf. figure 10(a). Interestingly, however, we observe that during the transient flow stages the flocs for $Co = 6 \times 10^{-8}$ have larger average diameters $\bar D_f$ and $\bar D_g$ than those for $Co = 1.2 \times 10^{-7}$, even though they contain fewer primary particles, cf. figures 10(b) and 10(c). The explanation for this finding is given by figure 10(d), which indicates that for $Co = 1.2 \times 10^{-7}$ the flocs have a higher average fractal dimension $\bar n_f$ and are more compact than those for $Co = 6 \times 10^{-8}$, which can be deformed more easily by turbulent stresses.

Figure 10. Temporal evolution of various floc size measures for different values of the cohesive number $Co$, with $St = 0.02$, $G = 0.29$, $\rho _s = 2.65$ and $\eta /D_p = 3.30$ (cases Flo1–5). (a) Average number of primary particles per floc $\bar N_p$; (b) average characteristic floc diameter $\bar D_f$; (c) average floc gyration diameter $\bar D_g$; and (d) average fractal dimension of flocs $\bar n_f$. Note that case Flo5 with $Co = 1.2 \times 10^{-7}$ has not yet reached the equilibrium stage by the end of the simulation. For higher $Co$-values, the equilibrium stage is characterized by larger flocs with more primary particles. During the transient stages, however, intermediate $Co$-values can give rise to flocs that are more elongated and hence larger than those at higher $Co$-values, in spite of having fewer primary particles.

In summary, as a general trend we observe that during the equilibrium stages weaker turbulence, lower Stokes numbers and higher cohesive numbers result in larger and more compact flocs.

4.4. Floc size distribution during the equilibrium stage

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 11(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, cf. figure 7. This finding is consistent with the experimental observations by Braithwaite et al. (Reference Braithwaite, Bowers, Nimmo Smith and Graham2012) in an energetic tidal channel. Corresponding results for different $St$-values display a similar trend (not shown).

Figure 11. 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$ (cases Flo6–9). (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$ (cases Flo1–4).

Figure 11(b) shows the size distributions for different values of the cohesive number. For larger values of $Co$, we find that the peak of the distribution decreases and shifts to larger flocs, while the exponential decay rate with increasing floc size is reduced.

4.5. Change in floc microstructure

In the following, we analyse the deformation in time of those flocs that maintain their identity over the time interval $\Delta T$, by keeping track of their characteristic diameter $D_f$. Accordingly, we distinguish between those flocs within the fraction $\theta _{id}$ whose value of $D_f$ increases or stays constant during $\Delta T$, $\theta _{id,gro}$, and those whose diameter decreases, $\theta _{id,shr}$

(4.7)\begin{equation} \theta_{id} = \theta_{id,gro} + \theta_{id,shr}. \end{equation}

Equations (4.2) and (4.7) thus yield

(4.8)\begin{equation} \theta_{br} + \theta_{id,gro} + \theta_{id,shr} + \theta_{ad} = 1 . \end{equation}

For the choice of $\Delta T = 3$, figure 12(a) displays the evolution of these fractions for the representative case Flo6. Interestingly, we find that $\theta _{id,gro}$ is consistently much larger than $\theta _{id,shr}$, which indicates that of those flocs who maintain their identity during $\Delta T$, many more see their value of $D_f$ increase than decrease. Hence, it is much more common for these flocs to deform from a compact shape to an elongated one than vice versa. This consistent difference between $\theta _{id,gro}$ and $\theta _{id,shr}$ can be maintained only if the elongated flocs eventually break. As a general trend, turbulent stresses thus stretch cohesive flocs before eventually breaking them. This confirms earlier numerical results by Nguyen et al. (Reference Nguyen, Rasmuson, Thalberg and Bjo2014) and Gunkelmann et al. (Reference Gunkelmann, Ringl and Urbassek2016), who employed conceptually simpler models with ‘sticky’ cohesive particles and observed that compact flocs have greater strength than elongated ones.

Figure 12. Evolution of the floc number fractions displaying different behaviors. (a) Of those flocs that maintain their identity during $\Delta T$, many more are being stretched than shrink, resulting in $\theta _{id,gro} \gg \theta _{id,shr}$ (case Flo6 with $G = 1.49$). (b) The fraction $\theta _{id,gro}$ that is being stretched increases for more intense turbulence. (c) The fraction $\theta _{id,shr}$ that shrinks decreases for stronger turbulence. For (b,c) the colour coding of the curves is identical, and the other parameter values are $Co = 1.2 \times 10^{-7}$ and $St = 0.06$ (cases Flo6–9).

The influence of the shear rate $G$ on the fractions $\theta _{id,gro}$ and $\theta _{id,shr}$ during equilibrium is displayed in figures 12(b) and 12(c), respectively. For larger values of $G$, the fraction $\theta _{id,gro}$ grows, while $\theta _{id,shr}$ is reduced, which reflects the fact that more intense turbulence tends to elongate the cohesive flocs more strongly. Figures 13(a) and 13(b) indicate that larger $St$-values also promote the stretching of those flocs that maintain their integrity, as they increase $\theta _{id,gro}$ and reduce $\theta _{id,shr}$. Figures 14(a) and 14(b) show that smaller $Co$-values result in the elongation of those flocs that maintain their identity, whereas stronger cohesive forces prompt the flocs to assume a more compact shape.

Figure 13. Evolution of floc number fractions for different values of $St$. (a) Of those flocs that maintain their identity during $\Delta T$, the fraction $\theta _{id,gro}$ that is stretched increases with $St$. (b) The fraction $\theta _{id,shr}$ whose diameter $D_f$ decreases is reduced for larger $St$. The other parameter values are $Co = 1.2 \times 10^{-7}$ and $G = 0.91$ (cases Flo10–13).

Figure 14. Evolution of floc number fractions for different values of $Co$. (a) Of those flocs that maintain their identity during $\Delta T$, the fraction $\theta _{id,gro}$ that is stretched increases for weaker cohesive forces. (b) The fraction $\theta _{id,shr}$ whose diameter $D_f$ decreases is reduced for weaker cohesive forces. The other parameter values are $St = 0.02$ and $G = 0.29$ (cases Flo1–5).

4.5.1. Orientation of elongated flocs

We now investigate the alignment of the elongated flocs with the principal strain directions of the turbulent velocity field. Towards this end, we define an Eulerian fluid velocity difference tensor $\boldsymbol A$ for each floc at time $t$ as

(4.9)\begin{equation} {\boldsymbol A}(m,n) = \frac{\boldsymbol u_{f,c}(n) - \boldsymbol u_{f,j}(n)}{\boldsymbol x_{c}(m) - \boldsymbol x_{p,j}(m)} , \end{equation}

where $m,n = 1, 2, 3$ represent the $x$-, $y$- and $z$-components, respectively, of the tensor and vectors. Also, ${\boldsymbol x}_c = (x_{c}, y_{c}, z_{c})^\textrm {T}$ denotes the location of the floc's centre of mass, and the fluid velocity averaged over the volume of the floc is written as $\boldsymbol {u}_{f,c} = \sum _{1}^{N_p}(\boldsymbol {u}_{f,i}) / N_p$. The location and fluid velocity at the centre of the primary particle $j$ that is located the farthest away from the floc's centre of mass are denoted as ${\boldsymbol x}_{p,j} = (x_{p,j}, y_{p,j}, z_{p,j})^\textrm {T}$ and $\boldsymbol {u}_{f,j} = (u_{f,j}, v_{f,j} w_{f,j})^\textrm {T}$. The orientation of the floc is defined as ${\boldsymbol x}_{f} = {\boldsymbol x}_{p,j} - {\boldsymbol x}_{c}$. Especially for large flocs, ${\boldsymbol x}_c$ and ${\boldsymbol x}_{p,j}$ can be multiple grid spacings apart from each other. We remark that $\boldsymbol A$ is defined by sampling the velocity difference at points separated along a line, and it thus represents a simplified approach for considering the influence of the fluid velocity gradients on the whole floc, compared with employing the full coarse-grained velocity gradient tensor (Pumir, Bodenschatz & Xu Reference Pumir, Bodenschatz and Xu2013). Hence, $\boldsymbol A$ differs from the standard, locally evaluated fluid velocity gradient tensor (Ashurst et al. Reference Ashurst, Kerstein, Kerr and Gibson1987; Pumir & Wilkinson Reference Pumir and Wilkinson2011; Voth & Soldati Reference Voth and Soldati2017).

We decompose this Eulerian velocity difference tensor $\boldsymbol A = \boldsymbol S + \boldsymbol Q$ into the symmetric velocity difference tensor $\boldsymbol S = {\boldsymbol S}^\textrm {T}$, which is similar but not identical to the strain rate tensor, and the anti-symmetric tensor $\boldsymbol Q = -{\boldsymbol Q}^\textrm {T}$. The three eigenvalues $r_m$ of the velocity difference tensor $\boldsymbol S$ are ordered as $r_1 > r_2 > r_3$. We remark that the intermediate eigenvalue $r_2$ is automatically zero, by nature of the definition of $\boldsymbol S$. With the three eigenvalues we associate three corresponding orthonormal eigenvectors ${\boldsymbol e}_m$

(4.10)\begin{equation} \boldsymbol S {\boldsymbol e}_m = r_m {\boldsymbol e}_m . \end{equation}

We define a modified vorticity vector $\boldsymbol \omega = \omega {\boldsymbol e}_{\omega }$ based on the anti-symmetric tensor $\boldsymbol Q$, with magnitude $\omega$ and unit direction vector ${\boldsymbol e}_{\omega }$ (Pumir & Wilkinson Reference Pumir and Wilkinson2011).

We furthermore define a modified deformation gradient tensor $\boldsymbol B$ that characterizes the Lagrangian deformation experienced by a fluid element extending from the floc's centre of mass to its primary particle $j$, over the time interval from $t$ to $(t + \Delta t)$, as

(4.11)\begin{equation} {\boldsymbol B}(m,n) = \frac{\boldsymbol x_{c}(m) - \boldsymbol x_{p,j}(m)}{[\boldsymbol x_{c}(n)+ \Delta t \boldsymbol u_{f,c}(n)] - [\boldsymbol x_{p,j}(n) + \Delta t \boldsymbol u_{f,j}(n)]} . \end{equation}

This modified deformation gradient tensor $\boldsymbol B$ provides a Lagrangian description of the fluid stretching (Parsa et al. Reference Parsa, Guasto, Kishore, Ouellette, Gollub and Voth2011; Ni, Ouellette & Voth Reference Ni, Ouellette and Voth2014). It differs from the standard locally evaluated deformation gradient tensor, for the same reasons mentioned earlier for the Eulerian velocity difference tensor $\boldsymbol A$.

The Lagrangian stretching tensor $\boldsymbol C = \boldsymbol B \boldsymbol B^T$, obtained from the two symmetric inner products of $\boldsymbol B$ with itself, is similar but not identical to the left Cauchy–Green strain tensor commonly used to define stretching in a Lagrangian basis (Chadwick Reference Chadwick2012). The three eigenvalues of the Lagrangian stretching tensor $\boldsymbol C$ are ordered as $r_{L1} > r_{L2} > r_{L3}$, and the three corresponding orthonormal eigenvectors are ${\boldsymbol e}_{Lm}$

(4.12)\begin{equation} \boldsymbol C {\boldsymbol e}_{Lm} = r_{Lm} {\boldsymbol e}_{Lm} . \end{equation}

In the following, we investigate the alignment of ${\boldsymbol x}_{f}$ and ${\boldsymbol e}_{\omega }$ with ${\boldsymbol e}_{m}$ and ${\boldsymbol e}_{Lm}$, respectively.

We focus on those elongated flocs with $n_f \leqslant 1.2$ and $N_{p} \geqslant 2$, and firstly analyse their alignment with the eigendirections ${\boldsymbol e}_m$ of the Eulerian velocity difference tensor and the vorticity vector ${\boldsymbol e}_{\omega }$ in terms of the magnitude of the angle $\alpha$ between them. We divide the elongated flocs into three different groups, according to the ratio of their characteristic diameter $D_f$ and the Kolmogorov length scale $\eta$. The alignment of small flocs with $D_f/\eta < 0.8$ and medium-size flocs with $0.8 \leqslant D_f/\eta \leqslant 1.2$ is indicated in figures 15(a) and 15(b), respectively. The alignment of large flocs with $D_f/\eta > 1.2$ is not shown. The results indicate that the modified vorticity vector ${\boldsymbol e}_{\omega }$ is always aligned with the intermediate eigenvector ${\boldsymbol e}_2$, which is consistent with the previous finding by Ashurst et al. (Reference Ashurst, Kerstein, Kerr and Gibson1987). We observe that medium-size flocs are strongly aligned with the intermediate eigenvector ${\boldsymbol e}_2$ and the vorticity vector ${\boldsymbol e}_{\omega }$, as shown in figure 15(b). This result is consistent with previous findings for microscopic axisymmetric rod-like particles in turbulence by Pumir & Wilkinson (Reference Pumir and Wilkinson2011), who noticed that the vortex stretching term ${\boldsymbol A} {\boldsymbol \omega }$ promotes, and the viscous term $\nabla ^2 {\boldsymbol \omega } / Re$ opposes, the alignment of ${\boldsymbol x_f}$ with ${\boldsymbol e}_{\omega }$. In contrast, figure 15(a) shows that small flocs tend to align themselves with the extensional strain direction ${\boldsymbol e}_1$. For large flocs, we did not observe preferential alignment of the flocs with any of the three eigendirections of the Eulerian velocity difference tensor (not shown).

Figure 15. Floc alignment with the principal directions of the symmetric Eulerian velocity difference tensor for the representative case Flo9. Results include both the flocculation and the equilibrium stages, for all elongated flocs with $n_f \leqslant 1.2$ and $N_{p,local} \geqslant 2$. The upper two frames show the alignment of the floc orientation ${\boldsymbol x}_f$ with the eigendirections ${\boldsymbol e}_m$ of the symmetric Eulerian velocity difference tensor, and with the vorticity vector ${\boldsymbol e}_{\omega }$: (a) small flocs with $D_f / \eta < 0.8$ and (b) medium-size flocs with $0.8 \leqslant D_f / \eta \leqslant 1.2$. Small flocs are preferentially aligned with the extensional strain direction, while medium-size flocs tend to align themselves with the intermediate strain direction. The lower two frames show the alignment with the eigendirections ${\boldsymbol e}_{Lm}$ of the Lagrangian deformation tensor: (c) the floc orientation ${\boldsymbol x}_f$ and (d) the vorticity vector ${\boldsymbol e}_{\omega }$. Both the flocs and the vorticity vector tend to be aligned with the strongest Lagrangian stretching direction.

The alignment of the elongated flocs and the modified vorticity vector with the eigendirections ${\boldsymbol e}_{Lm}$ of the Lagrangian stretching tensor is shown in figures 15(c) and 15(d), respectively. The results indicate that the elongated flocs are perfectly aligned, and the modified vorticity vector is strongly aligned with the direction corresponding to the largest eigenvalue ${\boldsymbol e}_{L1}$ of the Lagrangian stretching tensor $\boldsymbol C$. This alignment is consistent with, but even more pronounced than the corresponding previous findings by Parsa et al. (Reference Parsa, Guasto, Kishore, Ouellette, Gollub and Voth2011) and Ni et al. (Reference Ni, Ouellette and Voth2014), due to our definition of the modified deformation gradient tensor $\boldsymbol B$. The perfect alignment of ${\boldsymbol x}_{f}$ with ${\boldsymbol e}_{L1}$ suggests that the present Lagrangian stretching tensor $\boldsymbol C$ is well suited for analysing the instantaneous alignment of flocs in turbulent flows.

4.6. Floc size vs Kolmogorov length scale

Several authors have hypothesized that for sufficiently strong turbulence the median floc size should be of the same order as the smallest turbulent eddies (McCave Reference McCave1984; Fettweis et al. Reference Fettweis, Francken, Pison and Van den Eynde2006; Coufort et al. Reference Coufort, Dumas, Bouyer and Liné2008; Kuprenas et al. Reference Kuprenas, Tran and Strom2018). Others have suggested that even the largest flocs cannot exceed the Kolmogorov length scale (Verney et al. Reference Verney, Lafite, Brun-Cottan and Le Hir2011). In the following, we discuss data from the present simulations in order to explore this issue.

Figure 16 discusses case Flo9, with $\eta /D_p = 2.25$, $G = 0.62$, $St = 0.06$ and $Co = 1.2 \times 10^{-7}$. Figure 16(a) compares both the average and the maximum floc size to the Kolmogorov scale. It demonstrates that for all times the average floc diameter $\bar D_f$ is smaller than the Kolmogorov length scale $\eta$. However, at any given time the largest floc diameter $D_{f,max}$ is several times larger than $\eta$. We now define ‘big’ flocs as those whose diameter $D_f$ is larger than $\eta$, and we indicate their fraction as

(4.13)\begin{equation} \theta_{big} = \frac{N_{f,big}}{N_f} , \end{equation}

where $N_{f,big}$ is the number of big flocs at a given moment. Analogous to equation (4.8), we also define the fractions of big flocs that grow, break or maintain their identity, so that we have

(4.14)\begin{equation} \theta_{big,br} + \theta_{big,id,gro} + \theta_{big,id,shr} + \theta_{big,ad} = \theta_{big} . \end{equation}

Here the subscripts $br$, $ad$, $id$, $gro$ and $shr$ have the same meanings as in (4.8). Figure 16(b) demonstrates that $\theta _{big}$ plateaus around a value of 0.2, so that at any given time approximately 20 % of all flocs are larger than the Kolmogorov scale; $\theta _{big,id,shr}$ levels off around 0.1, which indicates that a substantial fraction of these big flocs deform towards a more compact shape while maintaining their identity over $\Delta T = 3$. Figure 16(c) shows that the ratio $\theta _{big,br}/\theta _{br}$ is stable around 0.6, so that about 60 % of those flocs that break are larger than the Kolmogorov scale $\eta$. The ratio $\theta _{big,id,gro}/\theta _{id,gro}$ levels off around 0.2, meaning that of those flocs that become elongated while maintaining their identity, only about 20 % are big. Hence we can conclude that most of the big flocs tend to either become more compact or to break, but that some continue to grow. This finding is consistent with previous experimental observations by Stricot et al. (Reference Stricot, Filali, Lesage, Spérandio and Cabassud2010), who found that the breakage of big flocs is not instantaneous and depends on the floc strength.

Figure 16. Constraint on the floc size by the Kolmogorov length scale, for case Flo9 with $\eta /D_p = 2.25$, $G = 0.62$, $St = 0.06$ and $Co = 1.2 \times 10^{-7}$. (a) Temporal evolution of the average and maximum floc diameters, $\bar D_f$ and $D_{f,max}$. The dashed horizontal line indicates the Kolmogorov length scale $\eta$. (b) The fraction $\theta _{big}$ of flocs that are larger than $\eta$, and the fraction $\theta _{big,id,shr}$ of big flocs maintaining their identity that become more compact. (c) The ratios $\theta _{big,br}/\theta _{br}$, $\theta _{big,id,gro}/\theta _{id,gro}$ and $\theta _{big,ad}/\theta _{ad}$. (d) Average time interval $\Delta t_{big,gro}$ over which big flocs exhibit continuous growth.

Figure 16(d) addresses the time scale over which big flocs grow. The duration of the continuous growth of the big flocs is denoted by $\Delta t_{big,gro}$. We remark that $\Delta t_{big,gro}$ is measured for all big flocs until their $D_f$ is smaller than $\eta$. The results indicate that, on average, flocs larger than the Kolmogorov scale keep growing only for the relatively short time period of $\Delta t_{big,gro} \approx 4.8$. This is consistent with previous observations for controls on floc growth in tidal cycle experiments by Braithwaite et al. (Reference Braithwaite, Bowers, Nimmo Smith and Graham2012), who found that big flocs cannot resist the turbulent stresses for long, and that they are torn apart quickly. This relatively quick breakage of large flocs in the simulations also agrees with our findings in § 4.5, which showed that flocs are being continually stretched until they break.

To summarize, while the size of an individual floc can be larger than the Kolmogorov length for a brief period of time, once $D_f$ becomes bigger than $\eta$, the floc tends to break relatively soon. Given that the physical parameter ranges listed in table 1 represent common fluid–particle systems in nature, our simulation data suggest that the average floc size $\bar D_f$ is effectively limited by the Kolmogorov length scale $\eta$ in such systems. We remark, however, that for other classes of primary particles with potentially much stronger bonds it may be possible, in principle, to form flocs that are significantly larger than the Kolmogorov scale.

For cases Flo14 and Flo15, figure 17 discusses corresponding results regarding the time scale over which big flocs grow. Flo14 employs an increased shear rate $G = 2.7$ along with $\eta /D_p = 1.08$, while Flo15 has $G = 7.4$ and $\eta /D_p = 0.65$. We remark that the ratio $\eta /D_p$ is widely used to classify the primary particles as either ‘small’ if $\eta /D_p > 1$, or as ‘finite size’ if $\eta /D_p \leqslant 1$ (Fiabane et al. Reference Fiabane, Zimmermann, Volk, Pinton and Bourgoin2012; Chouippe & Uhlmann Reference Chouippe and Uhlmann2015; Costa et al. Reference Costa, Boersma, Westerweel and Breugem2015). Hence Flo14 addresses the small particle scenario, while Flo15 considers finite-size particles. Interestingly, figure 17(a) shows that the time interval $\Delta t_{big,gro}$ over which big flocs grow for Flo14 is smaller than the corresponding value for Flo9 in figure 16(d). This observation indicates that the constraint on floc growth by the turbulent eddies becomes stronger for increasing shear rate $G$, which is consistent with experimental findings for small particles by Braithwaite et al. (Reference Braithwaite, Bowers, Nimmo Smith and Graham2012). Those authors had found that the time lag before big flocs break becomes shorter for larger $G$. However, a further increase of the shear rate to $G = 7.4$ in case Flo15, which means that the primary particles now fall into the finite-size category, yields a longer time lag $\Delta t_{big,gro} \approx 20.5$, as shown in figure 17(b). While the detailed reasons for this observation will require further investigation, we can conclude that the enhanced control on floc growth by the Kolmogorov length scale for stronger turbulent shear is seen to hold for small primary particles with $\eta /D_p > 1$, although it does not necessarily apply for finite-size primary particles with $\eta /D_p \leqslant 1$.

Figure 17. Average time interval $\Delta t_{big,gro}$ over which big flocs exhibit continuous growth: (a) $\eta /D_p = 1.08$, $Co = 1.2 \times 10^{-7}$, $St = 0.38$ and $G = 2.7$ (case Flo14); (b) $\eta /D_p = 0.65$, $Co = 1.2 \times 10^{-7}$, $St = 1.25$ and $G = 7.4$ (case Flo15).