Solution of population balances with growth and nucleation by high order moment-conserving method of classes

Chemical Engineering Science 62 (2007) 2277 – 2289 www.elsevier.com/locate/ces

Solution of population balances with growth and nucleation by high order moment-conserving method of classes Ville Alopaeus ∗ , Marko Laakkonen, Juhani Aittamaa Laboratory for Chemical Engineering and Plant Design, Helsinki University of Technology, P.O.B. 6100, FIN-02015 HUT, Espoo, Finland Received 6 July 2006; received in revised form 8 November 2006; accepted 22 December 2006 Available online 16 January 2007

Abstract The high order method of classes, developed in our earlier work [Alopaeus, V., Laakkonen, M., Aittamaa J., 2006a. Solution of population balances with breakage and agglomeration by high order moment-conserving method of classes. Chemical Engineering Science 61, 6732–6752] for solution of population balances (PBs), is extended to problems with growth and primary nucleation. The growth problem leads to a hyperbolic partial differential equation with fundamentally different numerical characteristics than the PB with breakage and agglomeration only. However, we show that the principle of moment conservation in the numerical solution can also be applied to this advection-type problem, leading to extremely accurate numerical solutions. The method is tested for two numerical cases. The first one is mass transfer induced particle growth, and the second one is primary nucleation with constant growth (similar to the Riemann advection problem). For mass transfer induced growth, we first analyze functional form of the growth rate from mass transfer correlation viewpoint, and derive a general analytical solution for the power-law growth. The numerical results from the moment conserving method are also compared to one well established high resolution numerical method for advection problems, namely the Lax–Wendroff method with van Leer flux limiter. It was shown that the present method is far superior by predicting the distribution moments with several order of magnitudes lower numerical error. For the Riemann problem with constant growth rate, the present method predicts the shock front location exactly without any numerical diffusion. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Population balance; Method of classes; High resolution methods; Mass transfer; Size distributions; Advection; Riemann problem

1. Introduction Population balance (PB) modeling is an important tool to analyze any particulate phase systems. The main idea behind PB modeling is to formulate balances for distributed variables along one or several internal coordinates. “Internal” refers here to those coordinates that are orthogonal to any spatial coordinates of the system, i.e., distributions are modeled at each point in space. During the last few decades, the PB framework has found applications in the fields of crystallization, agglomeration, gas–liquid and liquid–liquid systems, just to mention a few (Randolph and Larson, 1988; Ramkrishna, 2000). In order to close the general PB equation, separate ∗ Corresponding author. Tel.: +358 9451 2641.

E-mail addresses: Ville.Alopaeus@tkk.fi, Ville.Alopaeus@hut.fi (V. Alopaeus). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2006.12.070

models are required. These include, e.g. breakage and agglomeration rates of particles, daughter size distribution after breakage, growth or shrinkage rate due to mass transfer, and rate of primary nucleation. In this work, we focus on the numerical solution of PBs by the method of classes, where the continuous distribution is discretized into separate categories, each of which is represented by a pivot element. This pivot is a single variable on the internal coordinate for each category representing the particles belonging to that category. The particle number density is then represented as Dirac delta functions with respect to the internal coordinate instead of a continuous distribution. In principle, a numerical solution for pure growth problem is a relatively straightforward task, and can be found by using the method of lines. The category method proposed by Kumar and Ramkrishna (1997) was actually based on this idea, essentially setting the pivot elements to follow the particle growth.

2278

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

In their method, however, primary nucleation brings additional difficulties and relatively complicated algorithm is required to add new categories during the simulation. With their method, also determination of the breakage and agglomeration tables before simulation cannot be done since these tables depend on the location of pivot points. This is a highly undesirable feature if these models are to be used in process simulators or with CFD, since changing the number of categories requires the number of simulation variables as well as the breakage and agglomeration tables to be updated during simulation. This updating is the most time consuming task in the numerical solution, and may increase the computational time by an order of magnitude. As pointed out, e.g. by Gunawan et al. (2004), the growth problem is mathematically similar to the advection problem that has been examined quite thoroughly in the field of fluid dynamics. As these methods use fixed grids, they allow construction of breakage and agglomeration tables in advance. These methods have also been under considerable scrutiny in the past due to their wide-spread use. Therefore, we will compare the methods developed in this work mainly to one of these well established numerical methods for advection (the Lax–Wendroff scheme with van Leer flux limiter, as proposed by Gunawan et al., 2004). Recently, Qamar et al. (2006) compared several high order methods in problems related to crystal growth. They found that instead of fully upwind differencing (as in Lax–Wendroff scheme), a weighted average between upwind and central differencing gives slightly better results. All the weighted average methods are still only second order accurate in smooth regions and first order accurate near discontinuities. In the first part of this work (Alopaeus et al., 2006a), we developed a numerical method for solution of PB equations for breakage and agglomeration. It was based on conservation of arbitrary set of moments in the discretized particle size domain. The method is applicable to arbitrary internal coordinate discretization, and a further degree of freedom is introduced from the chosen set of conserved moments. If only zeroth and third moment (number and volume) were conserved in breakage or agglomeration, a low order method resulted in. This is the state of the art for all method of classes. If more than two moments are conserved, we end up with a high order method. In the first part of this paper, we showed that these high order moment conserving methods are several orders of magnitude more accurate than the state of the art low order methods with comparable discretization and computational cost. In this part, we extend the high order moment conserving method for growth and primary nucleation. Only one internal coordinate is considered, but the method can be extended to several internal coordinates in a quite straightforward manner. The simplest extension to multidimensional problems may be the method of dimensional splitting (see LeVeque, 2002 for further discussion). Validation of the method through numerical examples with multiple internal coordinates would, however, require quite extensive testing. Therefore, extension to multidimensional problems is left for the future.

