| Home | E-Submission | Sitemap | Contact Us |  
Environ Eng Res > Volume 19(2); 2014 > Article
Lee and Molz: Numerical Simulation of Turbulence-Induced Flocculation and Sedimentation in a Flocculant-Aided Sediment Retention Pond


A model combining multi-dimensional discretized population balance equations with a computational fluid dynamics simulation (CFD-DPBE model) was developed and applied to simulate turbulent flocculation and sedimentation processes in sediment retention basins. Computation fluid dynamics and the discretized population balance equations were solved to generate steady state flow field data and simulate flocculation and sedimentation processes in a sequential manner. Up-to-date numerical algorithms, such as operator splitting and LeVeque flux-corrected upwind schemes, were applied to cope with the computational demands caused by complexity and nonlinearity of the population balance equations and the instability caused by advection-dominated transport. In a modeling and simulation study with a two-dimensional simplified pond system, applicability of the CFD-DPBE model was demonstrated by tracking mass balances and floc size evolutions and by examining particle/floc size and solid concentration distributions. Thus, the CFD-DPBE model may be used as a valuable simulation tool for natural and engineered flocculation and sedimentation systems as well as for flocculant-aided sediment retention ponds.

1. Introduction

In recent years, various best management practices (BMPs) related to the control of sediments during storm events have been developed. Among these BMPs, several suggest that the removal of clay and other colloidal-sized particles by retention or detention ponds may be enhanced by the addition of flocculating agents. A few operators are now experimenting with the addition of such agents to sediment inflow, which can greatly improve the retention properties of ponds in some cases. Contemporary literature and sediment pond operators support the conclusion that flocculant-aided sediment retention ponds are going to be increasingly important in future as a means of minimizing the detrimental effects of erosion and nonpoint-source water pollution [13]. To date, use has been driven more by practicing engineers than by academics. However, the operation of such ponds is complex, involving turbulent flow of variable intensity, different pond geometries, particle growth due to flocculation, sedimentation of particle size classes at different rates, and various schemes for time-dependent flocculant additions. Most existing pond systems are not designed in a consistent manner based on fundamental principles. For example, many designs are based simply on an ad hoc rule, such as a set pond volume per hectare of drained area [4]. Hence, the entire field would benefit from a better understanding of flocculation and sedimentation processes and the availability of a realistic, physically based model for designing and optimizing the automated operation of sediment retention ponds. Therefore, there is a need for a realistic theory describing flocculation and non-homogeneous turbulent sedimentation in retention ponds, a practical method for solving the rather complex governing equations, and performance of the required small- and large-scale experiments necessary to characterize the parameters and functions that the theory contains. This paper primarily discusses the mathematical formulation and computation underlying flocculation and sedimentation processes in flocculant-aided sediment retention ponds.
One of the most realistic ways to simulate flocculation and non-homogeneous turbulent sedimentation in retention ponds is by applying population balance equations (PBEs) within a computational fluid dynamics (CFD) framework for solving the Navier–Stokes equations (mass and momentum conservation equations). PBEs have been used to simulate particle/floc aggregation phenomena for many scientific and engineering applications. Most PBEs derive from the Smoluchowski equation describing a simple particle/floc aggregation process. They are now generalized by incorporating various additional processes, such as particle/floc breakage models, shaping and growth strategies, chemical interaction models, etc. [516]. Application of PBEs, ranging from fundamental scientific research to advanced engineering applications, has become more practical as computational speed and capacity have increased. However, such applications are still at the forefront of engineering research.
In multi-dimensional simulation, such as simulation of retention pond dynamics, conventional PBEs continue to be computationally demanding even with modern computer technologies. For example, in conventional PBEs, particles/flocs aggregate in a sequential manner from singlet, to doublet, then to triplet, and so on. Thus, conventional PBEs require thousands to millions of particle/floc size classes and associated differential transport-reaction equations to simulate particle/floc growth from nano- or micro-sized constituent monomers to milli-sized aggregates that settle due to gravity. To overcome this computational difficulty in multi-dimensional simulations, the discretized PBE (DPBE) and quadrature method of moments (QMOM) have been proposed.
In the QMOM approach, moments of the particle/floc size distribution, instead of number concentrations of particles/flocs, are used as dependent variables in differential transport-reaction equations to reduce the computational overloads that occur in multi-dimensional applications. The lower-order moments then yield key monitoring indices, such as particle/floc sizes and solid concentrations, indirectly through the use of product-difference algorithms [10, 12, 14, 17, 18]. Thus, QMOM provides computational advantages but imposes difficulties for scientists or engineers attempting to understand results because of the more abstract formulation and resulting algorithms.
In the DPBE methodology, particle/floc number concentrations can be tracked as dependent variables in differential transport-reaction equations, similar to the method of conventional PBEs. However, the DPBE formulation differs from that of conventional PBEs because particles/flocs of DBPE are assumed to double their size from singlet to doublet and then to quadruplet, etc., in the flocculation process. Thus, with only dozens of defined particle/floc size classes, particles/flocs can grow to sizes susceptible to gravity-induced settling, which are thousands to millions of times larger than the size of monomers [7, 12, 19, 20]. In contrast to QMOM, DPBE directly tracks key indices, such as particle/floc sizes and solid concentrations, by simply integrating differential equations without additional data processing steps. Thus, with respect to clearness of results, the DPBE approach may be more intuitive and advantageous than QMOM. Therefore, in this research, the discretized particle transport-reaction model combined with a fluid dynamics model (CFD-DPBE model) was set up, and its applicability was tested in a model pond system. The mathematical formulation and application strategy of the CFD-DPBE model were studied in a two-dimension computational domain representing the vertical and flow-parallel cross section of a flocculant-aided sediment retention pond.

