Flocculation of suspended cohesive particles in homogeneous isotropic turbulence

Abstract We investigate the dynamics of cohesive particles in homogeneous isotropic turbulence, based on one-way coupled simulations that include Stokes drag, lubrication, cohesive and direct contact forces. We observe a transient flocculation phase, followed by a statistically steady equilibrium phase. We analyse the temporal evolution of floc size and shape due to aggregation, breakage and deformation. Larger turbulent shear and weaker cohesive forces yield smaller elongated flocs. Flocculation proceeds most rapidly when the fluid and particle time scales are balanced and a suitably defined Stokes number is $O(1)$. During the transient stage, cohesive forces of intermediate strength produce flocs of the largest size, as they are strong enough to cause aggregation, but not so strong as to pull the floc into a compact shape. Small Stokes numbers and weak turbulence delay the onset of the equilibrium stage. During equilibrium, stronger cohesive forces yield flocs of larger size. The equilibrium floc size distribution exhibits a preferred size that depends on the cohesive number. We observe that flocs are generally elongated by turbulent stresses before breakage. Flocs of size close to the Kolmogorov length scale preferentially align themselves with the intermediate strain direction and the vorticity vector. Flocs of smaller size tend to align themselves with the extensional strain direction. More generally, flocs are aligned with the strongest Lagrangian stretching direction. The Kolmogorov scale is seen to limit floc growth. We propose a new flocculation model with a variable fractal dimension that predicts the temporal evolution of the floc size and shape.


