Ind. Eng. Chem. Res. 2010, 49, 11633–11644
11633
Bivariate Extension of the Quadrature Method of Moments for Batch Crystallization Models Shamsul Qamar,*,†,‡ Saima Noor,‡ Qurrat ul Ain,§ and Andreas Seidel-Morgenstern† Max Planck Institute for Dynamics of Complex Technical Systems, Sandtorstrasse 1, 39106 Magdeburg, Germany, Department of Mathematics, COMSATS Institute of Information Technology, Park Road Chak Shahzad Islamabad, Pakistan, and Pakistan Institute of Science and Technology, Nilore Islamabad, Pakistan
This Article presents a bivariate extension of the quadrature method of moments for solving two-dimensional batch crystallization models involving crystals nucleation, size-dependent growths, aggregation, and dissolution of small nuclei below certain critical size in a dissolution unit. In this technique, orthogonal polynomials of lower order moments are used to find the quadrature abscissas (points) and weights. Several benchmark problems with different combinations of processes are considered in this Article. The accuracy and efficiency of the proposed method are validated against the analytical solutions and the high-resolution finite volume scheme. Excellent agreements were observed in all test problems. It was found that the current method is very efficient and accurate as compared to the high-resolution finite volume scheme. Introduction Crystallization is a commonly used separating and purifying technique in chemical, pharmaceutical, semiconductor, and food industries. A control of the crystal’s shape and size is very important for achieving the desired product and for improving the down stream processing. Simulation of the underlying process enables one to investigate the effects of different operating conditions. The observed data can be used for controlling and optimizing the quality of a product. The fines removal and their subsequent dissolution in an external dissolution unit was found useful for improving the quality of product crystal size distribution (CSD). It withdraws and dissolves excessive fines from the annular zone of crystallizer, which are generated during periods of high supersaturation. This effectively shifts the CSD toward right and often makes the distribution narrower. Fines dissolution is usually exercised in industrial applications of continuous crystallization.1 However, it is also practiced to control batch crystallizers.2-4 Population balance models (PBMs) have been widely used for modeling the dynamics of crystallization processes since the mid-1960s.5,6 These equations form hyperbolic partial differential equations and simulate a wide range of particulate processes including comminution, crystallization, granulation, flocculation, combustion, and polymerization. The major phenomena that influence these processes include growth, nucleation, aggregation, breakage, dissolution, and inlet and outlet streams. During the past decades, several numerical methods have been developed for solving population balance equations (PBEs). These methods are used to simulate either the evolution of moments of CSD or the CSD itself. The available methods include the method of moments,5,7-10 the method of characteristics,11,12 the method of weighted residual or orthogonal collocation,13 the Monte Carlo simulation,14,15 the finite difference schemes or discrete population balances,16 the spectral methods,17-19 and the highresolution finite volume schemes.20-22 Hulburt and Katz5 pioneered the application of the method of moments to PBMs of single property variable. Different * To whom correspondence should be addressed. E-mail: qamar@ mpi-magdeburg.mpg.de. † Max Planck Institute. ‡ COMSATS Institute. § PINSTECH.
methods were proposed for solving the closure problem raised by Hulburt and Katz,5 which are briefly discussed by Diemer and Olson.23 The quadrature method of moments (QMOM) was introduced by McGraw10 for modeling aerosol evolution. The method is based on the approximation of integrals involving the particle size distribution (PSD) through a quadrature method. In this method, the product difference (PD) algorithm24 was used for finding the quadrature abscissas and weights. Barrett and Webb25 have compared the method with other available approaches, such as Laguerre quadrature approximation and the finite element method, by solving the aerosol general dynamic equation. Afterward, the QMOM was extended for simultaneous aggregation and breakage problems.9 Moreover, Fan et al.26 have proposed the direct quadrature method of moments (DQMOM) as an alternative to the PD algorithm for finding the quadrature points and weights. In this method, variables appearing in the quadrature approximation are tracked directly by solving the convection equation for quadrature abscissas and weights. Later, the DQMOM was further improved by introducing an adaptive factor in the moment equations.27 A comparison of different QMOM and discussion about their limitations is also available in the literature.28 Recently, a new QMOM method was introduced for solving the PBE, which simultaneously solves the differential equations for the moments and a system of nonlinear equations resulting from the quadrature approximation as a differential algebraic equation system.29 Finally, Qamar et al.30 have used an alternative orthogonal polynomials-based QMOM method for solving one-dimensional batch crystallization models. Hulburt and Katz5 also outlined the bivariate extension of the method of moments by considering the evolution of two radii of curvature of ellipsoidal particles in a continuously fed batch reactor. The authors have briefly summarized the possible solutions and associated difficulties, such as the involvement of a large number of mixed moments and problems related to the reconstruction of bivariate distribution from bivariate moments. Their work did not get attention until the article by Wright et al.31 on the extension of QMOM method for modeling the dynamics of a population of inorganic nanoparticles undergoing simultaneous coagulation and particle sintering. For validation, the authors have compared the QMOM results with those obtained from the high resolution discrete method. They
10.1021/ie101108s 2010 American Chemical Society Published on Web 09/29/2010
11634
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
found that bivariate QMOM is very efficient and accurate as compared to the given discrete method. Later, the method was further elaborated by Rosner and Pyyko¨nen32 for the simulation of coagulating and sintering nanoparticles in flames. Recently, Zucca et al.33 have used DQMOM for solving bivariate PBEs. The accuracy of their method was assessed by a comparison with Monte Carlo method simulations. This Article extends the quadrature method of moment30 for solving two-dimensional batch crystallization models involving crystals nucleation, size-dependent growths, aggregation, and dissolution of small nuclei below certain critical size. In this technique, orthogonal polynomials of lower order moments are employed to find the quadrature abscissas (points) and weights. For the sake of better accuracy and optimal computational cost, a three-point quadrature method is derived, which needs a third-order orthogonal polynomial constructed from a matrix of 36 moments. The proposed method is not restricted to the calculation of specified number of moments; one can calculate as many moments as required. The numerical results of current QMOM are evaluated against the analytical solutions and the semidiscrete high-resolution finite volume scheme (FVS).22 In the case of finite volume scheme, the time evolution of selected moments is obtained by integrating the full bivariate distribution with Trapezoidal rule. In both methods, the resulting systems of ordinary differential equations (ODEs) were solved by the Runge-Kutta method of order four. It was observed that, even in the bivariate case, QMOM can be efficiently used to simulate the underlying physical process and to accurately predict its behavior. It is important to mention that the current QMOM follows an approach similar to that presented by Wright et al.31 The quadrature method of moments has low computational costs and better accuracy, but no information about the particle size distribution (PSD) is available. Instead, only a finite number of moments associated with the real distribution are calculated. There are several physical and optical properties of aerosol and clouds that can be easily estimated from the knowledge of lower order moments instead of the distribution itself.10,34 Complicated phenomena, such as nucleation, growth, and transport, can change the shape of aerosol size distribution. However, it is difficult to incorporate these phenomena in the models describing aerosol in multidimensional environments, for instance, models of atmospheric transport and turbulent jet flows. The modeling of size distribution not only produces numerical difficulties but also contains more information than actually needed. Therefore, the method of moments can be a good alternative in such situations. Moreover, QMOM is suitable when PBM is implemented in computational fluid dynamics (CFD) codes.34,35 In such scenarios, the external features like turbulent flow properties play an important role, and, hence, the distribution is a function of both internal and external coordinates, which amplifies the computational time. The moments of internal coordinates convert the given PBE to a system of transport equations that can be easily included in the given CFD code. Furthermore, several techniques are available for reconstructing PSD from finite number of moments.36 This Article is organized as follows. In section 2, a two-dimensional batch crystallization model for simultaneous processes is presented. In section 3, the proposed bivariate quadrature method of moments is derived. In section 4, numerical test problems are presented. Finally, section 5 gives conclusions and remarks. Batch Crystallization Model In this section, a mathematical model for an ideally mixed batch crystallizer equipped with a fines dissolution unit is introduced. Crystallization processes form a disperse system where the solid phase is dispersed in the continuous medium of the liquid phase. In the framework of PBEs, the state of an individual solid crystal is defined by internal coordinates representing its size. A population of crystals is characterized by its CSD, which is mathematically described by a number density function n(t,x,y) as a function of time t and size coordinates (x,y). This function represents the (average) number of crystals per crystal size. The rate of change in CSD is described by the population balance equation:5,6
n(0, x, y) ) n0(x, y)
(2)
Here, n0(x,y) ∈ R2 denotes the CSD of seed crystals added at the beginning of the batch process, G1(t,x) g 0 and G2(t,y) g 0 are the crystals growth rates along the characteristic directions x and y, B0(t) g 0 is the nucleation rate at minimum crystal size (x0,y0), and δ is the Dirac delta distribution. The third term on the right-hand side quantifies a reduction in the number of solid particles due to removal of fines from the crystallizer to the dissolution unit. The removal of fines is characterized by the death function h(x,y). In that term, Vcrz denotes the volume of the crystallizer, and V˙ represents the volumetric flow rate of the outgoing stream. A balance law for the liquid phase takes the form
(
)
∂Vc(x, y) ∂Vc(x, y) ∞ ∞ dm(t) + G2(t, y) )m ˙ in(t) - m ˙ out(t) - Fc 0 0 G1(t, x) n(t, x, y) dx dy (3) dt ∂x ∂y Here, Vc represents the volume of a single crystal, and Fc is a crystal density. Because of external fines dissolution, this equation
∫ ∫
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
has two mass fluxes (streams). The first one, m ˙ out(t), is a liquid stream containing fines from the crystallizer to the dissolution unit. The second one, m ˙ in(t), denotes the incoming pure and particle free liquid stream from the dissolution unit to the crystallizer. The third term on the right-hand side of eq 3 quantifies the effect of the growing crystals, which reduce the amount of solute in the liquid phase. The mass fluxes are defined as m ˙ out(t) ) w(t)Fsolu(T)V˙ FcV˙ ∞ ∞ m ˙ in(t) ) m ˙ out(t - tp) + V (x, y)h(x, y) × Vcrz 0 0 c n(t - tp, x, y) dx dy - FcVc(x0, y0)B0(t - tp)
(4)
∫ ∫
(5)
The second term on the right-hand side of eq 5 quantifies an increase in the solute mass due to external fines dissolution. The mass fraction w(t) is given as w(t) )
m(t) m(t) + msolv(T)
(6)
Here, msolv(T) is mass of the solvent, and Fsolu is density of the solution; both are depending on the temperature T. The temperature can be constant (isothermal case) or can be a decreasing function of time (nonisothermal case). Moreover, tp g 0 represents a residence time in the dissolution unit (the pipe). It is defined as Vp V˙
tp )
(7)
where Vp is the volume of the pipe. The size-dependent crystal growths can be defined as37,38 G1(t, x) ) kgx[S(t) - 1]gx(1 + R1x)R2 :) g1(t)g2(x)
(8)
G2(t, y) ) kgy[S(t) - 1]gy(1 + R3y)R4 :) g3(t)g4(y)
(9)
where kgx and kgy are growth rates constant. The exponents gx and gy denote growths orders, and R1, R2, R3, and R4 are constants. Moreover, S(t) denotes supersaturation of the dissolved component: S(t) )
w(t) wsat(T)
(10)
where wsat(T) represents saturated mass fraction. The nucleation rate is typically defined as1,39 B0(t) ) kb[S(t) - 1]b
∫ ∫ ∞
0
∞
0
Vcn(t, x, y) dx dy
(11)
where kb denotes the nucleation rate constant, and the exponent b represents the nucleation order. Note that the above model reduces to a batch crystallization model without fines dissolution when the third term on the righthand side of eq 2 and the first two terms on the right-hand side of eq 3 are zero. Equations 4 and 5 are then not needed. The klth moment of CSD is defined as µk,l(t) )
∫ ∫ ∞
0
∞ k l
0
x y n(t, x, y) dx dy, k, l ) 0, 1, 2, ...
(12) After the above moments transformation was applied, PBE 2 gives5,29
dµk,l(t) dt
∫ ∫ ∞
)-
0
∞
11635
[kxk-1ylG1(t, x)n(t, x, y) +
0
lxkyl-1G2(t, y)n(t, x, y)] dx dy
∫ ∫ x y h(x, y)n(t, x, y) dx dy + x y B (t) 1 + ∫ ∫ ∫ ∫ β(u, V, x', y')(u + x ) (V + y ) 2 V˙ Vcrz
-
∞
∫ ∫ ∞
0
∞ k l
∞
0
∞
0
3
3 k/3
3 l/3
3
0
∞ k l
0
k l 0 0 0
0
∞
0
-
∞
0
x y n(t, x, y)
∫ ∫ ∞
0
∞
0
×
n(t, u, V)n(t, x', y') dx' dy' du dV β(x, y, x', y')n(t, x', y') dx' dy' dx dy
(13)
The integral terms in the above equation cannot be replaced by moments due to size-dependent growth rates, due to death function h(x,y), and due to complicated integrands of the aggregation term. Moreover, the CSD n(t,x,y) is not available at each time step. Therefore, a standard quadrature approximation, such as the trapezoidal rule or Simpson’s rule, will not close the ODE-system 14. However, the quadrature method of moments (QMOM) can be applied to overcome this closure problem as explained below. Gaussian Quadrature Method for Bivariate Moments In the Gaussian quadrature method, a definite integral of a function is approximated as a weighted sum of function values at specified points within the domain of integration. Usually, a definite integral is approximated by summing the functional values at a set of equally spaced points, each value multiplied by a certain weight. However, the Gaussian quadrature rule gives a freedom for choosing both weights and points (abscissas) at which the function is evaluated. Moreover, the chosen abscissas may not be equally spaced. Let us consider an integral of the form ∫ba∫dcψ(x,y)f(x,y) dx dy, where ψ(x,y) denotes a non-negative weight function, and f(x,y) is a nonspecified function of x and y. Next, we can find a set of weights wi and abscissas (xi, yi) such that the approximation
∫ ∫ b
a
d
c
N
ψ(x, y)f(x, y) dx dy ≈
∑ w f(x , y ) i
i
i
(14)
i)1
is exact if f(x,y) is a sufficiently smooth function. Orthogonal polynomials are the backbones of the Gaussian quadrature method. To calculate a special orthogonal polynomial of order N, one has to construct a set of polynomials that includes exactly one polynomial of order i for each i ) 1, 2, ..., where all of them are mutually orthogonal over a specified weight function ψ(x,y).40 This can be explained by defining the scalar product of two functions r(x,y) and s(x,y) over a weight function ψ(x,y) as
〈r|s〉 )
∫ ∫ b
d
a
c
ψ(x, y)r(x, y)s(x, y) dx dy
(15)
Now, two functions are said to be orthogonal if their scalar product is zero. If there is no classical weight function ψ(x,y) given, as in our case, one needs additional information to obtain the weights and abscissas. In the current situation, one can use µk,l(t) )
∫ ∫ ∞
0
∞ k l x y n(t, x, y) 0
N
dx dy ≈
∑x yw
k l i i i
(16)
i)1
where xi and yi are the quadrature points in both characteristic directions, and wi represents quadrature weights. In this study, a threepoint bivariate quadrature formula is derived, that is, N ) 3. After the quadrature rule 16 was applied, the integral terms in eqs 3 and 13 can be approximated as
11636
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
(
)
N ∂V(xi, yi) ∂V(xi, yi) dm(t) G1(t, xi) + G2(t, yi) )m ˙ in(t) - m ˙ out(t) - Fc wi dt ∂x ∂y i)1
∑
(18)
Because of our assumptions, eqs 4 and 5 can be rewritten as
{
m ˙ out(t) ) w(t)Fsolu(T)V˙
(19)
m ˙ out(t),
m ˙ in(t) )
if t e tp
FcV˙ m ˙ out(t - tp) + Vcrz
N
∑ w V (x , y )h(x , y ) - F V (x y )B (t - t ), i c
i
i
i
i
c c
0 0
0
p
if t > tp
(20)
i)1
In the following, orthogonal polynomials will be derived for finding quadrature points and weights. Let us define z(x, y) :) x + xy + y
(21)
and a recursion relation of the form p-1 ) 0,
p0 ) 1,
pi ) (z - ai)pi-1 - bipi-2,
i ) 1, 2, ...
(22)
with ai )
〈zpi-1|pi-1〉 , i ) 1, 2, ... 〈pi-1|pi-1〉
(23)
bi )
〈pi-1|pi-1〉 , i ) 2, 3, ... 〈pi-2|pi-2〉
(24)
Because n(t,x,y) is used as a weight function ψ(x,y), therefore
〈pi|pi〉 )
∫ ∫ b
a
d
c
n(t, x, y)p2i dx dy
(25)
These definitions can be used to calculate orthogonal polynomials one after another until the required Nth-order polynomial is obtained. For understanding, the first two polynomials are calculated. The first-order polynomial is given as p1(z) ) (z - a1)p0 ) z - a1
(26)
By using eqs 21-23, a1 can be calculated as
〈zyp0|p0〉 a1 ) ) 〈p0|p0〉
∫ ∫ (x + xy + y)n(t, x, y)p ∫ ∫ n(t, x, y)p dx dy ∞
0
∞
2 0
0
∞
0
∞
0
dx dy
2 0
)
µ1,0 + µ1,1 + µ0,1 µ0,0
(27)
Thus, the first-order orthogonal polynomial is expressed as p1(z) ) z -
µ˜ 1,1 µ0,0
(28)
where µ˜ 1,1 ) µ1,0 + µ1,1 + µ0,1 Next, p2 can be defined as
(29)
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010 p2(z) ) (z - a2)p1 - b2p0
11637
(30)
where, according to eqs 23 and 28: 〈zp1|p1〉 a2 ) ) 〈p1|p1〉
∫ ∫ zn(t, x, y)p dx dy ∫ ∫ n(t, x, y)p dx dy ∞
∞
0
2 1
0
∞
∞
0
)
3 2 - 2µ0,0µ˜ 1,1µ˜ 2,2 + µ˜ 1,1 µ˜ 3,3µ0,0 2 2 µ˜ 2,2µ0,0 - µ0,0µ˜ 1,1
2 1
0
(31)
Here, µ˜ 1,1 is given by eq 29 and µ˜ 2,2 ) µ2,0 + µ0,2 + 2(µ2,1 + µ1,2 + µ1,1) + µ2,2
(32)
µ˜ 3,3 ) µ3,0 + µ0,3 + 3(µ3,1 + µ1,3 + µ2,1 + µ1,2 + µ3,2 + µ2,3 + 2µ2,2) + µ3,3
(33)
Similarly, eqs 24 and 28 give
∫ ∫ ∞
b2 )
〈p1|p1〉 ) 〈p0|p0〉
0
∞
0
(
n(t, x, y) z -
∫
∞
0
µ1,1 µ0,0
)
2
dx dy )
2 µ˜ 2,2µ0,0 - µ˜ 1,1
n(t, x, y) dx dy
2 µ0,0
Therefore, eq 30 for the second-order orthogonal polynomial becomes p2(z) )
2 2 z2(µ0,0µ˜ 2,2 - µ˜ 1,1 ) + z(µ˜ 1,1µ˜ 2,2 - µ0,0µ˜ 3,3) + µ˜ 1,1µ˜ 3,3 - µ˜ 2,2 2 µ0,0µ˜ 2,2 - µ˜ 1,1
(34)
In this manner, higher order polynomials can be calculated. A third-order polynomial is given as p3(z) ) z3 + η2z2 + η1z + η0
(35)
where η2 ) η1 ) η0 )
2 2 2 µ˜ 2,2µ˜ 4,4µ˜ 1,1 - µ0,0µ˜ 4,4µ˜ 3,3 + µ˜ 2,2µ0,0µ˜ 5,5 + µ˜ 3,3 µ˜ 1,1 - µ˜ 5,5µ˜ 1,1 - µ˜ 2,2 µ˜ 3,3 3 2 2 µ˜ 2,2 - µ˜ 2,2µ˜ 4,4µ0,0 - 2µ˜ 2,2µ˜ 3,3µ˜ 1,1 + µ˜ 3,3 µ0,0 + µ˜ 4,4µ˜ 1,1 2 2 2 µ˜ 2,2µ˜ 5,5µ˜ 1,1 + µ0,0µ˜ 4,4 - µ0,0µ˜ 5,5µ˜ 3,3 - µ˜ 4,4µ˜ 3,3µ˜ 1,1 - µ˜ 2,2 µ˜ 4,4 + µ˜ 2,2µ˜ 3,3 3 2 2 µ˜ 2,2 - µ˜ 2,2µ˜ 4,4µ0,0 - 2µ˜ 2,2µ˜ 3,3µ˜ 1,1 + µ˜ 3,3 µ0,0 + µ˜ 4,4µ˜ 1,1
(36)
2 3 2 2µ˜ 2,2µ˜ 4,4µ˜ 3,3 - µ˜ 2,2 µ˜ 5,5 - µ˜ 3,3 - µ˜ 4,4 µ˜ 1,1 + µ˜ 5,5µ˜ 3,3µ˜ 1,1 3 2 2 µ˜ 2,2 - µ˜ 2,2µ˜ 4,4µ0,0 - 2µ˜ 2,2µ˜ 3,3µ˜ 1,1 + µ˜ 3,3 µ0,0 + µ˜ 4,4µ˜ 1,1
Here, µ˜ 1,1, µ˜ 2,2, and µ˜ 3,3 are given by eqs 29, 32, and 33. Moreover, µ˜ 4,4 ) µ4,0 + µ0,4 + 6(µ2,2 + 2µ3,2 + 2µ2,3 + 2µ3,3 + µ4,2 + µ2,4) + 4(µ1,3 + µ3,1 + µ4,1 + µ1,4 + µ4,3 + µ3,4) + µ44
(37)
µ˜ 5,5 ) µ5,0 + µ0,5 + 5(µ4,1 + µ1,4 + µ5,1 + µ1,5 + µ5,4 + µ4,5) + µ55 + 10(µ3,2 + µ2,3 + 3µ3,3 + 3µ4,3 + 3µ3,4 + 2µ4,2 + 2µ2,4 + 2µ4,4 + µ5,2 + µ2,5 + µ5,3 + µ3,5)
(38)
The roots of this polynomial are zi, and the weights wi can be calculated as40 wi )
〈pN-1|pN-1〉 , pN-1(zi)p'N(zi)
i ) 1, 2, 3,
N)3
(39)
From the above equations, zi ) xi + xiyi + yi and wi are known. The next step is to find the abscissas xi and yi. For that purpose, three more equations are needed as given below: 3
µ˜ 2,1 )
∑xzw
) x1z1w1 + x2z2w2 + x3z3w3
(40)
2 i i i
) x1z21w1 + x2z22w2 + x3z23w3
(41)
∑xz w
) x1z31w1 + x2z32w2 + x3z33w3
(42)
i i i
i)1 3
µ˜ 3,2 )
∑xz w i)1 3
µ˜ 4,3 )
3 i i i
i)1
where
11638
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
Figure 1. Test problem 1 (case 1): Results of two-component aggregation with exponential distribution and constant kernel.
µ˜ 2,1 ) µ2,0 + µ2,1 + µ1,1
(43)
µ˜ 3,2 ) µ3,0 + µ1,2 + 2(µ3,1 + µ2,2 + µ2,1) + µ3,2
(44)
µ˜ 4,3 ) µ4,0 + µ1,3 + 3(µ2,2 + µ2,3 + µ3,1 + 2µ3,2 + µ3,3 + µ4,1 + µ4,2) + µ4,3 (45)
x2 )
µ˜ 4,3 - (z1 + z3)µ˜ 3,2 + z1z3µ˜ 2,1 z2w2(z3 - z2)(z1 - z2)
(47)
x3 )
µ˜ 4,3 - (z1 + z2)µ˜ 3,2 + z1z2µ˜ 2,1 z3w3(z3 - z1)(z3 - z2)
(48)
z1 - x1 z2 - x2 z3 - x3 , y2 ) , y3 ) 1 + x1 1 + x2 1 + x3
(49)
Next, The solution of eqs 40-42 gives
x1 )
µ˜ 4,3 - (z2 + z3)µ˜ 3,2 + z2z3µ˜ 2,1 z1w1(z1 - z2)(z1 - z3)
y1 )
(46)
Finally, the resulting system of ordinary differential equations (ODEs) in eqs 17 and 18 can be solved by any standard ODE-
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
11639
Figure 2. Test problem 1 (case 2): Results of two-component aggregation with Gaussian-type initial distribution and constant kernel.
solver. In this study, a Range-Kutta method of order four was used. Note that multiple three-point quadrature formulas are possible. For example, instead of z ) x + xy + y, one can use either z ) xy or z ) x + y to derive a third-order orthogonal polynomial (35) utilizing few moments out of 36 moments. Similarly, different combinations of moments can be used in eqs 40-42 to obtain xi or yi. The quadrature points and weights determined in this manner can be used in eqs 17 and 18 to propagate any number of moments from their initial values. However, the three-point quadrature methods, which do not employ all 36 moments, may give error in the calculation of those moments, which are not used in eqs 35 and 40-42. For that reason, the current three-point quadrature method is an optimum one.
To improve the accuracy of the quadrature formula, one can also derive the 12-point quadrature formula, which utilizes all 36 moments to obtain 12 quadrature points (xi,yi) and weights wi by using31 12
µk,l(t) )
∑x yw, k l i i i
k, l ) 1, 2, ..., 6
(50)
i)1
In this case, we have 36 equations for 36 unknowns. The resulting linear system can be solved for the required quadrature points and weight with the help of conjugate-gradient minimization algorithm.40 From the initial guess of quadrature points and weight, a 36-dimensional minimization can be used to locate values of the parameters that satisfy eq 50.31,40 However, the 36-dimensional minimization can be very difficult and compu-
11640
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
Figure 3. Test problem 1 (case 3): Results of two-component aggregation with exponential distribution and sum kernel. Symbols represent QMOM solution, and lines denote FVS solution.
Figure 4. Test problem 1 (case 4): Results of two-component aggregation with exponential distribution and Brownian-type kernel. Symbols represent QMOM solution, and lines denote FVS solution.
tationally expensive. A good initial guess can reduce the overall computational cost up to certain level. For example, the initial guess values can be obtained from the combination of four threepoint quadrature formulas, which can be grouped together to form a set of 12 quadrature points and weights. Numerical Test Problems The following three numerical test problems are considered for validation of the current QMOM. Test Problem 1: Two-Component Aggregation Problem. The analytical solutions are only available for twocomponent aggregation problems with constant kernel. A square mesh of 40 × 40 mesh elements and a geometric grid of the form xi ) 10-6 + 2(i-Nx)/3xmax,
yj ) 10-6 + 2(j-Ny)/3ymax, i, j ) 0, 1, 2, 3, ...
(51)
were considered. Here, Nx,Ny denote the maximum discretization points in the x- and y-directions. Moreover, xmax,ymax represent the maximum characteristic lengths in each direction. Case 1: Exponential Initial Distribution and Constant Kernel. The initial distribution is given as n(0, x, y) ) 9x2y2 exp(-x3 - y3) and the analytical solution is given as41
(52)
n(t, x, y) )
(
( )
4tx3y3 1/2 6xy 2 exp(-x3 - y3)I0(θ), with θ ) t+2 t+2 (53)
)
Here, I0 is the modified Bessel function of first kind of order zero. In Figure 1, the numerical results of QMOM and finite volume scheme (FVS)22 are plotted against the analytical solutions. The numerical results of both schemes are in good agreement with analytical solutions. However, one can see an overestimation in the result of FVS as compared to QMOM results. The plots of relative error show that QMOM Table 1. Parameters for Test Problem 2 description growth rate constant growth rate constant growth rate exponent growth rate exponent nucleation rate constant nucleation rate exponent density of crystals initial solute mass saturated mass fraction mass of seeds mass of solvent density of solution constant (eq 57) constant (eq 57) constant (eq 57) volume of the crystallizer volumetric flow rate volume of the pipe
symbols kgx kgy gx gy kb b Fc m (0) wsat mseeds msolv Fsolu σ jx jy Vcrz V˙ Vp
value
unit -5
0.68 × 10 1.37 × 10-5 0.73 0.73 3.42 × 107 2.35 1250 0.09915 0.090681 2.5 × 10-3 0.8017 1000 2.1 × 10-4 7.0 × 10-4 1.0 × 10-3 10-3 2 × 10-5 2.4 × 10-4
(m)/(min) (m)/(min) (1)/(m3 min) (kg)/(m3) kg kg kg (kg)/(m3) m m m m3 (m)/(min) m
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
11641
Figure 5. Test problem 2: Results for size-independent and size-dependent growth rates. Symbols represent QMOM solution, and lines are for FVS solution.
has produced less error in the solution as compared to the FVS. As expected, QMOM was found very efficient and accurate for solving bivariate aggregation problems. The program for QMOM is written in Matlab software under Linux operating system and was compiled on a computer with Intel(R) Core 2 Duo processor of speed 2 GHz and memory (RAM) 3.83 GB. The CPU time for QMOM was 0.2 s, which is much lesser than hours required by FVS if programmed in Matlab. Especially, FVS is very timeconsuming for solving aggregation problems due to the involvement of several nested loops. Case 2: Gaussian-type Initial Distribution and Again Constant Kernel. The initial distribution is given as
n(t, x, y) )
72x2y2
√t(t + 2)
The analytical solution is given as42
(54)
exp(-2x3 - 2y3)[I0(θ) - J0(θ)], with θ ) 4(xy)3/2
( t +t 2 )
1/4
(55)
Here, J0 and I0 are the Bessel function and the modified Bessel functions of first kind of order zero, respectively. The numerical results are shown in Figure 2. The results of QMOM and finite volume scheme (FVS) are plotted against the analytical solutions. A good agreement was found in both numerical schemes. Table 2. Test Problem 2: Errors in Mass Balances (Size-Independent Growth) description
n(0, x, y) ) 144x5y5 exp(-2x3 - 2y3)
3
CPU time absolute error relative error QMOM (s)
without fines dissolution 2.77 × 10-17 fines dissolution without delay 5.3 × 10-6 fines dissolution with delay 3.96 × 10-3
2.7 × 10-16 1.5 × 10-5 1.15 × 10-2
5.28 5.30 5.34
11642
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
Figure 6. Test problem 3: Results for nucleation and growth processes only.
Figure 7. Test problem 3: Results for nucleation, growth, and aggregation processes. Symbols are used for QMOM solution, and lines denote FVS solution.
However, one can see that QMOM gives better results as compared to the FVS scheme, which is also clear from the plots of relative error. Case 3: Exponential Initial Distribution and Sum Kernel. The initial distribution is the same as in case 1. However, instead of the constant kernel, a sum kernel β(t,x,y,x′,y′) ) β0(x3 + y3 + x′3 + y′3) is considered with β0 ) 1. In the two-component case, we have no exact solutions for other than a constant kernel. Therefore, the numerical results of QMOM and FVS are compared in Figure 3. The figure shows that QMOM can be efficiently applied to the aggregation problem with sum kernel. Once again, both schemes give comparable solutions. Case 4: Exponential Initial Distribution and Browniantype Kernel. The initial distribution is the same as in case 1. Here, a Brownian-type kernel:31 β(t, x, y, x', y') )
(x3 + x'3 + y3 + y'3)2 (x3 + x'3)(y3 + y'3)
(56)
is considered. The numerical results of QMOM and FVS are compared in Figure 4. The figure verifies that QMOM can be applied to aggregation problems of complicated kernels. The numerical results of both schemes are comparable. However, the QMOM has solved this problem more efficiently as compared to the FVS. Test Problem 2: Nucleation and Growth with Fines Dissolution. In this problem, a batch crystallizer is considered in which nucleation and growth are the dominant phenomena and is equipped with an external fines dissolution unit.30 The growth rate can be size-dependent, and a time-delay in the dissolution pipe is also included in the model. The aggregation process is neglected in this process; that is, the last two terms on the right-hand side of
(19) are absent. The initial data are taken as n(0, x, y) )
mseeds Vs(0)√2πσ
(
exp
-(x - jx)2 - (y - jy)2 2σ2
) (57)
where Vs(0) )
∫ ∫ ∞
0
∞
0
Vc
√2πσ
(
exp
)
-(x - jx)2 - (y - jy)2 dx dy 2σ2 (58)
In this problem, rectangular-shaped crystals of volume Vc ) x2y are assumed, where x is width and y is length of the crystal. Let (x0,y0) ) (0,0) and (xmax,ymax) ) (0.0025m, 0.005m). The interval [0,xmax] × [0,ymax] is subdivided into 200 × 400 grid points, and the final simulation time is 600 min. In the case of size-dependent growth rates, R1 ) 150 and R3 ) 300, while R2 ) 1 ) R4 in all cases. For the size-independent case, R1 ) 0 ) R3. The kinetic parameters and other constants are given in Table 1. The crystallizer was kept at a constant temperature of 33 °C. The following death function h(x,y) is assumed in this problem: h(x, y) )
[ (
)]
1.5 × 10-3 1 x2 y2 2.5 × 10-3 exp + 2 , σ) 2 3 σ σ √2πσ √2π (59)
In Figure 5, the plots of normalized moments for sizeindependent and size-dependent growth rates are presented. Analytical solutions are not available for this problem; therefore, results of current scheme are compared to those obtained from the FVS. Here, symbols are used for the numerical results of
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
QMOM, and lines represent the FVS results. Both schemes give overlapping results. It can be seen that fines dissolution diminishes the number of crystals; that is, µ0,0 reduces. However, the total volume, µ2,1(t), of crystals increases. It can also be observed that solute mass improves due to fines dissolution. Table 2 gives a comparison of the absolute and relative errors in mass balances of the QMOM for the size-independent growth rates. No significant changes were observed for the case of sizedependent growth rate. The table shows that QMOM needs about 5 s to solve this problem, which is much lesser than the time needed by FVS scheme for solving the same twodimensional problem. For a mesh of 200 × 200 grid points, the FVS scheme needs about 10 min when programmed in C++ software and needs hours when programmed in the Matlab software. Test Problem 3: Nucleation, Growth, and Aggregation. Suppose a stiff nucleation takes place at the minimum crystal size (x0,y0) ) (0,0) as a function of time, as given below: n(t, 0, 0) ) 0.5 + exp(-105t2)
(60)
The dimensionless particle size range is 0 e x,y e 2.0. The initial normal distribution is given as n(0, x, y) )
1
√2πσ2
(
exp -
(x - jx)2 + (y - jy)2 2σ2
)
(61)
where σ ) 0.125 and jx ) jy ) 0.4. The growth rates are chosen as G1 ) G2 ) 1.0, while the aggregation rate is considered as β ) 100. The numerical results are shown in Figures 6 and 7. The results of FVS are represented by lines, and symbols are used for the QMOM solutions. It can be observed that both schemes give comparable results. A decrease in the number of particles µ0,0 can be observed in the case of aggregation process. Conclusions A bivariate three-point QMOM was proposed for solving twodimensional batch crystallization models describing crystals nucleation, size-dependent growth, aggregation, and dissolution of small nuclei below certain critical size. In this technique, orthogonal polynomials were used to find the quadrature points and weights. In this study, a third-order orthogonal polynomial was used, which utilizes a matrix of 36-moments. The method has capability to calculate all important moments. Moreover, the proposed method is efficient due to the available analytical expression of orthogonal polynomial. Three numerical test problems with different combinations of processes were considered. The numerical results of QMOM were validated against the analytical solutions and the results of finite volume scheme. It was observed that, even in the bivariate case, QMOM can be used to model the underlying physical process and to predict its behavior efficiently and accurately. The main disadvantage of QMOM is the unavailability of complete PSD. However, in some physical processes, the PSD is not required, and lowerorder moments are sufficient to recover valuable quantities. In addition, it is also possible to reconstruct the PSD by tracking only a few lower-order moments. Acknowledgment This work was partially supported by the Higher Education Commission (HEC) of Pakistan through grant no. 1268. Literature Cited (1) Mersmann, A. Crystallization Technology Handbook, 2nd ed.; Marcel Dekker, Inc.: New York, 2001.
11643
(2) Jones, A. G.; Chianese, A. Fines destruction during batch crystallization. Chem. Eng. Commun. 1987, 62, 5–16. (3) Rohani, S.; Tavare, N.; Garside, J. Control of crystal size distribution in a batch cooling crystallizer. Can. J. Chem. Eng. 1990, 68, 260–267. (4) Zipp, G. L.; Randolph, A. D. Selective fines destruction in batch crystallization. Ind. Eng. Chem. Res. 1989, 28, 1446–1448. (5) Hulburt, H. M.; Katz, S. Some problems in particle technology. Chem. Eng. Sci. 1964, 19, 555–574. (6) Randolph, A.; Larson, M. A. Theory of Particulate Processes, 2nd ed.; Academic Press, Inc.: San Diego, CA, 1988. (7) Barrett, J. C.; Jheeta, J. S. Improving the accuracy of the moments method for solving the aerosol general dynamic equation. J. Aerosol Sci. 1996, 27, 1135–1142. (8) Madras, G.; McCoy, B. J. Reversible crystal growth-dissolution and aggregation breakage: Numerical and moment solutions for population balance equations. Powder Technol. 2004, 143-144, 297–307. (9) 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. (10) McGraw, R.; Nemesur, S.; Schwartz, S. E. Description of aerosol dynamics by the quadrature method of moments. Aerosol Sci. Technol. 1997, 27, 255–265. (11) Lim, Y. I.; Lann, J.-M. L.; Meyer, L. M.; Joulia, L.; Lee, G.; Yoon, E. S. On the solution of population balance equation (PBE) with accurate front tracking method in practical crystallization processes. Chem. Eng. Sci. 2002, 57, 3715–3732. (12) Qamar, S.; Warnecke, G.; Elsner, M. P. On the solution of population balances for nucleation, growth, aggregation and breakage processes. Chem. Eng. Sci. 2009, 64, 2088–2095. (13) Rawlings, J. B.; Witkowski, W. R.; Eaton, J. W. Modelling and control of crystallizers. Powder Technol. 1992, 69, 3–9. (14) Smith, M.; Matsoukas, T. Constant-number Monte Carlo simulation of population balances. Chem. Eng. Sci. 1998, 53, 1777–1786. (15) Tandon, P.; Rosner, D. E. Monte Carlo simulation of particle aggregation and simultaneous restructuring. J. Colloid Interface Sci. 1999, 213, 273–286. (16) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-I. A fixed pivot technique. Chem. Eng. Sci. 1996, 51, 1311–1332. (17) Dorao, C. A.; Jakobsen, H. A. Numerical calculation of the moments of the population balance equation. J. Comput. Appl. Math. 2006, 196, 619–633. (18) Dorao, C. A.; Jakobsen, H. A. A least squares method for the solution of population balance problems. Comput. Chem. Eng. 2006, 30, 535–547. (19) Dorao, C. A.; Lucas, D.; Jakobsen, H. A. Prediction of the evolution of the dispersed phase in bubbly flow problems. Appl. Math. Model. 2008, 32, 1813–1833. (20) Gunawan, R.; Fusman, I.; Braatz, R. D. High resolution algorithms for multidimensional population balance equations. AIChE J. 2004, 50, 2738–2749. (21) Qamar, S.; Elsner, M. P.; Angelov, I.; Warnecke, G.; Seidel-Morgenstern, A. A comparative study of high resolution schemes for solving population balances in crystallization. Comput. Chem. Eng. 2006, 30, 1119–1131. (22) Qamar, S.; Warnecke, G. Solving population balance equation for two-component aggregation by a finite volume scheme. Chem. Eng. Sci. 2006, 62, 679–693. (23) Diemer, R. B.; Olson, J. H. A moment methodology for coagulation and breakage problems: Part I - analytical solution of the steady-state population balance. Chem. Eng. Sci. 2002, 57, 2193–2209. (24) Gordon, R. G. Error bounds in equilibrium statistical mechanics. J. Math. Phys. 1968, 9, 655–663. (25) Barrett, J. C.; Webb, N. A. A comparison of some approximate methods for solving the aerosol general dynamic equation. J. Aerosol Sci. 1998, 29, 31–39. (26) 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. (27) 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. (28) Grosch, R.; Briesen, H.; Marquardt, W.; Wulkow, M. Generalization and numerical investigation of QMOM. AIChE J. 2006, 53, 207–227. (29) Gimbun, J.; Nagy, Z. K.; Rielly, C. D. Simultaneous quadrature method of moments for the solution of population balance equations, using a differential algebraic equation frameworks. Ind. Eng. Chem. Res. 2009, 48, 7798–7812. (30) Qamar, S.; Mukhtar, S.; Ali, Q.; Seidel-Morgenstern, A. A Gaussian quadrature method for solving batch crystallization models. AIChE J. 2010, accepted.
11644
Ind. Eng. Chem. Res., Vol. 49, No. 22, 2010
(31) 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. (32) Rosner, D. E.; Pyyko¨nen, J. J. Bivariate moment simulation of coagulating and sintering nanoparticles in flames. AIChE J. 2002, 48, 476–491. (33) Zucca, A.; Marchisio, D. L.; Vani, M.; Barresi, A. A. Validation of bivariate DQMOM for nanoparticle processes simulation. AIChE J. 2007, 53, 918–931. (34) McGraw, R. Properties and evolution of aerosol with size distributions having identical moments. J. Aerosol Sci. 1998, 29, 761–772. (35) Marchisio, D. L.; Pikturna, J. T.; Fox, R. O.; Vigil, R. D.; Barresi, A. A. Quadrature method of moments for population balance equations. AIChE J. 2003, 49, 1266–1276. ¨ ncu¨l, A. A.; The´venin, D. Technique for (36) John, V.; Angelov, I.; O the reconstruction of a distribution from a finite number of its moments. Chem. Eng. Sci. 2007, 62, 2890–2904. (37) Abegg, C. F.; Stevens, J. D.; Larson, M. A. Crystal size distributions in continuous crystallizers when growth rate is size dependent. AIChE J. 1968, 14, 118–122.
(38) Jancˇic´, S. J.; Grootscholten, P. A. M. Industrial Crystallization; Delft University Press: Delft, Holland, 1984. (39) Miller, S. M.; Rawlings, J. B.; Witkowski, W. R. Model identification and control of solution crystallization processes: A review. Ind. Eng. Chem. Res. 1993, 32, 1275–1296. (40) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes: The art of scientific computing, 3rd ed.; Cambridge University Press: New York, 2007. (41) Lushnikov, A. A. Evolution of coagulating systems III. Coagulating mixtures. J. Colloid Interface Sci. 1976, 54, 94–101. (42) Gelbard, F.; Seinfeld, J. H. Coagulation and growth of a multicomponent aerosol. J. Colloid Interface Sci. 1978, 63, 357–375.
ReceiVed for reView May 17, 2010 ReVised manuscript receiVed September 1, 2010 Accepted September 5, 2010 IE101108S