2. Material and Methods

2.1. Background and Mathematical Models

The CFD-DPBE model consists of 1) CFD software to obtain the Reynolds-averaged turbulent flow field and 2) multi-dimensional DPBE software, containing particle/floc aggregation and break-up kinetics to simulate transport, flocculation, and sedimentation within the pre-obtained flow field.

2.1.1. Computational fluid dynamics

The Reynolds-averaged continuity and Navier–Stokes (RANS) equations, containing a two-equation k – ɛ turbulence model, were solved using FLOW-3D software (Flow Science Inc., Santa Fe, NM, USA) to simulate turbulent fluid motion within a retention pond. In the CFD-DPBE model, particles/flocs are assumed to travel via fluid motion and to aggregate or disintegrate due to impact and shear forces or effects [14, 17, 21].
The RANS equations consist of mass and momentum conservation equations in the differential form given by Eqs. (1) and (2), respectively (we now use the summation convention in writing the 3D equations, wherein sums over the three spatial coordinates are understood when an index is repeated).
In Eqs. (1) and (2), i and j are indices, xi represents coordinate directions (i = 1, 2, 3 for x, y, z directions, respectively), 〈Ui〉 is the time-averaged velocity component (m/s), t represents time (s), ρ is the fluid density (kg/m3), p is the piezometric pressure (kg/m/s2), and ν is the kinematic viscosity of the fluid (m2/s). A symmetric second-order tensor 〈ui uj〉 representing Reynolds normal or shear stresses (m2/s2) is modeled with Eq. (3), and νT, the turbulent viscosity (N·s/m2), is specified by Eq. (4). In Eqs. (3) and (4), δij is Kroneker’s delta, Cμ is 0.09, a model constant, k represents turbulent kinetic energy (m2/s2), and ɛ is the turbulent energy dissipation rate (m2/s3). These two energy terms are obtained from Eqs. (5) and (6), which together represent the k–ɛ turbulence model. They are solved simultaneously with the RANS equations. The velocity gradient ( G=ɛ/v, /s), which is obtained from the two-equation k – ɛ turbulence model, causes particle/floc aggregation or break-up kinetics in the DPBE and thus serves as a coupling term between the turbulent flow field (CDF problem) and the DPBE [14]. In Eqs. (5) and (6), model fitting constants have been found to be σk = 1.0, σɛ = 1.3, Cɛ1 = 1.44, and Cɛ2 = 1.92 [17]. P, specified by Eq. (7), is the turbulent kinetic energy production term (m2/s3).

2.1.2 Multi-dimensional discretized population balance equations

