Simultaneous Quadrature Method of Moments for the Solution of

7798

Ind. Eng. Chem. Res. 2009, 48, 7798–7812

Simultaneous Quadrature Method of Moments for the Solution of Population Balance Equations, Using a Differential Algebraic Equation Framework Jolius Gimbun,†,‡ Zoltan K. Nagy,*,† and Chris D. Rielly† Department of Chemical Engineering, Loughborough UniVersity, Leicestershire, LE11 3TU, U.K., and Faculty of Chemical & Natural Resources Engineering, UniVersiti Malaysia Pahang, Lebuhraya Tun Razak, 26300 Gambang, Pahang, Malaysia

The quadrature method of moments (QMOM) is a recent technique of solving population balance equations for particle dynamics simulation. In this paper, an alternative solution for the QMOM is described and thoroughly tested, which is based on the formulation and simultaneous solution of a semi-explicit differential algebraic equation (DAE) system. The DAE system consists of the ordinary differential equations resulting from the application of the method of moments, as well as a system of nonlinear algebraic equations derived by applying the quadrature theory for the approximation of the moments. It is shown that the proposed approach provides an efficient procedure for evolving the quadrature abscissas and weights from the QMOM. The Jacobian matrix of the DAE system is provided analytically to make the solution more robust. The DAEQMOM method is compared to the well-established method for solving QMOM based on the product difference (PD) algorithm. The numerical results are compared to the analytical solutions in the case of breakage, aggregation, growth, and nucleation mechanisms. Excellent agreements are found on the moment evolution predicted by both methods. However, the DAE-QMOM method is found to be more accurate and robust than the PD-QMOM in some cases. Additionally, the DAE-QMOM is also capable of providing the solution significantly faster than the PD-QMOM method. 1. Introduction The population balance framework has been accepted for some time as the most fundamental approach for modeling particulate, droplet, or bubble dynamics in multiphase systems. A population balance is a powerful tool in evaluating the design of particulate related equipment such as crystallizers (e.g., Randolph and Larson1), particularly when coupled with computational fluid dynamics (CFD) software.2-6 Many solution techniques have been developed to solve the population balance equations, ranging from the simple standard method of moments (MOM), to the more advanced quadrature method of moments (QMOM), method of classes and direct numerical simulation approaches such as high resolution finite volume7-9 and finite element schemes,10 or kinetic Monte Carlo techniques.11 The standard method of moments is one of the most common methods of solving population balance equations (PBEs). In the MOM, the PBE for a spatially homogeneous systems is transformed into a set of ordinary differential equations (ODEs) by multiplying the population balance equation by powers of the characteristics size (Lk) in a length based PBE and integrating it, giving equations in terms of the moments.1 For nonhomogeneous systems, e.g., in the case of the CFD-PBM (population balance model), the MOM reduces the dimensionality of the problem thus making the solution simpler. The MOM is known to be an efficient method to solve the population balance equation, but it suffers from a closure problem for cases involving size-dependent growth, coalescence, or aggregation processes. Therefore, the MOM is not directly applicable for gas-liquid systems, where bubble breakage and coalescence are dominating mechanisms. The method of characteristics (MOCh) is an alternative technique of solving the population balance equation as * To whom correspondence should be addressed. E-mail: [email protected]. † Loughborough University. ‡ Universiti Malaysia Pahang.

demonstrated, e.g., by Lee et al.12 The MOCh for a first-order partial differential equation (PDE) determines the lines, called characteristic lines, along which the PDE reduces to a set of ODEs. Once the ODEs are found, they can be solved easily and transformed into a solution for the original PDE. The MOCh has the advantage of solving the PDE directly for the particle size distribution with any required level of resolution along with the moments; however, this method does not solve the closure problem. Therefore, the MOCh cannot be directly applied to model systems involving breakage and coalescence, such as in the case of gas-liquid dispersions. Other classic solution techniques, such as the method of classes (MOC), are capable of solving the population balance equation for growth, nucleation, breakage, coalescence, and aggregation processes. Although the terms “size class” or “sectional” had been used earlier by authors such as Gillette,13 Sutugin and Fuchs,14 Tolfo,15 and Gelbard et al.,16 the method only became readily available for implementation after a detailed derivation of the MOC for particle growth, breakage, and aggregation had been introduced by Marchal et al.17 The MOC has been applied to many problems involving breakage and coalescence or aggregation processes for gas-liquid dispersions,18,19 solid-liquid dispersions,20,21 and liquid-liquid dispersions.22-24 This method requires the whole particle size distribution (PSD) to be resolved into discrete size classes in order to obtain an accurate solution. The number of discretized equations increases with the number of size classes employed, and hence the solution for the MOC can be numerically expensive to obtain,25 especially for multiple coordinate systems, e.g., when coupled with three-dimensional CFD simulations. Most CFD-PBM studies26,27 only considered the Sauter mean diameter for calculations of the two-phase flow in CFD due to limitation in computational resources. However, solution with multiple particle/bubble size may become feasible in future as the computational power continues to increase. Furthermore, the PBE solution obtained from MOC with a small number of

10.1021/ie900548s CCC: $40.75  2009 American Chemical Society Published on Web 07/09/2009

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009 21

classes (e.g., Hounslow et al. ) introduces up to 5% error in the moments prediction, which may affect prediction of the particle size distribution. McGraw28 proposed an attractive method for solving the population balance equation, namely the quadrature method of moments (QMOM), which utilizes the quadrature theory to avoid the closure problem with standard MOM simulations. The QMOM by McGraw28 is based on the product difference (PD) algorithm of Gordon.29 Application of the QMOM has been extended into aggregation, coagulation, and breakage systems by Wright et al.,30 Rosner and Pyykonen,31 Fan et al.,32 and Marchisio et al.33 However, the PD algorithm is not always the best approach for computing the quadrature points from the moments of the particle size distribution (e.g., Lambin and Gaspard34) because, for a larger number of moments, the method is sensitive to small errors (e.g., Gautschi35). Therefore, the applicability of QMOM is limited to no more than six quadrature points28 and even fewer for more complex cases such as diffusion-controlled growth with secondary nucleation. Later, McGraw and Wright36 proposed a new method, namely the Jacobian matrix transformation (JMT) method, which provides an alternative approach to the PD algorithm. This method is also based on the projection of the solution into the power moments, yielding a set of equations for the evolution of the quadrature rule. The JMT method, however, still suffers from the problem of ill conditioning, requiring the number of moments to be small.37 The occurrence of singularities in the QMOM with PD or JMT is the consequence of using the power moments, thus making the calculation of population balance unfeasible for a higher number of moments and hence a higher number of quadrature points. Fan et al.32 proposed the direct quadrature method of moments (DQMOM) as another alternative to the PD algorithm. The DQMOM is based on the idea of tracking directly the variables appearing in the quadrature approximation, by solving the convection equation for quadrature points and weights directly. The authors claim that some shortcomings of the PD-QMOM are avoided by solving the weights and abscissas directly especially when concerning a CFD-PBM problem. As Marchisio and Fox38 pointed out, the DQMOM can still be ill-conditioned for a higher number of moments because of the intrinsic problem of finding the roots of high order polynomials.39 Su et al.40 suggested an improvement to DQMOM by introducing an adaptive factor in the moment equations which they claim to be able to improve the DQMOM accuracy for aggregation processes without adding a significant computational expense (only 4% higher with the adaptive factor). The adaptive DQMOM does not increase the original DQMOM accuracy for other systems such as growth and breakage, but it provides about 2% faster solution. Recently, Alopaeus et al.41 proposed another solution, namely the fixed quadrature method of moments (FQMOM). This method exploits the zeros of orthogonal polynomials to find the quadrature points. Unlike the PD-QMOM, which is restricted to an even number of moments (twice the number of the quadrature points), the FQMOM can be applied to any number of moments. The disadvantage of the FQMOM is that an optimized arbitrary constant must be established before it can be applied accurately for a specific problem. Finding that optimized constant might be easy for simple cases where analytical solution is available, but can be difficult for more complex systems. Most of the currently available QMOM solutions are restricted to using a small number of moments. Although a small number of moments may be sufficient to describe the particulate