Introduction
Individual cohesive particles suspended in liquid or gaseous fluid flows tend to form larger aggregates, due to attractive inter-particle forces that cause the primary particles to flocculate.This mechanism plays a dominant role in environmental processes such as sediment erosion and transport in rivers and oceans, or soil erosion by wind (Winterwerp 2002;Guo & He 2011;Wang et al. 2013;Tarpley et al. 2019).In planetary astrophysics, corresponding processes influence the coagulation of dust during the formation of protoplanetary disks (Ormel, Spaans & Tielens 2007;Schäfer, Speith & Kley 2007;Ormel et al. 2009Ormel et al. , 2011)).The emergence of large aggregates due to the flocculation of cohesive primary particles is also highly relevant in the context of a wide range of industrial processes, such as the ingestion of dust in gas turbine engines (Bons, Prenter & Whitaker 2017;Sacco et al. 2018), or the use of membrane separation technologies for wastewater treatment and the production of potable water (Bratskaya et al. 2006;Leiknes 2009;Moghaddam, Moghaddam & Arami 2010;Kang et al. 2012).Similarly, the operation of certain types of medical equipment, for example dry powder inhalers (Yang, Wu & Adams 2013a, 2015;Tong et al. 2013Tong et al. , 2016)), involves the formation of agglomerates or flocs.The flocculation process is strongly affected by the turbulent nature of the underlying fluid flow.Small-scale eddies modify the collision dynamics of the primary particles and hence the growth rate of the flocs, while turbulent stresses can result in the deformation and breakup of larger cohesive flocs.Hence the dynamic equilibrium between floc growth and breakup is governed by a complex and delicate balance of hydrodynamic and inter-particle forces.
A host of experimental studies have provided considerable insight into key aspects of the development of flocs in turbulent shear flows, such as their growth rate (Biggs & Lant 2000;Yu et al. 2006;Xiao et al. 2010;Kuprenas, Tran & Strom 2018), the equilibrium size distribution (Chaignon et al. 2002;Bouyer, Liné & Do-Quang 2004;Rahmani, Dabros & Masliyah 2004;Lee, Hyeong & Cho 2020) and the transient shape of the flocs (Maggi, Mietta & Winterwerp 2007;He et al. 2012;Guérin et al. 2017).Based on the early pioneering work by Levich (1962), several of these investigations have employed a population balance approach to formulate models for the temporal floc evolution (Winterwerp 1998;Maggi et al. 2007;Son & Hsu 2008, 2009;Shin, Son & Lee 2015).Alternative approaches based on the classical work by Smoluchowski (1918) propose statistical collision equations (Ives & Bhole 1973;Yang et al. 2013b;Klassen 2017).Most of the above approaches do not incorporate detailed information on the overall floc strength, which varies with the floc size and shape, and with the strength of the bonds between the primary cohesive particles (Dizaji, Marshall & Grant 2019).Moreno-Atanasio & Ghadiri (2006), on the other hand, consider the dependence of the overall floc strength on the number and strength of the bonds within the floc.Nguyen et al. (2014) and Gunkelmann, Ringl & Urbassek (2016) observe that loosely structured agglomerates fragment more easily during collisions than densely packed ones.
In recent years, highly resolved numerical simulations have begun to provide a promising new avenue for gaining insight into the interplay of hydrodynamic, inertial and inter-particle forces during the growth, deformation and breakup of aggregates (Marshall & Li 2014).The study by Zhao et al. (2020) focuses on a conceptually simple cellular model flow in order to explore the competition between inertial, drag and cohesive forces during the flocculation process.The authors find that floc growth proceeds most rapidly if the fluid and particle time scales are in equilibrium, so that a suitably defined Stokes number is of order unity.Based on simulations in a similar model flow, Ruan, Chen & Li (2020) suggest a criterion for the breakup of aggregates.Dizaji et al. (2019) investigate the dynamics, collision and fragmentation of flocs in shear flows, via two-way The paper is structured along the following lines.Section 2 briefly reviews the governing equations for the fluid flow and the particle motion, and it describes the computational approach.It identifies the governing dimensionless parameters and quantifies the range over which they will be varied in the present investigation.The properties of the turbulent flow fields are described in § 3, and their statistically stationary and isotropic nature is discussed.Starting from 10 000 randomly distributed individual particles, we then analyse the temporal evolution of the floc size and shape as a result of aggregation, deformation and breakage in § 4. Here, we will distinguish between the transient flocculation stage and the equilibrium stage, and we will discuss the underlying physical mechanisms.We will furthermore analyse the alignment of the flocs with regard to the principal strain directions of the turbulent velocity field, and we will focus on how the Kolmogorov scale affects the maximum floc size.Subsequently, we introduce the new flocculation model in § 5, and we compare its predictions to those obtained from existing models.Section 6 summarizes the main findings of the current investigation, and presents its key conclusions.

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 ρ f and kinematic viscosity ν is governed by where u f = (u f , v f , w f ) T denotes the fluid velocity vector and p indicates the hydrodynamic pressure.We employ the spectral approach of Eswaran & Pope (1988) to obtain the forcing term F tur , which generates and maintains statistically stationary turbulence, as implemented in Chouippe & Uhlmann (2015).Here, F tur is non-zero only in the low-wavenumber band where the wavenumber vector |κ| < κ f , with κ f = 2.3κ 0 and κ 0 = 2π/L 0 , with L 0 denoting the length of the physical domain.The origin κ = 0 is not forced.In addition to the cutoff wavenumber κ f , the random forcing process is governed by the dimensionless parameter D s = σ 2 T 0 L 4 0 /ν 3 , where σ 2 and T 0 indicate the variance and the time scale of the random process, respectively.Regarding the details of evaluating F tur from κ f and D s , we refer the reader to the original work by Eswaran & Pope (1988).
We approximate each primary suspended particle i as a sphere moving with translational velocity u p,i = (u p,i , v p,i , w p,i ) T and angular velocity ω p,i .These are obtained from the linear and angular momentum equations where the primary particle i moves in response to the Stokes drag force F d,i = −3πD p μ(u p,i − u f ,i ), and the particle-particle interaction force 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 μ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 2009;Biegert et al. 2017;Chen et al. 2019;Vowinckel et al. 2019).Here, u p,i indicates the particle velocity evaluated at the particle centre, u f ,i = N i 1 (φ i,k 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, u f ,k is the fluid velocity at the centre of the grid cell k and φ 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. 2019).Also, m p denotes the particle mass, μ 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 ρ p .The parameter F c,i accounts for the direct contact force F con,ij in both the normal and tangential directions, as well as for short-range normal and tangential forces due to lubrication F lub,ij and cohesion F coh,ij , where the subscript ij indicates the interaction between particles i and j.Also, I p = πρ p D 5 p /60 denotes the moment of inertia of the particle and T c,i represents the torque due to particle-particle interactions, where we distinguish between direct contact torque T con,ij and lubrication torque T lub,ij .Within a large floc, we account for all of the individual binary particle interactions.
The lubrication force F lub,ij is accounted for based on Cox & Brenner (1967) as implemented in Zhao et al. (2020).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. (2017), the collision force 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. (2017).The cohesive force 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. (2019), 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. (2020).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.

