On the representation of QMOM as a weighted-residual method—The dual-quadrature method of generalized moments

Computers and Chemical Engineering 35 (2011) 2186–2203

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

On the representation of QMOM as a weighted-residual method—The dual-quadrature method of generalized moments Paulo L.C. Lage ∗ Programa de Engenharia Química, COPPE, Universidade Federal do Rio de Janeiro, P.O. Box 68502, 21941-972 Rio de Janeiro, RJ, Brazil

a r t i c l e

i n f o

Article history: Received 14 February 2010 Received in revised form 24 March 2011 Accepted 31 May 2011 Available online 12 June 2011 Keywords: Population balance Simulation Multiphase flow Multiphase reactors Quadrature-based methods Weight-residual methods

a b s t r a c t Quadrature-closed moment-based methods for the solution of the population balance equation were analyzed. The QMOM was shown to be a particular case of a method based on generalized moments (QMoGeM). Then, a weighted-residual method (WRM) based on the generalized moments was derived (WRMoGeM). If it is closed by the Gauss–Christoffel quadrature, the resulting method (QWRMoGeM) was shown to be identical to the QMoGeM, giving the correct representation of QMOM as a WRM. The WRMoGeM formulation was used to derive a new method (DuQMoGeM) that employs two quadrature rules, one for discretizing the particulate system and other to accurately integrate the integrals in the equations for the generalized moments. A Galerkin version of this method was implemented and used to solve several examples with known analytical solutions. The DuQMoGeM solutions for breakage and aggregation problems were shown to be more accurate than the QMOM solutions. © 2011 Elsevier Ltd. All rights reserved.

1. Introduction Dispersed two-phase flows occur in a large number of engineering problems and a way to model particulate system behavior is through the population balance framework (Ramkrishna, 2000). There has been a considerable interest in the development of numerical methods to solve the population balance equation (PBE) in the past few years (Bove, Solberg, & Hjertager, 2005; Dorao & Jakobsen, 2006a, 2007; Grosch, Briesen, Marquardt, & Wulkow, 2007; Marchisio & Fox, 2005; McGraw, 1997; Su, Gu, Li, Feng, & Xu, 2008, among many others). Most of these works developed variants of the method of moments closed by quadrature (QMOM) and the others applied the least-squares methods (LSM) to solve population balance problems. A review of the existing methods is not in the aim of the present work, but some comments are given in the following. The method of moments closed by quadrature, QMOM, introduced the usage of an approximate closure of the moment equations by using the n-point Gauss–Christoffel quadrature, that is, a Gaussian quadrature whose weight function is the particle number distribution function (McGraw, 1997). The Gauss–Christoffel quadrature was calculated by the product-difference algorithm, PDA, developed by Gordon (1968). Afterwards, Marchisio and

∗ Tel.: +55 21 2562 8346; fax: +55 21 2562 8300. E-mail address: [email protected] 0098-1354/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2011.05.017

Fox (2005) developed the direct QMOM, DQMOM, in which the moments of the PBE are solved directly for the abscissas and weights of the Gauss–Christoffel quadrature, avoiding the computational cost of the quadrature calculation during the solution. The Gauss–Christoffel quadrature leads to a straightforward discretization of the internal coordinate domain. Considering the particle size as the internal coordinate, this discretization generates n particulate phases with different particle sizes in an Eulerian multiphase model (Silva & Lage, 2007) for a polydisperse two-phase flow. This makes this method quite amenable to be implemented in the existing CFD codes (Fan, Marchisio, & Fox, 2004; Marchisio, Vigil, & Fox, 2003; Silva, Damian, & Lage, 2008; Silva & Lage, 2007). However, QMOM/DQMOM have some drawbacks. First, they do not give a representation of the distribution function. Besides, the computation of the Gauss–Christoffel quadrature from the regular moments is usually an ill-conditioned problem (Gautschi, 2004), which limits the number of points in the quadrature rule that can be calculated with accuracy. But the Gauss–Christoffel quadrature with only a few points is not accurate enough to represent the terms that appears in the moments of the PBE. This can lead to very inaccurate solutions (Dorao & Jakobsen, 2006b, 2006c) and even loss of positiveness of the measure defined by the moments of the distribution function (Petitti et al., 2010). Moreover, they do not handle well problems where two abscissas converge to each other or when there is a particle efflux at a boundary, as in a dissolution or evaporation problem. Several modifications of QMOM or DQMOM were developed in the past years trying to overcome these difficulties. Alopaeus,

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

Nomenclature a A B b c c D f G g H L m P(x | x ) P q R S t x xmax y

aggregation frequency three-dimensional aggregation matrix particle birth term breakage frequency coefficient in polynomial approximation of f vector of f expansion coefficients particle death term particle number distribution function, weight function of measure d = f(x)dx growth matrix the particle growth rate the Heaviside step function breakage matrix vector of generalized moments probability distribution of daughter property values, x, upon breakage of a particle of property x matrix that transforms c into m moment of growth term PBE residual particle source term time particle property maximum value of particle property polynomial approximation

Greek letters ı(·) Dirac delta function ıij Kronecker delta ¯ a known measure with a positive definite inner d product ¯ weight function of measure d(x) = ω(x)dx ω  dimensionless zero order moment of solution given by Eq. (53)  test function or particle property number of daughters upon breakage ϑ Subscripts a related to aggregation b related to breakage Superscripts () relative to property  ¯ – relative to measure d a analytical solution P relative to function P

Laakkonen, and Aittamaa (2006a) used a quadrature with fixed abscissas and weights calculated from the moments of the distribution function to avoid the problems associated to the calculation of the Gauss–Christoffel quadrature, calling this modification as the fixed QMOM, FQMOM. Similarly, Su, Gu, and Xu (2009) reported the fixed-point QMOM, FPQMOM, a QMOM that uses fixed-point Gaussian quadratures to establish the position of the abscissas using also the moments to calculate the quadrature weights. Following Yoon and McGraw (2004), Attarakih, Drumm, and Bart (2009) used a two-point equal-weight quadrature with the abscissas and common weight determined from three moments of the distribution, avoiding the converging abscissa problem of the Gauss–Christoffel quadrature. Regarding the representation of the distribution function, there have been several attempts. Strumendo and Arastoopour (2008) used an expansion of the distribution function in a complete set

2187

of trial functions, that is, a basis of the space of functions defined in a finite domain on the real line was employed to represent the distribution. They called the method as FCMOM (Finite size domain Complete set of trial functions Method Of Moments). They assumed the PBE to be valid on a finite domain, which implies model modification in the case of aggregation. In order to cope with growth problems, they used moving boundaries to transform the PBE to the fixed domain for which the basis is chosen. They basically used the Legendre polynomial basis in the [−1,1] range. The moving boundary technique can largely improve the approximation of the distribution function for growth problems. The method solves for N regular moments of the transformed PBE and uses their values to obtain the coefficients of the distribution function expansion using the Legendre polynomials, which is used to approximate the aggregation and breakage terms present in the governing equations. Therefore, both the regular moments and the expansion coefficients are present in the final form of the governing equations. Several problems with growth were solved with good results due to the moving boundary technique. Problems with pure aggregation and simultaneous growth and aggregation were also solved with good representation of the distributions. Errors in the moments of the distribution were not analyzed. No breakage problem was solved. No discrete characterization of the distribution was sought. Moment-preserving methods give good approximations of the main characteristics of the distributions but they do not give a good representation of the distribution because only a few moments can accurately be solved for. On the other hand, the classical sectional discretization method (Kumar & Ramkrishna, 1996) preserves only two moments but give reasonably good representation of the distribution. So, it is easy to understand the development of discretization methods that preserve more than two moments (Alopaeus, Laakkonen, & Aittamaa, 2006b) and the derivation of a sectional form of QMOM by Attarakih et al. (2009), called SQMOM. In SQMOM, the internal variable domain is divided in sections which are represented by the so-called primary particles. A moment-preserving quadrature is calculated from the sectional moments and used to establish the characteristics of the secondary particles. The assumption of the preservation of two moments allows the determination of the primary particles characteristics from those of the secondary particles. Thus, the method solves for the regular sectional moments of the PBE in all sections and the breakage and aggregation terms in the governing equations are approximated by the moment-preserving quadrature used in each section. Examples involving breakage and aggregation in open and closed vessels were solved and the results were quite good for the first 4 moments (using a 2 point quadrature) but large numbers of sections were necessary to represent the distributions. No growth problem was solved. Recently, Massot, Laurent, Kah, and de Chaisemartin (2010) applied the theory of the canonical moments to solve the Hausdorff finite moment problem and reconstruct the distribution function using a maximum entropy principle. They were mainly interested in the pure growth PBE with a constant negative growth rate (vaporization). An integral form of the PBE was derived and used to develop a modified DQMOM algorithm using the method of characteristics. The moment effluxes at a given boundary were calculated from the reconstructed distribution. The method can be applied sectionally. The results for problems with constant and smooth-varying evaporation laws were quite good, with the first few moments showing relative errors of a few percent. The least-squares method is robust but computationally intensive, requiring iterative solution if a non-linear term is present, as the aggregation terms of the PBE. It uses a functional expansion of the number density function in terms of the internal coordinates of the PBE and it can also use time-space functional expansions (Dorao & Jakobsen, 2006a, 2007). The method relies on the min-

2188

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