7799

dynamics in simple cases, a larger number of moments might be required in other calculations, especially when an inversion technique is applied to discover the particle size distribution. Diemer and Ehrman,42 for example, needed at least 10 moments to obtain acceptable reconstructions of the complete particle size distribution. Therefore, a better solution technique for the population balance that is capable of solving for higher order moments is needed. This paper describes the possibility of solving the PBE via QMOM, without resorting to the PD algorithm. An alternative solution technique for the QMOM is thoroughly tested, based on the simultaneous solution of the moment equations and quadrature approximation as a semiexplicit differential algebraic equation (DAE) system. A similar DAE based formulation of the QMOM was also studied previously by Grosch et al.,43 although their work focused mainly on crystallization systems with growth and nucleation mechanisms. The result from the DAE-QMOM method is compared with the solution from the PD algorithm, as well as analytical solutions, for various mechanisms including growth, nucleation, aggregation, and breakage. 2. Population Balance Equation The dynamic population balance equation for a closed homogeneous system can be written with diameter as the internal coordinate as33,44

where β, a, G, B, b, L0, and δ are the aggregation kernel, breakage kernel, growth rate, nucleation rate, daughter particle size distribution, size of the nuclei, and Dirac delta function, respectively, whereas both L and λ are the particle characteristic length. The PBE in eq 1 can be further simplified using a moment transformation, where the kth moment of the distribution, µk, is given by µk )

∫

∞

0

n(L)Lk dL

(2)

The lower order moments (i.e., zeroth to third) are related to the physical description of the particle size distribution; i.e., µ0 is related to the total number of particles, µ1 is related to the total particle diameter, µ2 is related to the total particle surface area, and µ3 is related to the total particle volume. After the moment transformation, the PBE of eq 1 is represented by a set of ODEs in terms of the moments:45

The moment equations represented by eq 3 are solvable for growth (but not size-dependent growth) and nucleation problems using the standard method of moment technique; however, it is not possible to solve the breakage and coalescence terms due to the closure problem, since the integrations cannot be written in term of the moments. Therefore, eq 3 needs to be transformed

7800

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009

again into a quadrature method of moments formulation to eliminate the closure problem. The essence of the quadrature closure is to consider the number density n(L) as a general weight function and to approximate the integrals that appear during the transformation of the PBE to moment equations in terms of a set of abscissas and weights. The QMOM employs a quadrature approximation:28 µk )

∫

∞

0

For the case of N quadrature points, the DAE-QMOM method requires 2N moment ODEs and another 2N algebraic equations for the quadrature approximation. The first 2N moments are generated using eq 7, i.e. dµ0 )0 dt dµ1 G0 G0 G0 ) w + w + ... + w dt L1 1 L2 2 LN N l dµ2N-1 ) (2N - 1) dt

N

n(L)Lk dL ≈

∑wL , k i i

for k ) 0, 1, 2

(4)

i)1

where wi are the weights, Li are the abscissas, and N is the number of quadrature points. This quadrature approximation is exact if the function in eqs 4 are polynomials up to the order of 2N - 1. After applying the quadrature rule, the moment transformed PBE can be written as41

(8)

The other 2N algebraic equations are generated using the quadrature approximation given by eqs 4. 0 ) w1 + w2 + ... + wN - µ0 0 ) w1L1 + w2L2 + ... + wNLN - µ1 l 0 ) w1L12N-1 + w2L22N-1 + ... + wNLN2N-1 - µ2N-1 (9)

Now the closure problem has been eliminated, and hence the PBE in eq 5 is solvable by the means of the quadrature method of moments by following the evolution of wi and Li, as well as µk. The moments are nonlinearly related to the weights and abscissas by eqs 4.

Finally, the problem is represented by 4N equations containing 4N unknowns made up from 2N moments, N weights, and N abscissas. Equations 8 and 9 together represent a semi-explicit DAE system which can be solved using standard DAE solution techniques and software. In semi-explicit form the generic DAE can be written as

3. Numerical Techniques for QMOM The QMOM calculations require integration of the ODEs eq 5 for k ) 0, ... 2N - 1, alongside the nonlinear algebraic eqs 4. These equations may be numerically solved simultaneously as a set of coupled DAEs. However, most of the QMOM solutions described in the literature have been performed using the product difference algorithm of Gordon28 to solve eqs 4. Solutions for the population balance equation via the PD-QMOM are well established for growth, nucleation, breakage, and aggregation/ coalescence (e.g., McGraw;28 Marchisio et al.33). A variety of alternative methods are available to solve the QMOM, including FQMOM, JMT, and DQMOM (see section 1); however, those techniques are not widely used. Therefore, the solution from the DAE-QMOM is only compared with the most common technique for the QMOM, based on the PD algorithm. 3.1. DAE-QMOM. The DAE method is an attractive alternative method for solving the QMOM since it arises from the natural mathematical formulation of the QMOM approximation problem The ODE eq 5 is generated from the moment equations, whereas the algebraic eqs 4 are obtained from the quadrature rule. The DAEQMOM method is illustrated in this section using a diffusioncontrolled (size-dependent) growth model assuming the nucleation, breakage, and coalescence or aggregation rates are zero. In the diffusion-controlled growth model the growth rate is given by G0 G(L) ) L

(6)

where G0 is a constant. The population balance equation written in the form of the quadrature approximation is given by N dµk )k G(Li(t))Li(t)k-1wi(t) dt i)1

∑

(7)

f x˙(t) ) f(x b, b, z t)

(10)

0 ) g(x b, b, z t)

(11)

where b x(t) contains the differential variables and b(t) z contains the algebraic variables. The ODE eq 10 for b x(t) depends on additional algebraic variables b(t), z and the solution is forced to satisfy the algebraic constraints given by eq 11.46 In the case z ) [w1, w2, ..., wN, L1, of QMOM b x ) [µ0, µ1, ..., µ2N-1] and b L2, ..., LN]. The semi-explicit system eqs 10 and 11 can be written as b (t, b M y )y˙ ) b F(t, b y)

(12)

b is a mass matrix, and where M y)

()

F)

()

x z

and f g

The mass matrix for the DAE in eq 10 and 11 is given by

( )

b ) I 0 M 0 0

(13)