Non-dimensionalization
In order to render the above governing equations dimensionless, we consider primary particles with diameter D p = 5 μ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 × 10 −4 m.As time scale of the random turbulent forcing process we select T 0 = 7.81 × 10 −5 s.By choosing L 0 , T 0 and ρ f = 1000 kg m −3 as the characteristic length, time and density scales, we obtain the characteristic velocity scale U 0 = L 0 /T 0 = 8 m s −1 , which is similar to values employed in previous investigations (Chen et al. 2019;Chen & Li 2020).We employ L 0 and U 0 to define the turbulence Reynolds number Re = L 0 U 0 ρ f /μ.
The dimensionless continuity and momentum conservation equations can then be expressed as while the dimensionless equations of motion for the primary cohesive particles take the form Flocculation in homogeneous isotropic turbulence ρs = ρ p /ρ f .The dimensionless direct contact and lubrication forces, F con,ij and F lub,ij , are accounted for based on Zhao et al. (2020), while the dimensionless cohesive force F coh,ij is defined as (2.9) Here, ζmin = 0.0015 Dp and hco = 0.05 Dp represent the surface roughness of the particles and the range of the cohesive force, respectively.Also, n represents the outward-pointing normal on the particle surface, while ζn,ij is the normal surface distance between particles i and j.The cohesive number Co indicates the ratio of the maximum cohesive force F coh,ij at ζn,ij = hco /2 to the characteristic inertial force where the Hamaker constant A H is a function of the particle and fluid properties, and the characteristic distance ζ 0 = 0.00025D p .Vowinckel et al. (2019) 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. 2020).
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 Dp , the total number of particles N, the density ratio ρ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.

Computational set-up
The triply periodic computational domain Ω has a dimensionless size of L x × L y × L z = 1 × 1 × 1, with the number of grid cells N x × N y × N z = 128 × 128 × 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 2018), and it is in line with the earlier study of Chen et al. (2019).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 η = 1/(Re 3 ) 1/4 and u rms = (2k/3) 1/2 , respectively, where and k denote the domain-averaged dissipation rate and kinetic energy of the fluctuations.The Taylor Reynolds number Re λ = λu rms Re of the turbulence is based on the Taylor microscale λ = √ 15 u rms /(Re ) 1/2 .To provide a more complete quantitative description of the fluid shear, we define the vorticity fluctuation amplitude which can also be regarded as the turbulent shear rate.For additional details with regard to these quantities, we refer the reader to Pope (2001).