2. The PB In the first part of this work (Alopaeus et al., 2006a), the following discretized PB was considered NC

 dYi = Y˙i,in + (Li , Lj )Li g(Lj )Yj dt j =1

+

NC NC  

((Li , Lj , Lk )F (Lj , Lk )

j =1 k=1

+ A (Li , Lj , Lk )Li FB (Lj , Lk ))Yj Yk −Y˙i,out −g(Li )Yi − Yi

NC 

(F (Li , Lj ) + FB (Li , Lj ))Yj

j =1

+

NC  (Li , Lj ) j =1

Lj

Gj Yj + (Li )Li .

(1)

Note that the above equation is shown in a very general form, and contributions from any breakage (g), agglomeration (F ), growth (G) or primary nucleation () event may affect the size distribution in any category. In that respect, there are yet no assumptions regarding the above equation, except that length variable was taken to be the internal coordinate. The method presented here can, however, be extended to other internal coordinates without any further complication. An important feature of the method of classes adopted here, that should be kept in mind throughout the present development, is that the discretized size distribution is not density distribution with respect to the internal coordinate, but all particles at each category are assumed to be located exactly at the pivot element. For particle size density distributions this implies that the continuous variable has a dimension 1/m4 , where three dimensions refer to external (spatial) dimensions, and one to the internal coordinate (particle size). For the discretized distribution, the variable has a dimension of 1/m3 , i.e., it describes number density of particles with respect to external coordinates only. It is very important to understand that the discretized PB formulation is not the physical model describing evolution of the particle size distribution. Rather, it is a numerical method to solve the underlying physical model, which is an integropartial–differential equation. Since the method of classes is perhaps the most easily understood numerical method (besides Monte-Carlo simulations) for solution of PB models, the numerical solution procedure of counting each particle and adding it into appropriate “bin” is easily confused with the physical model. In this work, we first revise the moment equations for continuous and discrete distributions:  ∞ k = n(L)Lk dL (continuous distribution), (2) 0

k =

NC  i=1

Yi zik

(discrete distribution).

(3)

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

These will be used in the development of the method and also later in verification against analytical results. In the following chapters, the equivalence of moments and particle number concentrations in various categories is used in a truncated form k =

 (N)

Yi zik .

(4)

i=(1)

This allows a set of known moments for any particular event affecting particle size distribution to be distributed to a number of categories. These contributions can then be summed to cover the whole discretized distribution. In this way moments from all events (breakage, agglomeration, growth, primary nucleation, etc.) can be conserved. This is a fundamentally different viewpoint from methods previously presented for numerical solution of PB equations, since we do not try to find an explicit expression for time evolution of any particular category (except in the general sense of Eq. (1)), but to find contributions from each category to the neighboring categories. This leads to simpler formulation of the numerical problem, and allows more flexibility for the problem specification. Mathematically speaking, we use the property of superposition for partial differential equations with overlapping high resolution discrete elements. 3. Moment conserving numerical methods for growth In this chapter we develop a moment conserving numerical method for growth. In order to simplify the development, breakage and agglomeration are neglected here and only contribution from growth to the size distribution is considered. The method is, however, directly applicable to situations where other source terms such as breakage and agglomeration are simultaneously present by just summing contributions (time derivatives of particle number concentrations) from various sources. Primary nucleation is also neglected for now, but considered later in this paper in a separate chapter. For pure growth, the PB becomes d(n(L)) j(G(L)n(L)) =− . dt jL

In order to develop a moment-conserving method for growth, we first carry out a moment transformation for pure growth in category j dG(k)j dt

Note that this is a hyperbolic partial differential equation. It behaves fundamentally differently than the PB equation with pure breakage and agglomeration source terms. The most distinctive feature of the pure growth problem is that sharp fronts are conserved during the solution. Actually the growth problem is exactly similar to the advection equation. Numerical solutions of this problem are well known to be either low order (i.e., low accuracy) at least near discontinuities, or the solutions tend to oscillate. This is widely studied, e.g. in the fields of fluid dynamics and computational physics. Therefore, we expect that the best practices for numerical methods for growth and for other advection problems are mutually exchangeable. However, more detailed application of the methods developed in this work to other problems than solution of the PB equation with growth is outside the scope of the present work.

 =

L+ L−

(k)L(k)−1 G(L)n(L) dL,

(6)

where L− and L+ refer to the limits of j th category. Subscripts should be interpreted as G for growth, (k) is the order of the moment, and j refers to the category for which the moment transformation for growth is calculated. Since the particles in each category are assumed to be located exactly at the pivot elements as Dirac delta functions, the moments can be simply calculated as dG(k)j dt

(k)−1

= (k)Lj

Gj Yj .

(7)

In order to speed up the numerical solution, it is preferable to split the moment evolution equation into two parts, one of which depends only on the discretization and another that can generally change during the time integration. The growth table is preferably constructed so that it contains only grid dependent part. Then construction of this table can be done before the time integration in case of fixed pivot elements, and the table needs to be updated only if the grid is adjusted. It is also preferable to calculate the moments for each grid points with a shifted origin. This reduces condition numbers (all of them) for the transformation matrix, implying better accuracy for the matrix inversion especially for those pivot elements that are far from the true origin (zero particle size). This shifted origin should be close to the pivot element sizes into which the moments are distributed. The exact location is not very important, and very often it can be chosen to be half grid size smaller or larger than the smallest pivot element into which the moments are distributed. The time evolution of the moments due to growth can then be written as dG(k)j dt