where I is the identity matrix of dimension nx × nx and nx is the number of differential states. MATLAB is an efficient software tool for numerical computation and can be used readily for the solution of index 1 DAE problems of the type shown in eq 12. The differential index of DAE is defined as the number of differentiations needed to transform the DAE into an explicit ODE. MATLAB can solve the DAE of index 1 via the ode15s solver or ode23t solver.47 In this work, the ode15s solver is employed to solve the DAE-QMOM. The ode15s is based on a variant of the backward differentiation formula (BDF) called

Ind. Eng. Chem. Res., Vol. 48, No. 16, 2009 48

the numerical differentiation formula (NDF). According to Shampine and Reichelt,47 many tactics adopted for ode15s resemble those found in well-known codes such as LSODE49 and VODE.50 It was developed to integrate stiff ODEs and DAEs of the form of eq 12. The simplified Newton method is also implemented in ode15s to perform a correction to the current iteration (to satisfy the algebraic constraints), thus solving the DAE system. The requirement for the solution of the DAE system is the existence of a feasible initial condition for y0, that is, if there is a vector such that M(t0,y0)y˙ ) F(t0,y0). The ode15s is implemented to detect automatically the DAE system and then to perform an automatic computation of consistent initial conditions for a robust computation. The DAE system might lead to instability of the solution, especially when a more complicated function is considered, e.g., involving the breakage and coalescence kernels. To overcome such a problem, it is necessary to provide the analytical Jacobian matrix. The Jacobian of the DAE system is defined as follows:

[ ]

∂f1 ∂f1 ∂f1 ··· ∂y1 ∂y2 ∂y4N ∂f2 ∂f2 ∂f2 ··· ∂y4N J ) ∂y1 ∂y2 l l l l ∂f4N ∂f4N ∂f4N ··· ∂y1 ∂y2 ∂y4N

(14)

where f1 to f2N are the right-hand sides (RHS) of the differential moment equations, f2N+1 to f4N are the RHS of the algebraic equations, y1 to y2N represent the moments, y2N+1 to y3N are the weights, and y3N+1 to y4N are the abscissas. The Jacobian matrix for a given equation system can be generated using the symbolic computation tool in MATLAB. If for example, for N ) 2, for the diffusion-controlled growth only case, then the Jacobian matrix, J, becomes

7801

starting from the moments. The components in the first column of matrix P are Pi,1 ) δi1,

i ) 1, ..., 2N + 1

(16)

where δi1 is the Kronecker delta. The components in the second column of P are Pi,2 ) (-1)i-1µi-1,

i ) 1, ..., 2N + 1

(17)

The remaining components can be computed from an iterative equation as follows: Pi,j ) P1,j-1Pi+1,j-2 - P1,j-2Pi+1,j-1, j ) 3, ..., 2N + 1 and i ) 1, ..., 2N + 2 - j

[

For example, when N ) 2, the matrix P becomes 1 0 P) 0 0 0

µ0 -µ1 µ2 -µ3 0

µ1 -µ2 µ3 0 0

µ2 - µ12 -µ3 + µ2µ1 0 0 0

µ3µ1 - µ22 0 0 0 0

]

(18)

(19)

The coefficients of the continued fraction (Ri) are generated by setting the first element equal to zero (R1 ) 0) and computing the others according to the following recursive relationship: Ri )

P1,i+1 , P1,iP1,i-1

i ) 2, ..., 2N

(20)

A symmetric tridiagonal matrix is then obtained from sums and products of Ri: ai ) R2i + R2i-1,

i ) 1, ..., 2N - 1

(21)

and bi ) - √R2i+1R2i-1,

i ) 1, ..., 2N - 2

(22)

where ai and bi are the diagonal and subdiagonal of the Jacobian matrix. Once the tridiagonal matrix is determined, the abscissas and weights can be found by calculating the eigenvalues and eigenvectors of the tridiagonal matrix. The eigenvalues represent the abscissas, and the weights can be found from wj ) µ0υj12 The Jacobian matrix can be divided into three major blocks. The first block, J(1,1) to J(2N,2N), contains the derivatives of the right-hand sides of the moment equations with respect to the moments. This particular block will remain unchanged, unless the RHS includes the moments. The second block, J(1,2N+1) to J(2N,4N), contains the derivatives of the RHS of the differential equations, with respect to the weights and abscissas. Therefore this block always changes whenever the mechanism change. The third block, J(2N+1,1) to J(4N,4N), is derived from the algebraic equations of the quadrature approximation and is the independent of the mechanisms used in the PBE. 3.2. PD-QMOM. The solution of the QMOM via the PD algorithm was first introduced by McGraw28 and has been further developed for aggregation and breakage by Marchisio et al.33 The PD algorithm is employed to calculate the weights and the abscissas from the moments in the QMOM. The first step is the construction of a matrix P with components Pi,j

(23)

where υj1 is the first component of the jth eigenvector. It should be noted that the elements of the eigenvectors should be normalized such that the norm of each eigenvector is 1.0. The eigenvalues and eigenvectors can be calculated numerically using the eig function in MATLAB. 4. Validation of the QMOM Algorithms The performance of the DAE-QMOM is compared with the well-established PD-QMOM method and with a number of analytical solutions. Six cases were tested, namely (i) diffusioncontrolled growth (see section 3.1), (ii) constant growth and primary nucleation, (iii) power-law growth, (iv) breakage, (v) constant aggregation kernel, and (vi) sum aggregation kernel. In each case, except for the power-law growth, the initial size distribution was n0(L) ) 3L2

N0 -L3/V0 e V0

(24)

7802


solver (ode15s) is employed for both the DAE-QMOM and the PD-QMOM. All solutions in this paper were computed using three quadrature points (six moments). The error of the predicted moment is calculated as follows: % error )

Figure 1. Discontinuity of eq 29 at t ) 1 s. Upper right is the enlarged view of the discontinuity.

where N0 ) 1 m-3 and V0 ) 1 m3 or 0.001 m3. The particle size distribution in eq 24 is chosen because it has an analytical solution available from the literature. Unless otherwise stated, in ode15s the relative tolerance is set at 10-12 and the absolute tolerance is set at 10-10 for all cases tested in this work for both the DAE-QMOM and the PD-QMOM. The relative tolerance (RelTol) represents the error that applies to all components of the solution vector, whereas the absolute tolerance (AbsTol) applies to individual components of the solution vector. The error in each integration step err(i) for any state variable x(i) in the solution vectors satisfies the condition err(i) e max(RelTol · abs(x(i)),AbsTol(i)). A similar ODE/DAE

µanalytical - µcalculated × 100 µanalytical

(25)

4.1. Growth and Nucleation Examples. Case 1: DiffusionControlled Growth. The moment equations for the diffusioncontrolled growth have been described in section 3.1 (eqs 6-8). The initial distribution is given by eq 24 with N0 ) 1 m-3, V0 ) 1 m3, and G0 ) 0.01 m2/s. The zeroth moment is the number density of particles per unit volume and can be obtained by integrating the initial distribution in eq 24 from zero to infinity. For the growth only case, µ0 will remain constant. Exact solutions are available for zeroth and even moments only, given by McGraw28 as µ0 ) N0 ) constant

(26)

µ2 ) 2G0µ0t + µ2(0)