imization of a functional which always include the norm of the residual of the PBE. Therefore, it gives the “best” approximated solution for the considered functional expansion. Recently, Dorao and Jakobsen (2008) developed the least-square spectral method (LSSM) and Zhu, Dorao, Lucas, and Jakobsen (2009) applied it to solve the PBE coupled to an one-dimensional steady-state two-fluid model for a bubble column. Therefore, the application of LSSM to the PBE solution as part of a multifluid model requires that the whole model is solved by LSSM. Dorao and Jakobsen (2006b, 2006c) derived a relationship between the classical QMOM (McGraw, 1997) and the well-known method of weighted residuals (WRM). These MOM–WRM were then compared to QMOM (Dorao & Jakobsen, 2006c) and to LSM (Dorao & Jakobsen, 2006b) by solving integral-differential equations that resemble breakage problems. They stated that QMOM is a case of one of their MOM–MWR based on a Lagrange interpolant polynomial. Unfortunately, this is incorrect. In the present work, the correct relationship between the method of moments closed by quadrature and the method of weighted residuals are derived using generalized moments. The resulting weighted-residual method based on generalized moments can be further extended to a dual-quadrature method, which was called the dual-quadrature method of generalized moments, DuQMoGeM. In this method, different quadrature rules can be used for the two relevant tasks in such problems. A quadrature rule of high accuracy is used to calculate the integrals that are needed to close the system of equations for the generalized moments. The Gauss–Christoffel quadrature that uses the particle number distribution function as its own weight function is used for discretizing the particle density function, generating representative particulate “phases” that can be employed in PB-CFD simulations. Therefore, the method accuracy and the number of particulate “phases” are given by different quadrature rules. This article was organized as follows. The form of the population balance equation used in this work was given in Section 2. In Section 3, the quadrature method of generalized moments (QMoGeM) was derived and simplified to the classical QMOM. In Section 4, the weighted-residual method based on generalized moments (WRMoGeM) was derived and its equivalence to the QMoGeM was deducted. The dual-quadrature weighted-residual method of generalized moments (DuQMoGeM) was introduced in Section 5. The numerical solution of several problems that have analytical solutions were presented in Section 6. From the analysis of the results, conclusions were drawn in Section 7.

stituted by xmax without incurring in any approximation of the PBE. Further simplification was achieved by assuming binary aggregation and using a x property that is conserved during aggregation and breakage, as the particle mass. Although the usage of the method of characteristics or a moving boundary technique may improve the method developed in this work for growth problems, the particle growth phenomenon was included in the basic formulation given below. Thus, the resulting simplified form of the PBE is given as (Ramkrishna, 2000): ∂f (x, t) ∂ + [g(x, t)f (x, t)] = Ba (x, t) − Da (x, t) + Bb (x, t) ∂t ∂x − Db (x, t) + S(x, t)



1 2

Ba (x, t) =

x

a(x − x , x , t)f (x − x , t)f (x , t) dx



xmax

Da (x, t) =



(2)

0

0 xmax

Bb (x, t) =

a(x, x , t)f (x, t)f (x , t) dx

(3)

ϑ(x , t)b(x , t)P(x|x , t)f (x , t) dx

(4)

x

Db (x, t) = b(x, t)f (x, t)

(5)

In the above equations, 

P(x|x , t) = 0,



∀x > x



x



P(x | x ,

t) has the following properties:



xmax

P(x|x , t) dx = 0

(6) P(x|x , t) dx = 1,

(7)

0



x 

xP(x|x , t) dx = 0

xmax

xP(x|x , t) dx =

0

x ϑ(x , t)

(8)

Eq. (8) guarantees that property x is conserved during particle breakage. Defining the breakage and aggregation operators as: Lb f (x, t) ≡ Db (x, t) − Bb (x, t)

(9)

La f (x, t) ≡ Da (x, t) − Ba (x, t)

(10)

the residual of Eq. (1) can be written as: ∂f (x, t) ∂ + [g(x, t)f (x, t)] + Lb f (x, t) ∂t ∂x

R(x, t; f ) ≡ 2. The population balance equation A simple form of the population balance equation, PBE, was considered. It can be applied to an isolated population of particles that is homogeneous in space and characterized by a single property, which is usually the particle mass or volume. Although the range of values of this property is usually from zero to infinity, it has indeed a finite maximum because particle processes are usually physically limited. For instance, for a disperse twophase flow in a duct, aggregation or growth cannot produce particles larger than the duct diameter. These limitations must be embodied in the rates of the particle processes. However, for some theoretical problems, the range [0, ∞) must be kept for the particle property. In the following, a dimensionless coordinate, x, in a finite or semi-finite range will be employed. In order to cope with these two cases, the particle number distribution function, f(x, t) was defined for x ∈ [0, xmax ), xmax =∞ or a finite vale, and t ∈ [0, ∞). It should be noted that all derivations were performed for the [0, ∞) range and, when physical limitations exist, the upper limits of the integration were sub-

(1)

where g is the growth rate, S is a source function, for instance, due to nucleation, and the aggregation and breakage terms are given by:

+ La f (x, t) − S(x, t) = 0

(11)

It should be noted that the aggregation operator, Eq. (10), is nonlinear. 3. The quadrature method of generalized moments (QMoGeM) 3.1. The generalized moment of the PBE Similarly to any distribution function, only the integrals of the particle number distribution function, f(x, t), have some physical meaning for the particle population. Therefore, it is interesting to consider numerical methods that solve for some integrals of f instead of solving for f itself. Considering that j (x), j = 0, . . ., 2n − 1 are important particle properties, the corresponding important population properties are:



()

j (t) =

xmax

j (x)f (x, t) dx = (j , 1)d , 0

j = 0, . . . , 2n − 1 (12)

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

where the inner product used in the rightmost term is defined by:



xmax

(u, v)d =

u(x)v(x) d(x)

(13)

0

for any functions u and v, with the measure d(x) defined as d(x) = f(x, t) dx. () Evolution equations for j (t) can be directly obtained by operating Eq. (1) with to:

 xmax 0

j (x) · dx. After some manipulation, this led



()

dj

3.3. The QMOM g(x, t)˙ j (x)f (x, t) dx

+ [g(x, t)f (x, t)j (x)]0max −



xmax





xmax

+ 0



0

j (x) −

0

()

If j (x) = xj , j become the regular moments, j , of f. In this case, Eq. (18) become:



1  (x + x ) a(x, x , t)f (x, t)f (x , t) dx dx 2 j

xmax

+

known or assumed to be zero. For the upper size limit, the regularity condition requires that there is no particle flux, g(xmax , t)f(xmax , t) = 0. At the lower size boundary, the particle flux g(0, t)f(0, t) can be known for nucleation with particles of zero size. However, particles may disappear at the lower size boundary with a non-zero flux, that is, g(0, t) = / 0 and f(0, t) is unknown. This problem occurs for the case of constant negative growth rate, which cannot be solved by QMoGeM (or QMOM).

xmax

x

dt

dj dt



+

xmax

j (x)S(x, t) dx



xmax



+



j (x )P(x |x, t) dx =

0

x

j (x )P(x |x, t) dx

(15)





1 ˜ ( , a)d(x) , 1 2 j

() + (b, j − ϑ j ) d(x)

d(x )

= (j , S)dx

(16)



j

wi wk i −



1 j ( + k ) a( i , k , t) 2 i

n 

j

wi b( i , t)[ i − ϑ( i , t)Pj ( i , t)] = (xj , S)dx

(19)

where



xmax

Pj (x, t) =

 x P(x |x, t) dx =

0

x

x P(x |x, t) dx

j

j

(20)

0

4. The weighted-residual method based on generalized moments (WRMoGeM) A weighted residual method (WRM) consists of obtaining an approximation for a given equation by making its residual orthogonal to a given finite set of test functions, j :

˜ j = j (x + x ) and where  xmax

(u, v)dx =

j−1

wi g( i , t)j i

recovering the QMOM formulation (McGraw, 1997). Again, the xmax fluxes in the term [g(x, t)f (x, t)xj ]0 must be known or assumed zero.

x

+ (j , a)d(x) −



n

0

+ [g(x, t)f (x, t)j (x)]0max − (g, ˙ j )d(x)

dt

n 

i=1

Using again the inner product defined by Eq. (13), Eq. (14) can be rewritten as: d(j , 1)d

−

i=1 k=1

(14)

0

where ˙ j (x) = dj (x)/dx and

 n

0

j (x, t) =

xmax

+ [g(x, t)f (x, t)xj ]0

i=1



b(x, t)[j (x) − ϑ(x, t) j (x, t)]f (x, t) dx

=

2189

u(x)v(x) dx

(17)

0

Regardless of the actual properties taken into consideration, j , it is quite clear from Eq. (16) that the Gaussian quadrature approximation based on the measure d(x) = f(x, t) dx is well suited for approximating all integral terms except the source term. In the following, it is assumed that (j , S)dx is known by analytical or numerical calculation.

(R, j ) = 0,

j = 0, . . . , N − 1

(21)

N−1 using an adequate inner product. If the set {j }0

represents a basis of an N-dimensional subspace, then Eq. (21) guarantees that the residual of the equation is zero in this subspace. Of course, f must also be approximated, using a functional expansion in the same or another basis, to carry on the method. 4.1. The order of approximation in WRMoGeM

3.2. The QMoGeM The first 2n generalized moments can be used to determine the abscissas, k , and weights, wk , of the n-point Gaussian quadrature based on the measure d(x) = f(x, t) dx (Gautschi, 1994). This quadrature rule can be used to simplify Eq. (14) to ()

dj

x

+ [g(x, t)f (x, t)j (x)]0max −

dt

n 

wi g( i , t)˙ j ( i )

i=1

 n

+

n



wi wk j ( i ) −



1  ( + k ) a( i , k , t) 2 j i

i=1 k=1

+

n 



wi b( i , t)[j ( i ) − ϑ( i , t) j ( i , t)] = (j , S)dx

(18)

i=1

Since QMoGeM does not provide a functional approximation for the x distribution, the moment fluxes in [g(x, t)f (x, t)j (x)]0max must be

Dorao and Jakobsen (2006c) tried to formulate the QMOM as a weighted residual method. They chose a pure breakage problem because the PBE becomes linear in f, which simplifies the deduction. However, they made a serious mistake: they assume that the number density function, f(x), was approximated by the (n − 1)-order Lagrange interpolating polynomial based on the n abscissas of the Gauss–Christoffel quadrature. As the n-point Gaussian quadrature of a given integrand corresponds to its approximation by the (2n − 1)-order Hermite interpolating polynomial (see Theorem 1.48 of Gautschi, 2004), the order of the f expansion in the WRM version of QMOM should be 2n − 1. This can be easily proved by assuming that f(x) is expanded using the basis of polynomials, ϕn (x), n = 0, 1, . . ., which are orthogonal in the inner product defined for a known absolute continuous ¯ measure d(x) = ω(x)dx, that is:



xmax

ϕi (x)ϕj (x)ω(x) dx = ıij ||ϕi ||2 ¯

(ϕi , ϕj )d¯ ≡ 0

d

(22)

2190

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

The f expansion and its (2N − 1)-order approximation are then written as: ∞ 

f (x) = ω(x)

ci ϕi (x) ⇒ y2n−1 (x) = ω(x)

ϕi (x)j (x)ω(x) dx

