6992
Ind. Eng. Chem. Res. 2009, 48, 6992–7003
CFD Modeling of Nucleation, Growth, Aggregation, and Breakage in Continuous Precipitation of Barium Sulfate in a Stirred Tank Jingcai Cheng,† Chao Yang,*,†,‡ Zai-Sha Mao,† and Chengjun Zhao§ Key Laboratory of Green Process and Engineering, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China; Jiangsu Marine Resources DeVelopment Research Institute, Lianyungang 222005, China; and Shijiazhuang Chemical Fiber Corporation, Hebei 050032, China
In this work, the precipitation of barium sulfate (BaSO4) in a continuous stirred tank reactor (CSTR) is modeled. The flow field is obtained through solving the single-phase Reynolds averaged Navier-Stokes equations with a standard single-phase k-ε turbulence model. The population balance equation is solved through the standard method of moments (SMM) and the quadrature method of moments (QMOM) both with and without aggregation and breakage terms. In the cases of precipitation simulation without aggregation and breakage, the results predicted from 2-node QMOM, 3-node QMOM, and SMM are very close. Thus, 2-node QMOM could replace SMM and be well-incorporated into an in-house CFD code to simulate the precipitation in CSTR with acceptable accuracy. The predicted area-averaged crystal size d32 decreases almost linearly with increasing feed concentration, and the deviation from experimental data becomes significant at high feed concentration. Numerical simulation using 2-node QMOM with the Brownian motion and shear-induced aggregation kernels as well as a power-law breakage kernel indicates that the predicted d32 shows good qualitative agreement with experimental results, and the quantitative agreement is achieved when the appropriate breakage rate equation is adopted. 1. Introduction Precipitation or reactive crystallization of sparingly soluble salts is an important industrial operation that is widely used to produce fine and bulk chemicals, pharmaceuticals, biochemicals, catalysts, pigments, ceramics, etc. Precipitation is a very complex process, because it is complicated by several interacting phenomena, and for this reason, it has attracted much attraction. Large-scale production of particulate products is often carried out in a continuous process whereas batch or semibatch mode is usually for small-scale production. As a result of very low solubility of precipitated substances and for economic reasons, high reactant concentrations are commonly applied. Because of the high supersaturation levels resulting from high reactant concentrations, the involved mechanisms of primary nucleation, crystal growth, and aggregation proceed nearly simultaneously and eventually breakage takes place. The precipitation of barium sulfate (BaSO4) has been used as a model reaction and widely studied in the past few decades. Numerous theoretical and experimental studies have been reported in the literature, trying to understand important aspects of the precipitation process.1-8 The previous model-based experimental studies seldom took into account the local flow and turbulence in a precipitator. Nevertheless, computational fluid dynamics (CFD) has been applied to simulate the precipitation process in simple and complex geometries (e.g., stirred tanks) in recent years. This is usually implemented through solving the standard momentum and mass transport equations together with the population balance equation (PBE). Successful CFD simulations of precipitation in a stirred tank have been reported in many studies.9-14 By implementing the standard method of moments (SMM) without considering aggregation and breakage, Jaworski * To whom correspondence should be addressed. Tel.: +86-1062554558. Fax: +86-10-62561822. E-mail:
[email protected]. † Chinese Academy of Sciences. ‡ Jiangsu Marine Resources Development Research Institute. § Shijiazhuang Chemical Fiber Corporation.
and Nienow11 successfully simulated the precipitation of BaSO4 in a continuous stirred tank reactor (CSTR). They found that the residence time and the crystal shape factor significantly affected local supersaturation and volume-averaged crystal size, but in all cases the coefficients of crystal size variation obtained were very close to 1.0, which differed considerably from the experimental results in the range from 0.4 to 0.7.15 The influence of activity coefficient on the crystallization process of BaSO4 ¨ ncu¨l et al.,5 and the in the coaxial pipe mixer was studied by O PBE was modeled using SMM without aggregation and breakage. The precipitation of BaSO4 in a CSTR was also modeled by Wang et al.9 The effects of feeding location, feed concentration, impeller speed, and residence time were investigated, and the discrepancy between predictions of mean particle size and experimental results was found to increase with increasing feed concentrations, which was attributed to the neglect of aggregation. The influence of turbulent mixing on BaSO4 precipitation in a semibatch stirred tank was investigated by Vicum and Mazzotti,13 and the PBE was solved using SMM without aggregation and breakage. It was found that, at low concentrations, the predictions were quantitatively correct, but at high concentration, the overprediction was evident with the CFDbased mixing-precipitation model. Most of the above-discussed studies on the modeling of a precipitation process usually neglected aggregation and breakage, since the PBE was mainly solved using SMM, which has great difficulties in closing nonconstant aggregation and breakage kernels. However, aggregation was clearly evidenced in BaSO4 precipitation.2,16-19 Aggregation and breakage have an important influence on the quality of particulate products, especially for dense particle systems, and should not be neglected in the modeling of precipitation. In recent years, some works were reported on modeling BaSO4 precipitation with particle aggregation incorporated. The PBE is usually solved using SMM with a restrictive simplifying aggregation kernel or quadrature method of moments (QMOM) with a more complex aggregation term. Marchisio et al.16 studied the
10.1021/ie9004282 CCC: $40.75 2009 American Chemical Society Published on Web 07/08/2009
Ind. Eng. Chem. Res., Vol. 48, No. 15, 2009
turbulent precipitation of BaSO4 in a tubular reactor both experimentally and computationally and found that the aggregation and growth in an aqueous solution took place simultaneously with increasing concentrations. The PBE was solved using SMM with aggregation and without breakage. Although only an assumed constant laminar aggregation kernel was considered, the agreement with experimental results was much improved, suggesting that turbulent aggregation and a nonconstant kernel were necessary for better prediction. Gavi et al.20 studied numerically the nanoparticle precipitation of BaSO4 in a confined impinging jet reactor (CIJR). The PBE was solved using QMOM considering the Brownian motion induced aggregation, and the CFD predictions were found to be in good agreement with experimental data. Mixing (macro-, meso-, and micromixing) plays an important role in precipitation, especially precipitation of nanoparticles, and the role of micromixing in precipitation is generally recognized. However, according to Marchisio and Barresi,21 the role of micromixing in fast reactions varies depending on the operating conditions. In stirred tanks, nucleation takes place mainly at the feeding point, and according to Wijers et al.,22 the meso-mixing time criterion should be used for scale-up. Taguchi et al.23 found that macro-mixing was particularly important for stirred tank arrangements with micromixing being of lesser significance, which supports an earlier finding by Fitchett and Tarbell.24 Van-Leeuwen et al.25 showed that the effect of impeller speed on crystal size was limited and could be negligible for impeller speeds larger than 100 rpm. Wong et al.18 revealed that, in the concentration range applied, impeller speed had small effect on crystal size and morphology. Also, mixing sensitive reactions has been modeled for a number of cases via CFD without a micromixing model,26,27 and results showed good agreement with experimental data. In this work, the two feeding nozzles of the CSTR are placed comparatively far from each other in order that the chance of direct contact of both fresh reactants is reduced; thus, the effects of mixing are reduced. In the case of a double-feed-precipitation process, when the feeding nozzles are so close that both fresh reactants contact each other before any significant dilution with the bulk, mixing between fresh reactants becomes more important than a direct dilution with the bulk.12 Up to now, there still exists a lack of understanding concerning the relevance of the micromixing model in CFD simulations of fast reactions. In consideration of these, micromixing effects have not been included in the modeling at the moment. To our best knowledge, the full simulation of a precipitation process in a CSTR has not been reported yet. This work is an attempt to model the continuous precipitation process including nucleation, crystal growth, aggregation, and breakage in a stirred tank, so that the design and scale-up of the precipitation reactor could be better understood. Notwithstanding theoretical analysis and experimental determination are still needed to fully understand some fundamental aspects of precipitation, it is beneficial to find some reliable kinetic expressions for aggregation and breakage and for aggregation efficiency. The investigation of a precipitation process from a modeling perspective is appealing and promising to derive a fully predictive model. This work is organized as follows. First, the population balance equations and moment transformation methods (i.e., SMM and QMOM) are presented, and then aggregation and breakage kernels and aggregation efficiency expression are described. Finally, numerical simulations of a precipitation process in the CSTR both with and without aggregation and
6993
breakage are carried out, and comparison of model predictions with experimental data from the literature is performed. 2. Precipitation Models 2.1. General Population Balance. The Reynolds-averaged form of population balance equation (PBE) expressed in terms of number density function n(L;x,t) is28,29 ∂n(L;x, t) + ∇ · [un(L;x, t)] - ∇ · [Γef∇n(L;x, t)] ) ∂t ∂ [ ( ) ( G L n L;x, t)] + Ba(L;x, t) - Da(L;x, t) + ∂L Bb(L;x, t) - Db(L;x, t)
(1)
where G(L) is the growth rate and Ba(L;x,t), Da(L;x,t), Bb(L;x,t), and Db(L;x,t) represent the birth and death rates due to aggregation and breakage, respectively. In the case of the particle volume and particle length having the relationship υ ∝ L3, more precisely υ ) L3, the above length-based birth and death rates can be rigorously derived from the corresponding volume-based birth and death terms.30 The moments of the particle size distribution (PSD) are defined as follows: mk )
∫
+∞
n(L;x, t)LkdL
0
(2)
Applying the definition of moments to eq 1 results in ∂mk + ∇ · [umk] - ∇ · [Γef∇mk] ) (0)kJ(x, t) + ∂t +∞ kLk-1G(L)n(L;x, t)dL + Bak (x, t) - Dak (x, t) + Bbk (x, t) 0
∫
Dbk (x, t)
(3)
where J(x,t) is the nucleation rate and Bak (x, t) )
1 2
∫
+∞
0
n(λ;x, t)
∫
+∞
0
R(u, λ)β(u, λ)(u3 + λ3)k/3n × (u;x, t) du dλ
Dak (x, t) )
∫
+∞
0
Bbk (x, t) )
∫
+∞
Lkn(L;x, t)
∫
+∞
0
Lk
Dbk (x, t) )
0
(4)
R(L, λ)β(L, λ)n(λ;x, t) dλ dL (5)
∫
a(λ)b(L|λ)n(λ;x, t) dλ dL
(6)
∫
Lka(L)n(L;x, t) dL
(7)
+∞
0 +∞
0
where β(L,λ) is the collision kernel, R(L,λ) is the collision efficiency, a(λ) represents the breakage rate, and b(L|λ) is the fragmentation distribution function. In the case of size-independent growth without aggregation and breakage, the standard moment method (SMM) can be applied to compute the moments directly without requiring additional knowledge of the number density function. Then, the steady-state set of five equations used for computing five crystal size distribution moments, from the zeroth to the fourth, is expressed as28 ∇ · [umk - Γef∇mk] ) 0kJ(x) + jmk-1G (k ) 0-4) (8) For dealing with particle secondary processes such as aggregation and breakage, the quadrature method of moments (QMOM) is a good choice.30-33 The QMOM has been tested
6994
Ind. Eng. Chem. Res., Vol. 48, No. 15, 2009
and compared with other techniques, showing its great potential, and it can solve the PBE with good accuracy by requiring a very small number of equations.30,34 Therefore, this method is very suitable in computationally demanding CFD codes.32,35 The QMOM is based on the following Gaussian quadrature approximation: mk )
The nucleation rate J(x) adopted here is developed by Bałdyga et al.39 and used by other researchers.9,14,16 J(x) ) 2.83 × 1010(∆c)1.775(#m-3s-1), ∆c e 10mol m-3(heterogeneous)
{
2.53 × 10-3(∆c)15.0(#m-3s-1), ∆c > 10mol m-3(homogeneous) (16)
Nd
∫
+∞
n(L;x, t)Lk dL ≈
0
∑wL
k i i
(9)
i)1
where abscissas (Li) and weights (wi) are determined through the product-difference (PD) algorithm.31 By applying eq 9 to eq 3, the steady-state moment equation is ∇ · [umk] - ∇ · [Γef∇mk] ) (0) J(x) + k k
∑
Lk-1 i G(Li)wi
The local values of growth rate G are calculated from the two-step model,1 which includes the surface integration step and the molecular diffusion step. This model has been used in many studies for calculating the growth rate of BaSO4 precipitation.1,9,12,14,16 G ) kr(√cAscBs - √Ksp)2 ) kd(cA - cAs) ) kd(cB - cBs) (17)
+
i
1 2
∑ w ∑ w R β (L i
i
j ij ij
3 i
∑L w ∑wβ + ∑ a bj w - ∑ L a w
+ L3j )k/3 -
j
k i i
i (k) i i i
i
j ij
j
k i i i
(10)
i
where βij ) β(Li,Lj), Rij ) R(Li,Lj), ai ) a(Li) and bj(k) i )
∫
+∞
0
Lkb(L|Li) dL
(11)
2.2. Species Transport Equations. The transport equations for all chemical entities, in molar concentration ci (i ) Ba2+, Cl-, Na2+, SO42-, BaSO4), are applied for the steady-state operation in the following form: ∇ · [Fuci - Γef∇(Fci)] ) Si
(12)
The effective diffusion coefficient, Γef, is computed as a sum of the molecular diffusivity, Γ, and the turbulent diffusivity, Γt ) νt/Sct. The turbulent Schmidt number, Sct, is assumed to be ¨ ncu¨l et al.36 The source 0.7, as in Jaworski and Nienow11 and O term Si is taken to be equal to the specific crystal growth rate, Sg(eq 13), with a minus sign for Ba2+ and SO42-, a plus sign for BaSO4, and zero for other nonreacting ions (Cl- and Na+).36 The specific crystal growth rate Sg is related to the second moment of crystal size distribution, m2, and volumetric crystal shape factor, kV:11 Sg ) (
FBaSO4
Si ) (3m2G)kV F MBaSO4
cBa2+cSO42Ksp
∆c ) (Sa - 1)√Ksp
2kBT (Rc,i + Rc,j)2 3µ Rc,iRc,j
(18)
βFl(L(i), L(j)) ) ψGsh(Rc,i + Rc,j)3
(19)
βBr(L(i), L(j)) ) (13) and
2.3. Nucleation and Growth Kinetics. The supersaturation ratio Sa and supersaturation ∆c are defined as follows: Sa ) γac
where cAs and cBs are reactant concentrations on crystal surface, kr is a kinetic constant that is equal to 5.8 × 10-8 (m/s)(m3/mol)2,1,16 and kd is the mass transfer coefficient with a constant value of 10-7 (m/s)/(m3/mol) considered in this work.14,40 For given values of cA and cB, the growth rate G can be found by solving eq 17 using a Newton-Raphson method. 2.4. Collision Rate. In general, the process of aggregation comprises two steps. First, particles must be brought into close proximity by a transport mechanism, giving rise to a collision. Then, an aggregate will be formed if the net interparticle force is attractive and strong enough to win thermal agitation and hydrodynamic drag in order to make particles adhere or fuse.32,41 There are several mechanisms that can induce relative movements among particles and, hence, lead to collisions. Brownian motion induced collisions (perikinetic coagulation) are the controlling mechanism mainly for submicrometer particles (usually smaller than 1 µm in diameter). For particles with a diameter approximately in the range 1-50 µm, collisions are caused mostly by velocity gradients (orthokinetic coagulation). This shear-induced collision is the prevailing mechanism in particulate processing systems. Collisions from differential sedimentation or inertia (fluid acceleration) will become important for aggregates larger than ∼50 µm.42 The collision kernel functions for Brownian and shear-induced collision are expressed as43
(14) (15)
where γac is the activity coefficient, cBa2+ and cSO2are the ion 4 concentrations of Ba2+ and SO42-, respectively, and Ksp is the solubility product, which is 1.10 × 10-10 kmol2/m6 at 25 °C.26 The activity coefficient is computed according to Bromley’s method,37 which provides accurate data up to 6 M ionic strengths;38 all the ionic strengths considered in this work lie within this range. Secondary nucleation is negligible compared with primary nucleation; thus, only primary nucleation is considered here.
respectively, where L(i) is the characteristic particle size of an aggregate consisting of i primary particles (hereafter referred to as an i-sized aggregate), Rc,i is the collision radius of an i-sized aggregate, kB is the Boltzmann constant, T is the absolute temperature, Gsh is the characteristic velocity gradient (shear rate) of the flow field, and ψ is a numerical constant that depends on the type of flow. For a simple shear flow (e.g., laminar flow), ψ ) 4/3 and Gsh refers to the only nonzero component of velocity gradient tensor. In an isotropic turbulent flow and for particles smaller than the Kolmogorov length scale with negligible inertia, Saffman and Turner43 derived ψ ) 1.29, where the shear rate is related to local energy dissipation rate ε and the kinematic viscosity ν through Gsh )(ε/V)1/2. It is assumed that the two collision mechanisms superimpose and the overall collision rate function is expressed as
Ind. Eng. Chem. Res., Vol. 48, No. 15, 2009
β(L(i), L(j)) ) βBr(L(i), L(j)) + βFl(L(i), L(j))
(20)
The collision radius of an aggregate is proportional to the number of primary particles in the aggregate, i, according to the following relation Rc,i ) R0(i/kg)1/df, i > 3
(21)
where R0 is the radius of the primary particles forming the aggregate, df is the mass fractal dimension, and kg is a constant usually taken to be unity.44,45 According to Lattuada et al.,45 the scaling relation of eq 21 is considered to hold for the aggregates consisting of more than three primary particles, however, it is assumed that it can be extended to all aggregate sizes in order to simplify the calculations. Hence, the solid volume υi for a fractal aggregate is related to its characteristic size Rc,i by Jiang and Logan46 υi ) υ0(Rc,i /R0)df
(22)
where υ0 is the volume of the primary particles. 2.5. Collision Efficiency. Solid particles are subject to hydrodynamic interaction when approaching each other in a sheared fluid, as the viscous fluid layer between them produces resistance forces. Furthermore, a newly formed aggregate will be split into its original two parts if it does not have sufficient time for restructuring after collision, thus resulting in null aggregation. All these are reflected in the collision efficiency, R(L,λ), which is the ratio of the actual aggregation rate to the theoretical collision rate given by eq 20. Collision efficiency is a function of hydrodynamic interactions, interparticle forces, and the structure of aggregates (porosity and permeability). Several models have been developed in the literature44,47,48 to calculate the collision efficiency of porous aggregates (permeable floc models). Since the first step in the formation of clusters from monomers is a doublet formation, the collision efficiency for doublets can be estimated by investigating collisions between impermeable solid particles (impermeable flocs model).35,42,49,50 In a stirred precipitation system, aggregates are compact because they are subject to strong flow shear. Thus, the collision efficiency calculated from the impermeable flocs model is adopted in this study. The collision efficiency is then expressed as35 Rij ) kceFl-0.18, 10 < Fl < 105
(23)
where Fl is the flow number computed by Fl ) 12πµRc,iRc,j(Rc,i + Rc,j)2 /8AHR0
(24)
The prefactor kce ) 0.43 was obtained by fitting the particle numbers with time in the doublet formation stage.35 The Hamaker constant AH ranges approximately between 10-20 and 10-19 J for solid-liquid systems,51 and the value of 5.0 × 10-20 J is used in all simulations. 2.6. Breakage Kernel. The breakage of aggregate can occur by several mechanisms; however, only the hydrodynamic stress induced fragmentation, which is most commonly encountered in precipitation systems, is considered. Several expressions for the breakage rate kernel have been developed by different authors, and they are summarized by Marchisio et al.32 A semitheoretical expression that has found applications to a wide variety of fragmentation phenomena is the power-law breakage kernel:52,53
a(L) ) c1ν ε L
x y γ
6995
(25)
where c1 is a dimensionless empirical constant. By fitting the model to experimental data, Peng and Williams54 found that the exponent γ can be assumed to be between 1 and 3 and γ ) 2 is usually used. The fragment distribution function used in this work is the uniform distribution. Some other fragment distribution functions, such as symmetric fragmentation and erosion, have been used:30
{
2 3 b(L|λ) ) 6L /λ if 0 < L < λ 0 otherwise
(26)
3. Computational Details A continuous flow stirred tank reactor (CSTR) with basic dimensions identical to those assumed by Jaworski and Nienow11 is modeled in this study. Figure 1 shows the schematic diagram of the tank and the impeller. The cylindrical tank has a diameter of T ) 0.27 m with a flat bottom and is filled with a solution up to H ) 0.27 m. It is equipped with four standard baffles of width B ) 0.1T and a standard six-bladed Rushton turbine impeller of diameter D ) T/3, located at C ) T/2 from the bottom. The impeller speed is set to N ) 120 rpm. Two tubes are simulated, and through each a solution containing either BaCl2 or Na2SO4 is fed. Five concentrations are applied, i.e., 0.05, 0.1, 0.2, 0.3, and 0.4 mol/L. The modeling is always carried out with stoichiometric quantities of two reactants being introduced to the precipitator. Thus, the feed rates of the two solutions are always equal and set to 18 mL/s. With the net working volume of 15.4 L, the mean residence time, τ, is then ∼430 s. The feed tubes are located at the free surface, on the opposite sides of the shaft, midway between two neighboring baffles and roughly at the radial location of 5T/12. The reaction mixture exit is located in the tank bottom with the same radial location as that of the feed tubes. All the simulations undertaken in this work were programmed with FORTRAN language. The three-dimensional flow filed in the stirred tank is first obtained through solving the Reynolds averaged Navier-Stokes equations (RANS). Since the solid particles are smaller than 20-30 µm, the solid concentration is low, the solid particles follow closely the liquid, and the influence of the dilute solid phase on the flow field can be neglected, the single-phase RANS and a standard singlephase k-ε turbulence model can be applied.11,14,32 The form of the single-phase RANS in the cylindrical coordinates system and the numerical procedure have been detailed elsewhere.55 Figure 2 shows three predicted velocity components at different radial positions using four different grid numbers. These radial positions are 0.75H from the bottom and midway between two neighboring baffles. It is shown that the simulation results with 36 × 72 × 75 are very close to those with 36 × 72 × 90, and this can also be found in other positions of the tank. Thus, the grid of 36 × 72 × 90 (radial, azimuthal, and axial) is adopted here. The action of impeller is modeled using a modified “inner-outer” iterative procedure.55 Its main advantage is that the calculated flow parameters on the surface of the “inner” and “outer” regions are not averaged in the present procedure, thus preserving the pseudoperiodical turbulent properties. Wall functions are applied to all solid surfaces. After reaching convergence, with the sum of normalized residuals well below 10-4, the velocity and turbulence fields are saved and kept unchanged for the subsequent simulations of reactive precipitation. On the basis of the already known flow field, the concentration distribution of chemical species and lower-order moments
6996
Ind. Eng. Chem. Res., Vol. 48, No. 15, 2009
Figure 1. Sketch of the stirred tank precipitator and impeller.
are obtained through simultaneous solution of five equations for the chemical entity, eq 12, and the set of moment equations. For SMM, five moment equations, eq 8, are used, and for QMOM, the set of moment equations, eq 10, are applied. In the second stage of these simulations, a zero flux at the wall and the free surface are enforced for the variables such as ci and mk. The area-averaged crystal size d32, which is of great interest in precipitation processes, is computed by N
d32 )
N
∑m V/∑m V 3 i
i)1
2 i
(27)
i)1
where Vi is the volume of cell i and N is the total number of grid cells. It should be noted that both SMM and QMOM cannot track the PSD directly, which is also an important criterion for accessing the properties and qualities of crystallization products. However, starting from the first six moments, the PSD can be reconstructed, and first results are promising56 though it is not an easy task. Concerning the shape factor of BaSO4 crystals, kV, it is very hard to determine its value since a number of morphologies have been observed and many different values have been reported or used.2,7,11,16 According to Jones,57 kV ) π/6 for spherical particles and kV > π/6 for other shapes. As mentioned in Section 2.1, the length-based birth and death terms are originally derived from the volume-based PBE with the assumption of particle volume and particle length having the relationship υ ∝ L3.30 Actually, when defining υ ) L3 (i.e., kV) 1.0), the length-based PBE is rigorously deduced from the volume-based PBE.30,35 On the basis of these, kV) 1.0 was adopted in this paper. The first five moments of PSD, from m0 to m4, are computed with SMM. This method is presented and implemented here mainly for the purpose of comparing the results of QMOM without aggregation and breakage with those of SMM. The QMOM is implemented with either two nodes (Nd ) 2) or three nodes (Nd ) 3) when taking no account of aggregation and breakage processes (only the first two terms on the right of
Figure 2. Velocity components at different radial positions using four different grid numbers (radial positions are 0.75H from the bottom and midway between two neighboring baffles): (a) radial component, (b) azimuthal component, and (c) axial component.
eq 16 are used). By applying QMOM with Nd nodes, the first 2Nd moments will be tracked. For example, if Nd ) 3, the first six moments, from m0 to m5, are calculated. In the full simulations of BaSO4 precipitation with QMOM, some parameters in the expressions for collision rate, collision efficiency, and breakage rate should be determined first. The
Ind. Eng. Chem. Res., Vol. 48, No. 15, 2009
fractal dimension of the particles df ) 2.8 is used in this work. It should be noted here that the QMOM in the form presented in this work has assumed Euclidean (nonfractal) aggregates; hence, it will present an inconsistency when introducing a fractal dimension. However, in a stirred precipitation system, particles are subject to strong shear stresses and/or hydrodynamic breakage causes particles restructuring and ripening, which make particles become compact (df is very close to 3). Therefore, the application of the QMOM (assuming df ) 3) and the calculations of the collision rate and collision efficiency (assuming df ) 2.8) seem acceptable. This approach was also validated elsewhere.32,35,58 As already mentioned above, the exponent for the particle size in the power-law breakage kernel in eq 25 can assume values between 1 and 3, and 2 is usually taken.35,54 Therefore, by forcing a parabolic dependence on the particle size (γ ) 2), other exponents, x and y, can be found through a dimensional analysis. The resulting exponent for the turbulent dissipation rate is y ) 1 and that for the kinematic viscosity is x ) -2. As mentioned, c1 is a dimensionless constant and its value can be determined empirically. Hence, a best value for c1 among three adopted values will be recommended when best agreement between CFD simulations and experimental data is obtained. On the basis of the experimental disruption kernel of CaCO3 during precipitation obtained by Wo´jcik and Jones,59 the approximate order of magnitude of c1 or an approximate range of c1 was estimated prior to simulations. This was found to be very helpful to find a roughly reliable breakage rate when starting full simulations. The power-law breakage law used in this paper is semiempirical. Several other theoretical-based breakage laws have been developed and summarized by Marchisio et al.32 Since some parameters (e.g., interparticle force and distance between primary particles, etc.) for these theoretical expressions are hard to obtain or estimate, their applications and practicability have been limited compared with the semitheoretical power-law breakage kernel, although an empirical constant needs to be determined. Regarding the primary particle radius R0, the accurate estimation of its value is very difficult. An assumed value was thus adopted in this work. In precipitation, nuclei are usually of the size of several molecules. Nuclei will grow into crystallites rapidly. Since there are a huge number of nuclei and crystallites, collision and aggregation among these tiny particles would inevitably take place. Therefore, it seems reasonable to assume that the primary particle is of the size of several nuclei. Heineken et al.3 assumed the radius of barium sulfate nuclei being 0.5 × 10-9 m in their model based on the analysis of continuous precipitation process of BaSO4. Actually, we have conducted test simulations using different values of R0, i.e., 2.0 × 10-9, 1.0 × 10-8, and 5.0 × 10-8 m, in the case of aggregation and breakage, and both 2-node and 3-node QMOM were performed. Since the order of magnitude of c1 estimated from the experimental disruption kernel of CaCO3 during precipitation59 was approximately 100, for convenience, c1 ) 1.0 was applied for the test simulations. The feed concentration was 0.1 mol/L. As shown in Figures 6 and 7, the biggest error committed by three adopted R0 values is ∼4% (for m0). Therefore, R0) 2.0 × 10-9 m was adopted in simulations incorporating aggregation and breakage. It has been shown that the QMOM with 2 nodes can work with acceptable accuracy for aggregation and breakage problems.32,60 Also, it can be found in Figure 6 that the predictions from 2-node and 3-node QMOM are very close (the biggest error is
