- Email: [email protected]
A discretized model for tracer population balance equation: Improved accuracy and convergence
Computers and Chemical Engineering 30 (2006) 1278–1292
A discretized model for tracer population balance equation: Improved accuracy and convergence J. Kumar a,c,∗ , M. Peglow b , G. Warnecke a , S. Heinrich c , L. M¨orl c a
Institute for Analysis and Numerics, Otto-von-Guericke-University Magdeburg, Universit¨atsplatz 2, D-39106 Magdeburg, Germany Institute for Process Engineering, Otto-von-Guericke-University Magdeburg, Universit¨atsplatz 2, D-39106 Magdeburg, Germany c Institute of Process Equipment and Environmental Technology, Otto-von-Guericke-University Magdeburg, Universit¨ atsplatz 2, D-39106 Magdeburg, Germany
b
Received 3 June 2005; received in revised form 16 February 2006; accepted 27 February 2006 Available online 12 May 2006
Abstract A new discretization for tracer population balance equation is developed. It is compared to the modified discretized tracer population balance equation of Peglow et al. [Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2006). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineers, 52, 1326–1332]. The new formulation provides excellent prediction of the tracer mass distribution in all test cases. Furthermore, the new formulation is more efficient from a computational point of view, it takes less computational effort and is able to give a very good prediction for a coarser grid. Additionally, it is independent of the type of grid chosen for computation. For finer grids, both formulations tend to produce the same results. The performance of the new formulation is illustrated by the comparison with various analytically tractable problems. Moreover, the new formulation preserves all the advantages of the modified discretized tracer population balance equation of Peglow et al. [Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2005). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineering (published online: 29 December 2005)] and provides a significant improvement in predicting tracer mass distribution and tracer-weighted mean particle volume during aggregation process. © 2006 Elsevier Ltd. All rights reserved. Keywords: Population balance; Discretization; Aggregation; Particle; Tracer; Batch
1. Introduction Tracer studies have become an important tool in the experimental and computational study of disperse phase population dynamics. Ilievski and Hounslow (1995) estimated the agglomeration kinetics parameters from tracer data. Pearson, Hounslow, and Instone (2001) performed some tracer experiments and observed the breakage rate dependence on the size and the age of the granules. Further, Hounslow, Pearson, and Instone (2001) developed a two internal coordinate population balance equation (PBE). The twodimensional PBE has been reduced to two one-dimensional PBEs: one for granule size distribution (GSD) and the other for tracer mass distribution (TMD). They also presented a discretized PBE for tracer. In this work, we develop a new framework for the discretization of tracer population balance equation (TPBE). The discretized tracer population balance equation (DTPBE) presented in Hounslow et al. (2001) which was later modified by Peglow et al. (2006), predicts good results only for a size independent kernel. For sum and product kernels it leads to significant overprediction over exact Abbreviations: DTPBE, discretized tracer population balance equation; GSD, granule size distribution; PBE, population balance equation; TMD, tracer mass distribution; TPBE, tracer population balance equation ∗ Corresponding author. Tel.: +49 391 67 12329; fax: +49 391 67 12129. E-mail address: [email protected] (J. Kumar). 0098-1354/$ – see front matter © 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2006.02.021
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1279
Nomenclature a¯ B, D c Er H I Iagg K n N M q Sp t u, v v0 x
volume average of newborn particles birth and death terms volume of tracer in a granule relative error Heaviside step function total number of intervals the index of aggregation correction factor number density number of particles per unit volume mass density adjustable parameter determining the geometric discretization ratio a parameter defined by Eq. (14) time volume coordinate initial average volume representative volume an interval
Greek symbols β agglomeration rate for reduced one-dimensional problem βˆ agglomeration rate for complete two-dimensional problem δ Dirac-delta distribution η a function defined by Eq. (29) λ a function defined by Eq. (22) µ moments of the PSD τ dimensionless time Subscripts agg agglomeration i index for interval p primary particle T tracer Superscripts ana analytical in feed to a CST mod modified
solutions. Also, the presented model has limitations for the grid type and can be implemented only for special geometric grids. Now the objective is to get a discretized model which is more accurate and can be applied on general grids. 2. Theory A two-dimensional particle size distribution is defined as f (t, v, c), where v and c are two distinct properties, granule volume (size) and tracer volume, respectively. We assume that the density is constant and mass is replaced by volume. Thus, the total number in a domain D is given by D
f (t, v, c) dv dc.
(1)
It should be noted that granule volume v contains volume of tracer and volume of particles, that is, c ≤ v. Fig. 1 shows this two-dimensional particle space. The two-dimensional population balance equation can be obtained by extending the classical one-
1280
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 1. Two-dimensional space.
dimensional PBE (Hulburt & Katz, 1964) to two-dimensional space as ∂f (t, v, c) 1 v min(c,) ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d = β(t, ∂t 2 0 max(0,c−v+) ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d. − β(t, 0
(2)
0
The integral limits corresponding to are trivial. To fix the limits of γ we consider the following apparent relations c − γ ≤ v − , γ ≤ and obviously 0 < γ ≤ c. Clearly these relations provide the integral limits in the formulation (2). The limits of integration can also easily be perceived from Fig. 2. The above 2D PBE must be supplemented with an initial condition: f (0, v, c) = f0 (v, c).
(3)
The first and the second terms on the right-hand side of Eq. (2) account for the formation and the loss of the particles with the properties v and c. Let us assume that the aggregation kernel βˆ depends only on time and size of the granules, but not on the tracer contents within the granules, i.e. βˆ = β(t, v, ). Then the 2D PBE can be converted into two 1D PBEs corresponding to the conventional number density n(t, v) and mass of tracer within granules M(t, v). The number density n(t, v), may be obtained from f by integrating over all possible tracer masses: v n(t, v) = f (t, v, c) dc. (4) 0
The mass of tracer within the granules of size v is given analogously by v M(t, v) = cf (t, v, c) dc. 0
Fig. 2. Limits of integrals in 2D PBE.
(5)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1281
2.1. Reduction to 1D PBE for granule size distribution Integrating Eq. (2) over all possible values of c, we have ∂n(t, v) 1 v v min(c,) ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d dc = β(t, ∂t 2 0 0 max(0,c−v+) v ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d dc. − β(t, 0
0
Reversing the order of integration, we obtain ∂n(t, v) 1 v v−+γ ˆ = β(t, v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dc dγ d ∂t 2 0 0 γ ∞ v ˆ v, , c, γ)f (t, v, c)f (t, , γ) dc dγ d. − β(t, 0
(6)
0
0
(7)
0
ˆ v, , c, γ) = β(t, v, ), then by If we assume that the aggregation kernel depends only on time and the size of both granules, i.e. β(t, substituting c − γ = p and using the definition of density function n, we get ∞ ∂n(t, v) 1 v β(t, v − , )n(t, v − )n(t, ) d − n(t, v) β(t, v, )n(t, ) d. (8) = ∂t 2 0 0 2.2. Reduction to 1D PBE for tracer mass distribution Proceeding as before, multiplication of Eq. (2) by c and integrating over all possible values of c yields: v min(c,) ∂M(t, v) 1 v ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d dc β(t, c = ∂t 2 0 0 max(0,c−v+) v ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d dc. c − β(t, 0
0
(9)
0
Again, reversing the order of integration and some simple substitution gives v ∞ ∂M(t, v) = β(t, v − , )M(t, v − )n(t, ) d − M(t, v) β(t, v, )n(t, ) d. ∂t 0 0
(10)
Eqs. (8) and (10) are ordinary integro-differential equations which have to be solved numerically. Hounslow, Ryall, and Marshall (1988) developed a discretized method for solving GSD and Hounslow et al. (2001) extended the same idea for the TMD. The discretized PBE was applicable only for special geometric grid. In the next section we first briefly discuss the discretization of TMD and then propose a new idea for the better accuracy and implementation. The preceding population balance equation for the tracer mass can also directly be obtained from the classical one-dimensional population balance equation (8). Note that the tracer mass distribution M(t, v), which represent the total tracer mass contained in particles size v, can be defined as the product of the number density of particles size v with the mass of tracer c(t, v) in each particle: M(t, v) = c(t, v)n(t, v).
(11)
It follows that the rate of change of total tracer mass contained in particles of size v as a results of coagulation is therefore given by ∞ ∂M(t, v) 1 v = β(t, v − , )[c(t, v − ) + c(t, )]n(t, v − )n(t, ) d − n(t, v)c(t, v) β(t, v, )n(t, ) d. (12) ∂t 2 0 0 Using the definition of M(t, v) we obtain the required population balance equation (10). 3. Previous work Hounslow et al. (2001) proposed a discretized PBE for tracer which was based on discretized PBE for GSD (Hounslow et al., 1988). The discretized model could be applied only for fixed geometric grid of the type vi+1 = 2vi . Moreover, the proposed discretization for TMD was not consistent with granule mass. In a recent paper, Peglow et al. (2006) overcame this inconsistency in discretized model by introducing some correction factors and extended it for the adjustable grids of the type vi+1 = 21/q vi . Peglow et al. (2006) used the approach of Litster, Smit, and Hounslow (1995) and Wynn Ed (1996) for the extension to adjustable grids.
1282
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Noting Ni =
vi+1 vi
n(t, v) dv and Mi =
vi+1 vi
M(t, v) dv, the final extended and consistent model is given by
i−S 1 2(j−i+1)/q dMi = βi−1,j (Mi−1 Nj + Ni−1 Mj )K1 + βi−q,i−q Ni−q Mi−q dt 21/q − 1 j=1
+
q−1 i+1−S p+1 21/q − 2(j−i)/q − 2−p/q βi−p,j (Mi−p Nj + Ni−p Mj )K2 21/q − 1 p=1 j=i+1−Sp
+
q
i−Sp
p=2 j=i−Sp −1
+
i−S 1 +1 j=1
2(j−i+1)/q − 1 + 2−(p−1)/q βi−p,j (Mi−p Nj + Ni−p Mj )K3 21/q − 1
2(j−i)/q 1 − 1/q 2 −1
βi,j Ni Mj K4 −
i−S 1 +1 j=1
2(j−i)/q βi,j Mi Nj K5 − 21/q − 1
I
βi,j Mi Nj ,
(13)
j=i−S1 +2
where
q ln(1 − 2−p/q ) Sp = Int 1 − ln 2
(14)
and K1 =
2(i−j)/q , 1 + 2(i−j−1)/q
K2 = K3 =
(15)
2(i−j)/q , 1 + 2(i−j−p)/q
(16)
K4 =
1 , 1 + 2(j−i)/q
(17)
K5 =
21/q − 1 −2i/q + 2(2i−j)/q (21/q − 1) − . 2(j−i)/q 2i/q + 2j/q
(18)
The function Int[x] gives the integer part of x. Let us call this formulation as modified DTPBE. The corresponding discretization of the GSD equation (8), proposed by Litster et al. (1995) and corrected by Wynn Ed (1996), has been summarized as follows: q i−S 1 2(j−i+1)/q dNi = βi−1,j Ni−1 Nj + dt 21/q − 1 j=1
+
i−Sp
p=2 j=i−Sp −1
2(j−i+1)/q − 1 + 2−(p−1)/q 1 2 βi−p,j Ni−p Nj + βi−q,i−q Ni−q 1/q 2 −1 2
q−1 i+1−S i−S 1 +1 p+1 21/q − 2(j−i)/q − 2−p/q 2(j−i)/q − β N βi,j Ni Nj − N i−p,j i−p j 21/q − 1 21/q − 1 p=1 j=i+1−Sp
j=1
I
βi,j Ni Nj .
(19)
j=i−S1 +2
Eqs. (13) and (19) can be used to calculate Mi and Ni simultaneously. In the next section, we propose a new discretization for the computation of these quantities. 4. New formulation The new discretized model is based on the cell average technique by Kumar, Peglow, Warnecke, Heinrich, and M¨orl (2005). The entire size domain in this technique is divided into several small sections. The particles within a section are assumed to be concentrated at a representative size. The idea was to take the average of all new born particles within the section and then assign them to the neighboring nodes such that pre-chosen properties are exactly preserved. The final set of discretized equations is given by dNi = Bi−1 λ− ai−1 )H(¯ai−1 − xi−1 ) + Bi λ− ai )H(xi − a¯ i ) + Bi λ+ ai )H(¯ai − xi ) i (¯ i (¯ i (¯ dt + Bi+1 λ+ ai+1 )H(xi+1 − a¯ i+1 ) − Ni i (¯
I k=1
βi,k Nk .
(20)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Here our variant of the Heaviside step function H, as well as λi , Bi and a¯i are given as follows ⎧ 1, x > 0 ⎪ ⎪ ⎨ 1 H(x) = , x=0 ⎪ 2 ⎪ ⎩ 0, x < 0, a − xi±1 , xi − xi±1
λ± i (a) =
j,k vi ≤(xj +xk )
(21)
(22)
j≥k
Bi =
1283
1 1 − δj,k βj,k Nj Nk 2
(23)
and a¯ i =
j≥k
1 Bi
j,k vi ≤(xj +xk )
1 1 − δj,k βj,k Nj Nk (xj + xk ). 2
(24)
The same idea can be adapted for the discretization of TPBE. Since all the particles are assumed to be concentrated at representatives sizes xi ’s, the number and tracer mass density can be expressed as n(t, v) =
I
Nj (t)δ(v − xj )
(25)
j=1
and M(t, v) =
I
Mj (t)δ(v − xj ).
(26)
j=1
Birth term. Substitution of n and M in the birth term of Eq. (10) gives j≥k
Bi,T =
j,k vi ≤(xj +xk )
1 1 − δj,k βj,k (Mj Nk + Nj Mk ). 2
(27)
Analogous to particle birth term in the formulation (20), the assignment of a fraction of this tracer birth and collection of all birth contributions form neighboring cells, the modified tracer birth rate at xi is given by mod Bi,T = Bi−1,T λ− ai−1 )ηi (¯ai−1 )H(¯ai−1 − xi−1 ) + Bi,T λ− ai )ηi (¯ai )H(xi − a¯ i ) + Bi,T λ+ ai )ηi (¯ai )H(¯ai − xi ) i (¯ i (¯ i (¯
ai+1 )ηi (¯ai+1 )H(xi+1 − a¯ i+1 ), + Bi+1,T λ+ i (¯
(28)
xi . a
(29)
where ηi (a) =
The factor η assigns the tracer mass in the same ratio as granule mass. This factor can be illustrated by a simple example, consider the birth contribution of granules Bi−1 λ− ai−1 ) from the cell i − 1 to the cell i. The total volume (mass) fraction carried with these i (¯ particle is given as Bi−1 λ− ai−1 )¯ai−1 . On the other hand, corresponding to Bi−1 λ− ai−1 ) particles we assign Bi−1 λ− ai−1 )xi volume i (¯ i (¯ i (¯ to the cell i. Now we wish to distribute the tracer mass same as granule mass: assigned granule mass Bi−1 λ− ai−1 )xi assigned tracer mass i (¯ = . = − actual granule mass actual tracer mass Bi−1 λi (¯ai−1 )¯ai−1 It follows assigned tracer mass xi = := ηi (¯ai−1 ). actual tracer mass a¯ i−1 Similarly we can construct other two factors appearing in the formulation (28). A detailed discussion about this factor can also be found in Peglow et al. (2006).
1284
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Death term. Substituting n and M in the death term of Eq. (10), we obtain Di,T = Mi
I
β(xi , xj )Nj .
(30)
j=1
The detailed calculations of tracer birth and death terms term has been put in Appendix A. The final set of discrete equations for tracer mass is given as dMi mod = Bi,T − Di,T . dt
(31)
Now the discretized equations (20) and (31) will be used for the calculation of GSD and TMD. These new discretized equations will be compared with the Hounslow’s discretized model by application to several analytically tractable problems in the next section. 5. Numerical results In order to show the effectiveness of the new discretization of TPBE we consider two different types of problems, aggregation in a batch mode of operation and the problem which was considered in Hounslow et al. (2001) for a mixed-suspension mixed-product removal (MSMPR) system. In order to measure the extent of aggregation, the index of aggregation is defined in Hounslow (1990) as ⎧ µ0 (t) ⎪ ⎪ continuous systems, ⎨ 1 − in µ0 (t) Iagg = (32) ⎪ µ (t) ⎪ ⎩1 − 0 batch systems, µ0 (0) where µ0 (t) is the total number of particles at time t (zeroth moment) and µin 0 is the zeroth moment of the feed. 5.1. Aggregation in a batch system We consider two properties of the system: granule volume and the concentration of primary particles in a granule. For the sake of simplicity and the availability of analytical solutions let us take the discrete problem of mono-disperse charge particles with dimensionless size unity. Aldous (1999) provides the particle size distribution for size-independent kernel βi,j = 1 and N0 = 1 as i−1 t 4 . (33) Ni (t) = (t + 2)2 t + 2 Here Ni is the number concentration of clusters containing i primary particles. Here the size of a granule v is identified by the number of primary particles i in the granule. The total number of primary particles in clusters containing i primary particles, say Np,i which corresponds to M(t, v) in continuous setting, is trivial in this case and given by Np,i (t) = iNi (t). The mean volume size of the primary particle distribution in this case is computed by iNp,i = 1 + t, for all t. v¯ p (t) = i i Np,i
(34)
(35)
Now Ni ’s and Np,i have been calculated simultaneously using discretized PBE for GSD and TMD, respectively. Although the analytical solution of Np,i is trivial and there is no need to use DTPBE for this case, we just consider this case in order to check the ability of the new scheme. A comparison of primary particle number distribution has been depicted in Fig. 3. The numerical results are obtained for a coarse geometric grid of the type vi+1 = 2vi in this section since both the formulations produce the same results for a fine grid. The degree of aggregation is chosen to be 0.98 in this case. The modified DTPBE gives overprediction at large size range. This overprediction can clearly be seen once again in corresponding Fig. 4. The mean volume size of the primary particle is shown in this figure. The figures conclude that both the formulations differ significantly to each other. The ratio of primary particle mass to granule mass is plotted in Fig. 5. The ratio, as expected, is the same and constant for both the formulations. Since for a large degree of aggregation analytical solutions for a sum and a product kernel are difficult to calculate, we are plotting only mean volume size in these cases. For the sum kernel βi,j = i + j, the mean volume size takes the following form: iNp,i v¯ p (t) = i = e2t , for all t. (36) N i p,i
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 3. Comparison of primary particle number distribution for constant kernel, Iagg = 0.98.
Fig. 4. Progress of tracer-weighted mean particle volume for size-independent aggregation in a batch system, Iagg = 0.98.
Fig. 5. Final ratio of primary particles computed using the DTPBE and discretized PBE for aggregation in a batch system, Iagg = 0.98.
1285
1286
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 6. Progress of tracer-weighted mean particle volume for aggregation in a batch system (a) sum kernel, Iagg = 0.90 and (b) product kernel, Iagg = 0.35.
For the product kernel βi,j = ij, the mean volume is given by 1 iNp,i v¯ p (t) = i = , 0 ≤ t < 1. 1 − t N i p,i
(37)
The temporal change of the primary particle mean volume for the sum kernel βi,j = i + j, has been drawn in Fig. 6(a). The computation has been carried out for a degree of aggregation Iagg = 0.90. The figure shows that the prediction at very short times is the same for both formulations, but at later times prediction by the new formulation is considerably better than that predicted by the modified DTPBE. In Fig. 6(b), the same comparison is made for a product kernel βi,j = i × j. The results have been obtained for Iagg = 0.35. A similar conclusion can be reached for this case. Table 1 Computation time (s) for both techniques (batch problem) Case
Kernel
q
Hounslow’s technique
New technique
1
Constant
1 2 3
1.04 3.42 7.34
0.68 1.28 2.04
2
Sum
1 2 3
1.09 4.10 9.18
1.07 1.67 5.01
3
Product
1 2 3
0.92 4.12 8.34
0.60 1.57 2.65
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1287
Fig. 7. Decay of total tracer mass for size-independent aggregation in Ilievski and Hounslow’s problem, Iagg = 1/3.
Fig. 8. Progress of tracer-weighted mean particle volume for a size-independent aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/3, (a) q = 1 and (b) q = 3.
1288
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Table 2 Analytical solutions of tracer-weighted mean particle volume β
v¯ T /v0
β0
1+
(u + v)
2Iagg t τ (1 − Iagg )2 1 − Iagg 4Iagg − 2 t Iagg + exp 3Iagg − 1 3Iagg − 1 τ 1 − Iagg
(u · v)
exp
t 1− τ
1 − 8Iagg
2
The computations were carried out in the programming software MATLAB on a Pentium-M processor with 1.6 GHz and 1024 MB RAM. The set of ordinary differential equations resulting from the discretized formulation is solved by a Runge–Kutta fourth order method with adaptive step-size control. We have considered several cases in order to check the computation time for both formulations. A comparison of CPU time taken by both formulations for the same demand of accuracy and other conditions is concluded in Table 1. The table indicates that the computation times are comparable for coarse grids but the modified DTPBE takes considerably more time for fine grids in each case.
Fig. 9. Progress of tracer-weighted mean particle volume for a sum aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/4, (a) q = 1 and (b) q = 3.
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1289
5.2. Aggregation of a mono-disperse feed in MSMPR Hounslow et al. (2001) considered the following problem for the testing of their discretization: ∞ ∂n(t, v) 1 v = β(t, v − , )n(t, v − )n(t, ) d − n(t, v) β(t, v, )n(t, ) d ∂t 2 0 0 n(t, v) ∂n(t, v) + B0 δ(v) − , n(t, 0− ) = 0, = 0, t=0 τ ∂t v ∞ ∂M(t, v) = β(t, v − , )M(t, v − )n(t, ) d − M(t, v) β(t, v, )n(t, ) d ∂t 0 0
(38)
M(t, v) , M(t, 0− ) = 0, M(0, v) = δ(v − v0 ). (39) τ Ilievski and Hounslow (1995) provided the solutions of decay of total tracer mass and tracer-weighted mean particle volume v¯ T for the size-independent kernel β(u, v) = β0 , the sum kernel β(u, v) = β0 × (u + v), and the product kernel β(u, v) = β0 × u × v. The mean volume v¯ T is defined as ∞ vM(t, v) dv v¯ T (t) = 0 ∞ . (40) 0 M(t, v) dv −
Fig. 10. Progress of tracer-weighted mean particle volume for a product aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/9, (a) q = 1 and (b) q = 3.
1290
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
The analytical solution of v¯ T for these kernels have been collected in Table 2. They mentioned that v¯ T gives meaningful values up to one-third degree of aggregation for a sum kernel and is valid up to one-eighth degree of aggregation for a product kernel. The decay of total mass of tracer in each case takes the following form: MT = M0 exp(−t/τ).
(41)
The numerical computations have been made with B0 = τ = v0 = β0 = 1. The degree of aggregation in case of the size-independent constant kernel is taken to be 1/3. Fig. 7 shows that both formulations predict exactly the same total tracer mass and are in excellent agreement with analytical results. The prediction of the total tracer mass is similar and is in excellent agreement with analytical results for all other kernels. In Fig. 8(a), the numerical and analytical temporal change of v¯ T have been compared. The numerical solutions by the new technique are in very close agreement with analytical results while the modified DTPBE gives overprediction at large times. A similar comparison has been made in Fig. 8(b) for a fine grid vi+1 = 21/q vi , q = 3. We can conclude from the figure that both formulations predict almost the same results. Figs. 9 and 10 compare the prediction of tracer-weighted mean particle volume for the sum and the product kernels respectively. Similar findings as before have been observed in this case. The numerical results by the new formulation are in good agreement with analytical results. The prediction is very poor by the modified DTPBE even at small times and a diverging behavior of solution can be observed at large times. The prediction by both formulations get close to each other for fine grids, shown in Figs. 9(b) and 10(b). 6. Conclusions We developed a new discretization for the TPBE which can be applied to any type of grid. Moreover, we compared our new discretization with the modified DTPBE of Peglow et al. (2006) for several problems. It has been found in each case that the new discretization predicts better results than the modified DTPBE. It is more efficient on coarse grids. Furthermore, it has been shown that for a fine grid both formulations predict the same results. Additionally, a comparison of CPU time taken by both formulations concludes that the new technique takes less computational time. Acknowledgements This work was supported by the Graduiertenkolleg-828, ‘Micro–Macro-Interactions in Structured Media and Particle Systems’, Otto-von-Guericke-Universit¨at Magdeburg. The authors gratefully thank for funding through this PhD program. Appendix A All steps in the calculation of birth and death term follows Kumar et al. (2005). A.1. Birth term The birth rate in the ith interval is given as vi+1 v Bi,T = β(t, v − , )M(t, v − )n(t, ) d dv. vi
0
The parameter t is not important here, we omit it in our derivation. Assuming v1 = 0, this term can be rewritten as follows vi+1 vi+1 v i−1 vj+1 Bi,T = β(v − , )M(v − )n() d dv + β(v − , )M(v − )n() d dv. vi
j=1
vj
j=1
vj
vi
vi
(A.1)
Substituting n(v) = Ik=1 Nk δ(v − xk ) and M(v) = Ik=1 Mk δ(v − xk ) in Eq. (A.1), we obtain vi+1 i−1 vj+1 I I Bi,T = β(v − , ) [Mk δ(v − − xk )] [Nk δ( − xk )] d dv vi
+
vi+1
vi
v
vi
β(v − , )
k=1 I k=1
[Mk δ(v − − xk )]
k=1 I k=1
[Nk δ( − xk )] d dv.
(A.2)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1291
Using the definition of Dirac-delta distribution in the first term and changing the order of integration in the second term, we get vi+1 i−1 I Bi,T = β(v − xj , xj ) [Mk δ(v − xj − xk )]Nj dv vi
j=1
+
k=1
vi+1
vi
vi+1
β(v − , )
I
[Mk δ(v − − xk )]
k=1
I
[Nk δ( − xk )] dv d.
The foregoing equation can be further simplified to vi+1 i−1 I Nj Bi,T = β(v − xj , xj ) [Mk δ(v − xj − xk )] dv + vi
j=1
(A.3)
k=1
vi+1
xi
k=1
β(v − xi , xi )
I
[Mk δ(v − xi − xk )]Ni dv.
(A.4)
k=1
Again, applying the definition of Dirac-delta in both the terms, we get Bi,T =
i−1
Nj
j=1
β(xk , xj )Mk +
vi ≤(xj +xk )
β(xk , xi )Mk Ni .
(A.5)
(xi +xk )
Since each term Mj Nk in Eq. (A.5) have a corresponding term Mk Nj except for j = k, Eq. (A.5) can be rewritten in the following way: j≥k 1 1 − δj,k β(xj , xk )(Mj Nk + Mk Nj ). (A.6) Bi,T = 2 vi ≤(xj +xk )
A.2. Death term Analogous to Kumar et al. (2005): vi+1 ∞ Di,T = M(v) β(v, )N() d dv. vi
(A.7)
0
Eq. (A.7), truncated up to vI+1 , can be rewritten as follows vi+1 I vj+1 Di,T = M(v) β(v, )N() d dv. vi
j=1
(A.8)
vj
Again using the Dirac-delta distribution, we obtain vi+1 I I vj+1 I Di,T = [Mk δ(v − xk )] β(v, ) [Nk δ( − xk )] d dv =
=
vi
k=0
I
j=1
Nj
vi+1
vi
j=1 I k=0
vj
vi
k=0
[Mk δ(v − xk )]β(v, xj ) dv =
I
vi+1
β(xi , xj )Nj Mi .
I k=0
[Mk δ(v − xk )]
I
β(v, xj )Nj dv
j=1
(A.9)
j=1
References Aldous, D. J. (1999). Deterministic and stochastic models for coalescence (aggregation and coagulation): A review of the mean-field theory for probabilists. Bernoulli, 5, 3–48. Hounslow, M. J. (1990). A discretized population balance for continuous systems at steady state. American Institute of Chemical Engineering., 36, 106– 116. Hounslow, M. J., Pearson, J. M. K., & Instone, T. (2001). Tracer studies of high shear granulation. II. Population balance modelling. American Institute of Chemical Engineering, 47, 1984–1999. Hounslow, M. J., Ryall, R. L., & Marshall, V. R. (1988). A discretized population balance for nucleation, growth and aggregation. American Institute of Chemical Engineering, 38, 1821–1832. Hulburt, H. M., & Katz, S. (1964). Some problems in particle technology. A statistical mechanical formulation. Chemical Engineering Science, 19, 555– 578. Ilievski, D., & Hounslow, M. J. (1995). Agglomeration during precipitations. II. Mechanism deduction from tracer data. American Institute of Chemical Engineering, 41, 525–535. Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., & M¨orl, L. (2005). Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique. Chemical Engineering Science, 61, 3327–3342.
1292
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Litster, J. D., Smit, D. J., & Hounslow, M. J. (1995). Adjustable discretized population balance for growth and aggregation. American Institute of Chemical Engineering, 41, 591–603. Pearson, J. M. K., Hounslow, M. J., & Instone, T. (2001). Tracer studies of high shear granulation. I. Experimental results. American Institute of Chemical Engineering, 47, 1978–1983. Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2006). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineers, 52, 1326–1332. Wynn Ed, J. W. (1996). Improved accuracy and convergence of discretized population balance of Lister et al. American Institute of Chemical Engineering, 22, 2084–2086.
A discretized model for tracer population balance equation: Improved accuracy and convergence J. Kumar a,c,∗ , M. Peglow b , G. Warnecke a , S. Heinrich c , L. M¨orl c a
Institute for Analysis and Numerics, Otto-von-Guericke-University Magdeburg, Universit¨atsplatz 2, D-39106 Magdeburg, Germany Institute for Process Engineering, Otto-von-Guericke-University Magdeburg, Universit¨atsplatz 2, D-39106 Magdeburg, Germany c Institute of Process Equipment and Environmental Technology, Otto-von-Guericke-University Magdeburg, Universit¨ atsplatz 2, D-39106 Magdeburg, Germany
b
Received 3 June 2005; received in revised form 16 February 2006; accepted 27 February 2006 Available online 12 May 2006
Abstract A new discretization for tracer population balance equation is developed. It is compared to the modified discretized tracer population balance equation of Peglow et al. [Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2006). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineers, 52, 1326–1332]. The new formulation provides excellent prediction of the tracer mass distribution in all test cases. Furthermore, the new formulation is more efficient from a computational point of view, it takes less computational effort and is able to give a very good prediction for a coarser grid. Additionally, it is independent of the type of grid chosen for computation. For finer grids, both formulations tend to produce the same results. The performance of the new formulation is illustrated by the comparison with various analytically tractable problems. Moreover, the new formulation preserves all the advantages of the modified discretized tracer population balance equation of Peglow et al. [Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2005). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineering (published online: 29 December 2005)] and provides a significant improvement in predicting tracer mass distribution and tracer-weighted mean particle volume during aggregation process. © 2006 Elsevier Ltd. All rights reserved. Keywords: Population balance; Discretization; Aggregation; Particle; Tracer; Batch
1. Introduction Tracer studies have become an important tool in the experimental and computational study of disperse phase population dynamics. Ilievski and Hounslow (1995) estimated the agglomeration kinetics parameters from tracer data. Pearson, Hounslow, and Instone (2001) performed some tracer experiments and observed the breakage rate dependence on the size and the age of the granules. Further, Hounslow, Pearson, and Instone (2001) developed a two internal coordinate population balance equation (PBE). The twodimensional PBE has been reduced to two one-dimensional PBEs: one for granule size distribution (GSD) and the other for tracer mass distribution (TMD). They also presented a discretized PBE for tracer. In this work, we develop a new framework for the discretization of tracer population balance equation (TPBE). The discretized tracer population balance equation (DTPBE) presented in Hounslow et al. (2001) which was later modified by Peglow et al. (2006), predicts good results only for a size independent kernel. For sum and product kernels it leads to significant overprediction over exact Abbreviations: DTPBE, discretized tracer population balance equation; GSD, granule size distribution; PBE, population balance equation; TMD, tracer mass distribution; TPBE, tracer population balance equation ∗ Corresponding author. Tel.: +49 391 67 12329; fax: +49 391 67 12129. E-mail address: [email protected] (J. Kumar). 0098-1354/$ – see front matter © 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2006.02.021
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1279
Nomenclature a¯ B, D c Er H I Iagg K n N M q Sp t u, v v0 x
volume average of newborn particles birth and death terms volume of tracer in a granule relative error Heaviside step function total number of intervals the index of aggregation correction factor number density number of particles per unit volume mass density adjustable parameter determining the geometric discretization ratio a parameter defined by Eq. (14) time volume coordinate initial average volume representative volume an interval
Greek symbols β agglomeration rate for reduced one-dimensional problem βˆ agglomeration rate for complete two-dimensional problem δ Dirac-delta distribution η a function defined by Eq. (29) λ a function defined by Eq. (22) µ moments of the PSD τ dimensionless time Subscripts agg agglomeration i index for interval p primary particle T tracer Superscripts ana analytical in feed to a CST mod modified
solutions. Also, the presented model has limitations for the grid type and can be implemented only for special geometric grids. Now the objective is to get a discretized model which is more accurate and can be applied on general grids. 2. Theory A two-dimensional particle size distribution is defined as f (t, v, c), where v and c are two distinct properties, granule volume (size) and tracer volume, respectively. We assume that the density is constant and mass is replaced by volume. Thus, the total number in a domain D is given by D
f (t, v, c) dv dc.
(1)
It should be noted that granule volume v contains volume of tracer and volume of particles, that is, c ≤ v. Fig. 1 shows this two-dimensional particle space. The two-dimensional population balance equation can be obtained by extending the classical one-
1280
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 1. Two-dimensional space.
dimensional PBE (Hulburt & Katz, 1964) to two-dimensional space as ∂f (t, v, c) 1 v min(c,) ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d = β(t, ∂t 2 0 max(0,c−v+) ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d. − β(t, 0
(2)
0
The integral limits corresponding to are trivial. To fix the limits of γ we consider the following apparent relations c − γ ≤ v − , γ ≤ and obviously 0 < γ ≤ c. Clearly these relations provide the integral limits in the formulation (2). The limits of integration can also easily be perceived from Fig. 2. The above 2D PBE must be supplemented with an initial condition: f (0, v, c) = f0 (v, c).
(3)
The first and the second terms on the right-hand side of Eq. (2) account for the formation and the loss of the particles with the properties v and c. Let us assume that the aggregation kernel βˆ depends only on time and size of the granules, but not on the tracer contents within the granules, i.e. βˆ = β(t, v, ). Then the 2D PBE can be converted into two 1D PBEs corresponding to the conventional number density n(t, v) and mass of tracer within granules M(t, v). The number density n(t, v), may be obtained from f by integrating over all possible tracer masses: v n(t, v) = f (t, v, c) dc. (4) 0
The mass of tracer within the granules of size v is given analogously by v M(t, v) = cf (t, v, c) dc. 0
Fig. 2. Limits of integrals in 2D PBE.
(5)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1281
2.1. Reduction to 1D PBE for granule size distribution Integrating Eq. (2) over all possible values of c, we have ∂n(t, v) 1 v v min(c,) ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d dc = β(t, ∂t 2 0 0 max(0,c−v+) v ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d dc. − β(t, 0
0
Reversing the order of integration, we obtain ∂n(t, v) 1 v v−+γ ˆ = β(t, v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dc dγ d ∂t 2 0 0 γ ∞ v ˆ v, , c, γ)f (t, v, c)f (t, , γ) dc dγ d. − β(t, 0
(6)
0
0
(7)
0
ˆ v, , c, γ) = β(t, v, ), then by If we assume that the aggregation kernel depends only on time and the size of both granules, i.e. β(t, substituting c − γ = p and using the definition of density function n, we get ∞ ∂n(t, v) 1 v β(t, v − , )n(t, v − )n(t, ) d − n(t, v) β(t, v, )n(t, ) d. (8) = ∂t 2 0 0 2.2. Reduction to 1D PBE for tracer mass distribution Proceeding as before, multiplication of Eq. (2) by c and integrating over all possible values of c yields: v min(c,) ∂M(t, v) 1 v ˆ v − , , c − γ, γ)f (t, v − , c − γ)f (t, , γ) dγ d dc β(t, c = ∂t 2 0 0 max(0,c−v+) v ∞ ˆ v, , c, γ)f (t, v, c)f (t, , γ) dγ d dc. c − β(t, 0
0
(9)
0
Again, reversing the order of integration and some simple substitution gives v ∞ ∂M(t, v) = β(t, v − , )M(t, v − )n(t, ) d − M(t, v) β(t, v, )n(t, ) d. ∂t 0 0
(10)
Eqs. (8) and (10) are ordinary integro-differential equations which have to be solved numerically. Hounslow, Ryall, and Marshall (1988) developed a discretized method for solving GSD and Hounslow et al. (2001) extended the same idea for the TMD. The discretized PBE was applicable only for special geometric grid. In the next section we first briefly discuss the discretization of TMD and then propose a new idea for the better accuracy and implementation. The preceding population balance equation for the tracer mass can also directly be obtained from the classical one-dimensional population balance equation (8). Note that the tracer mass distribution M(t, v), which represent the total tracer mass contained in particles size v, can be defined as the product of the number density of particles size v with the mass of tracer c(t, v) in each particle: M(t, v) = c(t, v)n(t, v).
(11)
It follows that the rate of change of total tracer mass contained in particles of size v as a results of coagulation is therefore given by ∞ ∂M(t, v) 1 v = β(t, v − , )[c(t, v − ) + c(t, )]n(t, v − )n(t, ) d − n(t, v)c(t, v) β(t, v, )n(t, ) d. (12) ∂t 2 0 0 Using the definition of M(t, v) we obtain the required population balance equation (10). 3. Previous work Hounslow et al. (2001) proposed a discretized PBE for tracer which was based on discretized PBE for GSD (Hounslow et al., 1988). The discretized model could be applied only for fixed geometric grid of the type vi+1 = 2vi . Moreover, the proposed discretization for TMD was not consistent with granule mass. In a recent paper, Peglow et al. (2006) overcame this inconsistency in discretized model by introducing some correction factors and extended it for the adjustable grids of the type vi+1 = 21/q vi . Peglow et al. (2006) used the approach of Litster, Smit, and Hounslow (1995) and Wynn Ed (1996) for the extension to adjustable grids.
1282
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Noting Ni =
vi+1 vi
n(t, v) dv and Mi =
vi+1 vi
M(t, v) dv, the final extended and consistent model is given by
i−S 1 2(j−i+1)/q dMi = βi−1,j (Mi−1 Nj + Ni−1 Mj )K1 + βi−q,i−q Ni−q Mi−q dt 21/q − 1 j=1
+
q−1 i+1−S p+1 21/q − 2(j−i)/q − 2−p/q βi−p,j (Mi−p Nj + Ni−p Mj )K2 21/q − 1 p=1 j=i+1−Sp
+
q
i−Sp
p=2 j=i−Sp −1
+
i−S 1 +1 j=1
2(j−i+1)/q − 1 + 2−(p−1)/q βi−p,j (Mi−p Nj + Ni−p Mj )K3 21/q − 1
2(j−i)/q 1 − 1/q 2 −1
βi,j Ni Mj K4 −
i−S 1 +1 j=1
2(j−i)/q βi,j Mi Nj K5 − 21/q − 1
I
βi,j Mi Nj ,
(13)
j=i−S1 +2
where
q ln(1 − 2−p/q ) Sp = Int 1 − ln 2
(14)
and K1 =
2(i−j)/q , 1 + 2(i−j−1)/q
K2 = K3 =
(15)
2(i−j)/q , 1 + 2(i−j−p)/q
(16)
K4 =
1 , 1 + 2(j−i)/q
(17)
K5 =
21/q − 1 −2i/q + 2(2i−j)/q (21/q − 1) − . 2(j−i)/q 2i/q + 2j/q
(18)
The function Int[x] gives the integer part of x. Let us call this formulation as modified DTPBE. The corresponding discretization of the GSD equation (8), proposed by Litster et al. (1995) and corrected by Wynn Ed (1996), has been summarized as follows: q i−S 1 2(j−i+1)/q dNi = βi−1,j Ni−1 Nj + dt 21/q − 1 j=1
+
i−Sp
p=2 j=i−Sp −1
2(j−i+1)/q − 1 + 2−(p−1)/q 1 2 βi−p,j Ni−p Nj + βi−q,i−q Ni−q 1/q 2 −1 2
q−1 i+1−S i−S 1 +1 p+1 21/q − 2(j−i)/q − 2−p/q 2(j−i)/q − β N βi,j Ni Nj − N i−p,j i−p j 21/q − 1 21/q − 1 p=1 j=i+1−Sp
j=1
I
βi,j Ni Nj .
(19)
j=i−S1 +2
Eqs. (13) and (19) can be used to calculate Mi and Ni simultaneously. In the next section, we propose a new discretization for the computation of these quantities. 4. New formulation The new discretized model is based on the cell average technique by Kumar, Peglow, Warnecke, Heinrich, and M¨orl (2005). The entire size domain in this technique is divided into several small sections. The particles within a section are assumed to be concentrated at a representative size. The idea was to take the average of all new born particles within the section and then assign them to the neighboring nodes such that pre-chosen properties are exactly preserved. The final set of discretized equations is given by dNi = Bi−1 λ− ai−1 )H(¯ai−1 − xi−1 ) + Bi λ− ai )H(xi − a¯ i ) + Bi λ+ ai )H(¯ai − xi ) i (¯ i (¯ i (¯ dt + Bi+1 λ+ ai+1 )H(xi+1 − a¯ i+1 ) − Ni i (¯
I k=1
βi,k Nk .
(20)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Here our variant of the Heaviside step function H, as well as λi , Bi and a¯i are given as follows ⎧ 1, x > 0 ⎪ ⎪ ⎨ 1 H(x) = , x=0 ⎪ 2 ⎪ ⎩ 0, x < 0, a − xi±1 , xi − xi±1
λ± i (a) =
j,k vi ≤(xj +xk )
(21)
(22)
j≥k
Bi =
1283
1 1 − δj,k βj,k Nj Nk 2
(23)
and a¯ i =
j≥k
1 Bi
j,k vi ≤(xj +xk )
1 1 − δj,k βj,k Nj Nk (xj + xk ). 2
(24)
The same idea can be adapted for the discretization of TPBE. Since all the particles are assumed to be concentrated at representatives sizes xi ’s, the number and tracer mass density can be expressed as n(t, v) =
I
Nj (t)δ(v − xj )
(25)
j=1
and M(t, v) =
I
Mj (t)δ(v − xj ).
(26)
j=1
Birth term. Substitution of n and M in the birth term of Eq. (10) gives j≥k
Bi,T =
j,k vi ≤(xj +xk )
1 1 − δj,k βj,k (Mj Nk + Nj Mk ). 2
(27)
Analogous to particle birth term in the formulation (20), the assignment of a fraction of this tracer birth and collection of all birth contributions form neighboring cells, the modified tracer birth rate at xi is given by mod Bi,T = Bi−1,T λ− ai−1 )ηi (¯ai−1 )H(¯ai−1 − xi−1 ) + Bi,T λ− ai )ηi (¯ai )H(xi − a¯ i ) + Bi,T λ+ ai )ηi (¯ai )H(¯ai − xi ) i (¯ i (¯ i (¯
ai+1 )ηi (¯ai+1 )H(xi+1 − a¯ i+1 ), + Bi+1,T λ+ i (¯
(28)
xi . a
(29)
where ηi (a) =
The factor η assigns the tracer mass in the same ratio as granule mass. This factor can be illustrated by a simple example, consider the birth contribution of granules Bi−1 λ− ai−1 ) from the cell i − 1 to the cell i. The total volume (mass) fraction carried with these i (¯ particle is given as Bi−1 λ− ai−1 )¯ai−1 . On the other hand, corresponding to Bi−1 λ− ai−1 ) particles we assign Bi−1 λ− ai−1 )xi volume i (¯ i (¯ i (¯ to the cell i. Now we wish to distribute the tracer mass same as granule mass: assigned granule mass Bi−1 λ− ai−1 )xi assigned tracer mass i (¯ = . = − actual granule mass actual tracer mass Bi−1 λi (¯ai−1 )¯ai−1 It follows assigned tracer mass xi = := ηi (¯ai−1 ). actual tracer mass a¯ i−1 Similarly we can construct other two factors appearing in the formulation (28). A detailed discussion about this factor can also be found in Peglow et al. (2006).
1284
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Death term. Substituting n and M in the death term of Eq. (10), we obtain Di,T = Mi
I
β(xi , xj )Nj .
(30)
j=1
The detailed calculations of tracer birth and death terms term has been put in Appendix A. The final set of discrete equations for tracer mass is given as dMi mod = Bi,T − Di,T . dt
(31)
Now the discretized equations (20) and (31) will be used for the calculation of GSD and TMD. These new discretized equations will be compared with the Hounslow’s discretized model by application to several analytically tractable problems in the next section. 5. Numerical results In order to show the effectiveness of the new discretization of TPBE we consider two different types of problems, aggregation in a batch mode of operation and the problem which was considered in Hounslow et al. (2001) for a mixed-suspension mixed-product removal (MSMPR) system. In order to measure the extent of aggregation, the index of aggregation is defined in Hounslow (1990) as ⎧ µ0 (t) ⎪ ⎪ continuous systems, ⎨ 1 − in µ0 (t) Iagg = (32) ⎪ µ (t) ⎪ ⎩1 − 0 batch systems, µ0 (0) where µ0 (t) is the total number of particles at time t (zeroth moment) and µin 0 is the zeroth moment of the feed. 5.1. Aggregation in a batch system We consider two properties of the system: granule volume and the concentration of primary particles in a granule. For the sake of simplicity and the availability of analytical solutions let us take the discrete problem of mono-disperse charge particles with dimensionless size unity. Aldous (1999) provides the particle size distribution for size-independent kernel βi,j = 1 and N0 = 1 as i−1 t 4 . (33) Ni (t) = (t + 2)2 t + 2 Here Ni is the number concentration of clusters containing i primary particles. Here the size of a granule v is identified by the number of primary particles i in the granule. The total number of primary particles in clusters containing i primary particles, say Np,i which corresponds to M(t, v) in continuous setting, is trivial in this case and given by Np,i (t) = iNi (t). The mean volume size of the primary particle distribution in this case is computed by iNp,i = 1 + t, for all t. v¯ p (t) = i i Np,i
(34)
(35)
Now Ni ’s and Np,i have been calculated simultaneously using discretized PBE for GSD and TMD, respectively. Although the analytical solution of Np,i is trivial and there is no need to use DTPBE for this case, we just consider this case in order to check the ability of the new scheme. A comparison of primary particle number distribution has been depicted in Fig. 3. The numerical results are obtained for a coarse geometric grid of the type vi+1 = 2vi in this section since both the formulations produce the same results for a fine grid. The degree of aggregation is chosen to be 0.98 in this case. The modified DTPBE gives overprediction at large size range. This overprediction can clearly be seen once again in corresponding Fig. 4. The mean volume size of the primary particle is shown in this figure. The figures conclude that both the formulations differ significantly to each other. The ratio of primary particle mass to granule mass is plotted in Fig. 5. The ratio, as expected, is the same and constant for both the formulations. Since for a large degree of aggregation analytical solutions for a sum and a product kernel are difficult to calculate, we are plotting only mean volume size in these cases. For the sum kernel βi,j = i + j, the mean volume size takes the following form: iNp,i v¯ p (t) = i = e2t , for all t. (36) N i p,i
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 3. Comparison of primary particle number distribution for constant kernel, Iagg = 0.98.
Fig. 4. Progress of tracer-weighted mean particle volume for size-independent aggregation in a batch system, Iagg = 0.98.
Fig. 5. Final ratio of primary particles computed using the DTPBE and discretized PBE for aggregation in a batch system, Iagg = 0.98.
1285
1286
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Fig. 6. Progress of tracer-weighted mean particle volume for aggregation in a batch system (a) sum kernel, Iagg = 0.90 and (b) product kernel, Iagg = 0.35.
For the product kernel βi,j = ij, the mean volume is given by 1 iNp,i v¯ p (t) = i = , 0 ≤ t < 1. 1 − t N i p,i
(37)
The temporal change of the primary particle mean volume for the sum kernel βi,j = i + j, has been drawn in Fig. 6(a). The computation has been carried out for a degree of aggregation Iagg = 0.90. The figure shows that the prediction at very short times is the same for both formulations, but at later times prediction by the new formulation is considerably better than that predicted by the modified DTPBE. In Fig. 6(b), the same comparison is made for a product kernel βi,j = i × j. The results have been obtained for Iagg = 0.35. A similar conclusion can be reached for this case. Table 1 Computation time (s) for both techniques (batch problem) Case
Kernel
q
Hounslow’s technique
New technique
1
Constant
1 2 3
1.04 3.42 7.34
0.68 1.28 2.04
2
Sum
1 2 3
1.09 4.10 9.18
1.07 1.67 5.01
3
Product
1 2 3
0.92 4.12 8.34
0.60 1.57 2.65
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1287
Fig. 7. Decay of total tracer mass for size-independent aggregation in Ilievski and Hounslow’s problem, Iagg = 1/3.
Fig. 8. Progress of tracer-weighted mean particle volume for a size-independent aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/3, (a) q = 1 and (b) q = 3.
1288
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Table 2 Analytical solutions of tracer-weighted mean particle volume β
v¯ T /v0
β0
1+
(u + v)
2Iagg t τ (1 − Iagg )2 1 − Iagg 4Iagg − 2 t Iagg + exp 3Iagg − 1 3Iagg − 1 τ 1 − Iagg
(u · v)
exp
t 1− τ
1 − 8Iagg
2
The computations were carried out in the programming software MATLAB on a Pentium-M processor with 1.6 GHz and 1024 MB RAM. The set of ordinary differential equations resulting from the discretized formulation is solved by a Runge–Kutta fourth order method with adaptive step-size control. We have considered several cases in order to check the computation time for both formulations. A comparison of CPU time taken by both formulations for the same demand of accuracy and other conditions is concluded in Table 1. The table indicates that the computation times are comparable for coarse grids but the modified DTPBE takes considerably more time for fine grids in each case.
Fig. 9. Progress of tracer-weighted mean particle volume for a sum aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/4, (a) q = 1 and (b) q = 3.
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1289
5.2. Aggregation of a mono-disperse feed in MSMPR Hounslow et al. (2001) considered the following problem for the testing of their discretization: ∞ ∂n(t, v) 1 v = β(t, v − , )n(t, v − )n(t, ) d − n(t, v) β(t, v, )n(t, ) d ∂t 2 0 0 n(t, v) ∂n(t, v) + B0 δ(v) − , n(t, 0− ) = 0, = 0, t=0 τ ∂t v ∞ ∂M(t, v) = β(t, v − , )M(t, v − )n(t, ) d − M(t, v) β(t, v, )n(t, ) d ∂t 0 0
(38)
M(t, v) , M(t, 0− ) = 0, M(0, v) = δ(v − v0 ). (39) τ Ilievski and Hounslow (1995) provided the solutions of decay of total tracer mass and tracer-weighted mean particle volume v¯ T for the size-independent kernel β(u, v) = β0 , the sum kernel β(u, v) = β0 × (u + v), and the product kernel β(u, v) = β0 × u × v. The mean volume v¯ T is defined as ∞ vM(t, v) dv v¯ T (t) = 0 ∞ . (40) 0 M(t, v) dv −
Fig. 10. Progress of tracer-weighted mean particle volume for a product aggregation kernel in Ilievski and Hounslow’s problem, Iagg = 1/9, (a) q = 1 and (b) q = 3.
1290
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
The analytical solution of v¯ T for these kernels have been collected in Table 2. They mentioned that v¯ T gives meaningful values up to one-third degree of aggregation for a sum kernel and is valid up to one-eighth degree of aggregation for a product kernel. The decay of total mass of tracer in each case takes the following form: MT = M0 exp(−t/τ).
(41)
The numerical computations have been made with B0 = τ = v0 = β0 = 1. The degree of aggregation in case of the size-independent constant kernel is taken to be 1/3. Fig. 7 shows that both formulations predict exactly the same total tracer mass and are in excellent agreement with analytical results. The prediction of the total tracer mass is similar and is in excellent agreement with analytical results for all other kernels. In Fig. 8(a), the numerical and analytical temporal change of v¯ T have been compared. The numerical solutions by the new technique are in very close agreement with analytical results while the modified DTPBE gives overprediction at large times. A similar comparison has been made in Fig. 8(b) for a fine grid vi+1 = 21/q vi , q = 3. We can conclude from the figure that both formulations predict almost the same results. Figs. 9 and 10 compare the prediction of tracer-weighted mean particle volume for the sum and the product kernels respectively. Similar findings as before have been observed in this case. The numerical results by the new formulation are in good agreement with analytical results. The prediction is very poor by the modified DTPBE even at small times and a diverging behavior of solution can be observed at large times. The prediction by both formulations get close to each other for fine grids, shown in Figs. 9(b) and 10(b). 6. Conclusions We developed a new discretization for the TPBE which can be applied to any type of grid. Moreover, we compared our new discretization with the modified DTPBE of Peglow et al. (2006) for several problems. It has been found in each case that the new discretization predicts better results than the modified DTPBE. It is more efficient on coarse grids. Furthermore, it has been shown that for a fine grid both formulations predict the same results. Additionally, a comparison of CPU time taken by both formulations concludes that the new technique takes less computational time. Acknowledgements This work was supported by the Graduiertenkolleg-828, ‘Micro–Macro-Interactions in Structured Media and Particle Systems’, Otto-von-Guericke-Universit¨at Magdeburg. The authors gratefully thank for funding through this PhD program. Appendix A All steps in the calculation of birth and death term follows Kumar et al. (2005). A.1. Birth term The birth rate in the ith interval is given as vi+1 v Bi,T = β(t, v − , )M(t, v − )n(t, ) d dv. vi
0
The parameter t is not important here, we omit it in our derivation. Assuming v1 = 0, this term can be rewritten as follows vi+1 vi+1 v i−1 vj+1 Bi,T = β(v − , )M(v − )n() d dv + β(v − , )M(v − )n() d dv. vi
j=1
vj
j=1
vj
vi
vi
(A.1)
Substituting n(v) = Ik=1 Nk δ(v − xk ) and M(v) = Ik=1 Mk δ(v − xk ) in Eq. (A.1), we obtain vi+1 i−1 vj+1 I I Bi,T = β(v − , ) [Mk δ(v − − xk )] [Nk δ( − xk )] d dv vi
+
vi+1
vi
v
vi
β(v − , )
k=1 I k=1
[Mk δ(v − − xk )]
k=1 I k=1
[Nk δ( − xk )] d dv.
(A.2)
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
1291
Using the definition of Dirac-delta distribution in the first term and changing the order of integration in the second term, we get vi+1 i−1 I Bi,T = β(v − xj , xj ) [Mk δ(v − xj − xk )]Nj dv vi
j=1
+
k=1
vi+1
vi
vi+1
β(v − , )
I
[Mk δ(v − − xk )]
k=1
I
[Nk δ( − xk )] dv d.
The foregoing equation can be further simplified to vi+1 i−1 I Nj Bi,T = β(v − xj , xj ) [Mk δ(v − xj − xk )] dv + vi
j=1
(A.3)
k=1
vi+1
xi
k=1
β(v − xi , xi )
I
[Mk δ(v − xi − xk )]Ni dv.
(A.4)
k=1
Again, applying the definition of Dirac-delta in both the terms, we get Bi,T =
i−1
Nj
j=1
β(xk , xj )Mk +
vi ≤(xj +xk )
β(xk , xi )Mk Ni .
(A.5)
(xi +xk )
Since each term Mj Nk in Eq. (A.5) have a corresponding term Mk Nj except for j = k, Eq. (A.5) can be rewritten in the following way: j≥k 1 1 − δj,k β(xj , xk )(Mj Nk + Mk Nj ). (A.6) Bi,T = 2 vi ≤(xj +xk )
A.2. Death term Analogous to Kumar et al. (2005): vi+1 ∞ Di,T = M(v) β(v, )N() d dv. vi
(A.7)
0
Eq. (A.7), truncated up to vI+1 , can be rewritten as follows vi+1 I vj+1 Di,T = M(v) β(v, )N() d dv. vi
j=1
(A.8)
vj
Again using the Dirac-delta distribution, we obtain vi+1 I I vj+1 I Di,T = [Mk δ(v − xk )] β(v, ) [Nk δ( − xk )] d dv =
=
vi
k=0
I
j=1
Nj
vi+1
vi
j=1 I k=0
vj
vi
k=0
[Mk δ(v − xk )]β(v, xj ) dv =
I
vi+1
β(xi , xj )Nj Mi .
I k=0
[Mk δ(v − xk )]
I
β(v, xj )Nj dv
j=1
(A.9)
j=1
References Aldous, D. J. (1999). Deterministic and stochastic models for coalescence (aggregation and coagulation): A review of the mean-field theory for probabilists. Bernoulli, 5, 3–48. Hounslow, M. J. (1990). A discretized population balance for continuous systems at steady state. American Institute of Chemical Engineering., 36, 106– 116. Hounslow, M. J., Pearson, J. M. K., & Instone, T. (2001). Tracer studies of high shear granulation. II. Population balance modelling. American Institute of Chemical Engineering, 47, 1984–1999. Hounslow, M. J., Ryall, R. L., & Marshall, V. R. (1988). A discretized population balance for nucleation, growth and aggregation. American Institute of Chemical Engineering, 38, 1821–1832. Hulburt, H. M., & Katz, S. (1964). Some problems in particle technology. A statistical mechanical formulation. Chemical Engineering Science, 19, 555– 578. Ilievski, D., & Hounslow, M. J. (1995). Agglomeration during precipitations. II. Mechanism deduction from tracer data. American Institute of Chemical Engineering, 41, 525–535. Kumar, J., Peglow, M., Warnecke, G., Heinrich, S., & M¨orl, L. (2005). Improved accuracy and convergence of discretized population balance for aggregation: The cell average technique. Chemical Engineering Science, 61, 3327–3342.
1292
J. Kumar et al. / Computers and Chemical Engineering 30 (2006) 1278–1292
Litster, J. D., Smit, D. J., & Hounslow, M. J. (1995). Adjustable discretized population balance for growth and aggregation. American Institute of Chemical Engineering, 41, 591–603. Pearson, J. M. K., Hounslow, M. J., & Instone, T. (2001). Tracer studies of high shear granulation. I. Experimental results. American Institute of Chemical Engineering, 47, 1978–1983. Peglow, M., Kumar, J., Warnecke, G., Heinrich, S., M¨orl, L., Hounslow, M. J. (2006). Improved discretized tracer mass distribution of Hounslow et al. American Institute of Chemical Engineers, 52, 1326–1332. Wynn Ed, J. W. (1996). Improved accuracy and convergence of discretized population balance of Lister et al. American Institute of Chemical Engineering, 22, 2084–2086.