(27)

µ4 ) 4G02µ0t2 + 4G0µ2(0)t + µ4(0)

(28)

An analytical size distribution at any time t can be obtained from41

Figure 2. Error in the moment evolution for diffusion-controlled growth of DAE-QMOM (continuous lines) and PD-QMOM (circles).

n(L) )

(

(L - 2G0t) 3N0 L√L2 - 2G0t exp V0 V0 2

3/2

)


7803

(29)

which then can be integrated numerically using the function quadl in MATLAB to obtain the analytical solution for the odd moments. An analytical integration is not possible for eq 29 due to the presence of a discontinuity when L2 - 2G0t < 0. Figure 1 shows the plot of eq 29 at t ) 1 s. Equation 29 returns an imaginary value at L < (2G0t)1/2, which results in a discontinuity. A numerical integration via the quadl function at a very small value of tolerance (function tolerance 10-12) was employed to calculate the odd moments. The lower limit of the quadl integration was set at L ) (2G0t)1/2, the point where the discontinuity starts. The moments for the diffusion-controlled growth were predicted remarkably well by both the DAE and PD methods, with a small error (mostly less than 10-8, see Figure 2) for the zeroth and even moments where the exact analytical solution is available. There is no significant difference in the accuracy of the even moments predicted by the DAE-QMOM and PDQMOM in this particular case, except for a very small difference in the fourth moment. The error for the odd moments is significantly larger than that for the even moments, due to the fact that these are obtained from numerical integration of eq 29 and not from the exact solution. Evolutions of the first six moments are shown in Figure 3; each moment is normalized with their initial value (µk(t)/µk(0)). In this case where only growth is considered, the particle number density (zeroth moment) remains unchanged. The DAE-QMOM and PD-QMOM were also tested to solve the diffusion-controlled growth at a looser tolerance for ode15s (absolute tolerance ) relative tolerance ) 10-6) to further evaluate the robustness and accuracy of the proposed solution. Figure 4 shows the relative error of the moment evolution from both DAE-QMOM and PD-QMOM. For the case of a looser tolerance, prediction of the moment evolution by DAE-QMOM and PD-QMOM did not change much except for the fourth moment, which was already showing some difference even for the tight tolerance. The error in the fourth moment shows that the accuracy of moment prediction in this case depends only on the tolerance setting because the magnitude of errors reflects the absolute tolerance applied. Figure 5 shows the evolution of the weights from DAE-QMOM and PD-QMOM at the looser tolerance setting. For the case of the growth only problem, the weights should be constant because there is no change in the particle number density. The individual weights are conserved perfectly by DAE-QMOM at looser tolerance; however, a significant error is observed in the case of the PD-QMOM method.. There is no problem of conserving the weights at tight tolerance setting for both DAE-QMOM and PD-QMOM, and therefore the results for weight evolution at tight tolerance are not shown. Even for the looser tolerance, the total weights are well conserved, even if the individual weights are not. In some cases when the individual weights are not conserved, some of the weights actually become zero and thus they no longer contribute to the PD solution. In such cases, the PD algorithm has a singularity problem and is therefore no longer capable of solving the QMOM. The conservation of individual weights in QMOM is important to ensure that the solution is robust and accurate. The weights are a crude approximation of the actual particle size distribution. In the case of a normal distribution the abscissa of the middle weight is similar to the mean diameter and the other two weights have a similar value (see Figure 6). The individual weights in the QMOM methods represent a

Figure 3. Moment evolution for diffusion-controlled growth from DAEQMOM (circles) and analytical solution (continuous lines).

unique solution and therefore should be conserved perfectly in the growth only problem. This numerical evaluation has clearly demonstrated the robustness and accuracy of the DAE method compared to the PD algorithm in solving the population balance equations. Case 2: Constant Growth and Primary Nucleation. The next level of complication is to include nucleation, so that the number density of particle (µ0) increases with time. The evolution of the zeroth moment for this case is given by dµ0 ) B0 dt

(30)

This particular case of constant growth rate, G, does not involve any closure problem and can be solved easily using the standard method of moments. The moments equation for this case for MOM is given as follows: dµk dµ0 ) B0 dt dt ) kGµk-1,

kg1

(31)

where in this case study B0 ) 0.1 m-3s-1. The initial distribution is again given by eq 24 with N0 ) 1 m-3, V0 ) 1 m3, and G ) 0.01 m/s. The number of moments is not restricted for MOM and as many can be generated as are needed. A similar stiff ODE solver, ode15s from MATLAB, was employed to solve the MOM for a proper comparison with the QMOM method. Furthermore, the analytical solution for this case can be derived easily and the first three moments are given as follows: µ0(t) ) B0t + µ0(0)µ1(t) t2 ) GB0 + Gµ0(0)t + µ1(0)µ2(t) 2 t3 t2 2 ) 2G B0 + 2G2µ0(0) + 2Gµ1(0)t + µ2(0) 3 2 (32) The moment evolutions calculated using the two different QMOM algorithms and the MOM are compared to the analytical solutions of eqs 32. Figure 7 shows very small errors in the moment evolutions calculated using the two different QMOM algorithms and MOM. In this particular case the errors of the PD-QMOM and DAE-QMOM appear to be of similar magni-

7804


Figure 4. Error of moment evolution for diffusion-controlled growth from DAE-QMOM (continuous lines) and PD-QMOM (circles) to the exact solution.

tude. Prediction of the MOM was also found to be at an accuracy similar to that for both QMOM solutions. Case 3: Power-Law Growth. The power-law model is commonly applied for crystallization growth and mass transfer related problems. In the case evaluated in this work, a powerlaw growth model with the following relations without the nucleation term is considered. G(L) ) bLp

(33)

0.81-p - 0.31-p 100(1 - p)

(34)

where b)

and p ) 0.5 or -0.5. The case studied here is similar to that evaluated by Alopaeus et al.51 However, the initial distribution j) is assumed to be a normal distribution with a mean size of L 0.5 m and a standard deviation σ ) 0.05 m instead of eq 24 in Alopaeus et al.’s51 work. An analytical size distribution at any time t for this particular case of power-law growth is given by n(L) )

(

)

(1 - p)bt p/(1-p) × L1-p √2πσ j )2 ((L1-p - (1 - p)bt)1/(1-p) - L exp 2 2σ 1

1-

(

)

(35)

j and σ are the initial mean diameter and standard where L deviation, respectively. A detailed derivation of the analytical

solution for power-law growth is given in the Appendix. Similar results were reported by Alopaeus et al.51 for eq 24. The final analytical size distribution in eq 35 cannot be integrated analytically to obtain the exact moments due to a discontinuity in the function. However, eq 35 can be integrated numerically in MATLAB using the quadl solver to give an approximate value of the moments for comparison with the QMOM solution. The numerical integration employed a tight absolute tolerance of 10-12. The relative tolerance used in the ODE integrations is set to 10-12 for both cases (p ) 0.5 and p ) -0.5); however, it was necessary to reduce the absolute tolerance to 10-8 because the ODE integrator failed to solve the problem at a tighter tolerance. The comparison of relative errors in the moment evolutions obtained using the DAE-QMOM and PD-QMOM methods for power-law growth with p ) 0.5 and p ) -0.5 are presented in Figures 8 and 9, respectively. In both cases the DAE-QMOM described the moment evolutions with a significantly lower relative error in most of the moments. For a negative value of p, the smaller particles grow faster than the bigger ones causing the particle size distribution to evolve as a shock wave,51 resulting in a major difficulty in conserving the individual weights (no nucleation in this case) generating the multiplicity problem in the solution from the PD algorithm. Due to this singularity problem the correct solution cannot be found even for a tighter tolerance (see Figure 10). As mentioned in the section Case 1: Diffusion-Controlled Growth, the weights for a