(5)

2279

= mG(k)j Gj Yj ,

(8)

where mG(k)j = (k)(Lj − L0j )(k)−1 .

(9)

The rest of the moment evolution equation, namely Gj Yj , is calculated during the time integration. If Gj is only a function of length, it can be included in the mGkj part and one multiplication is spared during the time integration. However, growth rate in general is also a function of other system properties besides size, such as degree of supersaturation in crystallization. Therefore, we retain the more general form here. Next, we use the equivalence between evolution of moments and its contribution to the source terms at each category (Eq. (4)), and we can express the product table for growth as (mGj ) = [A]

(L , Lj ) , Lj

(10)

2280

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

where [A] is a linear operator transforming the discretized number distribution space into moment space. The growth table vector refers to the column with elements from vector , whereas row j is fixed, i.e., we distribute the moments from category j into categories specified by the  vector. The elements of the [A] matrix are simply Aab = (z(b) − L0j )(a) .

(11)

And the growth table can be found from a simple matrix inversion as (L , Lj ) = [A]−1 (mGj ). Lj

(12)

Note that the denominator Lj is retained only to make dimensions consistent and table  dimensionless. In numerical calculations, there is no need to carry out the division since the above form appears as such in the discretized PB equation. There remains two further factors that characterize the method. One is the span over which the moments are distributed. It appears that distributing the elements in consecutive categories with an upwind scheme appears to give the best results, i.e.,   M j (1) = j − floor +I (13) 2 and the last element is j (M) = j (1) + M − 1,

(14)

where I = 1 for positive velocity (G > 0) and I = −1 for negative velocity (G < 0). The function floor () rounds towards next smaller integer. For a general scheme, it is advisable to calculate separate tables for positive and negative velocities before time integration. For boundary elements where elements of  would be smaller than 1 or larger than NC, we distribute the moments to the first or last M categories, respectively. In this way the information of the chosen conserved moments are forced to remain in the discretized number density space. If j(G(L)n(L))/jL term is large at the discretization size limits (L = 0 or L = LNC + LNC /2), this may result in instabilities originating from the boundaries, especially the upper limit (LN C + LNC /2). In these cases it is advisable to reduce the number of conserved moments at the boundaries, leading to partially absorbing boundary conditions. The methodology for this is proposed in the first part of this work (Alopaeus et al., 2006a) for agglomeration. Another factor characterizing the method is the choice of conserved moments. For pure growth, zeroth moment (number of particles) is a natural conserved property, and should be included in the  vector. Other moments can be chosen quite freely, but a sequence of first integer moments is quite good according to our experience. For mass transfer processes, third moment should also be conserved, since that gives the total mass in the particulate phase, and this is the natural conserved property in mass transfer. Conservation of only zeroth and first moment leads to the first order upwind method, which is well

known to suffer from quite large numerical diffusion in the numerical solution. Note also that even and odd number of conserved moments gives quite different span of product particles, and the upwind characteristic may be quite different for these two. Odd number gives distribution that is symmetric about the growing particle category, whereas even numbers give upwind scheme. With the above limits for the particle span, even number of moments seem to be a better choice than odd. Since the number of particles is the natural conserved property in growth, the growth matrix is further scaled for exact conservation of particle number (each column in the matrix must sum up to one). This is important only when the numerical solution is practically indistinguishable from the analytic solution, and we still want to compare relative error between the numerical method and the analytical solution. In practical calculation, the difference between the scaled and unscaled growth tables is negligible. The above derived method is in principle quite close to the spectral methods and the methods of weighted residuals. Spectral methods and methods of weighted residuals are known to be extremely accurate for certain cases, but they also possess some deficiencies which have discouraged their widespread use in numerical fluid flow analysis. In the present method, there is no need to explicitly calculate trial function coefficients or to apply Fourier transforms inherent to the spectral methods. The high order elements in the present method (spaces of discrete number concentrations mapped by the  vectors) are completely overlapping in the global discretized internal coordinate space. This makes the present method very flexible and allows the method to be applicable to arbitrary grids despite the high order resolution. 4. Primary nucleation Primary nucleation is spontaneous formation of small particles from supersaturated fluid. The rate of primary nucleation can be distributed to the relevant categories (usually a few smallest categories) by using the same moment conserving method than for other properties. In this work we assume that the size of the primary nucleates is infinitely small, so that the time derivative of their moments is simply dBk = B0 dt dBk =0 dt

for k = 0, for k > 0.

(15)

Other primary nucleate sizes can be used without any further complication by simply calculating the exact moments for the assumed size. Also a primary nucleate size distribution can be used if a model for that is available, by calculating the distribution moments for nucleates belonging to each category and then distributing these similarly to the daughter size distribution in the first part of this work (Alopaeus et al., 2006a). These moments can then be distributed as earlier:   −1 dB ((Li )Li ) = [A] . (16) dt

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

Note again that multiplication with Li is retained merely for keeping the dimensions consistent, and the product L appears always as such. For small primary nucleates, the nucleation source is best distributed to M first categories where M is the number of conserved moments. Other source terms that can be expressed with the same functional form besides the primary nucleation can be distributed exactly the same way than primary nucleates. The source term is generally best distributed evenly on both sides of the actual source if possible, or in M first or last categories. 5. Diffusion