Turbulence properties for different Re λ
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 η.Since the particle diameter needs to be larger than the grid spacing, and the number of grid points is limited, suitable values of η require a relatively low Re λ .On the other hand, it is known that for Re λ ≤ O(50) the turbulence may not be fully developed and isotropic (Mansour & Wray 1994).Hence this section presents a more detailed discussion of the turbulence properties for Re λ ≤ O(50).
Figure 1 shows the time-dependent evolution of the box-averaged Kolmogorov length η, the root-mean-square velocity u rms , the Taylor Reynolds number Re λ 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 η/h = 0.6656, Chouippe & Uhlmann (2015) 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 3 shows the temporal evolution of the domain-averaged magnitude of the velocity components 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.We define the instantaneous kinetic energy components in Fourier space, E 11 (κ), E 22 (κ) and E 33 (κ), as  κη κη The spectra confirm that the statistically stationary flow fields are approximately isotropic.
where κ = |κ| denotes the wavenumber.Figure 4 shows the time-averaged one-dimensional energy spectra.Only the wavenumbers below the cutoff wavenumber (κ 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 (2015, p. 10) for higher values of Re λ ≈ 60.We conclude that the present forcing scheme yields statistically steady flow fields that are approximately isotropic for the current range of Re λ -values.Flocculation in homogeneous isotropic turbulence 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 (2006) 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.

Flocculation of cohesive particles
We adopt a multiscale time-stepping approach in which the fluid motion is calculated with a time step t based on the criterion that the Courant-Friedrichs-Lewy number CFL ≤ 0.5.The particle motion, on the other hand, is evaluated with a much smaller time step t p = t/15.Since the computational approach maintains a contact duration of T c = 10 t = 150 t p (Biegert et al. 2017), 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 where the Kolmogorov Reynolds number Re η = ηu rms Re.Since the particle diameter D p is constant throughout the present investigation, St depends on the density ratio ρ 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. (2019), 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 ρ s and the size ratio η/D p are implicitly accounted for by St and G.

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 × 10 −7 , St = 0.06, G = 0.62, ρ s = 2.65 and η/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. 2020).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.
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 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 θ id = N f ,id /N f , θ ad = N f ,ad /N f , and θ br = N f ,br /N f .It follows that We found that a value of 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 θ id , θ ad and θ 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 ≈ const., we observe that θ ad / = θ 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.

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 as well as its gyration diameter D g (Chen et al. 2019), (4.4) Here, x p,i denotes the position of the centre of primary particle i, and the floc's centre of mass is evaluated as x c = N p i=1 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 (2006a,b), we then calculate the fractal dimension n f of the floc as a measure of its compactness.A dense, nearly spherical floc has n f ≈ 3, while for a linear floc n f ≈ 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 2006a,b;Maggi et al. 2007;Son & Hsu 2009).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 t 210, hydrodynamic forces deform the floc so that it becomes more compact, which reduces D f and increases n f .Later on, near t ≈ 240, the floc is being stretched, which modifies D f and n f in the opposite directions.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 × 10 −7 , St = 0.06, ρ s = 2.65, and 2.24 η/D p 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 Np = 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 Df and the average gyration diameter Dg increase for smaller G.This is consistent with previous observations by other authors in both laboratory experiments (He et al. 2012;Guérin et al. 2017) and river estuaries (Manning & Dyer 2002;Manning, Langston & Jonas 2010).Both Dg and Df remain smaller than the Kolmogorov length scale 0.0179 η 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 nf ,lar from only those flocs with three or more particles.Figure 7(d) shows that nf ,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. (2012), whereas previous numerical work by Chen et al. (2019) reports a constant value nf ,lar = 1.64., 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. (2012), 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 Np (not shown) exhibits corresponding trends.
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 × 10 −7 , G = 0.91 and η/D p = 1.85.The value of St is varied by changing the density ratio ρ s .Figures 9(a 1) are seen to grow most rapidly, consistent with our earlier findings for two-dimensional cellular flows (Zhao et al. 2020).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 nf ,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 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, ρ s = 2.65 and η/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 × 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 Np 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 × 10 −8 have larger average diameters Df and Dg than those for Co = 1.2 × 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 × 10 −7 the flocs have a higher average fractal dimension nf and are more compact than those for Co = 6 × 10 −8 , which can be deformed more easily by turbulent stresses.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.

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 (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 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.

Change in floc microstructure
In the following, we analyse the deformation in time of those flocs that maintain their identity over the time interval T, by keeping track of their characteristic diameter D f .Accordingly, we distinguish between those flocs within the fraction θ id whose value of D f increases or stays constant during T, θ id,gro , and those whose diameter decreases, θ id,shr θ id = θ id,gro + θ id,shr . (4.7) Equations ( 4.2) and (4.7) thus yield For the choice of T = 3, figure 12(a) displays the evolution of these fractions for the representative case Flo6.Interestingly, we find that θ id,gro is consistently much larger than θ id,shr , which indicates that of those flocs who maintain their identity during 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 θ id,gro and θ 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. (2014) and Gunkelmann et al. (2016), who employed conceptually simpler models with 'sticky' cohesive particles and observed that compact flocs have greater strength than elongated ones.
The influence of the shear rate G on the fractions θ id,gro and θ id,shr during equilibrium is displayed in figures 12(b) and 12(c), respectively.For larger values of G, the fraction θ id,gro grows, while θ 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 θ id,gro and reduce θ 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.

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 A for each floc at time t as where m, n = 1, 2, 3 represent the x-, y-and z-components, respectively, of the tensor and vectors.Also, x c = (x c , y c , z c ) 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 u f ,c = N p 1 (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 x p,j = (x p,j , y p,j , z p,j ) T and u f ,j = (u f ,j , v f ,j w f ,j ) T .The orientation of the floc is defined as Especially for large flocs, x c and x p,j can be multiple grid spacings apart from each other.We remark that 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 2013).Hence, A differs from the standard, locally evaluated fluid velocity gradient tensor (Ashurst et al. 1987;Pumir & Wilkinson 2011;Voth & Soldati 2017).
We decompose this Eulerian velocity difference tensor A = S + Q into the symmetric velocity difference tensor S = S T , which is similar but not identical to the strain rate tensor, and the anti-symmetric tensor Q = −Q T .The three eigenvalues r m of the velocity difference tensor 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 S. With the three eigenvalues we associate three corresponding orthonormal eigenvectors e m Se m = r m e m . (4.10) We define a modified vorticity vector ω = ωe ω based on the anti-symmetric tensor Q, with magnitude ω and unit direction vector e ω (Pumir & Wilkinson 2011).
We furthermore define a modified deformation gradient tensor 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 + t), as This modified deformation gradient tensor B provides a Lagrangian description of the fluid stretching (Parsa et al. 2011;Ni, Ouellette & Voth 2014).It differs from the standard locally evaluated deformation gradient tensor, for the same reasons mentioned earlier for the Eulerian velocity difference tensor A.
The Lagrangian stretching tensor C = BB T , obtained from the two symmetric inner products of 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 2012).The three eigenvalues of the Lagrangian stretching tensor C are ordered as r L1 > r L2 > r L3 , and the three corresponding orthonormal eigenvectors are e Lm Ce Lm = r Lm e Lm . (4.12) In the following, we investigate the alignment of x f and e ω with e m and e Lm , respectively.We focus on those elongated flocs with n f 1.2 and N p 2, and firstly analyse their alignment with the eigendirections e m of the Eulerian velocity difference tensor and the 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 1.2 and N p,local 2. The upper two frames show the alignment of the floc orientation x f with the eigendirections e m of the symmetric Eulerian velocity difference tensor, and with the vorticity vector e ω : (a) small flocs with D f /η < 0.8 and (b) medium-size flocs with 0.8 D f /η 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 e Lm of the Lagrangian deformation tensor: (c) the floc orientation x f and (d) the vorticity vector e ω .Both the flocs and the vorticity vector tend to be aligned with the strongest Lagrangian stretching direction.
vorticity vector e ω in terms of the magnitude of the angle α 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 η.The alignment of small flocs with D f /η < 0.8 and medium-size flocs with 0.8 D f /η 1.2 is indicated in figures 15(a) and 15(b), respectively.The alignment of large flocs with D f /η > 1.2 is not shown.The results indicate that the modified vorticity vector e ω is always aligned with the intermediate eigenvector e 2 , which is consistent with the previous finding by Ashurst et al. (1987).We observe that medium-size flocs are strongly aligned with the intermediate eigenvector e 2 and the vorticity vector e ω , as shown in figure 15(b).This result is consistent with previous findings for microscopic axisymmetric rod-like particles in turbulence by Pumir & Wilkinson (2011), who noticed that the vortex stretching term Aω promotes, and the viscous term ∇ 2 ω/Re opposes, the alignment of x f with e ω .In contrast, figure 15(a) shows that small flocs tend to align themselves with the extensional strain direction 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).
The alignment of the elongated flocs and the modified vorticity vector with the eigendirections 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 Flocculation in homogeneous isotropic turbulence modified vorticity vector is strongly aligned with the direction corresponding to the largest eigenvalue e L1 of the Lagrangian stretching tensor C.This alignment is consistent with, but even more pronounced than the corresponding previous findings by Parsa et al. (2011) and Ni et al. (2014), due to our definition of the modified deformation gradient tensor B. The perfect alignment of x f with e L1 suggests that the present Lagrangian stretching tensor 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 1984;Fettweis et al. 2006;Coufort et al. 2008;Kuprenas et al. 2018).Others have suggested that even the largest flocs cannot exceed the Kolmogorov length scale (Verney et al. 2011).In the following, we discuss data from the present simulations in order to explore this issue.
Figure 16 discusses case Flo9, with η/D p = 2.25, G = 0.62, St = 0.06 and Co = 1.2 × 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 Df is smaller than the Kolmogorov length scale η.However, at any given time the largest floc diameter D f ,max is several times larger than η.We now define 'big' flocs as those whose diameter D f is larger than η, and we indicate their fraction as 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 θ big,br + θ big,id,gro + θ big,id,shr + θ big,ad = θ big .(4.14) Here the subscripts br, ad, id, gro and shr have the same meanings as in (4.8).
Figure 16(b) demonstrates that θ 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; θ 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 T = 3. Figure 16(c) shows that the ratio θ big,br /θ br is stable around 0.6, so that about 60 % of those flocs that break are larger than the Kolmogorov scale η.The ratio θ big,id,gro /θ 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. (2010), who found that the breakage of big flocs is not instantaneous and depends on the floc strength.
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 t big,gro .We remark that t big,gro is measured for all big flocs until their D f is smaller than η.The results indicate that, on average, flocs larger than the Kolmogorov scale keep growing only for the relatively short time period of t big,gro ≈ 4.8.This is consistent with previous observations for controls on floc growth in tidal cycle experiments by Braithwaite et al. (2012), 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 η, 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 Df is effectively limited by the Kolmogorov length scale η 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 η/D p = 1.08, while Flo15 has G = 7.4 and η/D p = 0.65.We remark that the ratio η/D p is widely used to classify the primary particles as either 'small' if η/D p > 1, or as 'finite size' if η/D p 1 (Fiabane et al. 2012;Chouippe & Uhlmann 2015;Costa et al. 2015).Hence Flo14 addresses the small particle scenario, while Flo15 considers finite-size particles.Interestingly, figure 17(a) shows that the time interval 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. (2012).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 t big,gro ≈ 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 η/D p > 1, although it does not necessarily apply for finite-size primary particles with η/D p 1.