Figure 5. Evolution of weights at looser tolerance: DAE-QMOM (continuous lines) and PD-QMOM (dashed lines).

growth only problem should be perfectly conserved since there is no change in the total number of particles in the system. The tests on using loose (relative tolerance of 10-6 and absolute tolerance of 10-6) and tight (relative tolerance of 10-12 and absolute tolerance of 10-8) tolerance settings were carried out to assess the robustness of the new solution method. Figures 10 and 11 show that the PD-QMOM fails to conserve the weights at tight tolerance and is even worse for the looser tolerance setting. In contrast, the DAE-QMOM conserves the weights at both looser and tight tolerance settings. For p ) -0.5 shown in Figure 10 the error in the PD-QMOM solution becomes larger as time increases. The moments are computed as a product between the weights and abscissas; thus any problem in predicting the weights could translate

7805

Figure 6. Weights for normal distribution function for three quadrature points (Lj ) 0.5 and σ ) 0.05).

into an error in the moment calculation. The moments may be predicted fairly well even if the weights are not predicted correctly because the PD-QMOM counteracts the change in the weights by altering the abscissas and forcing the moments to obey the quadrature rule. However, ill conditioning may occur when one or more weights become zero, because they no longer contribute to the PD solution, and when this happens the PD-QMOM stops providing any solution. Thus it is important to obtain a correct prediction of the weights to maintain the robustness of the solution process. 4.2. Aggregation Kernel Examples. The next stage of complication is aggregation (coalescence), which leads to changes in both the particle size distribution and the number

Figure 7. Comparison of error for the first six moments for primary nucleation and constant growth of DAE-QMOM (continuous lines), PD-QMOM (circles), and MOM (triangles) compared with the analytical solution of eq 32.

7806


Figure 8. Relative errors of the first six moments for power-law growth at p ) 0.5 of DAE-QMOM (continuous lines) and PD-QMOM (circles) to the analytical solution.

Figure 9. Relative errors of the first six moments for power-law growth at p ) -0.5 of DAE-QMOM (continuous lines) and PD-QMOM (circles) compared to the analytical solution.

density of particles. Two different aggregation kernels were tested, namely the constant and sum kernel. The initial distribu-

tion employed for both aggregation kernels was given by eq 24 with N0 ) 1 m-3 and V0 ) 0.001 m3.


7807

Figure 10. Evolution of the weights for power-law growth at p ) -0.5 for DAE (continuous lines) and PD (dashed lines). (A) Looser tolerance (relative tolerance of 10-6 and absolute tolerance of 10-6); (B) tighter tolerance (relative tolerance of 10-12 and absolute tolerance of 10-8).

Figure 11. Evolution of the weights for power-law growth at p ) 0.5 for DAE (continuous line) and PD (dashed line). (A) Looser tolerance (relative tolerance of 10-6 and absolute tolerance of 10-6); (B) tighter tolerance (relative tolerance of 10-12 and absolute tolerance of 10-8).

Case 4: Constant Aggregation Kernel. The constant aggregation kernel is given by

Case 5: Sum Aggregation Kernel. The sum aggregation kernel53,54 is given by

β(L, λ) ) 1

β(L, λ) ) L3 + λ3

(36)

and in this case the breakage, growth, and nucleation terms are set to zero. The analytical solutions for the moments for this case can be obtained from the following equation:52

(

µk ) µk0

2 2 + N0β(L, λ)t

)

1-((k-1)/3)

(37)

In Figure 12, the numerical solutions using the DAE-QMOM and PD-QMOM methods are compared to the analytical solution for the first six moments. In this case the errors of both DAEQMOM and PD-QMOM are similar. The third moment which is related to the total particle volume appears to be conserved very well in this case, as shown in Figure 13. Each of the moments in Figure 13 are normalized with their initial value (µk(t)/µk(0)). For aggregation only problems, the zeroth moment (representing the number of bubble per unit volume) decreases with time. The first moment, which represents the total particle diameter, and the second moment, which is related to the total particle surface area, also decrease due to aggregation. The total particle volume is represented by the third moment. For problems involving only aggregation and/or breakage, the third moment should be conserved, because there is no addition of new volume into the system. This is however not true for problems involving growth and/or nucleation, where the total volume of particle is changing.

(38)

For this case the analytical solution is only available for the zeroth moment, which is given by µ0 ) N0e-N0V0t

(39)

The third moment is also conserved in this case. The percentage error for the zeroth moment is significantly smaller in the case of the DAE-QMOM method than for the PD-QMOM as shown in Figure 14. Results from this example showed very little difference in prediction accuracy for both PD-QMOM and DAE-QMOM for aggregation kernel. Although there is a very small difference in the prediction for the zeroth moment, the difference is negligible (less than 10-6) considering the error for other moments is larger than 10-3 (except for the µ3, which is perfectly conserved). 4.3. Breakage Kernel. In this case the breakage kernel was assumed to be proportional to the particle volume: g(L) ) L3

(40)

The uniform breakage function is assumed, which gives uniform probability of all fragment sizes: b(L, λ) )

6L2 λ3

(41)

7808


Figure 12. Comparison of errors for the first six moments for constant aggregation kernel using DAE-QMOM (continuous lines) and PD-QMOM (circles) compared to the analytical solution of eqs 4-37. Table 1. Comparison of the CPU Time in seconds for Different Solution Methods for the Population Balance Equation PBE solution

diffusion-controlled growth

constant growth and primary nucleation

0.781 1.375

0.516 1.703 0.562

DAE PD MOM a

power-law growth p ) 0.5

p ) -0.5

aggregationa

breakageb

0.406 0.435

1.630 fails

0.609 3.188

1.875 3.547

Constant. b Proportional to volume.

Table 2. Comparison of CPU Time of DAE and PD at Various Number of Quadrature Points (CPU Speed 3.8 GHz, MATLAB 2006a) N

tDAE (s) tPD (s)

2

3

4

5

0.38 1.72

0.49 4.74

1.21 8.15

2.81 18.26

The daughter particle distribution function, b(L,λ), determines the number and size of the daughter particle, L, after the breakage event of particle size λ. For the uniform breakage function, binary breakage with similar particle sizes is assumed. In this case the aggregation, growth, and nucleation terms are set to zero. The initial particle distribution is given by eq 24 with N0 ) 1 m-3 and V0 ) 1 m3. In this case the analytical solution at any time t is given by (e.g., Hounslow et al.21) n(L) ) 3L2(1 + t)2e-L (1+t) 3

(42)

which then can be integrated analytically using the int function in MATLAB to obtain the exact moments. The analytical integration was performed with the limits from zero to infinity. As shown in Figure 15, the errors from both the DAE-QMOM