d(n(L)) j(j(Dn(L))/jL) = . dt jL

(17)

The diffusion contribution can be derived by applying the moment transformation, integrating by parts twice, and remembering that the distribution is expressed as Dirac delta functions with peaks located only at the pivot points. Then we end up with the following expression for time rate of change for the moments due to diffusion dt

(k)−2

= (k)((k) − 1)Lj

D j Yj .

(18)

The diffusion contribution table can be calculated quite similarly than the convection contributions. We need only one table for diffusion due to symmetry (there is no upwind direction), and the elements are preferably distributed evenly on both sides of the category from where the particles are diffusing. Then the first element of the distribution vector should be calculated from j (1) = j − ceil(M/2) + 1. The function ceil(M/2) rounds to the nearest integer greater than or equal to M/2. As for the convection, the diffusion table is preferably divided into a part depending on the pivot element locations only, and another part that depends also on temporally dependent properties: dD (k)j dt

= mD (k)j Dj Yj .

(19)

D (L , Lj ) = [A]−1 (mDj ). Lj

(22)

This matrix is appended to the discretized PB equation formally identically to the convection term, but diffusion coefficient Dj is used instead of the velocity Gj . Note that the dimension of the diffusion contribution table is 1/m2 , whereas the convection contribution is 1/m.

mD (k)j = (k)((k) − 1)(Lj − L0j )(k)−2 .

6.1. High resolution methods for advection equation Gunawan et al. (2004) suggested using high resolution finite volume methods adapted from the field of compressible gas dynamics, also well established in other various fields related to fluid flow analysis, for growth term in the PB equations. A detailed description and analysis of various methods can be found from LeVeque (2002). Of those methods Gunawan et al. analyzed, the following was found to give best results (Lax–Wendroff method with Van Leer flux limiter, referred as HR1 in their paper and in the following comparisons)  Yim+1 − Yim dYi  = dt growth k Gi Yi − Gi−1 Yi−1 = −   h  Gi kGi − 1− (Yi+1 − Yi ) i 2h h    Gi−1 kGi−1 (Yi −Yi−1 ) i−1 , (23) − 1− h 2h

(20)

Equivalence of the moment time derivatives and the moments distributed to the neighboring categories is similar than for the convection part: D (L , Lj ) . Lj

In order to verify the numerical behavior of the developed moment-conserving method, we test it with two cases. First, however, we revise shortly one well established numerical method recently proposed for solution of PBs with growth. Then we test the moment conserving method for mass transfer induced particle growth, and finally to a primary nucleation with constant growth rate (equivalent to the Riemann advection problem). In all cases we use uniform grid and the relevant particle size interval is [0, 1]. In the moment conserving method, six first moments are always conserved (moments from zero to five).

i =

Discretization of the only dependent part is then

(mDj ) = [A]

The elements of the [A] matrix are also exactly the same than in the convection case. The diffusion table is

6. Numerical examples

Diffusion was not considered separately in the PB equation (1), but can be accounted for in a straightforward manner. Although diffusion equation is parabolic in its nature, whereas convection is hyperbolic, these two terms can be treated quite similarly in the present numerical solution. Diffusive spreading of the number density distribution can be accounted for by using a moment transformation for the diffusion equation

dD (k)j

2281

(21)

i =

m Yim − Yi−1 m − Ym Yi+1 i

| i | + i . 1 + | i |

,

(24)

(25)

The above equations are formulated for growth (G > 0). k is the time step, h is the grid size, and is the flux limiter. This method is used as a reference numerical solution when numerical solution developed in this work is analyzed. Although

2282

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

there is also an analytical solution for all the numerical examples considered here, a reference to a well established numerical method gives further faith in the methods developed in this work. A problem arises since time step size appears in the righthand side (time derivative) of the above equation. Then standard high order variable step size integrators cannot be used as such. On the other hand, low order time integration schemes such as explicit Euler are highly inaccurate, and a very small time steps are required for accurate solution. In our comparisons, we used a modified second order Runge–Kutta scheme where explicit Euler is first applied with half and full step sizes   k k m+1/2 Yk/2 (26) = Y m + Y˙ t m , , Y m , 2 2   k˙ m k m+1/2 m+1 m+1/2 , (27) + Y t + k/2, , Y Yk/2 = Y 2 2 Ykm+1 = Y m + k Y˙ (t m , k, Y m )

(28)

and finally the values after the step are calculated by applying a first order Richardson extrapolation with half and full step sizes as m+1 − Ykm+1 . Y m+1 = 2Yk/2

6.2. Size dependent growth The first numerical example is related to size dependent growth without primary nucleation. In Appendix A, size dependent growth rate is analyzed based on mass transfer correlations, and corresponding analytical solution is derived for general power-law growth

Initial distribution is taken to be

L30 3N0 2 , L exp − n(L0 ) = v0 0 v0

n(L) =

(29)

If the time derivative Y˙ does not depend on the step size, the method reduces to the classical second order Runge–Kutta formula.

G = bLp .

the initial distribution may appear quite chaotic with negative elements present. However, since a high number of moments are conserved, the discretized distribution contains more information about its continuous counterpart than if only values at pivot points were used. Then the exact distribution could in principle be reconstructed with high accuracy from the discrete counterpart. In case only values at pivot points are used in the discretization, information is lost already at the initial condition specification stage. Then the continuous distribution cannot be exactly reconstructed even theoretically from its discretized counterpart. The problem of reconstructing continuous distribution from the information available in the solution of its discrete counterpart is not considered here. Grid adaptation and distribution reconstruction with high order moment conserving scheme is left for the future. The discussion above concerns also the final analytical size distribution. In Appendix A, we derive the following analytical expression for the size distribution with a general power-law growth