With a given flow field obtained from CFD software, the multi-dimensional DPBE is used to simulate particle/floc transport and flocculation in the ponds. Following Prat and Ducoste [14], a generic mathematical model for the DPBE may be written as
In Eq. (8)ni = n(x, y, z, Di, t) is number concentration of flocs (/m3) of linear class size Di (m) (i=1, 2, …imax; D1DiDmax; for all Di, ni is called the population density function); x, y, z, and t are position and time; 〈Ux, 〈Uy〉, and 〈Uz mean fluid velocity components in the x, y, and z directions (m/s); ρ is fluid density (kg/m3); k = k(x, y, z, t) is turbulent kinetic energy (m2/s2); ɛ = ɛ(x, y, z, t) is turbulent energy dissipation rate (m2/s3); Cμ = 0.09 is standard value of a CFD model constant (see Eq. (4)); and ugi is settlement velocity (m/s) of the ith floc class due to gravity. On the left-hand side of Eq. (8), the respective terms in brackets represent the storage change (I), particle/floc mean advection (II), and the turbulent dispersion of the particle/floc (III), while on the right-hand side, the source/sink terms (IV) represent the net effects of aggregation, break-up, and settling due to gravity [14]. The coefficients or functions in these later terms are largely empirical and must be determined by experiment. The quantities depending on the turbulent fluid variables (〈Ux〉, 〈Uy〉, 〈Uz〉, k, and ɛ) couple the DPBE equations (Eq. (8)) to the CFD equations (Eqs. (1)(7)). However, as currently formulated, the CDF equations are solved independent of the DPBE.
To track particle/floc fates with the DPBE, both Eulerian and Lagrangian tracking methods are applicable. However, in this research, the Eulerian method rather than the Lagrangian method, which tracks individual particles or flocs, was applied to observe the distribution of scalars within the entire computational domain (Eulerian [12, 14, 18, 22, 23] and Lagrangian [24]).
To obtain the particle/floc settling velocity in Eq. (8), Stokes equation was used in the context of fractal theory, which represents the structural characteristics of aggregating particles/flocs. Even though many complex and elaborate particle/floc settling equations have been developed, including those involving interaction or drag coefficients with ambient flow, the standard Stokes equation was applied as a prototype in this research [25, 26]. Fractal theory describes particle/floc packing or growth structure with constituent monomers in which particle/floc size follows a power-law function with respect to the number of monomers in a given floc size (Eq. (10)) [2732]. Stokes equation combined with fractal theory is given by Eq. (9) [16, 3335]. In Eqs. (9) and (10), Di represents floc diameter of size class i (m), D0 is monomer diameter (m), Df is fractal dimension, k is lacunarity (generally set as 1, which implies no lacunarity), ρs is particle density (kg/m3), ρw is fluid density (kg/m3), g is gravitational acceleration (9.81 m/s2), and η is fluid viscosity (kg/m/s). In Eq. (10), 2i-1 represents the number of monomers forming an ith particle/floc by following the discretized size classification strategy of the DPBE, which will be described in the following section.

2.1.3 Kinetics of particle/floc aggregation and breakage

The core part of a multi-dimensional DPBE (Eq. (8)) is the sink and source terms, which characterize the aggregation and break-up kinetics ((agg/break)i). These terms are written as a series of differential equations in Eq. (11). The particle/floc number concentration in a certain discrete size range (ni) is used as a dependent variable of a partial differential equation. Following the discretization scheme of the DPBE, each mean particle size class contains two times the number of constituent monomers in the previous smaller class. Thus if δ is the beginning (irreducible) particle size, class 1 will contain particles of size δ, class 2 will contain particles of size 2δ, class 3 will contain sizes 3δ and 4δ, class 4 will contain 5δ through 8δ, class 5 will contain 9δ through 16δ, and so on. Because the maximum particle size in class i increases as 2(i-2), 30 classes will contain particles sizes varying from δ to 228δ, which represents a growth factor of more than 268 million. Ignoring transport and settling for notational convenience, the partial differential equations with discrete particle/floc size sets may be written as
In Eq. (11), the processes indicated by the various Roman numerals are (I) i-sized particle/floc generation by collision with other smaller particle/floc classes, (II) generation by collision within the i-1 class, (III) disappearance by collision with smaller classes, (IV) disappearance by collision with equal or larger classes, (V) disappearance by fragmentation of the i class, and (VI) generation by fragmentation of larger classes. These mechanisms are illustrated in Fig. 1.
Several empirical or theoretical factors or functions (α, β, a, and b) are incorporated into the aggregation and break-up kinetics. The collision efficiency factor (α) represents the physicochemical properties of solid and liquid to cause inter-particle attachments (aggregation), whereas the collision frequency factor (β) represents the mechanical properties of fluids to induce inter-particle collisions. In experimental and modeling applications, the collision efficiency factor (α) is generally used as an application-specific fitting parameter, and the collision frequency factor (β) is generally applied as a fixed theoretical function correlated with shear rate (β(i, j) = (G/6)·(Di + Dj)3)[9]. Based on experiments, Ding et al. [36] recently modified these collision efficiency and frequency factors by incorporating the concept of critical size (Dc), which subdivides two different aggregation kinetic regions (the fast and slow aggregation regions) with respect to particle/floc size. In our preliminary modeling and simulation study, application of the critical size concept was found to prevent particle/floc overgrowth beyond the highest particle/floc size and consequently to minimize unexpected mass losses caused by mass escaping out of the defined size range. Thus, in this research, Ding’s critical size (Dc) was used as a limiter to prevent unrealistic particle/floc overgrowth in aggregation and break-up kinetics. Eqs. (12) and (13) represent collision efficiency and frequency functions, respectively. Di is the diameter of the ith class particle/floc, and Dc is the critical diameter at which 50% of the collisions are successful in forming aggregates [36].
With respect to particle/floc breakage kinetics, the size-dependent kinetic function shown in Eq. (14) has been commonly applied in previous studies [9, 30, 36, 37]. To simulate the fate of the broken fragments, the binary break-up distribution function was applied in our pond simulations because of its simplicity and robustness in computation. In the discretized PBE scheme, the binary break-up distribution function becomes 2, because one parent particle/floc is assumed to produce two equally sized daughter fragments in the break-up process (Eq. (15)) [30, 36]. In Eqs. (14) and (15), a0 is the selection rate constant, and Vi is volume of the ith class particle/floc size.

