Volume 61, Issue 1, January 2006, Pages 75–86

Advances in population balance modelling

Second International Conference on Population Balance Modelling

Modeling spatial distribution of floc size in turbulent processes using the quadrature method of moment and computational fluid dynamics

  • Department of Civil Engineering, North Carolina State University, 208 Mann Hall, Raleigh, NC 27695-7908, USA

Abstract

A study was performed that utilizes the quadrature method of moments (QMOM) to model the transient spatial evolution of the floc size in a heterogeneous turbulent stirred reactor. The QMOM approach was combined with a commercial computational fluid dynamics (CFD) code (PHOENICS), which was used to simulate the turbulent flow and transport of these aggregates in the reactor. The CFD/QMOM model was applied to a 28 l square reactor containing an axial flow impeller and 100 mg/l concentration of 1 μm nominal clay particles. Simulations were performed for different average characteristic velocity gradients (40,70,90, and 150 s-1). The average floc size and growth rate were compared with experimental measurements performed in the bulk region and the impeller discharge region. The CFD/QMOM results confirmed the experimentally measured spatial heterogeneity in the floc size and growth rate. In addition, the model predicts spatial variations in the aggregation and breakup rates. Finally, the model also predicts that the transport of flocs into the high shear impeller discharge zone was responsible for the transient evolution of the average floc size curve displaying a maximum before decreasing to a steady-state floc size.

Keywords

  • Quadrature method of moment;
  • Computational fluid dynamics;
  • Turbulence;
  • Aggregation;
  • Breakage;
  • Flocculation

1. Introduction

In drinking water treatment as well as the chemical process industry, flocculation bench scale and pilot plant tests are conducted frequently to gather information for use in designing or retrofitting full-scale flocculation processes. In order to accurately translate the results found during a bench- or pilot-scale to a full-scale operation, a good understanding of the flocculation process is required. A significant body of research involving both experimental and numerical analyses has been done to understand the complex dynamics of the flocculation process (Thomas et al., 1999, Ducoste, 2002 and Hopkins and Ducoste, 2003). These research studies have included the impact of solution chemistry, particle type, particle surface conditioning, chemical coagulant type, fractal configuration of floc particles, and reactor fluid mechanics on the resulting flocculation performance. Results of these studies have clearly shown that all of these conditions play a significant role in the final particle size distribution.

Over the years, several researchers have developed population balance models (PBMs) to describe the flocculation process (Lu and Spielman, 1985, Koh et al., 1987, Kusters, 1991, Spicer et al., 1996a, Flesch et al., 1999, Ducoste, 2002 and Selomulya et al., 2003). In many cases, the focus of the PBM research has been to better characterize the aggregation/breakage dynamics and the flocs’ structural evolution (Kusters et al., 1997, Flesch et al., 1999 and Selomulya et al., 2003). The modeling results in several of these research studies have shown reasonable agreement between the predicted and experimental floc size distribution and the transient evolution of the average floc size (Flesch et al., 1999, Patil et al., 2001, Ducoste, 2002 and Selomulya et al., 2003).

In most PBM research, the flow field characteristics were incorporated directly into the aggregation and breakup dynamics using tank averaged turbulent quantities (Thomas et al., 1999). A major reason for this simplification was due to the complexity in solving the PBM equations, which was also the subject of extensive and ongoing research (Ramkrishna and Mahoney, 2002). Consequently, these types of stand-alone PBMs do not completely account for the spatial heterogeneity of the turbulent flow field that may influence the aggregation/breakup mechanism. Some researchers have tried to incorporate the turbulence spatial heterogeneity into the PBM by using turbulent values from regions inside the reactor that may significantly impact either aggregation or breakup (Ducoste, 2002). Other research studies have taken a different approach to incorporating the turbulence spatial heterogeneity and tried to include the fraction of time particles spent in a limited number of regions with different levels of turbulence into the aggregation and breakup kinetics (Koh et al., 1987, Kusters, 1991 and Maggioris et al., 2000). However, both PBM modeling methods would predict the same aggregation/breakup dynamics, steady-state floc size distribution and the transient evolution of the average floc size throughout the reactor volume.

A recent study has shown that physical/chemical conditions may exist such that the average floc size may not be spatially homogenous at low mixing speeds (Hopkins and Ducoste, 2003). In Hopkins and Ducoste (2003), the data showed that floc growth, average floc size, and variance about the average floc size varied spatially in the reactor at low tank average energy dissipation rate with larger floc sizes and growth rates in the bulk region and a larger variance in the impeller discharge region. The spatial variation in these three parameters disappeared with increasing tank average energy dissipation rate. In addition, Hopkins and Ducoste found that the transient evolution of the average floc size displayed a maximum before decreasing to a steady-state floc size. This type of transient evolution confirmed other research results and was theorized to be caused by floc erosion and fracture due to more frequent exposure to higher shear forces in the impeller discharge zone (Spicer et al., 1996b). The transient behavior in the average floc size as well as the spatial heterogeneity seen by Hopkins and Ducoste cannot be predicted by only solving the PBM equations.