A new flocculation model with variable fractal dimension
As indicated by figures 7-10, the average characteristic floc diameter Df and the average fractal dimension nf both increase during the flocculation stage, and then remain constant during the equilibrium stage.This indicates that flocs of larger size generally have a more compact shape, and that it is difficult for elongated flocs to keep growing in turbulent shear without breaking.Closer inspection indicates that for all of the cases listed in table 3 the relationship between these two quantities can be approximated well by a power law of the form (5.1) The condition that nf = 1 for an individual primary particle requires that k 1 = 1, while the value of k 2 varies as a function of St, Co and G. Typical fitting results are shown in figure 18(a) for cases Flo4 and Flo5.This power law relationship allows us to obtain the average fractal dimension nf during flocculation as a function of the average floc diameter Df , rather than assuming a constant fractal dimension, as was done in earlier investigations (Winterwerp 1998;Kuprenas et al. 2018;Zhao et al. 2020).The power law (5.1) is closely related to the earlier study by Khelifa & Hill (2006a,b).However, those authors assumed that an individual primary particle has n f = 3, and consequently they set k 1 = 3.For the exponent k 2 they proposed an empirical correlation of the form where Df ,char denotes the characteristic floc size that exhibits the characteristic fractal dimension nf ,char .As a general rule, Df ,char and nf ,char should be evaluated from experiments before one can then determine k 2 from (5.2).Should that not be feasible, Khelifa & Hill (2006a,b) suggested assuming constant values of nf ,char = 2 and Df ,char = 2000 μm, which yields a constant value for k 2 that depends only on the primary particle size D p .As figure 18(a) indicates, however, k 2 should be a function of G, St and Co even for a constant D p , since nf ,char = 2 is associated with different average floc sizes Df ,char /D p in cases Flo4 and Flo5.Hence, even though (5.2) has been widely used to describe the fractal dimension of flocs (Maggi et al. 2007;Son & Hsu 2009;Klassen 2017), we will now try to refine this scaling law by accounting for the dependence of k 2 on St, Co and G.By fitting the simulation results for all of the cases Flo1-15, we obtain a relationship for k 2 of the form k 2 = 0.44St −0.018 Co 0.096 G −1.5 , (5.3) with a coefficient of determination value R-squared of 0.97.We remark that in a laboratory experiment or field investigation it may be challenging to evaluate the Stokes number St as defined in (4.1), if the root-mean-square velocity u rms is unknown.To overcome this difficulty, we follow the approach taken in our earlier work (Zhao et al. 2020), where we defined the characteristic Stokes number St char and cohesive number Co char by employing the characteristic fluid velocity u char = 0.25(G/Re) . (5.7) For the specific range of physical parameters listed in table 1, a 3 = 1 yields optimal agreement with a maximum deviation of ±30 % from the simulation data.As we will see below, this value of a 3 is not universally optimal, so that a 3 will have to be recalibrated for other parameter ranges.In the following, we will compare predictions for the fractal dimension by the new relation (5.7) with corresponding ones by the earlier relation of Khelifa & Hill ((5.1)-(5.2)).
Employing the approach of Maggi & Winterwerp (2004), Maggi et al. (2007) estimate the time evolution of the average fractal dimension nf and floc size Df in experiments with constant turbulent shear rates G = 5, 10, 20 and 40 s −1 , respectively.The suspended cohesive sediment in the experiments has a primary particle diameter D p = 5 μm, density ρ p = 2650 kg m −3 , and concentration c = 0.5 g L −1 .Since the authors assume nf = 3 for flocs with one particle while we set nf = 1 for that situation, we have to convert their original experimental data before we can compare them with the present simulation results.The details of the conversion are discussed in Appendix A, and the converted data are presented in figure 18(b).In addition, we set the dimensional Kolmogorov length η = [μ/(ρ f G)] 0.5 m and the Hamaker constant A H = 1.0 × 10 −20 J to obtain the characteristic values St char and Co char according to (5.4)-(5.5).Since the experimental shear rates G = 5 ∼ 40 s −1 are much smaller than the simulation values G = 3.7 × 10 3 ∼ 9.5 × 10 4 s −1 , we have to recalibrate the constant a 3 required for our model (5.7) from the experimental data.Based on the fact that the exponent k 2 should decrease for increasing G, we obtain a 3 = 4 × 10 −6 for the minimum experimental shear rate G = 5s −1 , and a 3 = 4 × 10 −5 for the maximum experimental shear rate G = 40 s −1 , respectively.Figure 18 In order to develop a variable fractal dimension model for the transient stages, we build on the approach taken in our recent investigation (Zhao et al. 2020).There we conducted cohesive sediment simulations for a steady, two-dimensional cellular flow model.Based on the simulation data, we proposed an analytical flocculation model of the form Df = ( Np ) 1/n f D p , (5.8a)  in the value of a 1 .In order to assess this robustness, we ran the present model for a 1 = 2 and a 1 = 32, instead of the optimal value a 1 = 8 that we had obtained earlier from the calibration.The results, shown in figure 19 as dashed and dotted red lines, indicate that the model predictions are reasonably robust with regard to uncertainties in the value of a 1 .
To summarize, our new fractal relation (5.7) no longer has the limitation associated with assuming a constant value for k 2 in (5.1) and (5.2) when predicting the variable fractal dimension nf .In addition, we observe that predictions of the floc size Df and the fractal dimension nf as functions of time by the present flocculation model (5.7) and (5.8) are fairly robust with respect to uncertainties that arise when calibrating the empirical coefficients by means of experimental data.