(A.8)

(A.9)

with N0 = 1 m−3 and v0 = 0.001. In this example, L can be considered as dimensionless diameter, where the physical size is divided by a reference size. A special care needs to be taken when initial continuous distribution is discretized. If only distribution function values at the pivot points are employed, especially higher moments suffer from sparse discretization. Therefore, at least equal number of moments should be conserved in the initial discretization as in the actual numerical scheme. With this approach, the initial discretized distribution is exact in terms of the chosen moments, but there is no guarantee that the number concentrations at each category remain positive. In fact, when a coarse discretization is used with a high number of conserved moments,

3N0 −p 1−p L (L − (1 − p)bt)(2+p)/(1−p) v0

(L1−p − (1 − p)bt)3/(1−p) , × exp − v0

L ((1 − p)bt)1/(1−p) , n(L) = 0,

L < ((1−p)bt)1/(1−p) .

(A.13)

In Fig. 1, initial and final analytical distributions are shown with 40 evenly distributed categories and six conserved moments. Final distribution is calculated with p = −1 (diffusion controlled growth). Note that in the figure, the discrete variables are connected with straight lines between the number densities at pivot points, but the discrete distribution is actually defined only at these discrete points. It can be seen that there are quite large negative number concentrations even in the discretized analytical solution. The size of the negative elements is reduced if the grid would be refined, but due to a steep gradient near the shock locations, overshoot remains a problem in any discretized form of the distribution. Next, we test the method for p = −1, −0.5, 0 and 0.5. These values cover the region that was shown to be relevant for mass transfer induced growth. Note that for positive values of p the solution is spreading as the initial particle size distribution grows, since larger particles grow faster than small ones. For negative values of p the solution behaves as a shock wave, shock forming at particle size of L = ((1 − p)bt)1/(1−p) . The reason for this behavior is that small particles grow faster than large ones, reaching them in size. In the examples, we focus on the size interval [0, 1] and analyze grid refinement characteristics with various values of the growth exponent p. Physical simulation time was 100 s for each case. Growth rate parameter b was calculated as b=

0.81−p − 0.31−p . (1 − p) · 100

(30)

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

2283

Fig. 1. Initial and final analytical discretized distributions with 40 categories and six conserved moments. Final distribution corresponds to diffusion limited growth (p = −1).

10

p p p p

1

1000

0 0.5 -0.5 -1

p p p p

100

0.01

ERR%

ERR%

0.1

= = = =

0.001 0.0001

= = = =

0 0.5 -0.5 -1

10

1

0.00001 0.000001

0

20

40 60 Number of categories

80

100

0.1

0

20

40 60 Number of categories

80

100

Fig. 2. Relative errors with moment conserving method of classes (left) and the HR1 method (right). Six moments were set to be conserved in the present moment conserving method.

With this growth rate parameter initial distribution at L = 0.3 will grow to L = 0.8 after 100 s of integration. Since initial distribution consisted of negligible amount of particles larger than 0.3, it remains in the region [0, 1] during the whole simulation period unless excessive numerical diffusion will spread the largest particles out from the region. In Fig. 2, numerical results are shown for the two methods as functions of number of categories. Six moments were set to be conserved in the present moment conserving method of classes. The relative error refers to the maximum of absolute errors in eight first moments, i.e.,    calc,i − anal,i   · 100. ERR% = max  (31)  i=0,...,7 anal,i In order to emphasize the differences between the two methods, the ratio of the errors with the two tested methods is shown in Fig. 3. The errors for both methods are relatively high in magnitude due to the high number of moments compared in the

error functional, highest errors in the tested cases being slightly less than 3% for the moment conserving method and about 400% for the HR1 method for six categories. For a smaller number of compared moments, relative errors are lower for both the tested methods. In fact, the moment conserving method developed in this work gives exact results for the same number of moments that are set to be conserved in case of linear growth rate (p = 0). This is not true for the HR1 method. The moment conserving method is at least about two orders of magnitude more accurate than HR1 in all the tested cases, the improvement being smallest when p > 0. The reason is that in these cases, both the initial distribution and the growth rate are increasing for small particles. The ratio of volumes in the numerical grid is relatively high in these regions, since we used constant diameter interval discretization. This result in high numerical errors for the largest moments considered (although still at least two orders of magnitude smaller

2284

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

1

p p p p

MCM ERR% / HR ERR%

0.1

= = = =

0 0.5 -0.5 -1

0.01 0.001 0.0001 0.00001 0.000001 0.0000001

0

20

40 60 80 Number of categories

100

Fig. 3. Error in the moment conserving method of classes divided by the error in the HR1 method.

than the HR1 method). A partial improvement would be to use non-uniform grid, but we have not considered that option here. Generally two important observations can be made: the errors with the moment conserving method are already at least two orders of magnitude more accurate with a small number of grid points (less than 15), and the accuracy is improved far more rapidly than with the HR1 method as the number of grid points is increased. Actually, with six conserved moments, the present method is fifth order accurate. If the inherently better accuracy of the present method compared to the HR1 is used for reducing the number of grid points for equal tolerance, two advantages are imminent. The first is that simulation time is reduced approximately directly proportional to the number of grid points. The second advantage is that with a sparser grid, also time integration step size can be increased while still fulfilling the CFL stability criterion. This speeds up the overall simulation time further. Oscillations inherent to the present high resolution method are most severe when p = −1. In that case the distribution gradient is steepest at the location of the shock, implying some overshoot. This was visible even in the analytical solution (Fig. 1). In Fig. 4, numerical solution after 100 s of growth with six conserved moments is shown for 40 categories and diffusion controlled growth rate (p = −1). It can be seen that the solution oscillates slightly more than the analytical solution (Fig. 1), but the oscillation is confined to locations where the analytical solution gradients are also steep, and no spurious oscillations are visible in regions where solution is smooth. In the next example, oscillatory behavior is analyzed with a more illustrative example (the Riemann problem). 6.3. Primary nucleation with constant growth (the Riemann problem) In the second example, we test the method for primary nucleation with constant growth rate G. Initially, there are no particles present, and at time t = 0, primary nucleation starts at