One possible approach to account for the spatial heterogeneity in turbulence is to combine a computational fluid dynamics (CFD) model that describes the convection and turbulent diffusion of the floc throughout the reactor volume with the PBM equations. However, PBM equations structured in the traditional particle size class approach require the solution of a large number of scalars (i.e., one for each particle size class) at every discrete location in the reactor volume. Consequently, the combined convection/turbulent diffusion PBM equations would be computationally expensive if implemented in a commercially available CFD software.

As an alternative to the traditional PBM class equations, researchers have developed PBMs based on the method of moments (McGraw, 1997). One particular method of moment approach is the quadrature method of moment (QMOM). In the QMOM approach, the floc size evolution is tracked by solving a system of differential equations for lower order moments and utilizes a quadrature approximation as closure for these moment equations (McGraw, 1997). The QMOM approach has been used by researchers to investigate aerosol dynamics (McGraw, 1997, Wright et al., 2001 and McGraw and Wright, 2003), floc aggregation–breakage kinetics (Marchisio et al., 2003a and Marchisio et al., 2003b) and has also been combined with CFD for modeling simultaneous aggregation and breakage in a Taylor–Couette device (Marchisio et al., 2003c). The QMOM approach may not accurately predict the steady-state floc size distribution due to the limited number of tracked size classes. However, the QMOM approach does provide enough statistical information to predict the transient evolution of the average floc size and the reactor floc population, which are important parameters that can be used to evaluate the flocculation process. Until this study, the CFD/QMOM approach had never been used to evaluate the flocculation process in a turbulent heterogeneous environment.

The purpose of this study was to evaluate the CFD/QMOM approach in modeling a turbulent heterogeneous flocculation process. In this study, the QMOM approach was combined with a commercial CFD software to predict the convection and turbulent dispersion of flocs undergoing aggregation–breakage in a stirred reactor containing an axial flow impeller. The predicted average floc size and growth rate were compared to experimental measurements in the bulk region and in the impeller discharge region. Simulations were performed at 4 average characteristic velocity gradients (40, 70, 90, and 150 s-1) in a 28 l square reactor containing 100 mg/l of 1 μm clay particles.

2. Numerical methods

2.1. Population balance modeling and quadrature method of moments

The general transport equation that includes the kinetics of aggregation and breakup for a floc particle is described in Eq. (1). In Eq. (1), (I) is the transient term, (II) is the convective term, (III) is the diffusive term and (IV) is the source term describing aggregation–breakup dynamics. In (IV), term (A) represents the generation of size (d) flocs due to the collisions of smaller particles and term (B) represents the reduction of size (d) flocs due to aggregation with other floc sizes. Term (C) represents the generation of size (d) flocs due to breakup of larger flocs and finally, term (D) represents the reduction of size (d) flocs due to breakage into smaller floc sizes. In Eq. (1), (α) is the collision efficiency function, (β) is the collision frequency function, (kb) is the breakup frequency function, (F) is the fragment distribution function, (n) is the number concentration of flocs, (ρ) is the density, (kε) is the turbulent kinetic energy, and (εp) is the turbulent energy dissipation rate. The population density function is a multivariable function n(x,d,t) where x=(xi,xj,xk), d is the particle linear size and t the time dependence.

equation1
View the MathML source

In this study, a QMOM was used to track the particle size evolution by solving a system of differential equations for lower order moments. Using this approach, the number concentration in the aggregation and breakage terms are converted to moment terms using a quadrature approximation displayed in Eq. (2) (McGraw, 1997):

equation2
View the MathML source
(mk) is the kth moment value, (Li) is the abscissa or dimensionless size indicator and (wi) is the weight or dimensionless number concentration indicator. The moment transformation using Eq. (2) was applied to Eq. (1) describing aggregation–breakup and is displayed in Eq. (3):
equation3
View the MathML source
where the source term (IV) is described as follows:
equation4
View the MathML source
(Δmka) and (Δmkb) represent the kth moment aggregation and breakup source, respectively. QMOM is based on the theory of orthogonal polynomial and Gaussian quadrature approximation ( Press and Teukolsky, 1990 and McGraw, 1997). Starting from a moment sequence consisting of Nd moments, the quadrature procedure begins by finding Nd abscissas (Li) and Nd weights (wi). Li and wi can be extracted from the moments using Gordon's (1968) product-difference (PD) algorithm (McGraw, 1997).