Conclusions
In the present investigation we have employed one-way coupled simulations to explore the dynamics of cohesive particles in homogeneous isotropic turbulence.The simulations account for the Stokes drag, as well as lubrication, cohesive and direct contact forces.They demonstrate the existence of a transient flocculation phase which is characterized by the growth of the average floc size.This flocculation phase is followed by a statistically steady equilibrium phase governed by a balance between floc growth and breakup.The simulations provide information about the temporal evolution of the floc size and shape, as a result of aggregation, breakage and deformation, and as a function of the governing parameters.In general, we find that larger turbulent shear and weaker cohesive forces limit the floc size and result in elongated floc shapes.Flocculation proceeds most rapidly during the transient stage when the Stokes number of the primary particles based on the Kolmogorov scales is of order unity.During the transient stage cohesive forces of intermediate strength yield the largest flocs.On one hand, these intermediate cohesive forces are strong enough to result in the rapid aggregation of primary particles, but on the other hand they are not so strong as to pull them into a compact shape.During the equilibrium stage, stronger cohesive forces produce larger flocs.Small Stokes numbers and weak turbulence typically lead to a later onset of the equilibrium stage.The equilibrium floc size distribution exhibits a preferred size as function of the cohesive number.This distribution decays exponentially for larger floc sizes.The simulation results indicate that flocs are generally elongated by turbulent stress before they eventually break.We observe that flocs close to the Kolmogorov scale in size preferentially align themselves with the intermediate strain direction and the vorticity vector.Flocs that are smaller than the Kolmogorov scale, on the other hand, tend to align themselves with the direction of extensional strain.The simulation results furthermore demonstrate that flocs generally align themselves with the strongest Lagrangian stretching direction.The simulations show that the average floc size is effectively limited by the Kolmogorov scale, and can at most exceed it marginally.However, individual flocs can grow larger than the Kolmogorov scale for a limited amount of time.Based on the simulation data we propose a novel flocculation model that allows for a variable fractal dimension, which enables us to predict the temporal evolution of the floc size and shape, as a function of the governing dimensionless parameters, after some limited calibration.Predictions by the new model are fairly robust and cover a broad range of parameters.
Declaration of interests.The authors report no conflict of interest.

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