and PD-QMOM in this case are small, on the order of ∼0.1% or less, and similar. The third moment is also well conserved, as shown in Figure 16. The profiles for the evolution of the normalized moments for the breakage kernel (Figure 16) are the reverse of those for the aggregation kernel in Figure 13. For this breakage only problem, the first three moments increase steadily as the large particles break into smaller particles. It is clear from the test for breakage only that both the DAE-QMOM and PD-QMOM methods yielded similar accuracy. 5. Comparison of Computational Effort for PD-QMOM and DAE-QMOM Apart from the prediction accuracy, the performance of the DAE-QMOM, in terms of CPU times is also compared to that of the existing PD-QMOM method. The test was carried out by employing a similar ODE solver (ode15s) and tolerance (similar to those applied for cases 1-6) for both PD and DAE algorithms. The comparison of the CPU times for the various PBE solutions is shown in Table 1. DAE-QMOM is more than 3 times faster than the existing solution via PD-QMOM for constant aggregation and growth with the primary nucleation problems. For the breakage and diffusion-controlled growth


7809

a

Table 3. Comparison of Average Error of Moment Calculated by DAE and PD at Various Numbers of Quadrature Points DAE

µ0 µ1 µ2 µ3 µ4 µ5 µ6 µ7 µ8 µ9 a

PD

N)5

N)4

N)3

N)2

N)2

N)3

N)4

N)5

6.7 × 10-7 3.6 × 10-4 1.0 × 10-4 0.0 4.9 × 10-5 1.6 × 10-5 1.1 × 10-15 1.1 × 10-4 3.5 × 10-5 5.5 × 10-11

3.0 × 10-7 1.5 × 10-3 1.2 × 10-3 0.0 6.1 × 10-4 1.8 × 10-3 2.4 × 10-13 1.1 × 10-3 5.7 × 10-2 0.29

6.8 × 10-7 4.0 × 10-3 8.2 × 10-3 0.0 3.2 × 10-2 9.3 × 10-3 0.47 1.8 4.5 8.8

1.7 × 10-6 0.5 0.84 0.0 3.6 10 20 31 43 55

1.1 × 10-6 0.5 0.84 0.0 3.6 10 20 31 43 55

1.0 × 10-6 4.0 × 10-3 8.2 × 10-3 0.0 3.2 × 10-2 9.3 × 10-3 0.47 1.8 4.5 8.8

2.7 × 10-7 1.5 × 10-3 1.2 × 10-3 0.0 6.1 × 10-4 1.8 × 10-3 1.9 × 10-13 1.1 × 10-3 5.7 × 10-2 0.29

3.4 × 10-7 3.5 × 10-4 1.0 × 10-4 0.0 4.9 × 10-5 1.6 × 10-5 7.5 × 10-14 1.1 × 10-4 3.5 × 10-5 5.4 × 10-11

Error from the extrapolated moment is marked by italic font.

Figure 13. Moment evolution for constant aggregation kernel obtained from the exact solution (solid lines) and the DAE-QMOM method (circles).

Figure 14. Comparison of errors for the zeroth moment for sum aggregation kernel using the DAE-QMOM (continuous line) and PD-QMOM methods (circles).

problems, the DAE-QMOM is shown to be at least 2 times as fast as the PD-QMOM. It is also interesting to note that DAEQMOM provides a slightly faster solution than MOM for constant growth with primary nucleation. This test confirms that it is possible to speed up the computational time of the PBE at least by 2 times with the new method. In the case of the powerlaw growth with p ) -0.5, the solution time using the PDQMOM appears to be faster than the DAE-QMOM method. However, in this case the PD-QMOM method yields consistently larger errors than DAE-QMOM and fails to conserve the

weights, even for a tight tolerance setting. Therefore, the computation time of the PD-QMOM in this test cannot be compared to the DAE-QMOM solution. For power-law growth with p ) 0.5, the DAE-QMOM provides a slightly faster solution than PD-QMOM; however, in this particular case the difference is not very significant. The performance of DAE and PD based QMOM is also evaluated for different numbers of quadrature points ranging from two to five. In this case the constant aggregation kernel (case 4) is employed due to the availability of analytical moment expression (eq 37). The higher order moments that are not resolved directly from the QMOM solution for a low number of quadrature points have been extrapolated using the quadrature theory (eq 4). The error of the calculated and extrapolated moments at t ) 100 s is calculated using eq 25. Table 2 shows the CPU times required by PD and DAE based QMOM for different numbers of quadrature points. The CPU time generally increases with the number of quadrature points used. The penalty of increasing the number of quadrature points is more serious for PD-QMOM than for DAE-QMOM. The CPU time for the PD increases over 900% by increasing the number of quadrature points from two to five compared to just over 600% in the DAE solution. The DAE-QMOM is significantly faster than the PD-QMOM, with an increasing difference for a greater number of quadrature points. According to Marchisio et al.33 the error of the QMOM predictions decreases as the number of quadrature points increases. Results from this work as shown in Table 3 are in good agreement with Marchisio et al.’s33 findings. At N ) 2, the maximum error of the predicted moment is close to 1%, and when N ) 5 the maximum error falls to about 10-4%. It is also found that there is no significant increase in prediction accuracy from three quadrature points onward. Therefore, it is not worth solving the QMOM beyond three quadrature points, because the computational effort will increase very significantly with only a marginal increase in prediction accuracy. In terms of solution accuracy, the PDQMOM and DAE-QMOM appear to be similar in this particular case. It is also interesting to note that the third moment, which is related to the total particle volume, has been predicted accurately for all cases, even when a small number of quadrature points was used. However, there is a large error from the extrapolated moments via quadrature theory. This is due to the fact that the abscissas and weights of QMOM are a unique solution to represent the particle size distribution at a particular number of quadrature points only, i.e., for the N quadrature points which they are originally assigned. This can be seen clearly when the quadrature theory is used to solve for the weights and abscissas of a normal

7810


Figure 15. Comparison of error for the first six moments for proportional to volume breakage kernel obtained using the DAE-QMOM method (continuous lines) and PD-QMOM method (circles).

6. Conclusions

Figure 16. Moment evolution for proportional to volume breakage kernel using the DAE-QMOM method (solid lines) and the exact solution (circles).

distribution function. For N ) 2, two abscissas of similar weights are produced, while at N ) 3, one abscissa of larger weights is placed exactly at Lmean and the other two on the left and right have the same weight, significantly lower than the weights for Lmean. Thus a significant error is present when the weights and abscissas obtained from N ) 2 are employed to extrapolate the higher order moments, i.e., fourth and fifth, because the abscissas and weights for N ) 2 are no longer a unique solution for N ) 3.