2.2. Numerical Simulation

In the first step of the CFD-DPBE simulation procedure, the commercial CFD code (FLOW-3D) was used to generate a steady state flow field in the model pond. RANS and the two-equation k–ɛ turbulence models were selected to simulate flow velocities and turbulence. This resulted in nodal values for (〈Ux〉, 〈Uy〉, 〈Uz〉, k, and ɛ) (Eqs. (1)(7)). Three different flow conditions that represent low, moderate, and high turbulence conditions were simulated with FLOW-3D, thus data resulting from three steady state flow fields were obtained and saved for the subsequent multi-dimensional DPBE simulation.
After the CFD simulation, the multi-dimensional DPBE was solved with an in-house program based on the finite-difference method and codified with MATLAB (MathWorks, Naticks, MA, USA). In these simulations, two significant numerical obstacles were identified and overcome in our preliminary research. First, the complexity and nonlinearity of a large number of coupled partial differential (DBPE) equations in an advection-dominated application resulted in a computational overload. To increase computational efficiency, we applied an operator-splitting algorithm in which particle/floc advection was split from particle/floc dispersion-reaction [3840] (Table 1). Thus, the advection terms and the dispersion-reaction terms were applied sequentially at each time step. Second, a standard central-differencing finite difference method was not optimal for simulating advection-dominated flow conditions with high Peclet numbers. Previous studies have shown that upwind differencing methods produce much improved results for a given node separation [4144]. Leveque’s flux-corrected upwind algorithm was applied to solve scalar transport equations in advection-dominated conditions [41, 45]. In this algorithm, particle/floc concentration (ni) of a computational cell was updated at each time step with the inflow and outflow, which are determined by velocities through cell interfaces and concentrations of neighbor computational cells at each time step. Table 1 outlines the numerical scheme used to solve the multi-dimensional DPBE with operator splitting and flux-corrected upwind algorithms. At each time step, particle/floc advection equations were solved with LeVeque flux-corrected upwind scheme and stepwise particle/floc dispersion-reaction equations were calculated implicitly with the Gauss–Seidel iterative method.
Fig. 2 shows schematic diagrams of a flocculant-aided sediment retention pond that consists of a turbulent mixing zone at the inlet and a subsequent sedimentation basin. This turbulent mixing zone may function as an effective flocculation basin with high fluid turbulence. Chemical flocculant is assumed to be injected at the inlet of the pond, so particles/flocs will begin aggregating immediately after entering the basin.
Fig. 3 shows the two-dimensional computational domain representing a simplified turbulent mixing zone with dimensions of 2 m × 10 m (height × length). The size of each computational cell was set to 0.2 m × 0.2 m. Both inlet and outlet were treated as continuous boundaries (Fluxin = Fluxout), while the water surface was treated as a closed boundary (Fluxout = 0). The bottom layer of the mixing zone was set as a closed boundary for fluid but an open boundary for settling particles/flocs. In other words, for simplification of the model pond system, settling particles/flocs were allowed to move through the bottom layer of the zone, thereby leaving the domain, while fluid remained in the computational domain. The volumetric influent flow rate was initially set to the fixed value of 8 m3/m/min, which is equivalent to 2.5 min of mean hydraulic residence time (tmean = Volume/FlowRate)within the model mixing zone. However, to create different levels of fluid turbulence and to compare the effects of turbulent intensity on flocculation efficiency, influent flow velocities were set at three different values (0.222, 0.334, and 0.667 m/s) by adjusting inlet width. Influent clay particles (monomers) were assumed to be spheres with a diameter of 1 μm and density of 2.65 g/L. The influent solid concentration was set to 2 g/L, which is equivalent to a particle number concentration of 1.47 × 1015/m3.
In the CFD-DPBE simulation, three empirical model constants (Df, Dc, and a0) were used for aggregation and break-up kinetics. The fractal dimension (Df) was selected as 2.5, which is an intermediate value in the dataset obtained from previous studies [33, 35, 46, 47]. A critical diameter (Dc) and breakage kinetic constant (a0) were chosen rather arbitrarily as 100 μm and 10/s following Ding’s recent flocculation theory [36]. However, these constants are previous site-specific values, so it is ultimately recommended that more applicable constants be measured with settling and kinetic experiments appropriate for retention pond applications.

3. Results and Discussion