Figure 3 .
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.

Figure 4 .
Figure 4. Time-averaged one-dimensional energy spectra.The vertical dashed lines indicate the respective cutoff wavenumber of the turbulence forcing scheme, κ f η = 2.3(2π/L x )η.(a) Case Tur1 and (b) case Tur8.The spectra confirm that the statistically stationary flow fields are approximately isotropic.

Figure 5 .
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 (θ id ), add primary particles (θ ad ) or undergo breakage (θ br ) over the time interval T = 3.All results are for case Flo9 with Co = 1.2 × 10 −7 , St = 0.06, G = 0.62, ρ s = 2.65, η/D p = 2.25.

Figure 6 .
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 × 10 −7 , St = 0.1, G = 0.91.

Figure 7 .
Figure 7. Temporal evolution of various floc size measures for different turbulent shear rates G, with Co = 1.2 × 10 −7 , St = 0.06, ρ s = 2.65 and 2.24 η/D p 2.28 (cases Flo6-9).(a) Average number of primary particles per floc Np .(b) Average characteristic floc diameter Df .(c) Average floc gyration diameter Dg .(d) Average fractal dimension nf ,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 Figure 8 .
Figure7shows 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 × 10 −7 , St = 0.06, ρ s = 2.65, and 2.24 η/D p 2.28.As can be seen from figure7(a), a smaller shear rate results in a longer transient phase before the average number of primary particles per floc Np = 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.Figures7(a) and 7(b) indicate that both the average characteristic diameter Df and the average gyration diameter Dg increase for smaller G.This is consistent with previous observations by other authors in both laboratory experiments(He et al. 2012;Guérin et al. 2017) and river estuaries(Manning & Dyer 2002;Manning, Langston & Jonas 2010).Both Dg and Df remain smaller than the Kolmogorov length scale 0.0179 η 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 nf ,lar from only those flocs with three or more particles.Figure7(d)shows that nf ,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 byHe et al. (2012), whereas previous numerical work byChen et al. (2019) reports a constant value nf ,lar = 1.64.Figure8discusses the floc growth during the very early flow stages, as a function of the turbulent shear rate G.As seen in figure8(a), the evolution of Df (t) can be closely 921 A17-15 ) and 9(b) indicate that the equilibrium values of both Np and Df 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 d Df /dt for different St-values, as shown in figure 9(c).Initially flocs with intermediate Stokes numbers of O(

921Figure 9 .
Figure 9. Temporal evolution of various floc size measures, for different Stokes number values St, with Co = 1.2 × 10 −7 , G = 0.91 and η/D p = 1.85 (cases Flo10-13).(a) Average number of primary particles per floc Np ; (b) average characteristic floc diameter Df ; (c) early-stage flocculation rate d( Df )/dt obtained from exponential fits of Df (t); and (d) average fractal dimension nf ,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 ≈ O(1) exhibit the fastest growth.

921Figure 10 .
Figure 10.Temporal evolution of various floc size measures for different values of the cohesive number Co, with St = 0.02, G = 0.29, ρ s = 2.65 and η/D p = 3.30 (cases Flo1-5).(a) Average number of primary particles per floc Np ; (b) average characteristic floc diameter Df ; (c) average floc gyration diameter Dg ; and (d) average fractal dimension of flocs nf .Note that case Flo5 with Co = 1.2 × 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.

921Figure 12 .Figure 13 .
Figure 12.Evolution of the floc number fractions displaying different behaviors.(a) Of those flocs that maintain their identity during T, many more are being stretched than shrink, resulting in θ id,gro θ id,shr (case Flo6 with G = 1.49).(b) The fraction θ id,gro that is being stretched increases for more intense turbulence.(c) The fraction θ 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 × 10 −7 and St = 0.06 (cases Flo6-9).

921Figure 14 .
Figure 14.Evolution of floc number fractions for different values of Co.(a) Of those flocs that maintain their identity during T, the fraction θ id,gro that is stretched increases for weaker cohesive forces.(b) The fraction θ 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).