Fig. 4. Final distribution after 100 s of integration. M =6, NC=40 and p=−1.

rate B0 = 1/m s. Nucleate size is assumed to be infinitely small (L0 = 0). The problem is formally similar to the well-known Riemann problem for advection, which is one of the most tested examples for numerical methods of this kind. The problem has a trivial analytical solution n(L) =

B0 G

n(L) = 0

for L Gt, for L > Gt.

(32)

From this, the exact moments can also be easily calculated k =

B0 (Gt)k+1 . G(k + 1)

(33)

It proves that the moment conserving method gives exact results for all those moments that are set to be conserved even for a small number of discretization points. Most notably, the location of the front, 21 /0 , and the variance, 2 /0 −(1 /0 )2 are exactly predicted. The smallest possible number of discretization points for exact conservation of the desired moments is obviously equal to the number of conserved moments, since we cannot have a bijective mapping into a functional space with a smaller dimension, i.e., M number of moments cannot be exactly preserved if mapped to NC < M grid points. If the location of the front and the distribution variance are the only variables of interest, we can solve the problem exactly with only three discretization points. However, in the numerical simulations we used at least six conserved moments and six grid points. Although the number densities associated with only six grid points give the shape of the front only vaguely, as can be seen in Fig. 5, it is basically possible to reconstruct the solution

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

Fig. 5. Numerical solution of the Riemann problem with six categories.

with high accuracy from the available information. For six conserved moments, there are six independent variables that can be used in the reconstruction process. This reconstruction process usually requires some assumed form of the distribution. For example, in the constant growth rate case, we could easily reconstruct the distribution by using just two first moments if we assume a distribution that is constant for the interval [0, L] and zero elsewhere. This assumed distribution is exact in this case since the present method predicts the specified moments exactly for constant growth rate. For equal number of discretization points and conserved moments, the present method is actually very close to the one we proposed earlier for moment transformed PB equations (Alopaeus et al., 2006b), except that the variables to be solved in the time integration are the number densities, not moments. In Fig. 6, the numerical solution is shown with 100 categories and six conserved moments. Some observations can be made. If we focus on the six first categories, we find a region of substantial oscillation. There are two reasons for that. Firstly, since the primary nucleate size does not coincide exactly to any pivot size, the nucleate needs to be distributed to six first categories, resulting in oscillations via the primary nucleation source term. Secondly, since we do not reduce the accuracy (number of conserved moments) for the smallest particles, the growth source term from the first categories needs to be distributed one-sidedly to the larger particle sizes. This enhances the oscillations inherent to the high order methods. Actually, the oscillations at the smallest categories could be smoothed out without considerable loss in accuracy by simply using a constant number density for the first categories, but we want to show the unsmoothed results here in order to give better illustration of the method performance.

2285

Fig. 6. Solution of the Riemann problem with 100 categories and six conserved moments.

The other region where oscillations are visible is near the discontinuity (near the shock front). This overshooting, closely related to the Gibbs overshooting phenomenon, was quite expected for the present high order method. The overshooting is present only near the discontinuity, and these oscillations do not spread from the region near the discontinuity, as is the case with unlimited Lax–Wendroff scheme. The HR1 method removes oscillations, but slightly overpredicts the front location and distribution variance (essentially smoothing the solution to some extent). Especially with a small number of categories, the HR1 method smears out the solution quite notably. HR1 method performance with 6, 20, and 100 categories are shown in Fig. 7. As the PB methods are usually applied as a part of a complex overall process model, even with multiple external coordinates such as in the CFD calculations, the number of categories should be kept at minimum. In these cases, 20–40 categories are the practical upper limit with present computer capacity, if the solution with respect to external coordinates are truly three dimensional. 6.4. Diffusion Although diffusion is not usually included in the PB equation, we present a numerical comparison of the diffusion equation solution just to show the potential of the present method for diverse problems. There are some examples where a diffusive term has been included in the PB model (Randolph and Larson, 1988), but this term rises usually from improper spatial averaging of the system leading to variable growth rate in the region that was assumed to be homogeneous.

2286

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

10 1 0.1

1

10

100

1000

ERR %

0.01 0.001 0.0001 0.00001 0.000001 0.0000001 0.00000001

number of categories

Fig. 8. Maximum of the relative errors in eight first integer moments in the diffusion problem.