In the CFD simulation with the commercial FLOW-3D code, three steady state flow fields were obtained for the model-mixing zone. These flow fields are shown in Fig. 4 with case 1, low; case 2, moderate; and case 3, high turbulence conditions, which were induced by the different influent flow velocities of 0.222, 0.334, and 0.667 m/min, respectively. The arrows and contours in Fig. 4 represent mean flow velocity vectors (〈U〉) and shear rate distributions (G = (ɛ/ν)1/2), respectively. In the low turbulence condition (case 1), velocity vectors were uniformly directed from the inlet to the outlet, and shear rates were limited to a low level with a maximum shear rate of 13.5/s (Fig. 4(a)). However, in the high turbulence condition (case 3), a swirling zone was identified above the inlet, and high shear rates, with a maximum of 79.3/s, were observed near the inlet (Fig. 4(c)). The moderate turbulent flow condition (case 2) showed flow characteristics intermediate between the two extreme cases (Fig. 4(b)). Later in this paper, we will illustrate the theoretical effects of the velocity and shear rate distributions on flocculation efficiencies.
With steady state flow field data obtained from the CFD simulation, solutions to the multi-dimensional DPBE were obtained with an in-house program. At the beginning, the consistency and stability of the developed numerical algorithms were tested by monitoring solid mass balances and particle/floc size evolution.
Mass balances were calculated with Eq. (16) and monitored as shown in Fig. 5(a). In Eq. (16), Massin,acc, Massout,acc, Massdeposit,acc, and Massretained represent time-integrated masses caused by in-flux at the inlet, outflux at the outlet, deposition on the bottom, and retention in the pond, respectively, with time progression. Theoretically, Massin,acc should be equal to the sum of Massout,acc, Massdeposit,acc, and Massretained, that is, the mass balance calculated from Eq. (16) should be 100%. Contrary to our expectation, mass balances for low, moderate, and high turbulence conditions were all below, or slightly below, 100% at steady state conditions (99.7%, 97.8%, and 96.1%, respectively). However, these balances were in an acceptable error range, considering the approximating nature and complexity of the numerical methods. In Fig. 5(a), the mass fractions by particle/floc deposition on the bottom (Massdeposit,acc/Massin,acc) are also shown for the three different turbulence conditions. The mass fraction by deposition in the high turbulence condition (case 3) was found to be much higher than that in the low turbulence condition (case 1) because high turbulence enhanced flocculation and the subsequent sedimentation processes. These mass fractions and balances became stabilized as the mixing zone systems approached steady state conditions.
Similarly, the mass-weighted mean particle/floc diameter (Dmm), defined by Eq. (17) [48], flowing through the outlet was tracked with time progression to check numerical stability. In Eq. (17)mi represents the mass of all particles in the ith particle/floc size class, and M represents the total mass for all particle/floc size classes.
In Fig. 5(b), after the fastest growing phase, mass mean particle/floc diameters (Dmm) oscillated and then appeared to stabilize gradually. As mentioned in the previous section, the CFD-DPBE model consists of highly coupled nonlinear equations, which may produce fluctuating results. Strogatz [49] discussed the tendency of such nonlinear equations to produce oscillatory behavior in numerical simulations. A variety of phenomena can contribute to this behavior, including ‘chaotic’ behavior. Thus, the observed oscillatory behaviors shown in Fig. 5(b) were ascribed to the nonlinear nature of the CFD-DPBE model. Such behavior should be examined closely in future experimental and modeling studies.
After examining the consistency and stability of the CFD-PBE model and the mass-weighted mean particle/floc size (Dmm), solid concentration distributions at steady state conditions were investigated in the model-mixing zone. Figs. 6 and 7 show the distributions of mass-weighted mean particle/floc size and solid concentration, respectively, in the three different turbulent flow fields. In case 1, with low turbulence, mass mean particle/floc sizes were less than 27 μm, and solid concentrations were near-homogeneously distributed without particle/floc deposition. However, in case 3 with high turbulence, mass-weighted mean particle/floc sizes grew to 195 μm, which is a sufficient size for the particles/flocs to escape the computational system by settling and being deposited on the bottom. Thus, a longitudinal gradient of solid concentrations was observed in the computational domain due to particle/floc sedimentation. The moderate turbulent flow condition produced results approximately midway between the two extremes. Further, the swirling zones above the inlet in cases 2 and 3 were found to work as small flocculation compartments. Particles/flocs traveling through these swirling zones are more exposed to flocculation and thus tend to grow larger than those passing through the other zones. For example, in case 3, particles/flocs in the swirling zone grew to about 200 μm, while those right next to the swirling zone remained below 50μm.
Table 2 summarizes results from the CFD-PBE simulations after reaching steady state. Mass-weighted mean particle/floc size (Dmm) and deposited mass fraction (Massdeposit,acc/Massin,acc) in case 3, with the highest influent flow velocity and shear rate, were up to 7.5 and 12.1 times higher than those in case 1, with the lowest influent flow velocity and shear rate, respectively. As expected, turbulence in sediment retention ponds enhance flocculation efficiency in the mixing zone, at least to a certain extent. In Fig. 8, the cumulative mass distributions of particles/flocs flowing through the outlet are shown for the three different turbulence conditions studied. As expected, for the low turbulence condition (case 1), the particle/floc size distribution was more weighted toward the small size range than that for the moderate and high turbulence conditions (cases 2 and 3). Thus, in case 1, raw clay particles coming through the inlet are not aggregated properly in the turbulent mixing zone, so a large fraction of particles/flocs may not settle appropriately in the subsequent sedimentation basin. In conclusion, considering the results in Table 2 and Fig. 8 from the steady state CFD-DPBE simulations, the turbulence conditions in a turbulent mixing zone were found to have important effects on both flocculation and subsequent sedimentation efficiencies. Optimization of this situation is an important topic for both experimental and theoretical future studies.