A new QMOM method for solving the population balance equation has been proposed, which solves simultaneously the differential equations for the moments and the system of nonlinear equations resulting from the quadrature approximation as a differential algebraic equation system. The solutions from the proposed method are compared to the product difference based QMOM. Both methods have been validated against the analytical solutions for growth, nucleation, breakage, and coalescence. The results indicate that both methods are capable of predicting accurately the moment evolutions in all cases tested in this study, except for one case of power-law growth. The new numerical solution for the QMOM via DAE method has been demonstrated to be more robust and accurate than the PD algorithm in some cases, especially for growth only problems, characteristic of a large number of applications in particular in the field of crystallization. The DAE-QMOM did not improve the PD-QMOM predictions for the breakage and aggregation kernels evaluated in this work; however, there is a measurable improvement observed for the sum aggregation case. The DAEQMOM has a mathematically simpler formulation eliminating the additional complexity related to the product difference algorithm, providing a more robust solution framework. The only downside of the DAE-QMOM at present is the lack of an automatic Jacobian generator in the DAE solver. At present the Jacobian elements are computed symbolically in software such as MATLAB or MAPLE and then included manually within the DAE-QMOM algorithm. However, the use of DAE solvers with automatic differentiation would make the DAE-QMOM much simpler to implement than the PD-QMOM. Apart from


being able to provide a more accurate solution, the DAE-QMOM also provides a faster solution than both PD-QMOM and MOM for all cases evaluated in this work, assessed from the CPU times required to perform the calculation, which is a significant advantage of the method especially when the PBE has to be solved in the context of CFD simulations.

L ) (L0

1-p

+ (1 - p)bt)

(A.12)

L0 ) (L1-p - (1 - p)bt)1/(1-p)

(A.13)

Integrating eq A.10:

∫

Acknowledgment

n

n0

J.G. is grateful to the scholarship from Ministry of Higher Education, Malaysia, and Universiti Malaysia Pahang.

ln

dn ) -bp n

n ) -bp n0

∫L

PBE for the growth only problem is given by

∂n ∂n ∂G +G +n )0 ∂t ∂L ∂L

(A.2)

p-1

1-p 0

dt

(A.14)

1 dt + (1 - p)bt

(

(

]

)

)

After rearranging

(

n ) n0 1 +

The growth law is given by G ) bLp

t

0

t

0

[

(A.1)

∫L

t 1 ) -bp ln(L01-p + (1 - p)bt) (1 - p)b 0 1-p L + (1 p)bt 0 p ln ) p-1 L01-p (1 - p)bt p/(p-1) ) ln 1 + (A.15) L01-p

Appendix: Derivation of Analytical Solution for Power-Law Growth

∂(Gn) ∂n + )0 ∂t ∂L

7811

1/(1-p)

(A.3)

(1 - p)bt L01-p

)

p/(p-1)

(A.16)

Replace L0 with eq A.13: After differentiation, the growth law is

(

∂G ) bpLp-1 ∂L

( (

(A.4)

Replacing eq A.4 into eq A.2: ∂n ∂n +G ) -nbpLp-1 ∂t ∂L

(A.5)

From the method of characteristics, the characteristic lines in which the PDE reduces to ODEs can be represented as follows: n(L, t) ) n(L(s), t(s))

(A.8)

dL ) bLp dt

(A.9)

dn ) -nbpLp-1 dt

(A.10)

Integrating eq A.9 with respect to L and t: L

L0

dL )b Lp

∫ dt; t

0

)

(A.18)

(

)

j )2 ((L1-p + (1 - p)bt)1/(1-p) - L 1 exp × 2σ2 σ√2π (1 - p)bt p/(1-p) 1(A.19) L1-p

(

)

Literature Cited

Since t - t0 ) s - s0 assuming t0 ) s0 ) 0, thus t ) s. Equations A.8 become

∫

(

j )2 1 (L - L exp 2σ2 σ√2π

The final distribution at any time t is given by n(L) )

dt dL )1 ds ds dn )G ds ) -nbpLp-1

)

When the initial distribution is a normal distribution given

Applying the chain rule

Comparing eq A.7 to A.5 results in

)

)

by n0(L) )

(A.7)

1-p

(

(A.6)

dL ∂n dt ∂n ∂n ) + ∂s ds ∂L ds ∂L

)

p/(p-1) (1 - p)bt L - (1 - p)bt L1-p - (1 - p)bt + (1 - p)bt p/(p-1) ) n0(L0) L1-p - (1 - p)bt 1-p L - (1 - p)bt p/(1-p) ) n0(L0) L1-p (1 - p)bt p/(1-p) ) n0(L0) 1 (A.17) L1-p

n ) n0(L0) 1 +

[( 1 -1 p )L ]

1-p L L0

) bt

(A.11)

(1) Randolph, A. D.; Larson, M. A. Theory of particulate processes: Analysis and techniques of continuous crystallization; Academic Press: New York, 1971. (2) Gerstlauer, A.; Motz, S.; Mitrovic, A.; Gilles, E. Development, analysis and validation of population models for continuous and batch crystallizers. Chem. Eng. Sci. 2002, 57 (20), 4311–4327. (3) Jaworski, Z.; Nienow, A. W. CFD modelling of continuous precipitation of barium sulphate in a stirred tank. Chem. Eng. J. 2003, 91, 167–174. (4) Wei, H.; Zhou, W.; Garside, J. Computational fluid dynamics modeling of the precipitation process in a semibatch crystallizer. Ind. Eng. Chem. Res. 2001, 40, 5255–5261. (5) Woo, X. Y.; Tan, R. B. H.; Chow, P. S.; Braatz, R. D. Simulation of mixing effects in antisolvent crystallization using a coupled CFD-PDFPBE approach. Cryst. Growth Des. 2006, 6, 1291–1303. (6) Woo, X. Y.; Tan, R. B. H.; Braatz, R. D. Modeling and computational fluid dynamics-population balance equation-micromixing simulation of impinging jet crystallizers. Cryst. Growth Des. 2009, 9, 156–164.

7812