actly of size L = 1. Note that with this initial distribution also negative elements are generally present. These negative elements are rapidly diffusing away as the true distribution spreads, and the distribution starts to look like the actual Gaussian distribution. Size distribution is discretized with even discretization on the interval L = [0, 2]. Diffusion coefficient is D = 0.001 m2 /s and simulation time is 10 s. During this simulation time, only negligible amount of particles reach the discretized size interval limits, so that the problem of specifying appropriate boundary conditions is avoided. With the present formulation, diffusion out from the selected size interval is neglected (particles are forced to remain in the chosen region). Six moments are conserved in the numerical scheme, although due to the symmetric nature of diffusion, one row for each column of D (L , Lj )/Lj does not receive particles (five remaining active). This results in a fourth order accurate solution, which can also be seen from the slope of Fig. 8. The analytical solution for the diffusion case for t > 0 is available as (Crank, 1975)   1 L2 n(L) = √ . (35) exp − 4DG t 2 DG t

Fig. 7. Solution of the Riemann problem with HR1 method. (a) Six categories; (b) 20 categories and (c) 100 categories.

In this example, we consider a pure diffusion equation with a constant diffusion coefficient jn(L) j2 n(L) = DG , jt jL2

(34)

where one unit of particles is located originally at L=1 as Dirac delta pulse. The moments of this initial distribution (k = 1 for all k) are spread into six categories closest to L = 1 similarly than primary nucleates in case there was no category ex-

From this, the moments can be easily calculated. The same measure for relative error was used than in the previous cases, i.e., maximum of absolute values of relative errors in eight first integer moments. Results are shown in Fig. 8. The error in the three first integer moments specifying the average location and the variance of the distribution are already below 10−8 % even with the smallest number of categories (six). This shows that the distribution is diffused with the correct diffusion coefficient without almost any error from the numerical method. To summarize the findings in the numerical examples above, we believe that the best practice is to calculate the solution with a truly high resolution method, such as the moment conserving method developed in this work, and then smooth the final solution for visualization purposes if necessary. If some overall distribution properties are of interest, such as the mean size, total surface area or location of the front, the smoothing procedure should not be used for best possible accuracy.

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

Similar numerical behavior was found in the present examples as those studied in the first part of the work (Alopaeus et al., 2006a), namely that the system of equations seems to be less stiff if a somewhat higher number of categories than the number of moments was used. The solution accuracy cannot be improved freely by increasing both number of categories and number of conserved moments, due to numerical instabilities. This is also a serious limitation in all QMOM methods. On the contrary, accurate solutions are obtained by using a reasonably high number of moments (six in the examples in this paper), and a moderate number of categories. Usually 15–30 gives already extremely accurate results without any numerical instabilities. 7. Conclusion A high order method of classes for growth and primary nucleation was developed. The method is based on conservation of arbitrary set of moments for the growth rate and for the primary nucleate size. This method is an extension to our earlier work (Alopaeus et al., 2006a) for solution of PBs for breakage and agglomeration, but due to the fundamentally different nature of the growth problem (being purely hyperbolic partial differential equation), the numerical solution requires a completely new derivation and testing with appropriate numerical examples. The first such test example is mass transfer induced particle growth. Firstly, we analyzed the system from the mass transfer rate viewpoint and derived an analytical solution for the general power-law growth rate. Numerical solutions with the method developed in this work and also with the Lax–Wendroff method with van Leer flux limiter (referred as the HR1 method) were compared to the analytical solutions with relevant growth rates. The numerical examples showed that the present method gives results that are several orders of magnitude more accurate than the results calculated with the HR1 method. The second example was primary nucleation with constant growth rate. This leads to a problem that behaves fundamentally similarly than the well-known Riemann advection problem. It is one of the standard numerical examples used for evaluation of various numerical methods relevant in the field of fluid dynamics. It was shown that the present method solves the Riemann problem exactly in terms of the distribution moments even with a small number of categories. For large number of categories, the Gibbs oscillations inherent to high order methods are present, but they are limited to regions close to the shock fronts or other steep gradients. Since the distribution moments are solved with high accuracy, we expect that also a non-oscillatory solution can be reconstructed from the available information very accurately.

[A] b1 , . . . , b3

FB (L, ) g G h I k k L Lmax m mGkj M n n(L) N N0 NC p q t v v0 Vm Y x zi

parameters in mass transfer coefficient correlations linear operator between growth distribution space and growth moment space parameters in mass transfer coefficient correlations

growth rate parameter, m1−p /s primary nucleation rate, 1/m4 s total concentration, mol/m3 diffusion coefficient, m2 /s agglomeration rate between particles of size and L that does not result in subsequent immediate breakage, m3 /s agglomeration rate between particles of size and L that results in subsequent immediate breakage, m3 /s breakage rate, 1/s growth rate, m/s grid size, m index for flow direction, dimensionless time step size, s mass transfer coefficient, m/s internal coordinate, length in this paper, m maximum particle size; size of the largest size category, m parameter in mass transfer coefficient correlations size dependent part of the growth table, mk−1 number of moments to be conserved in the discretization parameter in mass transfer coefficient correlations density of particle size distribution, 1/m4 mass transfer flux, mol/m2 s initial total particle number density, m−3 total number of categories growth rate exponent, dimensionless size discretization parameter, dimensionless time, s terminal velocity, m/s initial total particle volume, m3 partial molar volume of the transferring component, m3 /mol particle number concentration, 1/m3 driving force for mass transfer, dimensionless pivot element, m

Greek letters (Li , Lj ) A (Li , Lj , Lk )

Notation a1 , . . . , a3

b B0 ct D F (L, )

2287

 k

daughter particle size distribution; probability that a particle of size Li is born when Lj breaks, 1/m agglomeration product distribution; probability that a particle of size Li is born when particles Lj and Lk agglomerate and subsequently break down, 1/m ratio of two consecutive gradients particle size, m viscosity, kg/m s kth moment of density distribution, mk−3