4. Conclusions

The main purpose of this research was to estimate the applicability of a novel CFD-DPBE combined model in simulating flocculation and sedimentation processes in the turbulent mixing zone of a sediment retention pond. In this modeling and simulation study, the following important findings were identified and discussed:
  1. The employed CFD software (FLOW-3D) was a useful tool for generating steady state flow field data, such as flow velocities and shear rates, which were used in subsequent multi-dimensional DPBE simulations.

  2. As an alternative to QMOM, DPBE formulation was applied to simulate a multi-dimensional flocculation/sedimentation process. Solution of the multi-dimensional DPBE provides more readily understandable results for engineers and scientists with slightly more computations than QMOM but well within the capabilities of modern personal computers for two-dimensional flow fields.

  3. A standard, central-differencing, finite difference approach was judged inadequate for simulating the flocculation and sedimentation processes in sediment retention ponds due to computational instability caused by nonlinearity, advection dominance, and complexity of the DPBE model. Thus, operator-splitting and LeVeque flux-corrected algorithms were applied to overcome the computational difficulties. The detailed numerical model is available from the authors upon request.

  4. In application of the CFD-DPBE model, increased turbulence was found to enhance flocculation and sedimentation efficiencies. However, methodology for optimizing this effect requires further study.

This research was limited to pure simulation work without experimental validation. Thus, in future research, batch kinetic experiments and bench- or full-scale pond tests are required to calibrate, validate, and fully understand the CFD-DPBE model. In addition, the irregular behavior shown in Fig. 5(b) requires further investigation.
In summary, the CFD-DPBE model was successfully applied to generate steady state flow field data and numerically simulate flocculation and sedimentation processes in the turbulent mixing zone of a sediment retention pond. Thus, the CFD-DPBE model was shown to be a promising simulator of flocculant-aided sediment retention ponds. Furthermore, the model may be applied to flocculation and sedimentation occurring in various natural and engineering systems, such as water/wastewater treatment, nano-material synthesis, or sediment-depositing estuary systems [6, 24, 36, 5052].


Primary funding for this study was provided by the Natural Resources Conservation Service of the US Department of Agriculture (NRCS-69-4639-1-0010) through the Changing Land Use and Environment (CLUE) Project at Clemson University. Additional support was provided by the Cooperative State Research, Education, and Extension Service of the US Department of Any opinions, findings, conclusions or recommendations expressed in this article are solely those of the authors and do not necessarily reflect the views of the USDA/NRCS.


1. Gowdy W, Iwinski SR, Woodstock G. Removal efficiencies of polymer enhanced dewatering systems. In : Proceedings of the 9th Biennial Conference on Stormwater Research & Watershed Management; 2007 May 2–3; Orlando, FL.

2. Harper HH. Current research and trends in alum treatment of stormwater runoff. In : Proceedings of the 9th Biennial Conference on Stormwater Research & Watershed Management; 2007 May 2–3; Orlando, FL.

3. Kang JH, Li Y, Lau SL, Kayhanian M, Stenstrom MK. Particle destabilization in highway runoff to optimize pollutant removal. J. Environ. Eng. 2007;133:426–434.

4. Akan AO, Houghtalen RJ. Urban hydrology, hydraulics, and stormwater quality. Hoboken: John Wiley & Sons; 2003.

5. Smoluchowski MV. Versuch einer mathematischen theorie der koagulationskinetik kolloider Losungen. Z. Phys. Chem. 1917;92:129–168.

6. Lawler DF, Wilkes DR. Flocculation model testing; particle sizes in a softening plant. J. Am. Water Works Assoc. 1984;76:90–97.

7. Hounslow MJ, Ryall RL, Marshall VR. A discretized population balance for nucleation, growth, and aggregation. AIChE J. 1988;34:1821–1832.

8. Spicer PT, Pratsinis SE. Shear-induced flocculation: the evolution of floc structure and the shape of the size distribution at steady state. Water Res. 1996;30:1049–1056.

9. Spicer PT, Pratsinis SE. Coagulation and fragmentation: universal steady-state particle-size distribution. AIChE J. 1996;42:1612–1620.

10. McGraw R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997;27:255–265.

11. Somasundaran P, Runkana V. Modeling flocculation of colloidal mineral suspensions using population balances. Int. J. Miner. Process. 2003;72:33–55.

12. Marchisio DL, Vigil RD, Fox RO. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003;258:322–334.
crossref pmid