(7) Gunawan, R.; Fusman, I.; Braatz, R. D. High resolution finite volume methods for simulating multidimensional population balance equations with nucleation and size dependent growth. AIChE J. 2004, 50, 2738–2749. (8) Gunawan, R.; Fusman, I.; Braatz, R. D. Parallel high resolution simulation of particulate processes with nucleation, growth, and aggregation. AIChE J. 2008, 54, 1449–1458. (9) Zhang, T.; Tade´, M. O.; Tian, Y. C.; Zang, H. High-resolution method for numerically solving PDEs in process engineering. Comput. Chem. Eng. 2008, 32, 2403–2408. (10) Ulbert, Z.; Lakatos, B. G. Dynamic simulation of crystallization processes: Adaptive finite element collocation method. AIChE J. 2007, 53, 3089–3107. (11) Rosner, D. E.; McGraw, R. L.; Tandon, P. Multi-variate population balances via moment- and Monte Carlo simulation methods. Ind. Eng. Chem. Res. 2003, 42, 2699–2711. (12) Lee, K.; Lee, J. H.; Yang, D. R.; Mahoney, A. W. Integrated runto-run and on-line model-based control of particle size distribution for a semi-batch precipitation reactor. Comput. Chem. Eng. 2002, 26, 1117–1131. (13) Gillette, D. A. A study of aging of lead aerosolssII: A numerical model simulating coagulation and sedimentation of a leaded aerosol in the presence of an unleaded background aerosol. Atmos. EnViron. 1970, 6, 451– 462. (14) Sutugin, A. G.; Fuchs, N. A. Formation of condensation aerosols under rapidly changing environmental conditions. Theory and method of calculation. J. Aerosol Sci. 1970, 1, 287–293. (15) Tolfo, F. A simplified model of aerosol coagulation. J. Aerosol Sci. 1977, 8, 9–19. (16) Gelbard, F.; Tambour, Y.; Seinfeld, J. H. Sectional representations for simulating aerosol dynamics. J. Colloid Interface Sci. 1980, 76, 541– 556. (17) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Crystallization and precipitation engineeringsI: An efficient method for solving population balance in crystallization with agglomeration. Chem. Eng. Sci. 1988, 43, 59–67. (18) Chen, P.; Sanyal, J.; Dudukovic, M. P. CFD modeling of bubble columns flows: Implementation of population balance. Chem. Eng. Sci. 2004, 59, 5201–5207. (19) Dhanasekharan, K. M.; Sanyal, J.; Jain, A.; Haidari, A. A generalized approach to model oxygen transfer in bioreactors using population balances and computational fluid dynamics. Chem. Eng. Sci. 2005, 60, 213–218. (20) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretizationsIII. Nucleation, growth and aggregation of particles. Chem. Eng. Sci. 1997, 52, 4659–4679. (21) Hounslow, M. J.; Pearson, J. M. K.; Instone, T. Tracer studies of high-shear granulation: II. Population balance modeling. AIChE J. 2001, 47, 1984–1999. (22) Chen, Z.; Pruss, J.; Warnecke, H.-J. A population balance model for disperse systems: Drop size distribution in emulsion. Chem. Eng. Sci. 1998, 53, 1059–1066. (23) Alopaeus, V.; Koskinen, J.; Keskinen, K. I. Utilization of population balances in simulation of liquid-liquid systems in mixed tanks. Chem. Eng. Commun. 2003, 190, 1468–1484. (24) Floury, J.; Legrand, J.; Desrumaux, A. Analysis of a new type of high pressure homogeniser. Part B. study of droplet break-up and recoalescence phenomena. Chem. Eng. Sci. 2004, 59, 1285–1294. (25) Costa, C. B. B.; Maciel, M. R. W.; Filho, R. M. Considerations on the crystallization modeling: Population balance solution. Comput. Chem. Eng. 2007, 31, 206–218. (26) Gimbun, J.; Rielly, C. D.; Nagy, Z. K. Modelling of mass transfer in gas-liquid stirred tanks agitated by Rushton turbine and CD-6 impeller: A scale-up study. Chem. Eng. Res. Des. 2009, 87, 437–451. (27) Moilanen, P.; Laakkonen, M.; Visuri, O.; Alopaeus, V.; Aittamaa, J. Modelling mass transfer in an aerated 0.2 m3 vessel agitated by Rushton, Phasejet and Combijet impellers. Chem. Eng. J. 2008, 142, 95–108. (28) McGraw, R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255–265. (29) Gordon, R. G. Error bounds in equilibirium statistical mechanics. J. Math. Phys. 1968, 9, 655–663.

(30) Wright, D. L.; McGraw, R.; Rosner, D. E. Bivariate extension of the quadrature method of moments for modeling simultaneous coagulation and sintering of particle populations. J. Colloid Interface Sci. 2001, 236, 242–251. (31) Rosner, D. E.; Pyykonen, J. J. Bivariate moment simulation of coagulating and sintering nanoparticles in flames. AIChE J. 2002, 48, 476– 491. (32) Fan, R.; Marchisio, D. L.; Fox, R. O. Application of the direct quadrature method of moments to polydisperse gas-solid fluidized beds. Powder Technol. 2004, 139, 7–20. (33) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322–334. (34) Lambin, P.; Gaspard, J. P. Continued-fraction technique for tightbinding systems. A generalized-moments method. Phys. ReV. B 1982, 26, 4356–4368. (35) Gautschi, W. Algorithm 726: ORTHPOL-a package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Trans. Math. Software 1994, 20, 21–62. (36) McGraw, R.; Wright, D. L. Chemically resolved aerosol dynamics for internal mixtures by the quadrature method of moments. J. Aerosol Sci. 2003, 34, 189–209. (37) Dorao, C. A.; Jakobsen, H. A. The quadrature method of moments and its relationship with the method of weighted residuals. Chem. Eng. Sci. 2006, 61, 7795–7804. (38) Marchisio, D. L.; Fox, R. O. Solution of population balance equations using the direct quadrature method of moments. J. Aerosol Sci. 2005, 36, 43–73. (39) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes in C: The Art of Scientific Computing, 2nd ed.; Cambridge University Press: Cambridge, 1992. (40) Su, J.; Gu, Z.; Li, Y.; Feng, S.; Xu, X. Y. An adaptive direct quadrature method of moment for population balance equations. AIChE J. 2008, 54, 2872–2887. (41) Alopaeus, V.; Laakkonen, M.; Aittamaa, J. Numerical solution of moment-transformed population balance equation with fixed quadrature points. Chem. Eng. Sci. 2006, 61, 4919–4929. (42) Diemer, R. B.; Ehrman, S. H. Pipeline agglomerator design as a model test case. Powder Technol. 2005, 156, 129–145. (43) Grosch, R.; Briesen, H.; Marquardt, W.; Wulkow, M. Generalization and numerical investigation of QMOM. AIChE J. 2007, 53, 207–227. (44) Rod, V.; Misˇek, T. Stochastic modelling of dispersion formation in agitated liquid-liquid systems. Trans. Inst. Chem. Eng. 1982, 60, 48–53. (45) Hulburt, H. M.; Katz, S. Some problems in particle technology. Chem. Eng. Sci. 1964, 19, 555–574. (46) Ascher, U. M.; Petzold, L. R. Computer Methods for Ordinary Differential Equations and Differential Algebraic Equations, 1st ed.; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1998. (47) Shampine, L. F.; Reichelt, M. W. MATLAB ODE suite. SIAM J. Sci. Comput. 1997, 18, 1–22. (48) Klopfenstein, R. W. Numerical differentiation formulas for stiff systems of ordinary differential equations. RCA ReV. 1971, 32, 447–462. (49) Hindmarsh, A. C. LSODE and LSODI, two new initial value ordinary differential equation solvers. ACM SIGNUM Newslett. 1980, 15, 10–11. (50) Brown, P. N.; Byrne, G. D.; Hindmarsh, A. C. VODE: A variablecoefficient ODE solver. SIAM J. Sci. Comput. 1989, 10, 1038–1051. (51) Alopaeus, V.; Laakkonen, M.; Aittamaa, J. Solution of population balances with growth and nucleation by high order moment-conserving method of classes. Chem. Eng. Sci. 2007, 62, 2277–2289. (52) Gelbard, F.; Seinfeld, J. H. Numerical solution of the dynamic equation for particulate systems. J. Comput. Phys. 1978, 28, 357–375. (53) Melzak, Z. A. The effects of coalescence in certain collision processes. Q. J. Appl. Math. 1953, 11, 231. (54) Scott, W. T. Analytic studies on cloud droplet coalescence. J. Atmos. Sci. 1968, 25, 54.

ReceiVed for reView April 3, 2009 ReVised manuscript receiVed June 16, 2009 Accepted June 21, 2009 IE900548S

Simultaneous Quadrature Method of Moments for the Solution of

Recommend Documents