921Figure 16 .
Figure 16.Constraint on the floc size by the Kolmogorov length scale, for case Flo9 with η/D p = 2.25, G = 0.62, St = 0.06 and Co = 1.2 × 10 −7 .(a) Temporal evolution of the average and maximum floc diameters, Df and D f ,max .The dashed horizontal line indicates the Kolmogorov length scale η.(b) The fraction θ big of flocs that are larger than η, and the fraction θ big,id,shr of big flocs maintaining their identity that become more compact.(c) The ratios θ big,br /θ br , θ big,id,gro /θ id,gro and θ big,ad /θ ad .(d) Average time interval t big,gro over which big flocs exhibit continuous growth.
(b) demonstrates that the present relation successfully reproduces the range of experimental data for different G-values, whereas Khelifa & Hill's relation does not account for variations in G.At the same time, we do need to keep in mind that the present model does require a recalibration of a 3 for different experimental parameter ranges.

Figure 19 .
Figure 19.Comparisons between the numerical data and predictions by the present model and the combined model listed in table 4, simulation data of case Flo4 with Co = 6.0 × 10 −8 , St = 0.02 and G = 0.29 is selected.(a) Calibration predictions for the temporal evolution of the average floc sizes Df , the calibrated coefficients in the models are a 2 = 0.5 and a 3 = 1, the constant k 1 = 1.(b) Comparisons for the temporal evolution of the average fractal dimension nf .

Table 2 .
Physical parameters of the single-phase turbulence simulations.As input parameters we specify the fluid Reynolds number Re = L 0 U 0 /ν and the characteristic parameter of the random turbulent forcing process D s = σ 2 T 0 L 4 0 Re 3 .The simulation then yields the Taylor Reynolds number Re λ = λu rms Re, the Kolmogorov scale η, the average root-mean-square velocity u rms and the shear rate G = 1/(Re η 2 ).All of these output quantities are obtained by averaging over space and time, after a statistically stationary state has evolved.

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 ρ s and η/D p are implicitly accounted for by St and G.