13. Rahmani NH, Dabros T, Masliyah JH. Evolution of asphaltene floc size distribution in organic solvents under shear. Chem. Eng. Sci. 2004;59:685–697.

14. Prat OP, Ducoste JJ. Modeling spatial distribution of floc size in turbulent processes using the quadrature method of moment and computational fluid dynamics. Chem. Eng. Sci. 2006;61:75–86.

15. Runkana V, Somasundaran P, Kapur PC. A population balance model for flocculation of colloidal suspensions by polymer bridging. Chem. Eng. Sci. 2006;61:182–191.

16. Lee DG, Bonner JS, Garton LS, Ernest AN, Autenrieth RL. Modeling coagulation kinetics incorporating fractal theories: a fractal rectilinear approach. Water Res. 2000;34:1987–2000.
crossref pmid

17. Fox RO. Computational models for turbulent reacting flows. Cambridge: Cambridge University Press; 2003.

18. Marchisio DL, Vigil RD, Fox RO. Implementation of the quadrature method of moments in CFD codes for aggregation–breakage problems. Chem. Eng. Sci. 2003;58:3337–3351.

19. Kumar S, Ramkrishna D. On the solution of population balance equations by discretization: I. A fixed pivot technique. Chem. Eng. Sci. 1996;51:1311–1332.

20. Ramkrishna D, Mahoney AW. Population balance modeling: promise for the future. Chem. Eng. Sci. 2002;57:595–606.

21. White FM. Viscous fluid flow. 2nd edNew York: McGraw-Hill; 1991.

22. Heath AR, Koh PT. Combined population balance and CFD modelling of particle aggregation by polymeric flocculant. In : Proceedings of the 3rd International Conference on CFD in the Minerals and Process Industries; 2003 Dec 10–12; Melbourne, Australia.

23. Lian G, Moore S, Heeney L. Population balance and computational fluid dynamics modelling of ice crystallisation in a scraped surface freezer. Chem. Eng. Sci. 2006;61:7819–7826.

24. Schwarzer HC, Schwertfirm F, Manhart M, Schmid HJ, Peukert W. Predictive simulation of nanoparticle precipitation based on the population balance equation. Chem. Eng. Sci. 2006;61:167–181.

25. Stokes GG. Mathematical and physical papers. 1:Cambridge: Cambridge University Press; 1880.

26. Brown PP, Lawler DF. Sphere drag and settling velocity revisited. J. Environ. Eng. 2003;129:222–231.

27. Jiang Q, Logan BE. Fractal dimensions of aggregates determined from steady-state size distributions. Environ. Sci. Technol. 1991;25:2031–2038.

28. Johnson CP, Li X, Logan BE. Settling velocities of fractal aggregates. Environ. Sci. Technol. 1996;30:1911–1918.

29. Spicer PT, Pratsinis SE, Raper J, Amal R, Bushell G, Meesters G. Effect of shear schedule on particle size, density, and structure during flocculation in stirred tanks. Powder Technol. 1998;97:26–34.

30. Flesch JC, Spicer PT, Pratsinis SE. Laminar and turbulent shear-induced flocculation of fractal aggregates. AIChE J. 1999;45:1114–1124.

31. Chakraborti RK, Atkinson JF, Van Benschoten JE. Characterization of alum floc by image analysis. Environ. Sci. Technol. 2000;34:3969–3976.

32. Chakraborti RK, Gardner KH, Atkinson JF, Van Benschoten JE. Changes in fractal dimension during aggregation. Water Res. 2003;37:873–883.
crossref pmid

33. Adachi Y. Dynamic aspects of coagulation and flocculation. Adv. Colloid Interface Sci. 1995;56:1–31.

34. Miyahara K, Adachi Y, Nakaishi K, Ohtsubo M. Settling velocity of a sodium montmorillonite floc under high ionic strength. Colloids Surf. A Physicochem. Eng. Asp. 2002;196:87–91.

35. Sterling MC, Bonner JS, Ernest AN, Page CA, Autenrieth RL. Application of fractal flocculation and vertical transport model to aquatic sol-sediment systems. Water Res. 2005;39:1818–1830.
crossref pmid

36. Ding A, Hounslow MJ, Biggs CA. Population balance modelling of activated sludge flocculation: investigating the size dependence of aggregation, breakage and collision efficiency. Chem. Eng. Sci. 2006;61:63–74.

37. Parker DS, Kaufman WJ, Jenkins D. Floc breakup in turbulent flocculation processes. J. Sanit. Eng. Div. 1972;98:79–99.

38. Langseth JO, Tveito A, Winther R. On the convergence of operator splitting applied to conservation laws with source terms. SIAM J. Numer. Anal. 1996;33:843–863.

39. Aro CJ, Rodrigue GH, Rotman DA. A high performance chemical kinetics algorithm for 3-D atmospheric models. Int. J. High Perform. Comput. Appl. 1999;13:3–15.