The abscissas and weights extraction from the moments can also be performed by using a modified sequence of moments, according to Prat and Ducoste (2004) and combined with Wheeler's algorithm (Press and Teukolsky, 1990). In both cases, the quadrature procedure was simplified to calculate only eigenvalues, which are the abscissas and eigenvectors where the weights are determined from the first component of each eigenvector of a (Nd×Nd) symmetric matrix. In this study, the second approach was used. In addition, a three-point (Nd=3) QMOM was used in the present work since it has been shown to sufficiently track the lower order moments (m0, m1, m2, m3, m4, m5) with reasonable accuracy (Marchisio et al., 2003c). This was confirmed in the present study by performing the QMOM approach with Nd=3 and 4.

2.2. Flocculation dynamics

In this study, a collision frequency function (β) (Eq. (5)) based on orthokinetic flocculation assuming turbulent shear (Saffman and Turner, 1956) was used.

equation5
βij=1/6.18(Li+Lj)3(εp/ν)0.5.
The agglomeration part of the source term includes a collision efficiency function (α) that accounts for unsuccessful particle collisions due to electrostatic repulsion or hydrodynamic retardation. The collision efficiency was calculated using a modified version of Adler's flow number equation (Eq. (6)), which represents the ratio of hydrodynamic shear forces to van der Waals forces between colliding floc particles (Adler, 1981; Kusters et al., 1997). The collision efficiency function using terms from the QMOM approach is described in Eq. (6):
equation6
View the MathML source
In Eq. (6), C1 represents an empirical constant that accounts for any reduction in collision efficiency due to residual electrostatic repulsion forces. Terms (μ), (ν), (d0) and (A) are, respectively, dynamic and kinematic viscosity, initial mean particle diameter (in meters) and Hamaker constant View the MathML source.

The breakup frequency function (kb) that describes fracture of particles in the viscous dissipation subrange was used and is shown in Eq. (7) (Kusters, 1991):

equation7
kbi=[4/(15π)](εp/ν)0.5exp[-C2(Liεp)].
The empirical constant, C2, represents the strength of floc particles and their resistance to shear forces. Both aggregation–breakup empirical constants (C1, C2) are functions of a number of parameters such as particle type, water and coagulant chemistry, and temperature. C1 and C2 were determined by fitting the CFD/QMOM model results with 40 s-1 experimental growth rate and steady-state average floc size as described in Section 3.

For the fragment distribution (View the MathML source), the CFD/QMOM assumes that floc breakage results in the formation of two fragments with mass ratio 1:4. (View the MathML source) is shown in Eq. (8):

equation8
View the MathML source
Kusters (1991) observed experimentally that parent flocs tend to breakup into uneven fragments with at least one daughter particle with a radius of about 70–80% of the parent floc radius. The fragment distribution function in Eq. (8) is consistent with Kusters’ experimental observation. In the present model, particles were allowed to break, once they reached a sufficient size such that the corresponding size indicator (abscissa) of the smallest daughter particle was equal to the smallest abscissa (L1) value of the initial floc size distribution. Based on Eqs. (3)(8), a CFD/QMOM model was developed with the commercial CFD code PHOENICS. The transport equations for the first six lower-order moments were implemented in the CFD code via user defined subroutines. At each time step of the transient simulation, moments were updated and new abscissas and weights were extracted from the moment sequence using turn-key routines (ORTHOG, GAUCOF) from Press et al. (1992). A schematic of the QMOM numerical procedure is shown in Fig. 1. The ORTHOG routine computes the coefficients of the Jacobi matrix by using the Wheeler's algorithm (Press and Teukolsky, 1990) and the GAUCOF routine computes the abscissas and the weight of the Gaussian quadrature formula. GAUCOF calls two other routines; TQLI that determines the eigenvalues and the eigenvectors of the symmetric, tridiagonal Jacobi matrix and EIGRST that sorts the eigenvalues in descending order (Press et al., 1992).

Full-size image (28 K)
Fig. 1. 

Flowchart for the determination of abscissas (Li) and weights (wi) from the moments (mi).

3. Numerical study

In this work, a three-dimensional simulation of a 0.114 m diameter A310 foil impeller in a square tank with dimensions of View the MathML source was performed using PHOENICS. The flocculation simulations were performed in two steps. First, the flow field was solved for four tank average characteristic velocity gradients (Gm=[εpmean/ν]0.5)=40,70,90, and 150 s-1 by changing the impeller velocity boundary condition. The fluid flow field within the reactor was determined by solving the equations of motion (White, 1991). The equations of motion include the Reynolds-averaged continuity and Navier Stokes equations. In this study, the two-equation k-ε model was used to characterize the turbulence in the reactor (White, 1991). They are called two-equation models because, in addition to the Reynolds-averaged continuity and Navier Stokes equations, the two-equation turbulence model requires writing balance equations for the turbulent kinetic energy, k, and the turbulent energy dissipation rate, ε. (White, 1991). For each case, calculations were performed until a steady-state solution was reached. Fig. 2 displays velocity vectors and the turbulent energy dissipation for an average characteristic velocity gradient View the MathML source in the mid-plane perpendicular to the reactor walls.