2288

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289



vector that contains information of the categories into which agglomerate product is distributed growth distribution: a table that determines the distribution to category of size Li due to growth of particle of size Lj density, kg/m3 flux limiter agglomeration product distribution; a table that determines the fraction of agglomerate from particles Lj and Lk to category of size Li distribution of primary nucleates vector of conserved moment orders

(Li , Lj )

(Li , Lj , Lk )

(Li ) 

Appendix A. Analysis of mass transfer induced growth In this analysis, we first consider mass transfer correlations of type Sh = a1 + b1 · Rem Scn .

(A.1)

This form is quite general and widely used in analysis of mass transfer in dispersed phase systems. In order to simplify the present analysis, we limit to single component transfer, and neglect convective mass transfer flux (drift flux). The analysis could be extended to true multicomponent systems in a straightforward manner if necessary. When the above coefficient correlation is written explicitly, we find that k=

m−n v m Lm−1

a1 · D + b1 m−n n−1 , L  D

(A.2)

where v is the terminal or slip velocity between particulate and continuous phase. The physical properties are not functions of size, whereas the terminal velocity is. If we assume that velocity can be expressed with a power law v ∝ Lp1

(A.3)

and lump all size-independent parameters together, we end up with a2 k= (A.4) + b 2 Lp . L Mass transfer induced growth can be expressed as G = NV m ,

(A.5)

where N is the mass transfer flux (in mol/m2 s) of the transferring component and Vm is its partial molar volume. In absence of convective flux across the interface between dispersed and continuous phases, the mass transfer flux is N = ct kx.

(A.6)

Combining these equations, the growth rate can be expressed as a3 G= + b 3 Lp . L

(A.7)

In order to derive an analytical expression for growth, we consider the following simplified form G = bLp .

(A.8)

This is also the limiting value for two relevant cases, namely growth by pure diffusion ((a3 /L)  b3 Lp ) and growth due to relative velocity induced mass transfer ((a3 /L)  b3 Lp ). The lumped parameter p can be expressed by the original ones as p =mp1 +m−1. In most mass transfer coefficient correlations, m takes values close to 0.5. The parameter p1 takes values from 2 (Stoke’s law for small particles) to 0.5 for relatively large particles where the drag coefficient is constant. For even larger particles, the drag coefficient experiences sudden bumps and it becomes increasing function of the Reynolds numbers, but this region is not relevant to this work (Clift et al., 1978). Therefore, the relevant limits for the parameter p are ranging from −0.25 to 0.5 for relative velocity induced mass transfer and −1 for purely diffusion controlled mass transfer, so that the total interval of interest ranges approximately from −1 to 0.5. Note that the special case of size independent mass transfer (also called as McCabe’s L law) is a special case of the present development with p = 0. Next, we derive an analytical size distribution for the general power-law growth. The initial particle number density distribution is assumed to be

L30 3N0 2 n(L0 ) = , (A.9) L exp − v0 0 v0 where L0 refers to the size variable at time t = 0. The derivation of the analytical size distribution can be done by first noting that the growth rate is G=

dL = bLp . dt

(A.10)

Integrating, we obtain  L  t −p d = b d, L0

(A.11)

0

L0 = (L1−p − (1 − p)bt)1/(1−p) L0 = L exp(−bt)

p  = 1,

p = 1.

(A.12)

Inserting this into the initial size distribution and noting the change of variable in the integration (dL0 → dL), we end up with n(L) =

n(L)=

3N0 −p 1−p L (L − (1 − p)bt)(2+p)/(1−p) v0

(L1−p − (1 − p)bt)3/(1−p) , × exp − v0

  3N0 2 L3 exp(−3bt) , L exp −3bt− v0 v0

p  = 1,

p=1. (A.13)

If p < 0, the number density is zero for particles smaller than L < ((1 − p)bt)1/(1−p) .

(A.14)

V. Alopaeus et al. / Chemical Engineering Science 62 (2007) 2277 – 2289

In these cases, the distribution follows a shock-wave kind of behavior with a shock front forming at the above limiting particle size for large times. References Alopaeus, V., Laakkonen, M., Aittamaa, J., 2006a. Solution of population balances with breakage and agglomeration by high order momentconserving method of classes. Chemical Engineering Science 61, 6732–6752. Alopaeus, V., Laakkonen, M., Aittamaa, J., 2006b. Numerical solution of moment-transformed population balance equation with fixed quadrature points. Chemical Engineering Science 61, 4919–4929. Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops, and Particles. Academic Press, New York. Crank, J., 1975. Mathematics of Diffusion. second ed. Clarendon Press, Oxford.

2289

Gunawan, R., Fusman, I., Braatz, R.D., 2004. High resolution algorithms for multidimensional population balance equations. A.I.Ch.E. Journal 50, 2738–2749. Kumar, S., Ramkrishna, D., 1997. On the solution of population balance equations by discretization—III. Nucleation, growth and aggregation of particles. Chemical Engineering Science 52, 4659–4679. LeVeque, R.J., 2002. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, Cambridge. Qamar, S., Elsner, M.P., Angelov, I.A., Warnecke, G., Seidel-Morgenstern, A., 2006. A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering 30, 1119–1131. Ramkrishna, D., 2000. Population Balances, Theory and Applications to Particulate Systems in Engineering. Academic Press, San Diego. Randolph, A.D., Larson, M.A., 1988. Theory of Particulate Processes. second ed. Academic Press, New York.