(27)

ci ϕi (x)

(23) x Gji = [g(x, t)j (x)ω(x)ϕi (x)]0max − (g ˙ j , ϕi )d¯ x

(ϕi , f )dx

(ϕi , 1)d ||ϕi ||2 ¯



1

=

||ϕi ||2 ¯ d

||ϕi ||2 ¯ d

= [g(x, t)j (x)ω(x)ϕi (x)]0max



xmax

ϕi (x)f (x) dx 0

g(x, t)˙ j (x)ϕi (x)ω(x) dx

(28)

0

i

(24)

||ϕi ||2 ¯



Lji = (j − j , bϕi )



d

are the generalized moments obtained by using the ¯ basis of polynomials that are orthogonal in the measure d(x). Due to the degree of exactness of 2n − 1 of the n-point Gaussian quadrature for the measure d(x) = f(x) dx, Eq. (24) can be exactly calculated by n 

||ϕi ||2 ¯



[j (x) − (x) j (x)]b(x)ϕi (x)ω(x) dx

˜ Ajik = (j (a, ϕk )d(x ¯  ) − 0.5(j a, ϕk )d(x ¯  ) , ϕi )



xmax



0

(25)

xmax

• for a finite interval, x ∈ [0, xmax ], ω(x) = 1, leading to the Legendre polynomials shifted to the given interval, • for a semi-infinite interval, x ∈ [0, ∞), ω(x) = xa e−x , which implies the generalized Laguerre polynomials, and • for an infinite interval, x ∈ (− ∞ , ∞), ω(x) = e−a2 x2 , which results in the Hermite polynomials.

¯ d(x)

[j (x) − 0.5j (x + x )]

0

× a(x, x )ϕi (x)ϕk (x )ω(x)ω(x )dx dx

d j=1

for i = 0, 1, . . ., 2n − 1. Therefore, the first 2n terms of the series expansion of f(x) can be exactly calculated using the 2n values of wj and j , j = 1, . . ., n of the n-point Gaussian quadrature, which can be calculated by the first 2n generalized moments of f(x) (Gautschi, 2004). Thus, the WRM equivalent to QMOM must use a polynomial approximation of (2n − 1) order. It is also clear that the best measure to be used in Eqs. (22) and (23) is the one for each ω(x) best approximates f(x, t), at least ¯ for a given time. However, as the measure d(x) = ω(x)dx must be known to allow the computation of the Gauss quadrature for the measure d(x), the usual choices in Eq. (22) are:

(29)

0

= wj ϕi ( j )

¯ d

xmax

=

(ϕ)

where i

1

xmax

−

(ϕ)

=

d

ci =

xmax

0

i=0

where

=



Pji = (j , ϕi )d¯ =

2n−1

i=0

ci =



where

and



(30)

xmax

sj = (j , S)dx =

j (x)S(x, t) dx

(31)

0

The first term of Eq. (26) was written using the results of the application of (j , ·)dx to Eq. (23):



2n−1 ()

j

= (j , 1)d = (j , f )dx =

Pji ci

(32)

i=0

or, in matrix form m = Pc

(33)

where m =

() [j ],

being

() j

the generalized moment using the

test function j . Therefore, Eq. (26) can be written as: ()

dj

dt



2n−12n−1

+

(Lji + Gji )Pil−1 l

()

l=0 i=0

It should be pointed out that the WRM using Eqs. (22) and (23) gives the best polynomial approximation of f(x)/ω(x) in the sense of the least-squares, that is, the one that gives the minimum of ||f (x)/ω(x) − yn (x)/ω(x)||2dx (Hildebrand, 1987). It should also be clear that the inner product that has to be used in Eq. (21) is the one defined for the measure dx.



2n−12n−12n−12n−1

+

−1 −1 Ajik Pkr Pil r l ()

()

= sj

(34)

r=0 l=0 i=0 k=0

for j = 0, 1, . . ., 2n − 1. If A = [Ajik ] is taken as a three-dimensional matrix 2n × 2n × 2n, Eq. (34) can be written in the following matrix form:

4.2. The WRMoGeM

dm + (L + G)P−1 m + {[(AP−1 )P−1 ]m}m = s dt

In this section, the weighted residual method based on the generalized moments was derived. Then, its reduction to the conventional QMOM was deducted. Differently from Dorao and Jakobsen (2006c), the PBE for simultaneous growth, breakage and aggregation, as given by Eq. (1), was considered. Eq. (14) can be obtained from Eq. (11), using the inner product (j (x), ·)dx . Then, the following equation can be deducted by substituting f by its approximation y2n−1 , after some manipulation:

which represents the weighted residual method based on the generalized moments. The application of the WRM to the solution of the PBE is not new. Ramkrishna (2000) described the earliest works. Usually, it is given by Eq. (26) which has to be solved for the coefficients of the series expansion of f. One of the drawbacks of the WRM method is the need of calculating the integrals present in the definition of Gji , Lji and Ajik at each time step when the growth rate and the breakage and aggregation functions depend on time, usually through a time-varying continuous-phase property. On the other hand, the computation of Pji is made only once if all j and ϕk are not functions of time.

d dt

2n−1 

Pji ci

i=0





2n−1

+

i=0

Gji ci +

 i=0



2n−12n−1

2n−1

Lji ci +

i=0 k=0

Ajik ck ci = sj

(26)

(35)

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

2191

4.3. The QWRMoGeM

, or on the choice of test the specific polynomial basis used, {ϕk }2n−1 0

The quadrature-closed weighted residual method based on the generalized moments consists of calculating the Gaussian quadrature for the measure d using the 2n generalized moments (Gautschi, 1994) and its use to evaluate the integrals in Eq. (34):

functions, {j }0 . If j = xj , the usual form of QMOM, Eq. (19), is obtained. In order to compare the present formulation with that given by Dorao and Jakobsen (2006c), consider Eq. (26) without aggregation or growth, with Pji time independent and in matrix form:

2n−1





2n−12n−1

2n−1

Gji Pil−1 l

()

=

l=0 i=0



2n−1

Gji ci =

l=0

dc + P−1 Lc = P−1 s dt

x

{[g(x, t)j (x)ω(x)ϕi (x)]0max

l=0

or using Eq. (33):

x

− (g ˙ j , ϕi )d¯ }ci = [g(x, t)j (x)f (x, t)]0max



2n−1

−

l=0

−

n 

ω x (g ˙ j , ϕi ) ci = [g(x, t)j (x)f (x, t)]0max f d ω( k )  ci ϕi ( k ) f ( k , t) 2n−1

wk g( k , t)˙ j ( k )

k=1



The above equations should be compared to Eqs. (33) and (43) of Dorao and Jakobsen (2006c) keeping in mind that the above vectors and square matrices are 2n-dimensional.

Eq. (34) is simplified if the choice j = ϕj , ∀ j is made. In this case ()

wk g( k , t)˙ j ( k )

(36)

j

dt





2n−1

Lji Pil−1 l

()

=

l=0 i=0



Lji ci =



(j − j , bϕi )

¯ d

ci

n



(j − j , b

l=0

k=1

ω( k )  ci ϕi ( k ) f ( k , t) 2n−1



− ( k ) j ( k ))b( k )

||ϕi ||2 ¯

n



wk (j ( k ) − ( k ) j ( k ))b( k )



Ajik

i=0 k=0

||ϕi ||2 ¯ ||ϕk ||2 ¯

2n−12n−1 (ϕ)

i +

d

d

d

(ϕ)

(ϕ)

k i

= sj (41)

d

||ϕj ||2 ¯

d

dcj dt





2n−12n−1

2n−1

+

(Lji + Gji )ci +

i=0

Ajik ck ci = sj

(42)

i=0 k=0 (ϕ)

l=0



 Lji + Gji

= ci ||ϕi ||2 ¯ , that leads to:

for j = 0, 1, . . ., 2n − 1. This system of equations can also be expressed in terms of the expansion coefficients ci , in the following form:

 ω ϕi ) ci = wk (j ( k ) f d

2n−1

=

d

i=0

l=0



(ϕ)

and Pji = ıji ||ϕi ||2 ¯ , giving i

2n−1

+

2n−1

l=0

=

(ϕ)

= j

(ϕ)

dj

k=1 2n−12n−1

(40)

4.4. The Galerkin form of WRMoGeM

l=0

n

−

dm + LP−1 m = s dt

x [g(x, t)j (x)f (x, t)]0max

=

(39)

(37)

The ci usually have magnitudes much smaller than those of i , making the numerical integration of Eq. (42) more accurate than that of Eq. (41). This form is quite suitable for generating a dual quadrature method, as explained in the next section.

k=1

5. The DuQMoGeM

and





2n−12n−12n−12n−1

2n−12n−1

−1 −1 Ajik Pkr Pil r l ()

()

=

r=0 l=0 i=0 k=0



i=0 k=0

2n−12n−1

=

The dual-quadrature weighted-residual method of generalized moments consists of using the (2n − 1)-order polynomial expansion for f, given by Eq. (23), and the WRMoGeM form of the PBE, Eq. (41) or (42), to solve for the first 2n gener-

Ajik ck ci

(j (a, ϕk )d(x ¯ ) −

1 ˜ ( a, ϕk )d(x ck ci ¯  ) , ϕi ) ¯ 2 j d(x)

(ϕ) 2n−1

i=0 k=0



2n−12n−1

=

(j (a,

i=0 k=0

=



n n  

ω ω 1 ˜ ω ϕ ) − ( a, ϕ ) ,ϕ ) c c f k d(x ) 2 j f k d(x ) i f d(x) k i

 wp wr

ω( r )  ϕk ( r )ck f ( r , t) 2n−1

n n  

 wp wr

(38)

1 j ( p ) − j ( p + r ) a( p , r ) 2

p=1 r=1

k=0

=





ω( p )  ϕi ( p )ci f ( p , t) 2n−1

×

, of the polynomial basis {ϕk }2n−1 , 0 ¯ for which the recursion coeffiorthogonal in the measure d ¯ k }∞ are known analytically and easily computed (for ¯ k, ˇ cients {˛ 0 instance, by the drecur routine of the ORTHPOL package, Gautschi, ¯ k }N−1 is all ¯ k, ˇ 1994). It should be noted that knowledge of {˛ 0 ¯ k , and abscissa, ¯ k , of that is needed to compute the weights, w ¯ measure (ORTHPOL the N-point Gaussian quadrature for the d alized moments, {k }0

i=0



1 j ( p ) − j ( p + r ) a( p , r ) 2

p=1 r=1

where the (2n − 1)-order approximation of f, Eq. (23), was used. If the results of Eqs. (36), (37) and (38) are substituted into Eq. (34), Eq. (18) is obtained. This proves that the QWRMoGeM should give identical results to the QMoGeM as long as both methods can be applied. It is interesting to note that this result does not depend on

dgauss routine, Gautschi, 1994) and also to compute the value of {ϕk }N . 0 (ϕ) 2n−1

From {k }0

2n−2

¯ k} and {˛ ¯ k, ˇ 0

, the recursion coefficients

n−1 {˛j , ˇj }0

can be calculated using the modified Chebyshev algorithm (ORTHPOL dcheb routine, Gautschi, 1994) from which the n n-point Gaussian quadrature for the measure d, {wk , k }1 , is easily computed (ORTHPOL dgauss routine, Gautschi, 1994). This gives a ¯ k }2n−1 are discretization of the particle density function. Besides, {ˇ 0

needed to calculate the norms ||ϕi ||2 ¯ , i = 0, . . . , 2n − 1. Therefore, d

for a discretization of the particle system into n phases, N ≥ 2n. From ¯ k }N−1 , the N-point Gaussian quadrature for the the values of {˛ ¯ k, ˇ 0 ¯ measure can also be easily computed. d

2192

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

a

10-11

-10

a

10

Pure Growth - Case 1 DuQMoGeM (n=2, N=4)

-11

-12

10

1 2 3

10-14

(φ),a

10-13

k

10-15

10-16 0.0

10-12

|μk(φ)-μ(φ),a |/|μk k

a a |μk -μk | / | μk |

|

10

-13

4.0

6.0

8.0

10.0

12.0

14.0

1 2 3

-14

10

Pure Growth - Case 1 QMOM ( n =2) 2.0

-15

10

0.0

16.0

2.0

4.0

6.0

8.0

0

b ξ1 ξ2 w1 w2

ξk(t) or w k(t)

10-1

10-2

10

10-2

4.0

6.0

8.0

14.0

16.0

10-3

Pure Growth - Case 1 DuQMoGeM - (n=2, N=4)

10

-4

0.0

2.0

4.0

6.0

8.0

-4

2.0

16.0

ξ1 ξ2 w1 w2

10-1

-3

0.0

14.0

0

Pure Growth - Case 1 QMOM (n=2) 10

12.0

10

ξk(t) or w k(t)

10

10.0

t

t

b

k

10

10.0

12.0

14.0

16.0

t Fig. 1. Case 1 with pure growth. QMOM solution for n = 2: (a) relative errors in the regular moments and (b) abscissas and weights of the Gauss–Christoffel quadrature.

10.0

12.0

t Fig. 2. Case 1 with pure growth. DuQMoGeM solution for n = 2, N = 4: (a) relative errors in the generalized moments and (b) abscissas and weights of the Gauss–Christoffel quadrature. ϕ

Lji = (ϕj − j , bϕi ) The integrals present in Gij , Lji and Ajik are inner products in ¯ measure that can, in principle, be calculated to any desired the d accuracy by the corresponding Gaussian quadrature by choosing a convenient value for N. Thus, n controls the number of moments solved for (2n) and the accuracy of the representation of the distribution function (2n term expansion) used in the weighted-residual formulation whereas N (N ≥ 2n) controls the accuracy of the resulting integrals in the moment equations. The moment flux at the boundaries of the internal variable domain can be evaluated by the approximation of the distribution function given by Eq. (23), if they are not zero. This is the only dependency of this method on the functional approximation of the distribution function which enables the DuQMoGeM to solve problems with non-zero moment efflux that cannot be solved by QMOM. Except for this case, the DuQMoGeM is a generalized moment method that does not depend on the approximation of the distribution function given by Eq. (23), which is an additional feature that can represent sufficiently smooth distribution functions. The approximation converges in the mean to the distribution function if it belongs to L2 [0, xmax ], and convergence is uniform if it belongs to C[0, xmax ]. Therefore, Gji = [g(x, t)ϕj (x)ω(x)ϕi (x)]x0max − (g ϕ˙ j , ϕi )d¯ x

= [g(x, t)j (x)ω(x)ϕi (x)]0max −

N  k=1

¯ k g( ¯ k , t)ϕ˙ j ( ¯ k )ϕi ( ¯ k ) w

(43)

=

N 

¯ d

ϕ ¯ k [ϕj ( ¯ k ) − ( ¯ k , t) j ( ¯ k , t)]b( ¯ k , t)ϕi ( ¯ k ) w

(44)

k=1

Ajik =

=



ϕj (a, ϕk )d(x ¯ ) −

N N  

1 (ϕ ˜ a, ϕk )d(x ¯  ) , ϕi 2 j



 ¯ d(x)



1 ¯ lw ¯ r ϕj ( ¯ l ) − ϕj ( ¯ l + ¯ r ) w 2

l=1 r=1

× ϕi ( ¯ l )ϕk ( ¯ r )a( ¯ l , ¯ r , t)

(45)

for i, j, k = 0, 1, . . ., 2n − 1. The evaluation of sj can also use the same quadrature as follows: sj =

ϕ

j

ω

 ,S

¯ d

=

N  k=1

¯k w

ϕj ( ¯ k ) S( ¯ k , t) ω( ¯ k )

(46)

Even though it is convenient to use the Gaussian quadrature in ¯ for the above integrations, an adaptive quadrature the measure d scheme with error control can be used in its place. It should be pointed out that if the distribution function does not belong to the L2 space (as in the case of distributions represented by Dirac delta functions), the functional expansion is lost. However, the method still solves for the scaled generalized moments and it is able to obtain the Gauss–Christoffel quadrature as a discrete representation of the distribution.

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

102

a

0

(φ),a

Pure Growth - Case 2 DuQMoGeM (n=4, N=8)

0.7

10-2

(φ),a (φ)

10

10

k

-8

10

-10

10

-12

0.0

ξk(t)

|/|μk

0.6 -4

10-6

|μk -μk

ξ1 ξ2 ξ3 ξ4 analytical numerical

0.9 0.8

|

10

c

2193

Pure Growth - Case 2 DuQMoGeM (n=4, N=8) 0.2

0.4

0.6

0.8

1.0

0 1 2 3 4 5 6 7 1.2

0.5 0.4 0.3 0.2 0.1 0.0

1.4

0.0

1.6

0.2

0.4

t 2.5

Pure Growth - Case 2 2.0

f(x,t)

1.5 1.0

d

t 0.0 0.4 0.8 1.2 1.6 numerical analytical

0.3 0.2

0.0

0.1

DuQMoGeM (n=4, N=8) 0.2

0.4

0.6

1.0

w1 w2 w3 w4 analytical numerical

0.4

0.5

-0.5 0.0

0.8

Pure Growth - Case 2 DuQMoGeM (n=4, N=8)

0.5

w k(t)

b

0.6

t

0.0 0.8

1.0

x

0.0

0.2

0.4

0.6

0.8

1.0

t

Fig. 3. Case 2 with pure growth. DuQMoGeM solution for n = 4, N = 8: (a) relative errors in the moments, (b) distribution function evolution and (c) abscissas and (d) weights of the Gauss–Christoffel quadrature.

The DuQMoGeM has also some similarities to the FCMOM (Strumendo & Arastoopour, 2008) because both use functional expansion of the distribution function and solve for the moments of the PBE. However, DuQMoGeM does not, a priori, truncate the domain and it solves for the scaled generalized moments that are the coefficients of the functional expansion. 6. Results The results of the application of the methods discussed above for several examples are presented in the following. All examples are transient problems with known analytical solutions. They were numerically integrated in time using the DASSL routine (Petzold, 1982) with given absolute ( a ) and relative ( r ) tolerances. Typically, a = 10−14 − 10−30 and r = 10−11 − 10−13 . Polynomial recursion coefficients and quadrature rules were calculated using routines of the ORTHPOL package Gautschi (1994). For the DuQMoGeM solutions, the N-point Gaussian quadrature for ¯ measure was used for convenience, choosing N as large as the d needed. All computations were performed in standard double precision using FORTRAN (GNU Fortran 4.1.2 in a Linux Intel Q6600 platform). Some CPU times are reported in the following for typical runs that generated the output files for analysis. The reported values were evaluated by averaging the CPU times of 10 runs using the same values of a and r in DASSL integration for both

methods. However, the programs implementing the QMOM and DuQMoGeM were not fully optimized and the DuQMoGeM program was written for time-independent growth, breakage and aggregation. Besides, the calculations of the recursion coefficients of the orthogonal polynomials, the corresponding quadrature rule and the polynomials values at the fixed quadrature abscissas that are executed only once and stored in memory were not timed. Therefore, the comparison of CPU times must be taken with caution.

6.1. Pure growth problems The present form of the DuQMoGeM was developed for breakage and aggregation problems. In Section 5, it has been formally extended to include growth problems. However, this is not the best way to perform this extension because weighted-residual methods are best suited to solve non-hyperbolic problems. A better extension of DuQMoGeM to growth problem should use the method of characteristics. Nonetheless, it was important to verify the method ability to handle the growth term in the present form. In Section 5, the growth term in Eq. (16) with j = ϕj , given as x

qj = [g(x, t)f (x, t)ϕj (x)]0max − (g, ϕ˙ j )d(x)

(47)

2194

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

102

a

c

ξ1 ξ2 ξ3 ξ4 analytical numerical

Pure Growth - Case 3 DuQMoGeM (n=4, N=8)

1.0

|

100

(φ)

|μk -μk

10

-2

ξk(t)

(φ),a

|/|μk

(φ),a

0.8

k 10

0 1 2 3 4 5 6 7

-4

Pure Growth - Case 3 10-6

10

DuQMoGeM (n=4, N=8)

-8

0.0

0.2

0.4

0.6

0.6

0.4

0.2

0.0 0.8

0.0

1.0

0.2

0.4

d

1.2

numerical analytical

1.0

0.30

w k(t )

f(x,t)

t=1

0.6

0.2

Pure Growth Case 3

0.25 0.20 0.15 0.10

DuQMoGeM (n=4, N=8)

Pure Growth - Case 3

0.05

DuQMoGeM (n=4, N=8)

0.0 0.0

1.0

w1 w2 w3 w4 analytical numerical

0.40 0.35

0.8

0.4

0.8

t

t b

0.6

0.00 0.2

0.4

0.6

0.8

1.0

x

0.0

0.2

0.4

0.6

0.8

1.0

t

Fig. 4. Case 3 with pure growth. DuQMoGeM solution for n = 4, N = 8: (a) relative errors in the moments, (b) distribution function and (c) abscissas and (d) weights of the Gauss–Christoffel quadrature.

was approximated in Eq. (42) using the fixed-point Gaussian quadrature and the distribution function expansion by



2n−1

qj

Gji ci

(48)

i=0

where the growth matrix was given by Eq. (43). However, as a dualquadrature method, DuQMoGeM can also use the Gauss–Christoffel quadrature for the second integral, giving: qj

[g(x, t)ϕj (x)ω(x)ϕi (x)]x0max

−

n 

wk g( k , t)ϕ˙ j ( k )

(49)

k=1

The usage of the alternative form given by Eq. (49) has a larger computational cost because the Gauss–Christoffel quadrature must be evaluated within the ODE integrator. However, it has been verified that it usually gives a better approximation of (g, ˙ j )d(x) for problems in [0, ∞) domain (positive growth rates) than the corresponding approximation using the Gauss–Laguerre quadrature given by Eqs. (48) and (43). 6.1.1. Negative growth rate problems The most interesting examples of the application of the present extension of DuQMoGeM to growth problems are those involving negative growth rates. For these problems, a finite domain can be used, enabling a better functional approximation of the distribu-

tion function. Therefore, the test cases with negative growth rate were analyzed in the [0,1] range using shifted Legendre polynomials as the orthogonal basis and using Eq. (48) for the growth term approximation. The following problems were analyzed: 1. Case 1: g(x) = − 0.5x with f(0, t) = 0.3ı(x − 0.3) + 0.7ı(x − 0.7), 2. Case 2: g(x) = − 0.5 with f(0, t) = 60x2 (1 − x)3 , 3. Case 3: g(x) = − 0.5 with f(0, t) = 1. The QMOM was also used to solve the first case because the moment fluxes at the domain boundaries are zero for the linear negative growth rate but it cannot be used for the other two cases with constant negative growth rate. The first growth case illustrates the behavior of the QMOM and DuQMoGeM when both methods can solve for the first moments even though the DuQMoGeM functional expansion is lost because the distribution function does not belongs to the L2 [0, 1] space. Due to the given initial distribution, both methods were applied with n = 2. Since Eq. (43) gives the exact value of Gij for a Gauss–Legendre quadrature for N ≥ 4, the minimum value of N = 4 was used in DuQMoGeM solutions. Figs. 1 and 2 show the results for the QMOM and DuQMoGeM solutions up to t = 15. Figs. 1(a) and 2(a) show the evolution of the relative errors in the regular and generalized moments, respectively. Since both methods calculated the zero-order moments without any error, their errors were not represented in these figures. Figs. 1(b) and 2(b) show the evolution

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203 4

a

10

Pure Growth - Case 4 DuQMoGeM (n=3, N=6)

3

ξk(t )

10

102

ξ1 ξ2 ξ3 analytical numerical

101

100

ξ1 ξ2 analytical eq.(48) eq.(49)

1.2

1.0

ξk(t )

a

2195

0.8

Pure Growth - Case 5 DuQMoGeM (n=2, N=4)

0.6

0.4

0.2 0.0

2.0

4.0

6.0

8.0

0.0

10.0

0.2

0.4

0.6

0.8

1.0

0.8

1.0

t

t

b

b

0.50

0.8 0.7

ξ1 ξ2 ξ3 analytical numerical

0.30

0.6 0.5

w k(t )

w k(t )

0.40

0.20

Pure Growth - Case 5 DuQMoGeM (n=2, N=4)

0.4 0.3

Pure Growth - Case 4 DuQMoGeM (n=3, N=6)

0.10

ξ1 ξ2 analytical eq.(48) eq.(49)

0.2 0.1

0.00 0.0

2.0

4.0

6.0

8.0

10.0

t Fig. 5. DuQMoGeM solution for Case 4 with pure growth using n = 3: (a) abscissas and (b) weights of the Gauss–Christoffel quadrature.

of the abscissas and weights of the Gauss–Christoffel quadrature for the QMOM and DuQMoGeM solutions, respectively. From these figures, it can be seen that both methods obtained the same solution. The mean CPU times for QMOM and DuQMoGeM solutions were 38 and 29 ms, respectively. The reason why the DuQMoGeM integration was faster for this example can be explained by the fact that all generalized moments tend to steady values after about t = 10. This occurs because the distribution function tends to a Dirac delta function at x = 0 as t→ ∞ and the generalized moments tend to the finite values of the Legendre polynomials at x = 0, which have magnitudes between 0.05 and 0.5 for k = 1, 2,3. On the other hand, the QMOM integrates the regular moments, k that continue to decrease without bound for k > 0. Therefore, the DASSL routine had almost no effort to perform the integration after t = 10 for the DuQMoGeM solution but had the same effort for the QMOM solution for all t. However, this comes with a drawback because the moments were integrated using a given relative tolerance. For the DuQMoGeM, this implies in the loss of significant digits as the generalized moments approach the limiting values, while the values of the regular moments integrated by the QMOM solution keep the number of accurate digits close to −log 10 ( r ) as they decrease to zero as t→ ∞. In case 2 of pure growth, the initial distribution function is exactly represented by the polynomial basis using the first 6 terms. However, for t > 0 more terms are needed. The DuQMoGeM solution for n = 4 and N = 8 is shown in Fig. 3. Fig. 3(a) shows the relative errors of the generalized moments. Due to the error in evaluating the moment fluxes at x = 0, the errors in the general-

0.0 0.0

0.2

0.4

0.6

t Fig. 6. DuQMoGeM solution for Case 5 with pure growth: (a) abscissas and (b) weights of the Gauss–Christoffel quadrature.

ized moments increases sharply up to t = 0.2 and reach over 100% for some moments after t = 0.6. Nonetheless, the distribution function is very well represented up to the time when it is almost null, as can be seen in Fig. 3(b). This happens because the values of the generalized moments are quite small and tend to zero for t ≥ 2. Therefore, their relative errors are large but their absolute errors are small. However, the error accumulation in the generalized moments eventually reached a point where the abscissas and weights of the Gauss–Christoffel quadrature could not be accurately calculated, as shown in Fig. 3(c) and (d), respectively, for t values larger than about 1. The distribution function of case 3 is initially represented by only the first term in the polynomial expansion. However, for t > 0, the solution is the step function f(x, t) = H(1 − x − 0.5t), x ∈ [0, 1]. Although the approximation of the distribution function is not accurate for low values of n, a solution could be obtained by the DuQMoGeM. Similarly to case 2, the solution could be carried out until the error accumulated by the inexact approximation of the moment fluxes was sufficiently large to impede the correct determination of the Gauss–Christoffel quadrature. The DuQMoGeM solution for n = 4 and N = 8 is shown in Fig. 4. Fig. 4(a) shows the relative errors of the generalized moments, which increase very quickly. Fig. 4(b) shows the distribution function and its approximation at t = 1, which shows the expected Gibbs phenomenon at the point of discontinuity. Fig. 4(c) and (d) shows, respectively, the abscissas and weights of the Gauss–Christoffel quadrature, whose simulated values begin to deviate from the analytical solution around t = 0.6.

2196

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

a

a 2.0

1.3 1.2 1.1

0 1 2 3 analytical, all k

1.6

1.4

0.9 0.8 0.7 0.6 0.5

1.2

0.4

Pure breakage - Case 6 QMOM (n =2)

1.0 0.0

2.0

4.0

Pure breakage - Case 6 QMOM (n=2)

1.0

k

ξk(t) or w k(t )

μk (t )/μk (0)

1.8

ξ1 ξ2 w1 w2

0.3

6.0

8.0

0.2 0.0

10.0

2.0

4.0

6.0

t

b b

8.0

10.0

t 1.0

2.0

0.8

Pure breakage - Case 6 QMOM (n=4)

1.6

0 1 2 3

1.4

ξk(t) or w k(t )

μk (t )/μk (0)

1.8

k

1.2

1.0 0.0

ξ1 ξ2 ξ2 ξ2 w1 w1 w1 w2

0.6

0.4

0.2

Pure breakage - Case 6 QMOM (n =4) 2.0

4.0

6.0

8.0

10.0

t Fig. 7. Case 6 with pure breakage. QMOM solution for k with (a) n = 2 and (b) n = 4.

The previous results for cases 1–3 were also obtained using the approximation given by Eq. (49). The results were basically identical to those obtained with Eq. (49) because both the Gauss–Christoffel and the fixed-point Gauss–Legendre quadratures employed in each case can obtain qj without any error for the considered growth rates. It should be pointed out that the accuracy of the DuQMoGeM solutions for these problems are greatly enhanced if the method is applied using a moving boundary technique as in the FCMOM (Strumendo & Arastoopour, 2008). This has been tested as part of the ongoing research in the extensions of the DuQMoGeM. 6.1.2. Positive growth rate problems Two examples of pure growth with positive growth rates were chosen to illustrate the DuQMoGeM behavior: • Case 4: g(x) = 0.5x with f(0, t) = x3 exp (− x), and √ • Case 5: g(x) = 0.5 x with f(0, t) = 0.3ı(x − 0.3) + 0.7ı(x − 0.7). Both were solved in the [0, ∞) domain using Laguerre polynomials. Case 4 is an example where the initial distribution function has an exact representation by Eq. (23) with ω(x) = x3 exp (− x) for n ≥ 1 and the Gauss–Laguerre quadrature calculates exactly the Gji matrix. This case could be solved using Eq. (48) with high accuracy, even for values of t when the approximation of the distribution

0.0 0.0

2.0

4.0

6.0

8.0

10.0

12.0

t Fig. 8. Case 6 with pure breakage. QMOM solution for k and wk with (a) n = 2 and (b) n = 4.

function given by Eq. (23) is very poor. Fig. 5 shows some results for the DuQMoGeM solution with n = 3 (N = 6). Fig. 5(a) shows the abscissas and Fig. 5(b) shows the weights of the Gauss–Christoffel quadrature obtained from the generalized moments. It can be seen that the agreement with the analytical solution is perfect up to t = 10, when the zeroth order moment of the distribution has increased almost 150 times its initial value. The QMOM solution with n = 3 is also exact for this problem and it is not shown here. Mean CPU times for these QMOM and DuQMoGeM solutions were 69 and 51 ms. Case 5 illustrates the behavior of both methods when the growth terms given by Eq. (47) involve the integral of a nonpolynomial function. Besides, the distribution is given by two Dirac delta functions. In this case, there is no functional approximation of the distribution, which compromises the approximation of qj given by Eq. (48). Fig. 6(a) and (b) shows, respectively, the abscissas and weights of the Gauss–Christoffel quadrature obtained using DuQMoGeM with Eqs. (48) and (49). It is clear that the usage of Eq. (49) gives the correct solution. This occurs because the Gauss–Christoffel quadrature calculates exactly the integral in Eq. (49) for the considered distribution function. The solution given by QMOM is exactly the same as that obtained by DuQMoGeM with Eq. (49) as both methods become similar. The only difference between them is the set of moments that are solved for. Mean CPU times for the solution with QMOM and DuQMoGeM using Eq. (49) were 22 and 24 ms.

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

10-9

a

a -10

10

2197

Pure Breakage - Case 6

-10

10

10-11

-11

10-12

|μk-μk |

(φ),a

10-12

(φ)

10

-15

10

10

DuQMoGeM (n=1, N=3)

10

0.0

10-18 0.0

2.0

2.0

4.0

6.0

10

10-15 -16 -17

-16

10

-14

Pure Breakage - Case 6

(φ)

-14

10-13

k 0 1 2 3 DuQMoGeM (n=2, N=5)

-13

10

(φ),a

(φ),a

|μ0-μ0 |/| μ0 |

10

8.0

10.0

4.0

t

b

6.0

8.0

10.0

8.0

10.0

t b

2.0

1.0

0.8

ξk(t) or w k(t )

ξk(t) or w k(t )

Pure Breakage - Case 6

1.5

ξ w 11

1.0

0.5

ξ1 ξ2 w1 w2

0.6

0.4

Pure Breakage - Case 6

0.2

DuQMoGeM (n=1, N=3) 0.0 0.0

2.0

4.0

6.0

DuQMoGeM (n=2, N=5) 8.0

0.0 0.0

10.0

2.0

t

4.0

6.0

t ()

Fig. 9. Case 6 with pure breakage. DuQMoGeM solution for n = 1 and N = 3: (a) k n

relative error, (b) {wk , k }1 .

6.2. Pure breakage problems in a finite domain At first, the examples given by Dorao and Jakobsen (2006b, 2006c) were considered as possible test cases, but they were discarded due to two reasons. First, the daughter particle breakage distribution chosen by them does not represent any particle breakage process as it does not satisfy the conditions given by Eqs. (6)–(8). Besides, their examples with known analytical solutions, except one, were presumably wrongly reported because the given values for f(x, t), b(x), h(x | x ) = (x )P(x | x ) and S(x, t) do not satisfy the assumed PBE. Therefore, two new examples were similarly derived for the pure breakage problem in a finite domain. A third example was taken from literature. Dorao and Jakobsen (2006b, 2006c) chose problems for each f(x, t) is not a function of x, which intensifies the effect of quadrature errors in the solution of QMOM. This is also done here, assuming the same function, f(x, t) = 2 − e−t , x ∈ [0, 1]. However, a binary breakage, (x) = 2, ∀ x, with an uniform daughter particle breakage distribution, P(x | x ) = H(x − x)/x , was considered. Two cases were analyzed according to the choice of the breakage frequency: 1. Case 6: b(x) = x2 , which implies that all integrands are in some finite-dimensional polynomial space, and 2. Case 7: b(x) = x1/3 , which leads to integrands that do not belong to any finite-dimensional polynomial space. For these two choices of b(x), the corresponding source terms that were deducted in order to the PBE be satisfied are:

()

Fig. 10. Case 6 with pure breakage. DuQMoGeM solution for n = 2 and N = 5: (a) k n

relative error, (b) {wk , k }1 .

1. Case 6: S(x, t) = 2x2 (2 − e−t ) − 2(1 − e−t ) and 2. Case 7: S(x, t) = 7e−t − 12 + 7(2 − e−t )x1/3 , For the QMOM simulations, both Pj (x, t) and (xj , S)dx were calculated analytically to let the quadrature error be only the one related to the breakage operator. For the DuQMoGeM simulations, the source terms were also calculated by the N-point Gaussian ¯ measure. For Case 7, this has an additional quadrature for the d negative impact in the DuQMoGeM accuracy because the source term cannot be exactly integrated by this Gauss quadrature. The behavior of the first four normalized moments for the QMOM solution is shown in Fig. 7. The analytical solution gives k (t)/k (0) = 2 − e−t , ∀ k, which were plotted in Fig. 7(a) but not in Fig. 7(b) for a better visibility. Fig. 7(a) shows the results for n = 2, where it can be seen that all moments but 0 deviate from the analytical solution. The behavior of the moments for this case resembles that for the first example of Dorao and Jakobsen (2006c). Fig. 7(b) shows the results for n = 4, where it is clear that the 4 first normalized moments converged to the analytical solution in the graph scale. Considering all 2n moments that are solved for by QMOM, the maximum absolute error at t = 10 reduces from 0.087 for n = 2 (and 0 ) to 0.0002 for n = 4 (and 4 ). However, the QMOM solution fails as a discretization procedure n both for n = 2 and 4. Fig. 8 shows the transient behavior of {wk , k }1 for both cases. For a time-independent particle number distribution function, all k must be time invariant. From the analytical solution, it is easy to realize that wk (t)/wk (0) = 2 − e−t , ∀k. These analytical solutions were not plotted in Fig. 8 for a better visibility.

2198

a

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

100

a

-2

10-6

(φ),a

10-8 10-10

10-10

(φ)

Pure breakage - Case 7 QMOM (n =2)

-12

10

10-14 0.0

10-12

DuQMoGeM (n=1) 10

-14

Pure Breakage - Case 7 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

10-16

0.9

0.0

2.0

4.0

t

b

N 100 500 2000 8000

(φ),a

a

a

|μk -μk | / | μk |

k

|μ0 -μ0 |/| μ0 |

0 1 2 3

10-4

10-8

-2

10-4

10

10-6

10

b

7.0 6.0

6.0

8.0

k

N 100 500 2000 8000

10.0

t

Pure breakage - Case 7 QMOM (n=2)

10

-2

10

-4

Pure Breakage - Case 7

10-6 10

ξ1 ξ2 w1 w2

3.0

(φ),a (φ)

4.0

|μk-μk |

ξk(t) or w k(t)

5.0 -8

10-10 10-12

2.0

10-14

1.0

10-16

0.0 0.0

10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0 1 2 3 DuQMoGeM (n=2)

-18

0.0

0.9

2.0

4.0

6.0

8.0

10.0

t

t ()

Therefore, it is clear that the transient variations of all abscissa values shown in these figures are inadequate. Besides, it is easy to see that all weights fail to follow its analytical solution for n = 2 or 4. Thus, this example clearly shows that QMOM can obtain accurate moments but they are not precise enough to generate an accurate quadrature for the discretization of the particle number distribution function. Therefore, its use for PB-CFD simulations is not recommended. Let now consider the solution of this case using the DuQMoGeM for n = 1 and 2, leading to a linear and a cubic polynomial expansion of f in Eq. (23). Since x ∈ [0, 1], the shifted Legendre polynomials are used in this expansion. For n = 1, the largest polynomial that is being numerically integrated for Lji determination is of 5th order. Therefore, N is chosen to be 3 in this case. For n = 2, the quadrature needs to be exact for a 9th order polynomial, and then N = 5 was chosen. Since their accuracy is very good, only the generalized moment errors were plotted in Figs. 9(a) and 10(a) for the cases with n = 1 and () 2, respectively. For the solution with n = 1, 1 = 0 identically. Fig. ()

10(a) shows that the largest error is for 0 . The transient behavn ior of {wk , k }1 for the solutions with n = 1 and 2 are shown in Figs. 9(b) and 10(b), respectively. Their behavior clearly follows the analytical solution showing that the moments were determined with enough accuracy for a precise discretization of the particle number distribution function. It is clear from Figs. 9 and 10 that both solutions are accurate to the machine precision at t = 0 and that the accuracy decreases along time integration due to the given integration tolerance

Fig. 12. Case 7 with pure breakage. k and (b) n = 2 and several N.

error in DuQMoGeM solution for (a) n = 1

employed at each time step in the DASSL routine. It should be pointed out that, for these two solutions, the accuracy cannot be made better by increasing the N value. As a matter of fact, a large N value accumulates more truncation error. Thus, the accuracy tends to decrease slowly as N increases for large N values. The difference in accuracy of the discretization of the particle density function between QMOM and DuQMoGeM can be fully 101

Pure Breakage - Case 8

ξk(t) or w k(t)

Fig. 11. Case 7 with pure breakage. QMOM solution with n = 2: (a) k relative error n and (b) {wk , k }1 .

ξ1 w1 analytical DuQMoGeM, N=2n QMOM

0

10

n=1 10-1

0

2

4

6

8

10

t Fig. 13. Case 8 with pure breakage. Gauss–Christoffel quadrature obtained by QMOM and DuQMoGeM with n = 1.

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

10

a

4

a

102

2199

104

Pure Aggregation - Case 9 DuQMoGeM (n=2, N=4)

103

(φ),a |μ(φ) | / |μ(φ),a | k -μk k

100 102

ξk(t)

10-2 10-4

0 1 2 3 4 5

Pure Breakage Case 8

-6

10

-8

10

DuQMoGeM (n=3, N=6)

-10

10

101

k

0

2

4

6

8

10

10

Pure Breakage Case 8

0.6

0.5

1.0

1.5

2.0

2.5

3.0

100

-1

10

w k(t)

ξk(t)

0.8

0.0

t

b

ξ1 ξ2 ξ3 analytical numerical

1.0

-1

10

t b

ξ1 ξ2 analytical numerical

0

ξ1 ξ2 analytical numerical

-2

10

0.4 -3

10

DuQMoGeM (n=3, N=6)

0.2

Pure Aggregation - Case 9 DuQMoGeM (n=2, N=4) -4

0.0 0

2

4

6

8

10

t

c

0.0

0.5

1.0

1.5

2.0

2.5

3.0

t

3.0 Fig. 15. DuQMoGeM solution of Case 9 with pure aggregation for n = 2: (a) abscissas and (b) weights of the Gauss–Christoffel quadrature.

Pure Breakage - Case 8

2.5

DuQMoGeM (n=3, N=6)

2.0

w k(t)

10

ξ1 ξ2 ξ3 analytical numerical

1.5 1.0 0.5 0.0 0

2

4

6

8

10

t Fig. 14. DuQMoGeM solution of Case 8 with pure breakage for n = 3: (a) relative () errors of k , (b) abscissas and (c) weights of the Gauss–Christoffel quadrature.

appreciated by comparing Figs. 8(a) and 10(b). These results show how striking differently are the results for QMOM and DuQMoGeM and that accuracy and number of particle classes in f discretization are independent for DuQMoGeM. The computational costs of the QMOM (n = 2) and the DuQMoGeM (n = 2, N = 5) were equal within the standard deviations of the mean CPU times, which were 32 and 31 ms, respectively. Case 7 with pure breakage is much more difficult than Case 6 because the chosen breakage frequency leads to integrands that do not belong to any finite-dimensional space of polynomials. However, the analytical solution for f is the same and it still belongs to ℘0 . The QMOM solution for any n ≥ 2 could not be performed up to the desired time (t = 10) because the error associated with

the Gauss–Christoffel quadrature is large for the present case and it accumulates up to the point when the Gauss–Christoffel quadrature cannot be further calculated. This can be seen in Fig. 11(a) and (b) that shows the relative error of the moments and the abscissas and weights of the Gauss–Christoffel quadrature, respectively. For n = 1, a QMOM solution up to t = 10 could be calculated but it embodied a large relative error in 0 that reached 66%. It is interesting to keep n = 1 or n = 2 fixed and analyze the behavior of the DuQMoGeM solutions for increasing N. Fig. 12(a) and (b) shows the evolution of the errors of the generalized moments for n = 1 and n = 2, respectively. For all solutions shown in Fig. 12(a), () 1 = 0 identically. It is clear that the quadrature error is the main source of the solution error for both cases, showing that the use of an adaptive quadrature routine would be advantageous for this case. It can also be seen that the moment integration errors built up more slowly for the n = 2 solution in comparison to the n = 1 solution. For large N (say 500), both solutions give sufficiently accurate () k to enable the calculation of wk and k with high accuracy. Since these results are quite similar to those shown in Figs. 9(b) and 10(b), they are not shown here. The final example of pure breakage comes from Ziff and McGrady (1985). It is similar to Case 6 but the initial distribution function is a Dirac delta: • Case 8: b(x) = x2 , P(x | x ) = H(x − x)/x , (x) = 2, S(x, t) = 0, f(x, 0) = ı(x − 1)

2200

a

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

a

4

10

10-14

Pure Aggregation - Case 10 DuQMoGeM (n=2, N=4) 10-15

ξ1 ξ2 analytical numerical

10-16

k

0 1 2 3

(φ)

(φ),a

2

10

|μk-μk |

ξk(t)

103

1

10

10

-17

10

-18

Case 11, Π(∞) = 1 DuQMoGeM (n=2, N=4)

100 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

0.0

t

b

b

10

-4

10

-5

Pure Aggregation - Case 10 DuQMoGeM (n=2, N=4) 0.0

0.5

1.0

1.5

2.0

8.0

10.0

12.0

10

-2

10

-4

10

-6

10.0

12.0

a

-3

ξ1 ξ2 analytical numerical

|μk-μk |

w k(t)

10

6.0

100

10-1

-2

4.0

t

100

10

2.0

2.5

10-8 10

-10

10

-12

0 1 2 3

Case 11, Π(∞) = 1 QMOM (n=2)

k

10-14

3.0

3.5

4.0

10

-16

t

0.0

2.0

4.0

6.0

8.0

t

Fig. 16. DuQMoGeM solution of Case 10 with pure aggregation for n = 2: (a) abscissas and (b) weights of the Gauss–Christoffel quadrature.

()

Fig. 17. Case 11 for (∞) = 1. Absolute error in k and k for (a) DuQMoGeM and (b) QMOM solutions with n = 2.

6.3. Pure aggregation in semi-infinite domain whose analytical solution is f (x, t) = exp(−tx2 )[ı(x − 1) + 2tH(1 − x)]

(50)

This example is quite interesting because QMOM can only solve it for n = 1 whereas DuQMoGeM can solve it for any value of n. The Gauss–Christoffel quadrature obtained by both methods for n = 1 are shown in Fig. 13. It is clear that the quadrature used by DuQMoGeM gave more accurate results than that used by QMOM. The relative error of 0 at t = 10 is 18% for QMOM and 9.8% for DuQMoGeM. An increase in N in DuQMoGeM does not lead to a better accuracy for this case because the integration error is associated to the presence of a Dirac delta distribution in the solution. These errors in 0 show that the discretization using only one size class is not adequated for this problem. The DuQMoGeM results for n = 3 (N = 6) are shown in Fig. 14. The relative errors in the generalized moments are plotted in Fig. 14(a). It can be seen that the relative error in the zeroth order moment, which is identical to 0 , was always below 0.1%. Therefore, the representation of the distribution by three size classes seems to be quite reasonable. The accuracy of the Gauss–Christoffel quadrature is shown in Fig. 14(b) and (c), which depicted its abscissas and weights, respectively. Although there is no perfect match to the analytically calculated solution, the errors are quite small.

The cases considered here use the sum kernel, a(x, x ) = x + x , with two different initial distributions: • Case 9: f(x, 0) = exp (− x), whose solution is (Gelbard & Seinfeld, 1978): f (x, t) =

 exp(−t − 2x + xe−t ) I1 (2x 1 − e−t ) √ −t x 1−e

(51)

• Case 10: f(x, 0) = ı(x − 1), whose solution is (Ramkrishna, 2000): ∞ i−1  [x(1 − e−t )]

f (x, t) = exp(t − x + xe−t )

i!

ı(x − i)

(52)

i=1

Both problems were solved by DuQMoGeM using ω(x) = exp(− x). Therefore, the initial distribution of Case 9 is exactly represented by Eq. (23) for n ≥ 1. Fig. 15(a) and (b) shows the abscissas and weights of the Gauss–Christoffel quadrature up to t = 3 obtained from DuQMoGeM solution with n = 2 (N = 4). The agreement with the analytical solution is so good that the curves coincide in these figures. The same solution was obtained with QMOM. Mean CPU times for the QMOM and DuQMoGeM solutions were 259 and 443 ms, respectively. The usage of the four-point Gauss–Laguerre quadrature instead of the two-point Gauss–Christoffel quadrature in the double integral present in the aggregation terms made DuQMoGeM be 70% slower than QMOM for this case.

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

a

4.0

Case 11, Π(∞) = 2 DuQMoGeM (n=4, N=8)

3.5 3.0

t 0.0 0.5 1.0 2.0 analytical numerical

2.5

f(x,t)

Case 10 is another example for which a QMOM solution can only be obtained for n = 1 whereas DuQMoGeM does not have this limitation. For instance, Fig. 16(a) and (b) shows the abscissas and weights of the Gauss–Christoffel quadrature up to t = 4 obtained from DuQMoGeM solution with n = 2 (N = 4) together with the one computed from the analytical solution. The agreement between the numerical and analytical values is excellent. It can be seen in Fig. 16(b) that the w2 value increase from zero at t = 0 to a maximum value and then decreases as the aggregation proceeds. Fig. 16(a) shows that 1 (0) = 1, representing the initial distribution, and that 2 > 1 . Therefore, the second class appears to give a better representation of the large particles generated by aggregation.

2201

2.0 1.5 1.0 0.5

6.4. Simultaneous breakage and aggregation in semi-infinite domain

0.0 0.0

0.5

1.0

1.5

2.0

x McCoy and Madras (2003) gave an excellent example to test the DuQMoGeM behavior for a simultaneous breakage and aggregation problem in the [0, ∞) domain:

b

101

Case 11, Π(∞) = 2, n=4

The solution behavior depend on the value of√B, or more specifically, on the value of (∞) = 0 (∞)/0 (0) = 2B, which is the ratio of the particle number density at the steady-state condition to its value at the initial condition. The analytical solution was given by McCoy and Madras (2003) as: 2

f (x, t) = [(t)] exp[−x(t)]

ξk(t)

• Case 11: (x) = 2, ∀ x, P(x | x ) = H(x − x)/x , b(x) = Bx, a = 1 and f(x, 0) = exp (− x).

0

10

(53)

where

-1

10

0.0

(54)

If (∞) = 1, (t) = 1, ∀ t and the solution is time invariant, that is, f(x, t) = exp (− x). For (∞) < 1 aggregation dominates and for (∞) > 1 breakage dominates. For a Laguerre polynomial expansion, given by Eq. (23), ω(x) = exp (− x) = f(x, 0) and the initial condition is given by c0 = 1 and ci = 0, ∀ i > 0. Therefore, for (∞) = 1, the numerical solution must be time invariant if the quadrature used for the breakage and aggregation terms is accurate. This makes this problem an interesting example. For the given breakage and aggregation functions, the integrands are, at most, of order 4n − 1, which should be exactly integrated by the N-point Laguerre quadrature. Therefore, N ≥ 2n, and the minimum value of N = 2n was chosen. The errors in the moments in DuQMoGeM and QMOM solution for n = 2 are shown in Fig. 17(a) and (b), respectively. In these figures, when the error is exactly zero the corresponding curve cannot be seen. For the DuQMoGeM solution, all errors are smaller than the required absolute tolerance of integration in DASSL (10−14 ). Although not shown n here, the numerical values of {wk , k }1 were constant to machine precision along the simulation. On the other hand, Fig. 17(b) shows that the second and third moments in QMOM solution have large errors. 2 has a error that does not decrease with time. Accordingly, n the calculated values of {wk , k }1 by the QMOM solution were not invariant. Finally, a case of dominant breakage with (∞) = 2 was considered. It has been found that the DuQMoGeM solution with n = 4 is very good as can be seen in Fig. 18(a), which shows the representation of the distribution for several values of t. The agreement between the analytical solution and the functional expansion given by Eq. (23) is very good. Fig. 18(b) and (c) shows, respectively, the abscissas and weights of the Gauss–Christoffel quadrature up to t = 4 obtained from the DuQMoGeM, QMOM and analytical

1.0

2.0

3.0

4.0

5.0

t

c

1

10

Case 11, Π(∞) = 2, n=4

0

10

ξ1 ξ2 ξ3 ξ4 analytical DuQMoGeM (N=2n) QMOM

-1

w k(t)

1 + (∞) tanh((∞)t/2) (t) = (∞) (∞) + tanh((∞)t/2)

ξ1 ξ2 ξ3 ξ4 analytical DuQMoGeM (N=2n) QMOM

10

-2

10

10-3

-4

10

0.0

1.0

2.0

3.0

4.0

5.0

t Fig. 18. DuQMoGeM solution of Case 11 for (∞) = 2 and n = 4: (a) representation of f(x, t) for x ∈ [0, 2] at several values of t, (b) abscissas and (c) weights of the Gauss–Christoffel quadrature.

solutions. It can be seen that QMOM gives the Gauss–Christoffel quadrature with a reasonable error while the DuQMoGeM results are quite close to the analytical values. Mean CPU times for these QMOM and DuQMoGeM solutions were 79 and 47 ms, respectively. In this case, DuQMoGeM was 66% faster than QMOM. 7. Conclusions This work considered the solution of the monovariate population balance equation with growth, breakage and aggregation. The

2202

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203

unique internal variable is an additive particle property (mass or volume) that is conserved during breakage and aggregation. First, the QMOM was expressed as a particular case of a more general quadrature-closed method based on generalized moments of the particle number distribution function (QMoGeM). Then, a weighted-residual method based on the generalized moments (WRMoGeM) was derived. Further, it was shown that, if WRMoGeM is closed by the Gauss–Christoffel quadrature based on the distribution function, the resulting method (QWRMoGeM) is identical to the QMoGeM. This way, the correct representation of QMOM as a weighted-residual method was then achieved, correcting the derivation present in the literature (Dorao & Jakobsen, 2006b, 2006c). The WRMoGeM formulation was used to derive a new moment-based method, that employs two quadrature rules. The Gauss–Christoffel quadrature based on the particle number distribution function is mainly used for discretizing the particulate system in several “phases”, although sometimes it should be preferred for integrating the growth term in the moment equations. Another quadrature rule is responsible to the accurate integration of the breakage and aggregation terms in the equations for the generalized moments. This method was called DuQMoGeM. In this work, the DuQMoGeM was implemented in a Galerkin formulation. Thus, f was represented by a product of a weight function and a series of polynomials that are orthogonal in the inner product defined by this weight function and the same polynomials were used to define the generalized moments. Therefore, the second quadrature used by DuQMoGeM was the easily available variable-order Gaussian quadrature based on that weight function. It should be pointed out that the f expansion is a bonus in this moment-based method. If the distribution is a continuous function, the f expansion will uniformly converge to it. If the distribution does not belong to the L2 space, the f expansion is lost but the method still solves for the generalized moments in a scaled form. Results for pure growth, pure breakage, pure aggregation and breakage-aggregation problems which have analytical solutions were quite promising, as the DuQMoGeM was able to obtain their solutions with good accuracy. However, the treatment of problems with growth should be improved. The DuQMoGeM clearly separates the problems of: • characterization of the particle population in several “phases”: it uses the generalized moments based on a convenient polynomial basis which allows the computation of the Gauss–Christoffel quadrature based on the particle number distribution function using a better conditioned procedure, • approximation of the particle number distribution function: it may use a basis in L2 , as the polynomial basis used in the present work, or even finite-element or spectral discrete approximations, • integration of the growth, breakage and aggregation terms: it may use a different Gaussian quadrature rule or the Gauss–Christoffel quadrature, as in the present work, or a set of rules related to the particle number distribution function approximation or even an automatic quadrature routine. All of these benefits came from the WRMoGeM formulation, which seems to be a good framework to develop new numerical methods for the accurate solution of the PBE that are also able to represent a particle system by an optimal and small set of particulate “phases”. Acknowledgments The author acknowledges the financial support from CNPq (grants no. 301672/2008-3, 485817/2007-1 and 476268/2009-5).

References Alopaeus, V., Laakkonen, M., & Aittamaa, J. (2006a). Numerical solution of momenttransformed population balance equation with fixed quadrature newblock. Chemical Engineering Science, 61, 4919–4929. Alopaeus, V., Laakkonen, M., & Aittamaa, J. (2006b). Solution of population balances with breakage and agglomeration by high-order moment-conserving method of classes. Chemical Engineering Science, 61, 6732–6752. Attarakih, M. M., Drumm, C., & Bart, H.-J. (2009). Solution of the population balance equation using the sectional quadrature method of moments (SQMOM). Chemical Engineering Science, 64, 742–752. Bove, S., Solberg, T., & Hjertager, B. H. (2005). A novel algorithm for solving population balance equations: The parallel parent and daughter classes. Derivation, analysis and testing, Chemical Engineering Science, 60, 1449. Dorao, C. A., & Jakobsen, H. A. (2006a). A least squares method for the solution of population balance problems. Computers and Chemical Engineering, 30, 535–547. Dorao, C. A., & Jakobsen, H. A. (2006b). Numerical calculation of the moments of the population balance equation. Journal of Computational and Applied Mathematics, 196, 619–633. Dorao, C. A., & Jakobsen, H. A. (2006c). The quadrature method of moments and its relationship with the weighted residuals. Chemical Engineering Science, 61, 7795–7804. Dorao, C. A., & Jakobsen, H. A. (2007). Time-space-property least squares spectral method for population balance problems. Chemical Engineering Science, 62, 1323–1333. Dorao, C. A., & Jakobsen, H. A. (2008). hp-adaptive least squares spectral element method for population balance equations. Applied Numerical Mathematics, 58, 563–576. Fan, R., Marchisio, D. L., & Fox, R. O. (2004). Application of the direct quadrature method of moments to polydisperse, gas–solid fluidized beds. Powder Technology, 139, 7–20. Gautschi, W. (1994). Algorithm 726: ORTHPOL—A package of routines for generating orthogonal polynomials and Gauss-type quadrature rules. ACM Transactions on Mathematical Software, 20, 21–62. Gautschi, W. (2004). Orthogonal polynomials: Computation and approximation. In Numerical mathematics and scientific computation. Oxford, Great Britain: Oxford University Press. Gelbard, F., & Seinfeld, J. H. (1978). Numerical solution of the dynamic equation for particulate systems. Journal of Computational Physics, 28, 357–375. Gordon, R. (1968). Error bounds in equilibrium statistical mechanics. Journal of Mathematical Physics, 9, 655. Grosch, R., Briesen, H., Marquardt, W., & Wulkow, M. (2007). Generalization and numerical investigation of QMOM. AIChE Journal, 53, 207–227. Hildebrand, F. B. (1987). Introduction to numerical analysis (revised 2nd ed.). Dove Publications Inc. Kumar, S., & Ramkrishna, D. (1996). On the solution of population balance equations by discretization—I. A fixed pivot technique. Chemical Engineering Science, 51, 1311–1332. Marchisio, D. L., & Fox, R. O. (2005). Solution of population balance equations using the direct quadrature method of moments. Journal of Aerosol Science, 36, 43–73. Marchisio, D. L., Vigil, R. D., & Fox, R. O. (2003). Implementation of the quadrature method of moments in CFD codes for aggregation-breakage problems. Chemical Engineering Science, 58, 3337–3351. Massot, M., Laurent, F., Kah, D., & de Chaisemartin, S. (2010). A robust moment method for evaluation of the disappearance rate of evaporating sprays. SIAM Journal of Applied Mathematics, 70, 3203–3234. McCoy, B., & Madras, G. (2003). Analytical solution for a population balance equation with aggregation and fragmentation. Chemical Engineering Science, 58, 3049–3051. McGraw, R. (1997). Description of aerosol dynamics by the quadrature method of moments. Aerosol Science and Technology, 27, 255. Petitti, M., Nasuti, A., Marchisio, D. L., Vanni, M., Baldi, G., Mancini, N., et al. (2010). Bubble size distribution modeling in stirred gas–liquid reactors with QMOM augmented by a new correction algorithm. AIChE Journal, 56, 36–53. Petzold, L. (1982). A description of DASSL: A differential/algebraic system solver (Technical report). Sandia National Laboratories. Ramkrishna, D. (2000). Population balances—Theory and applications to particulate systems in engineering. San Diego: Academic Press. Silva, L. F. L. R., Damian, R. B., & Lage, P. L. C. (2008). Implementation and analysis of numerical solution of the population balance equation in CFD packages. Computers and Chemical Engineering, 32, 2933–2945. Silva, L. F. L. R., & Lage, P. L. C. (2007). Implementation of an Eulerian multi-phase model in OpenFOAM and its application to polydisperse two-phase flows. In Proceedings of the 1st OpenFOAM international conference, volume CDROM. London, England: ICON, 16 pp. Strumendo, M., & Arastoopour, H. (2008). Solution of PBE by MOM in finite size domains. Chemical Engineering Science, 63, 2624–2640. Su, J., Gu, Z., Li, Y., Feng, S., & Xu, X. Y. (2008). An adaptive direct quadrature method of moment for population balance equations. AIChE Journal, 54, 2872–2887.

P.L.C. Lage / Computers and Chemical Engineering 35 (2011) 2186–2203 Su, J., Gu, Z., & Xu, X. Y. (2009). Advances in numerical methods for the solution of population balance equations for disperse phase systems. Science in China Series B: Chemistry, 52, 1063–1079. Yoon, C., & McGraw, R. (2004). Representation of generally mixed multivariate aerosols by the quadrature method of moments: I. Statistical foundation. Aerosol Science, 35, 561–576.

2203

Zhu, Z. J., Dorao, C. A., Lucas, D., & Jakobsen, H. A. (2009). On the coupled solution of a combined population balance model using the least-squares spectral element method. Industrial & Engineering Chemistry Research, 48, 7994–8006. Ziff, R. M., & McGrady, E. D. (1985). The kinetics of cluster fragmentation and depolymerization. Journal of Physics A: Mathematical and General, 18, 3027–3037.