Full-size image (66 K)
Fig. 2. 

Velocity vectors and local turbulent energy dissipation rate for an average characteristic velocity gradient View the MathML source.

In Fig. 2, two major zones of interest in the stirred tank have been outlined: (a) the bulk zone, located above the impeller, and (b) the discharge zone, located below the impeller. In this plane, seven locations have been selected to investigate the local behavior of the flocculation process. At these seven locations, the transient evolution of the normalized average floc size (d43n), abscissas (Li), weights (wi), moments (mk) and moments’ rate (Δmk) were observed. Point (1) was located in the center of a recirculation zone (RZ) just outside the impeller discharge zone. Points (2) and (3) were located in the bulk zone (BZ), far from the impeller region, while points (4)(7) were located at various points in the shear zone (SZ) below the impeller.

Once the flow fields were determined for each (Gm), Eq. (3) was solved for each moment (m0, m1, m2, m3, m4, m5) using velocities and turbulent values (i.e., turbulent energy dissipation rate and turbulent kinetic energy) as initial conditions. Initially, a uniform floc distribution was set for the entire computational domain. In the flocculation experiments, the initial floc distribution was not strictly mono-dispersed and ranged from 0.5 to 5 μm (Ducoste and Clark, 1998). As in Ducoste (2002), an average mean diameter of about 2 μm was adopted for the initial average floc size in this study. The initial condition for the floc size distribution was set by selecting three abscissas (L1=1; L2=2; L3=3) and three weights (w1=0.9; w2=0.08; w3=0.02), giving the initial values of the six moments (m0=1; m1=1.12; m2=1.40; m3=2.08; m4=3.80; m5=8.32) and the initial average particle size (d43i=1.83). This three-class floc size distribution assumes that, respectively, 90%, 8%, and 2% of the total number of particles have a diameter of 1, 2, and 3 μm. The values of empirical constants (C1) and (C2) were determined by performing a least-squares fit between the predicted steady-state average floc size (d43) and the experimentally measured average floc size. Adjustments were performed using the average characteristic velocity gradient (View the MathML source data set (Ducoste and Clark, 1998 and Hopkins and Ducoste, 2003). The first empirical constant (C1), which represents the aggregation efficiency, was determined with the slope of theaverage floc size evolution (d43) during the linear growth period at the beginning of the flocculation process where aggregation exceeds breakup. The other empirical constant (C2), which represents the strength of floc particles and their resistance to shear forces, was determined with the experimental steady-state value reached by the average floc size at the end of the flocculation process. Table 1 summarizes empirical constants (C1 and C2) and values of initial abscissas (Li) and weights (wi) used in the CFD/QMOM model.

Table 1.

Values for empirical constant and initial value of the PSD used in the CFD/QMOM model

Parameters
C1C2w1w2w3L1L2L3
242.50.900.080.02123
Full-size table

In the CFD/QMOM approach, the average floc size was calculated as (d43=m4/m3). m3 represents the total mass of the floc particles in the reactor and was used as an indicator of mass conservation throughout the simulation. In addition to local spatial values for (d43n=d43/d43i), a normalized volume-weighted average floc size (d43n) and average floc density (m0) were computed for the reactor volume as

equation9
View the MathML source
equation10
View the MathML source
where (Vj) is the volume of each (j) computational cell. The slopes of the initial growth region and the volume-weighted variance about the steady-state average floc size were defined as
equation11
Slope=Δ〈d43n〉/Δtime,
equation12
View the MathML source
Slope, d43n, and variance were also computed at points 1–7 in Fig. 1 in order to highlight the spatial behavior of the flocculation mechanism.

As discussed earlier, numerical results were compared to experimental measurements performed in a previous study (Hopkins and Ducoste, 2003). A detailed description of the experimental setup and procedure is provided in their study. Hopkins and Ducoste utilized a photometric dispersion analyzer (PDA). The PDA is an optical technique that measures light fluctuations in a flowing suspension. The light fluctuations are generated from random variations of the suspended particle size that is assumed to follow a Poisson distribution. Gregory (1985) showed that the root mean square (VRMS) value of the light fluctuations is a function of the square root of the particle concentration and particle size. By measuring the VRMS of a fluctuating signal and the average light intensity (VDC), the change in the flocculation performance can be determined through the ratio of VRMS/VDC. While this ratio does not directly indicate the particle size, it has been shown to describe the degree of flocculation (Gregory and Nelson, 1986, Gregory and Chung, 1995, Kang and Cleasby, 1995 and Matsui et al., 1998). More importantly, larger ratio values imply larger aggregate sizes for a given suspension (Burgess et al., 2002). In this study, the ratio values were normalized with the initial ratio value at the start of the flocculation process. This dimensionless ratio value will be compared to d43n (Eq. (9)).

4. Results and discussion

Fig. 3 displays the mean slope of the floc size growth region, normalized volume weighted-mean floc size (d43n), and the transient evolution of (d43n) for (Gm)=40 and 70 s-1. In Figs. 3a and b, the results confirm experimental observation that the floc growth rate increases with increasing (Gm) while average floc size decreases with increasing (Gm) due to higher breakup rates (Parker et al., 1972, Tambo and Wanatabe, 1979 and Oles, 1992; Spicer and Pratsinis, 1996a and Spicer and Pratsinis, 1996b; Selomulya et al., 2001 and Hopkins and Ducoste, 2003). In Fig. 3c, the transient evolution of d43n confirmed experimental observations that the floc size goes through an initial growth phase due to the aggregation rate exceeding the breakup rate (De Boer et al., 1989, Spicer et al., 1996a and Kusters et al., 1997). After this initial growth phase, the floc growth rate decreases due to an increasing floc breakup rate until a steady-state size is reached. This transient floc size evolution was apparent for the (View the MathML source mixing condition over the 30 min time period as measured by Hopkins and Ducoste (2003). Moreover, for the higher mixing conditions (i.e., (View the MathML source) (Figs. 3d and 4), the model predicts a peak region in the floc size evolution curve before decreasing to a final steady-state value. The lower steady-state value was typically 10–15% smaller than the peak value and was achieved in a shorter time period with increasing (Gm) (Fig. 4).

Full-size image (53 K)
Fig. 3. 

Growth rate slope vs. (Gm) value, (b) steady-state average floc size vs. (Gm), (c) average floc size evolution (〈d43n〉) for View the MathML source, (d) average floc size evolution (〈d43n〉) for View the MathML source.

Full-size image (32 K)
Fig. 4. 

Simulated tank averaged floc size evolution (〈d43n〉) for different values of (Gm).

In Hopkins and Ducoste (2003), the dimensionless ratio value displayed a higher peak region than the predicted d43n. The dimensionless ratio values were 20–30% higher than the dimensionless ratio at steady-state. The higher dimensionless ratio values relative to the final steady-state value was likely due to the structure of the floc aggregate that is typically non-spherical and is often described using fractals (Flesch et al., 1999, Ducoste, 2002 and Selomulya et al., 2003). In addition, the dimensionless ratio value also displayed an earlier decline to the steady-state value than the model results. This earlier decline in the experimental transient evolution is likely due to the erosion or restructuring of primary particles within the floc, which is an additional breakup mechanism not simulated in the current CFD/QMOM model.

The predicted transient evolution in d43n at the higher mixing conditions has also been observed experimentally by others (Spicer et al., 1996b, Bouyer et al., 2001 and Selomulya et al., 2001). Bouyer et al. (2001) experimentally observed a peak in the floc size evolution in a 2 l square reactor using a radial impeller for (View the MathML source. Selomulya et al. (2001) performed experiments in a stirred tank using an axial flow impeller (three-blade fluid foil Lightnin A310 impeller) and observed similar transient behavior for shear rates varying from 40 to 100 s-1 over a 100 min flocculation period. Spicer et al. (1996b) also noted a peak in the transient floc size evolution curve at (View the MathML source for flocculation times of about 3 h. Results of Bouyer et al., and Spicer et al. suggest that the lower mixing conditions could also exhibit a peak followed by a lower steady-state floc size at significantly longer flocculation mixing periods. Although flocculation simulations performed in this study were typically carried out for 30 min to match the standard duration time for water treatment flocculation processes, simulations performed for an equivalent 4-h period also displayed a peak followed by a lower steady-state floc size. The transition to the lower steady-state size at (View the MathML source occurred at 150 min; confirming the experimental results seen by others.

Researchers have theorized that the reduction in the steady-state floc size following a peak value was the result of primary particle erosion from the floc structure or fracture of the floc due to prolonged shearing in the high turbulence zones near the impeller (Williams et al., 1992, Spicer et al., 1996b, Selomulya et al., 2001, Bouyer et al., 2001 and Hopkins and Ducoste, 2003). Hopkins and Ducoste further suggested that the transient floc size evolution shown in Figs. 3c and 4 were the result of spatial heterogeneity in the floc size due to spatial variations in the aggregation and breakup rates. Larger flocs that are found in the bulk region and continuously formed due to lower breakup rates are transported to the impeller region where they are fractured to smaller floc sizes due to higher breakup rates. The degree of spatial variation in floc size would increase as (Gm) decreases. Hopkins and Ducoste measured the local variance at a point in the bulk region and the impeller discharge region and showed that the variance does decrease with increasing (Gm). However, only two points were measured. Numerical simulations were performed to investigate the degree of spatial variation in the floc size and are reported in Figs. 5 and 6.

Full-size image (137 K)
Fig. 5. 

Contour plot of particle size (d43n): (a) View the MathML source 1 min time period, (b) View the MathML source 5 min time period, (c) View the MathML source 15 min time period, (d) View the MathML source 1 min time period, (e) View the MathML source 5 min time period, (f) View the MathML source 15 min time period.

Full-size image (52 K)
Fig. 6. 

Tank-averaged variance vs. (Gm); (b) floc size evolution for View the MathML source at different locations in the reactor; (c) floc size evolution for View the MathML source at different locations in the reactor.

Fig. 5 displays, in the vertical symmetry plane, the spatial distribution of the normalized average floc size (d43n) for (View the MathML source and (View the MathML source at 1, 5 and 15 min into flocculation. Results show that regardless of (Gm), larger flocs can be found in the center of recirculation. In Figs. 5a–c for (View the MathML source, particles continue to grow throughout the 30 min simulated flocculation time. For higher shear rate, (View the MathML source (Figs. 5d–f), contour plots for 1,5 and 15 min show that flocs evolve everywhere through a maximum size before decreasing to a lower steady-state size.

Fig. 6a displays the volume-weighted variance about the steady-state average floc size as a function of (Gm). In Fig. 6a, the variance was found to decrease with increasing (Gm), confirming the limited spatial measurements of Hopkins and Ducoste (2003). Figs. 6b and c display the transient average floc size evolution for (Gm)=40 and 90 s-1, respectively, at different points shown in Fig. 2. As can be seen in Figs. 6b and c, the transient evolution of the local mean floc size at the different locations all have the same shape as the transient evolution of the tank-average mean floc size. However, the value of the local mean floc size does depend on the location in the reactor. For the steady-state value, Figs. 6b and c also show that the average floc size was larger in the center of the recirculation zone (pt 1), than in the bulk region (pt 3) and the impeller discharge zone (pt 5). The model also predicts that flocs in the bulk zone achieve a larger steady-state value than in the impeller discharge zone. While no experimental measurements were made in the center of the recirculation zone, Hopkins and Ducoste (2003) also found a larger steady-state PDA ratio value in the bulk region than in the impeller discharge zone.

In Figs. 6b and c, the model results also show that the fastest growth rates were found in the center of the circulation (pt 1). It was theorized that in the center of this recirculation zone, the flow field tends to prolong the contact between particles without significant floc breakup due to a lower local turbulence energy dissipation rate than the impeller discharge zone. In comparison, growth rates in bulk (pt 3) and impeller discharge zones (pt 5) were slightly smaller than the tank averaged value with the growth rate at (pt 3) slightly higher than the growth rate at (pt 5). The results in Figs. 6b and c seem to suggest that the recirculation zone (pt 1) plays a key role on the dynamics of the initial growth region for the entire reactor. In addition, the numerical results confirm experimental observations of the spatial variation in floc growth (Hopkins and Ducoste, 2003).

Ducoste (2002) suggested that cluster–cluster collisions could be used to characterize the type of aggregation in the bulk region or regions outside the impeller discharge zone at low (Gm) since the bulk region may contain turbulent zones that favor these types of collisions (Danielson et al., 1991 and Hansen and Ottino, 1996). Consequently, the onset of cluster–cluster aggregation in the bulk region would make the growth rate seem larger. In addition, the lower turbulent energy dissipation rate in the bulk region would allow the formation of larger steady-state average floc sizes. These larger flocs would then be fractured as a result of the higher turbulent environment in the impeller discharge region, which would consequently produce a lower steady-state floc size (Ducoste, 2002). In the impeller region, particle collisions may be mostly dominated by primary–primary or primary–cluster aggregation. The resulting growth rate in the impeller zone would therefore be lower due to these types of particle collisions albeit under higher local shear rates. However, the numerical results in this study displayed only a slight difference between the growth rate in the impeller discharge zone (pt 5) and the bulk zone (pt 3). One possible explanation is that the model does not account for the porous floc structure and how the structure may change due to the heterogeneous turbulent environment (Selomulya et al., 2003). In their study, Selomulya et al. showed that by modeling changes in the floc structure through changes in the fractal dimension, the transient evolution of the average floc size displayed a peak followed by a lower steady-state value. However, the results of this study showed that changes in floc structure were not the only reason for this type of transient floc size evolution. The results in Figs. 5 and 6 suggest that the transport of larger bulk region floc sizes into the high turbulent zone does contribute to this transient average floc size evolution. To better understand the production of the average floc size transient evolution curve, a closer investigation of the moment 4 rate equation (Δm4) was performed.

Fig. 7 displays (Δm4) transient evolution during the flocculation process. In Figs. 7a and b, the tank average (Δm4) was displayed along with aggregation (Δm4a) and breakage (Δm4b) portions for (View the MathML source and (View the MathML source, respectively. For (View the MathML source (Fig. 7a), results indicate that tank average Δm4 increases for the first 2 min since the breakage rate was negligible and was only influenced by the aggregation rate. This phase corresponds to the linear growth rate region of the transient average floc size evolution. When flocs reach a certain size, the breakage rate begins to increase and slow down the growth of flocs. Finally, after 10 min into flocculation, both the aggregation and breakage rates were equal resulting in a steady-state (m4) value. For (View the MathML source, the tank average (Δm4) always remains positive until it falls to zero within the 30 min of flocculation. As a result, the transient average floc size evolution results in a growth region followed by a single steady-state region as shown in Fig. 3c within that 30 min time frame.

Full-size image (75 K)
Fig. 7. 

Tank average moment 4 m4) rate equation with aggregation m4a) and breakage m4b) contributions for View the MathML source. (b) Tank average moment 4 m4) rate equation with aggregation m4a) and breakage m4b) contributions for View the MathML source. (c) Local variation of m4)View the MathML source. (d) Local variation of m4)View the MathML source.

In cases showing a peak followed by a lower steady-state value in the transient average floc size evolution such as for (View the MathML source (Fig. 6c), the modeling results show that the tank average (Δm4) stays positive during the initial floc growth stage and becomes negative once the breakage rate (Δm4b) exceeds the aggregation rate (Δm4a). This negative breakage rate was responsible for the lower steady-state floc size following a peak value (Fig. 3d) in the transient average floc size evolution. While the transient evolution of the tank-average (Δm4) would intuitively explain the dynamics of the average floc size evolution, it does not completely explain the occurrence of a larger negative breakage rate at some period during the flocculation process. A deeper understanding of this floc size evolution will come from analyzing the local moment rate (Δm4).

Figs. 7c and d, display the local moment rate (Δm4) for (View the MathML source and (View the MathML source at different tank locations. In Figs. 7c and d, the results clearly show that the local moment rate (Δm4) varies spatially in the flocculation reactor with certain regions maintaining a positive (Δm4) and other regions maintaining a negative (Δm4). As an example, (Δm4) was positive in the center of recirculation as well as in the bulk region while strongly negative in the impeller discharge zone. Fig. 8 displays contour plots of the local moment rate (Δm4) in the vertical symmetry plane for (View the MathML source and (View the MathML source. The results in Fig. 8 shows that a majority of the tank volume (85%) promotes favorable conditions for floc growth while regions in the impeller discharge zone (15% of the tank volume) promote conditions for floc breakup. The division between regions promoting growth versus breakup did not change significantly with increasing (Gm). However, the model does predict a higher breakup rate within the breakup promoting region with increasing (Gm).

Full-size image (97 K)
Fig. 8. 

Contour plot of moment 4 m4) rate equation (a) View the MathML sourceView the MathML source, (b) 40 s-1View the MathML source, (c) View the MathML sourceView the MathML source, (d) 150 s-1View the MathML source.

The results in Figs. 7c, d, and 8 seem to support the concept that floc particles, which were continuously growing in the growth promoting region, were being transported to regions of intense turbulence (or breakup promoting regions) such as the impeller discharge zone where they were fractured to smaller floc sizes. Moreover, the results in Figs. 7 and 8 confirm what researchers have theorized, that transport through the high shear zone of the impeller discharge region (Williams et al., 1992, Spicer et al., 1996b, Selomulya et al., 2001, Bouyer et al., 2001 and Hopkins and Ducoste, 2003) may result in a peak followed by a lower steady-state average floc size due to spatial heterogeneity in the floc size, aggregation rate, and breakup rate. However, final confirmation of the CFD/QMOM spatial distribution of the floc size will have to come from experimental techniques that can measure the in situ floc size at many points in a mixing reactor, simultaneously.

5. Conclusions

A study has been performed to model the transient spatial evolution of the average floc size (d43n) in a turbulent stirred reactor using a three-point quadrature method of moments (QMOM) for modeling the aggregation–breakage dynamics and computational fluid dynamics (CFD) to simulate the turbulent flow and transport of these aggregates in the stirred reactor. In the three-point QMOM, three abscissas and three weights were used to simulate the particle size distribution. Six moments (m0, m1, m2, m3, m4, and m5) were used to simulate the statistics of the particle size distribution. Simulations were performed at various shear rates (Gm=40,70,90, and 150 s-1). The CFD/QMOM flocculation model results have confirmed experimental observations that average floc size and variance about the average floc size will decrease with increasing (Gm) while the aggregate growth rate will increase with increasing (Gm).

The CFD/QMOM model results also showed that spatial heterogeneity in the turbulence energy dissipation rate promoted spatial heterogeneity in the average steady-state floc size and floc growth rate. The model confirmed experimental observation that the local steady-state average floc size and the floc growth rate were larger in the bulk regions and smaller in the impeller discharge zone. In addition, the CFD/QMOM model predicted that the maximum average steady-state floc size and floc growth rate were found in the center of the recirculation zone where prolonged contact between particles were favored at lower local energy dissipation rates that minimized floc breakup.

The modeling results also predicted a peak followed by a lower steady-state value in the transient evolution of the average floc size for simulated shear rates (Gm=70,90, 150 s-1) while there was no peak value in the transient floc size evolution curve for (View the MathML source during a 30 min simulated flocculation time. This unique transient evolution in the average floc size at (Gm)=70,90, and 150 s-1 was not observed with aggregation–breakup models using only traditional class size PBM approaches or stand-alone QMOM approaches. An explanation for the development of the transient floc size evolution curve at (Gm)=70,90, and 150 s-1 was developed by evaluating the spatial evolution of the moment 4 rate equation (Δm4). The modeling results revealed that the local aggregation and breakup rates varied spatially in the reactor volume. A large fraction of the reactor volume encompassing regions outside the impeller discharge zone was found to promote floc aggregation (85% of the tank volume) with (Δm4) remaining positive throughout the simulated flocculation period. Regions inside the impeller discharge zone and near the impeller (15% of the tank volume) were found to promote floc breakup with (Δm4) remaining negative beyond the initial growth phase. The CFD/QMOM results confirmed that the peak followed by a lower steady-state value in the transient evolution of the average floc size for simulated shear rates (Gm=70,90, 150 s-1) was primarily caused by transport of floc particles in a turbulent heterogeneous environment that promoted spatial heterogeneity not only in the floc size but also in the aggregation and breakup rates.

Notation

AHamaker constant View the MathML source
C1empirical constant in collision efficiency function
C2empirical constant in breakup frequency function
dparticle diameter
d0initial mean particle diameter
d43average floc size (d43=m4/m3)
d43iinitial average floc size
d43nnormalized average floc size (d43n=d43/d43i)
d43nvolume-weighted average floc size for the entire tank
Ffragment distribution function
Glocal characteristic velocity gradient
Gmaverage characteristic velocity gradient
kbsbreakup frequency function
kεturbulent kinetic energy
Liabscissas (or dimensionless size indicator) of the quadrature approximation
mkkth moment of the particle size distribution
m0floc density (moment 0th)
m0volume-weighted average floc density for the entire tank
n(x,d,t)number concentration of flocs
Ndnumber of nodes of the quadrature approximation
Vjvolume of the jth computational cell
wiweight (or dimensionless number concentration) of the quadrature approximation
Greek letters
αcollision efficiency function
βcollision frequency function
Δmkkth moment rate variation (Δmkmkamkb)
Δmkakth moment aggregation contribution
Δmkbkth moment breakage contribution
εpturbulent energy dissipation rate
μdynamic viscosity of fluid
νkinematic viscosity of fluid
ρdensity of fluid
Full-size table

References

    • Hansen and Ottino, 1996
    • J.M. Hansen, J. Ottino
    • Aggregation and cluster size evolution in nonhomogeneous flows

    • Journal of Colloid and Interface Science, 179 (1996), pp. 89–103

    • Kusters, 1991
    • Kusters, K.A., 1991. The influence of turbulence on aggregation of small particles in agitated vessels. Ph.D. Dissertation, University of Eindhoven, Netherlands.
    • Press et al., 1992
    • W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery
    • Numerical Recipes in Fortran; The Art of Scientific Computing

    • (second ed)Cambridge University Press, Cambridge (1992)

    • Prat and Ducoste, 2004
    • Sack, R.A., Donovan, A.F., 1972. Algorithm for gaussian quadrature given modified moments. Numerische Mathematik 18 (5), 465.
    • White, 1991
    • F.M. White
    • Viscous Fluid Flow

    • McGraw Hill, New York, NY (1991)

Corresponding author contact information
Corresponding author. Tel.: +1 919 515 8150; fax: +1 919 515 7908.