40. Badrot-Nico F, Brissaud F, Guinot V. A finite volume upwind scheme for the solution of the linear advection–diffusion equation with sharp gradients in multiple dimensions. Adv. Water Resour. 2007;30:2002–2025.

41. Durran DR. Numerical methods for wave equations in geophysical fluid dynamics. New York: Springer; 1999.

42. Rogers SE, Kwak D. An upwind differencing scheme for the incompressible Navier–Strokes equations. Washington, DC: National Aeronautics and Space Administration; 1988.

43. Alhumaizi K. Comparison of finite difference methods for the numerical simulation of reacting flow. Comput. Chem. Eng. 2004;28:1759–1769.

44. Timin T, Esmail MN. A comparative study of central and upwind difference schemes using the primitive variables. Int. J. Numer. Methods Fluids. 1983;3:295–305.

45. LeVeque RJ. High-resolution conservative algorithms for advection in incompressible flow. SIAM J. Numer. Anal. 1996;33:627–665.

46. Bushell GC, Yan YD, Woodfield D, Raper J, Amal R. On techniques for the measurement of the mass fractal dimension of aggregates. Adv. Colloid Interface Sci. 2002;95:1–50.
crossref pmid

47. Turchiuli C, Fargues C. Influence of structural properties of alum and ferric flocs on sludge dewaterability. Chem Eng J. 2004;103:123–131.

48. Hinds WC. Aerosol technology: properties, behavior, and measurement of airborne particles. 2nd edNew York: John Wiley & Sons Inc; 1999.

49. Strogatz SH. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Reading: Perseus Books Publishing; 1994.

50. Bungartz H, Wanner SC. Significance of particle interaction to the modelling of cohesive sediment transport in rivers. Hydrol. Process. 2004;18:1685–1702.

51. Winterwerp JC. On the flocculation and settling velocity of estuarine mud. Cont.Shelf Res. 2002;22:1339–1360.

52. Maggi F, Mietta F, Winterwerp JC. Effect of variable fractal dimension on the floc size distribution of suspended cohesive sediment. J. Hydrol. 2007;343:43–55.

Fig. 1
Diagram of aggregation and break-up processes for the i = 3 particle/floc size class in the discretized population balance equation.
Fig. 2
Schematic diagram of a flocculant-aided sediment retention pond with a turbulent mixing zone, sedimentation basin, and discharge drain.
Fig. 3
Schematic diagram of the computational domain representing a simplified turbulent mixing zone in a sediment retention pond.
Fig. 4
Steady state flow field profiles from computational fluid dynamics simulation for (a) case 1, low turbulence, (b) case 2, moderate turbulence, and (c) case 3, high turbulence. Arrows and colors represent flow velocities and shear rates, respectively.
Fig. 5
(a) Mass fractions and balances and (b) mass-weighted mean floc diameter (Dmm) with respect to dimensionless residence time, which is normalized by dividing real fluid residence time by theoretical mean residence time.
Fig. 6
Mass-weighted mean floc diameter (Dmm) distributions in the computational domain. The distributions of case 1 (a), case 2 (b), and case 3 (c) are listed from the top.
Fig. 7
Solid concentration distributions in the computational domain. The distributions of case 1 (a), case 2 (b), and case 3 (c) are listed from the top.
Fig. 8
Cumulative mass distribution of particle/floc sizes at the outlet of the model basin.
Table 1
The simplified numerical algorithm for solving the CFD-DPBE model
  • Initialization

    • - Supporting data (flow field data from CFD, solid and liquid properties)

    • - Computational system layout (dimensions, mesh)

  • DPBE calculation (operator splitting algorithm)

  • Post-processing

    • - Mass balance, particle/floc diameters, solid concentrations, etc.

[i] CFD-DPBE: computational fluid dynamics–discretized population balance equation, FDM: finite difference method.

Table 2
Flow field characteristics and flocculation/sedimentation efficiencies for three different turbulent conditions in the mixing zone
Flow turbulence level Flow field characteristic Flocculation/sedimentation efficiency

vin (m/s) G a (/s) Dmmb (μm) Massdeposit,accMassin,acc(%)
Case 1 Low 0.222 13.5 24.59 1.204 %
Case 2 Moderate 0.334 28.3 105.2 4.787 %
Case 3 High 0.667 79.3 183.2 14.54 %

a Maximum values in the computational domain.

b Averaged values along the outlet.

PDF Links  PDF Links
PubReader  PubReader
Full text via DOI  Full text via DOI
Download Citation  Download Citation
Editorial Office
464 Cheongpa-ro, #726, Jung-gu, Seoul 04510, Republic of Korea
TEL : +82-2-383-9697   FAX : +82-2-383-9654   E-mail : eer@kosenv.or.kr

Copyright© Korean Society of Environmental Engineers. All rights reserved.        Developed in M2community
About |  Browse Articles |  Current Issue |  For Authors and Reviewers