- Email: [email protected]
Numerical treatment of the population balance equation using a Spline-Galerkin method
Comprrtcls chmt. Engng Vol. 18, No. 9. pp. 775-783, 1994 Copyright @ 1994 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0098-1354/94 $7.00 + 0.00 009%1354(94)EOO@7-A
NUMERICAL TREATMENT OF THE POPULATION BALANCE EQUATION USING A SPLINE-GALERKIN METHOD L. D.
I Department
ERASMUS,’
D. EYRE’ and R. C.
EVERSON’
of Mathematics and Applied Mathematics and ‘Department Potchefstroom University for CHE, Potchefstroom, Republic
(Received
31 August
of Chemical Engineering, of South Africa
1993; final revision received 24 January 1994; received for publication I February 1994)
Abstract-This paper describes a numerical technique for solving the Lifshitz-Slyozov equation of continuity which applies to certain mass transfer proceses. The Lifshitz-Slyozov equation which describes a mechanism involving the transfer of atoms (or undissociated molecules) from smaller particles to larger particles dispersed in a supersaturated medium (Ostwald ripening) is also considered together with collision and subsequent coalescence of particles (ripening plus collection). Special attention is given to the case in which the net growth rate appearing in the Lifshitz-Slyozov equation is a nonsmooth function of the type I(v) 0~v 8, where v is the particle volume and 0
INTRODUCTION
Lifshitz and Slyozov (LS) (1961) equation is used to describe a number of physical and chemical processes involving the decay and growth of particles due to the transfer of atoms or undissociated molecules in a supersaturated environment, known as Ostwald ripening. Examples of these processes include precipitation of solutes from supersaturated solutions and the sintering of supported metal catalysts on an inert substate. Although more detailed microscopic descriptions of the equations of continuity have been derived for Ostwald ripening (Dadyburjor and Ruckenstein, 1977) these equations reduce to the macroscopic LS equation after approximation. Ruckenstein and Dadyburjor (1977) also proposed a combination of Ostwald ripening plus collision and subsequent coalescence that occurs during sintering of supported catalyst particles. Collision and subsequent coalescence will be referred to as collection and the combination as ripening plus collection. Let n(u, t) du be the number of particles with volume between D and u+ dv at time t. The LS equation is a continuity equation: The
775
where Z(u) = in is the growth rate of a particle of volume u, and S( u, C,n) a source term for particles of volume U. In general the source term can be a nonlinear function of the dependent variable n(n, t), as, for example, in the case where particles can be formed by collection. It will be assumed that equation (1) satisfies the usual boundary conditions, namely: n(o, 0) = W(U), n(0, t) = 0.
(2) (3)
The interpretation of these conditions is that there are no particles of zero size, and that the initial distribution of particles wO(u) is given. The total number density MO(t) and volume density M,(t) is given by the moments: MO(t) =
M,(t)=
- n(u, t) du, I 0
(4)
-n(v,r)udu,
(5)
I
0
respectively. If collection does not occur (S = 0) then the total number density of particles is conserved, so that A&(t) =A$, where No denotes the total initial number density _
L.D.
776
ERASMUS
Collection can be included in the LS equation by considering a source term given by the nonlinear expression :
1
”
2I 0 -n(u,
n(u’, t)n(u - u’, t)K(u’,
1)
00 n(u’,
t)K(u’,
u -u’)
u) du’.
I 0
du’
(6)
Here the collection kernel K(u’, u) describes the geometry and dynamics of collection mechanism. The equations (1) and (6) have been extensively investigated by many authors. A detailed discussion of the effect of coupling between collection and the LS equation can be found, for example, in the work of Ramabhadran et al. (1976) and Peterson et al. (1978). Of interest in the present study is a growth law appearing in the LS equation that has the form: I(u) = o/3v@ where O
1,
(7)
where as is the growth constant. A negative value of ua in equation (1) would correspond to decay rather than growth. The exponent p is determined from physical considerations. For example, the case fl= l/3 corresponds to the well known continuum diffusion growth law based on particle diameter. This can be seen by changing the independent variable in the equation of continuity from a volume v to the particle diameter D (for spherical particles). In this new variable the growth law can be written as:
et al.
easily be found and one must resort to numerical methods. Ramkrishna (1985) provides a survey of numerical methods for solving the population balance equation. A collocation method using a finite element basis of cubic splines has been described by Gelbard and Seinfeld (1978), and by Eyre et al. (1988). Viljoen et al. (1990) use the splinecollocation method to obtain numerical solutions of equations that describe decay plus collection. In the case of Langmuir (1918) evaporation, however, Viljoen et al. (1990) reported numerical instability of the spline-collocation method. This result suggests that the collocation approach is not suited to cases in which the growth law has a nonsmooth behaviour; as in the case of equation (7) with O
Before applying the projection method it is first convenient to introduce dimensionless variables and to map the volume variable v onto a finite interval. For this purpose it is convenient to introduce a nonlinear mapping:
(8) The case of 0 C/l < 1 is interesting from a numerical point of view because it is an example of a nonsmooth function (the function has an infinite derivative at v = 0). This behaviour in the growth law leads to a solution of the LS equation that is also a nonsmooth function. In the case of zero collection the exact solution of the LS equation for t>O is:
XW,
[( v
l-(1-&$
This solution exists provided that: v’-fl==(l
-@)a&
u=5
(IO)
When coupled to the particulate balance equation, an exact solution of the LS equation for realistic collection kernels, similar to solution (9), cannot
1+x l-x
( > ’
(11)
where x E [ - 1, 11. (An alternative mapping such as one based on a log or tan function is also possible and will lead to essentially the same results as those presented in this paper.) It is convenient also to introduce a dimensionless time t = A$$, where x is a characteristic rate constant of the process. The dependent variable n(u, t) is now transformed into a function of dimensionless arguments:
141-B) > 1. (9)
METHOD
?i(x,r)=--
5
2
No (1 -x)2
n(u, 0.
The initial function wO(v) is similar transformed so that:
2
iG”C-4
=G
Cl
_x)2
w(v)-
(13)
The parameter c is chosen so that the function fi(x, t) is reasonably distibuted over the interval
Treatment of the population balance equation
777
[ - 1, 11. This parameter remains fixed throughout the calculation. The transformed function fi(x, t) satisfies the LS equation:
The (unknown) solution of the LS equation is approximated by the linear combination
(14)
Substituting this approximation for the solution in the LS equation leads to a residual function:
where the dimensionless growth rate a(x) by:
is given
(22)
r(x, r) = m+l
a(x) = &*5@-‘(1 For the zero-collection ized so that:
+x)B(1 -X)*-8.
(15)
case rZ(x, t) is now normal-
C {Bj(x)fi(r) + [Bi(x)a’(x) + Bl(xMx)lJ(r)) j=O
i=U
I I --I
ii(x,t)dX=l.
(16)
I
ii@‘, T)ti(U, r)y(x’, U) --1
z)
I
l-u l_x
2
( >
dx’
1
--I
A(x’,
r)y(x’,
x) dx’,
(17)
where
r(n’,x)=
K(v’, 0)
08)
x
and u= u(x’, x) is the transformed ence , namely:
volume differ-
2(x-X’)-(l-X)(1-x’) (19)
U=2(x-x’)+(1-x)(l-x’)~
Discretization of equations (14) and (17) is carried out on a space of cubic splines. Cubic B-splines on x E [ - 1, l] are constructed as follows: the interval is partitioned by m nodal points {xi}K1 where - 1 = Xl-=ZX,< - . -
B!“(X)
=
o
Xj-lCXsXi,
(20)
otherwise.
7
The B-splines of order I (degree I- 1) are generated by the stable iterative method of Cox (1972) and deBoor (1972): *c,‘(x)
=
,
(x
--+,m~_,“w - (xi--wY-‘Yx)~
(21)
Xi-X;_/
Note that the index on the cubic B-splines runs from i=2,. . . , m + 3. A convenient notation is Bi (x) = B&(x) to denote the cubic B-spline which is nonzero over the interval (x~_~,_x~+~)so that the index now runs from i=O,. . . ,m+l.
2
j=(J
X
x
Bi(“)f,
(t)Bj(x’lfi(r)Y(x’,
u>
-I
m+l
2
l-_x
da?
(
x
-A@,
1
f-u x
The transformed collection term is
S(X, .,+;
xx-j
m+lm+l
-
>
+c
m+l c
idI
Ix
Bj(X')fi(t)y(X',
~icax~)
j=O
X)
dX’.
W)
--L
Note that the coefficients vi(Z)} are functions of the dimensionless time z. Where no confusion should arise this explicit dependence on time will be suppressed. Clearly the residual function r(x, r) should be small if the spline function (22) is to be a useful approximation to the solution of the original integro-differential equation. Let (q, r) denote the usual inner product: (+, r) =
1 Il(x)e
I -1
b,
w
where O(X) is some appropriately chosen weight factor. Furthermore, let {v*} be a chosen set of (m+2) test functions. The coefficients vi} satisfy the (m+ 1) system of ordinary differential equations: (7&,r)=O;
k=l,.
. . ,(m+l),
(25)
and the boundary condition (3). The approach usually taken when using spline methods for solving population balance equations is a collocation method. In the collocation method the test functions are chosen to be delta functions, ~~=6(x - u,), where uk E [ - 1, 11. An alternative approach, and the one taken in the present paper, is to choose the test functions to be cubic B-splines, t&k= Bk. This is the Galerkin method. Thus the residual is chosen to be zero only in the sense of an integral over the cubic B-splines. The problem now reduces to that of solving the system of ordinary differential equations: Lf(r)=-ivfr)+f(r)Cf(r),
(26)
L. D. ERASMUSet al.
778 7
2.0
6
AQ? 1.8
5
1.6
1.4
-0.6
-0.8
-0.4
0.4
0.2
0
-0.2
0.6
x Fig. 1.
Computed solution of the
LS
equation using the spline-Galerkin method. Since:
for the expansion coefficients f = (f&f,, . . . ,f,,,+t)T. The matrices in equation (26) are given by:
a’(x) Li.j=
'
-I
I
NjJ=
Bi(~)Bj(x)w(x)~9
=
%&_
B
[
(27)
(1-x)2-p (1 +x)*-B - (2 -/I)(1
1
I
-I
+
' I
+x)B(l
Bi(x)B;(n)a(x)w(x)dx
Bi(X)Bj(X)U'(X)CO(X)dX
(28)
-1
and
-x)‘-8
19
(30)
the integrand in the expression for Ni.j contains a weak singularity at x = - 1. This singularity can be removed by chasing an appropriate weight factor w(x) = (1 +x)‘-B. It should be remarked that this Table 1. Exact and computed moments MO(r) and M,(s) for the LS equation
' I
Bk(X)W(X)dX
-I
x I
Bi(U)Bj(X')y(X',
U)
MU
-1
5
0.2 0.4 0.6 0.8 1.0
l.orm 1.000 1.000 1.000 1.000 1.000
:.: 1:6 1.8 2.0
1.ooo l.OOU l.ooo l.oM) l.OOcl
0.0
X
l
--I
Bj(X’)y(X’,
X) dX’.
(29)
The integrals in equations (27-29) are evaluated numerically using a standard Gauss-Legendre quadrature formula.
(exact)
MU
(computed) l.OW l.ooo 1.000 l.uoO l.lXKl i.ooo l.olM l.ocKI l.ooO l.ooo l.tmO
M,
M,
(exact)
(computed)
2.ooo 3.224 4.656 6.265 8.036 9.951 12.002 14.180 16.478 18.887 21.406
2.000 3.067 4.412 5.979 7.733 9.651 11.714 13.912 16.240 18.674 21.219
Treatment of the population balance equation
779
0.8
0.6
-0.2
-0.4
-0.6
-0.8 __ -1
.I)
I
I
I
I
I
I
I
-0.8
I
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
x
Fig. 2. Error in the computed solution, ii(computed) - rl(exact).
2.5
2.0
1.5 G 5 IC 1.0
0.5
0 -“.t.
-‘.U
-“.4
-0.2
0
0.2
0.4
0.6
I Fig. 3. Computed
solution of the collection equation using the spline-Galerkin
method.
0.8
L. D. EWMUS et al.
780
Table2. Exact and computed moments M&) and M,(s) for the collection equation
T
MI (computed)
MO (exact)
MI (exact)
Ml (computed)
0.913 0.839
2.000 2.oiKl 2.000
2.000 2.aXl
0.723 0.675 0.634 0.597 0.564
?E 2:Oo0 2.Olm 2.wm 2.MKJ
0.508 0.535
f:E
where
1 Ai.j=
(33)
dx9
Bi(x)Bj(x) -I
0.0
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
1.ooo
1.000 0.909 0.833 0.769 0.714 0.667 0.625 0.588 0.556 0.527 0.500
0.7n
f:E 2.000 2.ooo 2.aKl 2.000 2.aiw 2.ooo 2.oLul
approach is used to remove unwanted singularities in the integrand but that the integrand still remains a non-smooth function with possible singularities in the derivative. The system of equations (26) is subject to the initial condition:
gi=
Next we turn to the question of temporal discretization. We introduce the time step Ar. Time is now expressed in terms of an index r = r/At, while the coefficient fi, , denotes an approximation to J(t) , A classical one-step e-method gives the (m + 1) algebraic equations: m+l
m+l
NO)
-0.8
Fig. 4. Computed
-0.6
-0.4
-0.2
#PI+,
-I
L~.jfj.r+(l-8)Ar i=O
(31)
= g,
r
m+l
m+l
Any standard method, such as collocation, could be used to determine the coefficients (f;(O)}. However, in keeping with the Galerkin approach, the coefficients (f,(O)} are obtained by solving the matrix equation:
Bi (X)WO(X)dX*
-1
=c
~fit;-(")Bj(x)=~OSx)-
Il
X
=l i-0
Nk.i-
2 j=O
(35)
fi,r,
Ck,i.jf;..r 1
where k=l,. . . ,m+l. The value of B can be varied from 0= 0 (explicit scheme) to 0 = 1 (fully implicit scheme). The value of 8= l/2 gives the Crank-Nicolson method.
0
0.2
0.4
0.6
0.8
solution of the combined ripening plus collection equation using the spline-Galerkin method.
1.0
Treatment
of the population
'r
781
equation
chosen. A condition that satisfies the boundary conditions (2) and (3) is:
6-
/J ‘1
SF &
4-
I’ I :
1c
3-
:
2-
\
\
2v 4& v exp -, wo( u) = vo> ua (
\ \ \ \ \
where
0.2
0.3
0.4
0.5
0.6
0.7
0.1
I 0.9
Fig. 5. Computed solution after a time t = 2. The solid line shows the double distribution for the combined ripening plus collection equation. The broken line shows the result for ripening alone.
Supplementing equation boundary conditions:
(35) with the discrete
(36)
~fi,r+tsj(_1)=09
j=O
yields an algebraic system for (m +2) expansion coefficients fi, r+l. (Note that only four terms are summed in equation (36) because each subinterval is spanned by only four nonzero basis functions.) In the case of collection this is a nonlinear system which is solved using a standard quasi-Newton solver (ZSPOW from the IMSL package). The set fi_, is used as an initial estimate.
3. NUMERICAL
EXAMPLES
A number of simple test problems will be considered. For this purpose a simple initial condition is
Table 3. Computed moments M,,(r) and M,(r) for the combination of ripening plus collection equation
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
is
the
average
particle
We begin with a numerical solution of the LS equation (ripening alone). The LS equation with growth law given by (7) has the solution:
x
0.0
M,,(O) = No and u.
(37)
volume.
1O0.1
balance
Ez
0:829 0.748 0.684 0.635 0.589 0.547 0.515 0.492 0.472
2.ooo 3.181 4.412 5.686 7.028 8.462 10.005 11.657 13.423 15.303 17.311
(1 +Bw -8)
n(u,t)=~
4Nol-(1-#3)$$ [ 1
For numerical calculations we have chosen the following parameters; the initial number of particles N,= 1 and the initial average volume of particles u0= 2. The initial condition now becomes: we(u) = u e-“.
(39)
We take /3= * and u,,~ = 1.547. Figure 1 shows the number of particles ii as a function of x for several values of the dimensionless time t. The curves represent a numerical solution of the LS equation using a basis of 19 cubic B-splines. The mapping parameter is 5 = 10 and the time step is At = 0.01. In this example there is a rapid growth of the small particles. In fact the left-hand-side of the particle spectrum is finite only for particle sizes which satisfy inequality (10). The larger particles grow more slowly as indicated by the small change on the right-hand-side of the particle spectrum. It should be remarked that only the positive part of the function is shown in the diagram. In practice a small negative part of the function appears on the lefthand-side of the particle spectrum. This negative value is only a few percent of the peak value and is the result of representing a nonsmooth function by a smooth cubic spline approximation. The negative values are set to zero for calculations at subsequent time intervals. The calculated moment using 19 spline basis functions and the exact moment MI are shown in Table 1. The calculated moment shows a small error (< 10%) which appears initially in the first step of the calculation and then remains constant throughout subsequent time steps. This error is reflected also in the difference between the numerical solution and the exact solution [given by equation (38)] as shown in Fig. 2. Again this error could be attributed to the smoothness properties of the cubic spline approximation. Next we turn to the collection equation. For the purpose of demonstrating the performance of the
L. D. ERASMUS er (II.
782
algorithm we consider the case of a constant kernel, K(u’, v) =$_ In this case the collection equation (Scott, 1968) has the solution:
I-_-
(W
where
MO=-
2No 2+gN,,t
and
M1=NovO.
(41)
Notice that in the case of collection alone the total number of particles is not conserved but that the total volume of particles is conserved. The solution of the collection equation should be well represented by a cubic spline approximation. We choose 5 = 1. Figure 2 shows the computed solution using 19 spline basis functions for several values of the time t. These curves are indistinguishable from the exact solution_ In this example there is little growth of the smaller particles, but rapid growth of the larger particles producing a long tail in the large particle end of the spectrum. It is clear that the effect of collection, while leading to growth in the average particle size, has a much different effect on the particle spectrum that does the growth due to ripening (LS equation). Table 2 shows the computed moments Mu and MI. These results agree extremely well with the exact values and illustrate the accuracy of the spline-Galerkin method for solving the collection equation. The numerical solution of the above two examples has been used to test the numerical procedures. An exact solution of the ripening plus collection equation for the growth law (6) and a constant (nonzero) collection kernel is not known. Thus one must rely on numerical approximations to obtain a solution of the integro-differential equation. A dimensionless measure of the collection/ripening rate is given by the ratio A = ~No~‘-~lu~. The relationship between g= 1 and cLj3 is such that A = 1 for a mapping parameter c= 10. The case A = 1 corresponds to particle growth from equal rates due to ripening and collection. Figure 4 shows the computed particle number spectrum. The rapid growth of the smaller particles can be attributed almost entirely to the effect of ripening. The rapid growth of the larger particles is an enhancement of that seen for the case of collection. An interesting feature of this spectrum is the formation of a double bump. Figure 5 shows the particle spectrum at time t = 2. The solution for both ripening and coltection (solid line) and the result for
ripening alone (broken line) are shown. The first peak in the solid line corresponds to the peak in the ripening spectrum. Table 3 shows the computed moments M,(t) and Ml(r). Although the exact solution for this last example is not known, for a combination of ripening and collection there is an expected decrease in the total number of particles MO and an increase in the total volume of particles Ml with increasing time. It can be remarked that calculations using more realistic kernels that describe sintering, such as those given by Ruckenstein and Pulvermacher (1973), also produce a double bump when combined with the growth law in equation (6). 4. CONCLUSIONS
Numerical solutions of the LS equation have been obtained in a case where the growth rate is a nonsmooth function of the particle volume. Numerical solutions for the combination of ripening plus collection have also been obtained from a numerical method that combines cubic splines and a Galerkin technique to determine the spline expansion coefficients. A novel feature of the ripening plus collection equation is the emergence of a double distribution in the particle size spectrum. NOMENCLATURE
Ai.i = Defined by equation (33) Bi = Cubic B-splines C *, r.i = Defined by equation (29) D = Particle diameter J(z) = Spline expansion coefficient f;. , = Expansion coefficient at discrete time step g, = Defined by equation (34) I( U) = Growth rate K(u’, U) = Collection kernel Li.i = Defined by equation (27) m = Number of nodal points MO = Number density defined by equation (5) MI = Volume density defined by equation (6) n(u, t) du = Average number of particles per unit volume NO = Initial average number of particles per unit volume Ni, j = Defined by equation (28) r = Residual function defined by equation (23) S = Source term defined by equation (6) t = Time u = Defined by equation (19) u, = Collocation points v = Particle volume B = Average particle volume u0= Initial average particle volume w,, = Initial volume distribution x = Dimensionless variable X, = Nodal points Greek letters a = Dimensionless growth rate /3= Exponent in growth term y = Dimensionless collection kernel 5 = Volume parameter
Treatment 8= A = 5= a, = T= x= vi = w =
of the population
Parameter in the e-method on rate Dimensionless ripening/collecti Constant collection kernel Growth parameter Dimensionless time Characteristic rate constant Test function Weight factor in equation (24) REFERENCES
Cox M. G.. The numerical evaluation of B-solines. J. Inst. _ Math. Abpf. 10, 134-149 (1972). Dadyburjor D. B. and E. Ruckenstein. Kinetics of Ostwald ripening. J. Cryst. Growth 40. 279-290 (1977). with B-splines. 1. Approx. deBoor C., On calculating Theory 6, 50-62 (1972). Evre D.. C. J. Wright and G. Reuter. Soline-collocation *with adaptive m&h grading for sob&g the stochastic collection equation. J. Comput. Phys. 78, 288-304 (1988). Gelbard F. and J. H. Seinfeld, Numerical solution of the dynamic equations for particulate systems. J. Comput. Phys. 23, 357-375 (1978).
balance
equation
783
Langmuir I., The evaporation of small drops. Phys. Reo. 12, 368 (1918). Lifshitz I. M. and V. V. Slyozov. The kinetics of precipitation from supersaturated solid solutions. J. Phys. Chem. Solids 19, 35-50 (l%l). Peterson T. W., F. Gelbard and J. H. Seinfeld, Dynamics of source-reinforced, coagulating, and condensing aerosols. J. Colloid Interface Sci. 63, 426-445 (1978). Ramabhadran T. E., T. W. Peterson and J. H. Seinfeld, Dynamics of aerosol coagulation and condensation. AIChE JI 22, 840-851 (1976). Ramkrishna D., The status of population balances. Reu. Chem. Engng 3,49-95 (1985). Ruckenstein E. and D. B. Dadyburjor, Mechanisms of aging of supported metal catalysts. J. Catalysis 48,73-80 (1977). Ruckenstein E. and B. Pulvermacher, Kinetics of crystallite sintering during heat treatment of supported metal catalysts. AZChE Jl 19, 356-364 (1973). Scott W. T., Analytic studies of cloud droplet coalescence I. J. Atmos. Sci. 35, 54-65 (1968). Viljoen H. J., D. Eyre and C.-J. W-right, Solving dynamic equations for the collection and evaporation of an aerosol. Can. J. Chem. Engng 68,938-942 (1990).
NUMERICAL TREATMENT OF THE POPULATION BALANCE EQUATION USING A SPLINE-GALERKIN METHOD L. D.
I Department
ERASMUS,’
D. EYRE’ and R. C.
EVERSON’
of Mathematics and Applied Mathematics and ‘Department Potchefstroom University for CHE, Potchefstroom, Republic
(Received
31 August
of Chemical Engineering, of South Africa
1993; final revision received 24 January 1994; received for publication I February 1994)
Abstract-This paper describes a numerical technique for solving the Lifshitz-Slyozov equation of continuity which applies to certain mass transfer proceses. The Lifshitz-Slyozov equation which describes a mechanism involving the transfer of atoms (or undissociated molecules) from smaller particles to larger particles dispersed in a supersaturated medium (Ostwald ripening) is also considered together with collision and subsequent coalescence of particles (ripening plus collection). Special attention is given to the case in which the net growth rate appearing in the Lifshitz-Slyozov equation is a nonsmooth function of the type I(v) 0~v 8, where v is the particle volume and 0
INTRODUCTION
Lifshitz and Slyozov (LS) (1961) equation is used to describe a number of physical and chemical processes involving the decay and growth of particles due to the transfer of atoms or undissociated molecules in a supersaturated environment, known as Ostwald ripening. Examples of these processes include precipitation of solutes from supersaturated solutions and the sintering of supported metal catalysts on an inert substate. Although more detailed microscopic descriptions of the equations of continuity have been derived for Ostwald ripening (Dadyburjor and Ruckenstein, 1977) these equations reduce to the macroscopic LS equation after approximation. Ruckenstein and Dadyburjor (1977) also proposed a combination of Ostwald ripening plus collision and subsequent coalescence that occurs during sintering of supported catalyst particles. Collision and subsequent coalescence will be referred to as collection and the combination as ripening plus collection. Let n(u, t) du be the number of particles with volume between D and u+ dv at time t. The LS equation is a continuity equation: The
775
where Z(u) = in is the growth rate of a particle of volume u, and S( u, C,n) a source term for particles of volume U. In general the source term can be a nonlinear function of the dependent variable n(n, t), as, for example, in the case where particles can be formed by collection. It will be assumed that equation (1) satisfies the usual boundary conditions, namely: n(o, 0) = W(U), n(0, t) = 0.
(2) (3)
The interpretation of these conditions is that there are no particles of zero size, and that the initial distribution of particles wO(u) is given. The total number density MO(t) and volume density M,(t) is given by the moments: MO(t) =
M,(t)=
- n(u, t) du, I 0
(4)
-n(v,r)udu,
(5)
I
0
respectively. If collection does not occur (S = 0) then the total number density of particles is conserved, so that A&(t) =A$, where No denotes the total initial number density _
L.D.
776
ERASMUS
Collection can be included in the LS equation by considering a source term given by the nonlinear expression :
1
”
2I 0 -n(u,
n(u’, t)n(u - u’, t)K(u’,
1)
00 n(u’,
t)K(u’,
u -u’)
u) du’.
I 0
du’
(6)
Here the collection kernel K(u’, u) describes the geometry and dynamics of collection mechanism. The equations (1) and (6) have been extensively investigated by many authors. A detailed discussion of the effect of coupling between collection and the LS equation can be found, for example, in the work of Ramabhadran et al. (1976) and Peterson et al. (1978). Of interest in the present study is a growth law appearing in the LS equation that has the form: I(u) = o/3v@ where O
1,
(7)
where as is the growth constant. A negative value of ua in equation (1) would correspond to decay rather than growth. The exponent p is determined from physical considerations. For example, the case fl= l/3 corresponds to the well known continuum diffusion growth law based on particle diameter. This can be seen by changing the independent variable in the equation of continuity from a volume v to the particle diameter D (for spherical particles). In this new variable the growth law can be written as:
et al.
easily be found and one must resort to numerical methods. Ramkrishna (1985) provides a survey of numerical methods for solving the population balance equation. A collocation method using a finite element basis of cubic splines has been described by Gelbard and Seinfeld (1978), and by Eyre et al. (1988). Viljoen et al. (1990) use the splinecollocation method to obtain numerical solutions of equations that describe decay plus collection. In the case of Langmuir (1918) evaporation, however, Viljoen et al. (1990) reported numerical instability of the spline-collocation method. This result suggests that the collocation approach is not suited to cases in which the growth law has a nonsmooth behaviour; as in the case of equation (7) with O
Before applying the projection method it is first convenient to introduce dimensionless variables and to map the volume variable v onto a finite interval. For this purpose it is convenient to introduce a nonlinear mapping:
(8) The case of 0 C/l < 1 is interesting from a numerical point of view because it is an example of a nonsmooth function (the function has an infinite derivative at v = 0). This behaviour in the growth law leads to a solution of the LS equation that is also a nonsmooth function. In the case of zero collection the exact solution of the LS equation for t>O is:
XW,
[( v
l-(1-&$
This solution exists provided that: v’-fl==(l
-@)a&
u=5
(IO)
When coupled to the particulate balance equation, an exact solution of the LS equation for realistic collection kernels, similar to solution (9), cannot
1+x l-x
( > ’
(11)
where x E [ - 1, 11. (An alternative mapping such as one based on a log or tan function is also possible and will lead to essentially the same results as those presented in this paper.) It is convenient also to introduce a dimensionless time t = A$$, where x is a characteristic rate constant of the process. The dependent variable n(u, t) is now transformed into a function of dimensionless arguments:
141-B) > 1. (9)
METHOD
?i(x,r)=--
5
2
No (1 -x)2
n(u, 0.
The initial function wO(v) is similar transformed so that:
2
iG”C-4
=G
Cl
_x)2
w(v)-
(13)
The parameter c is chosen so that the function fi(x, t) is reasonably distibuted over the interval
Treatment of the population balance equation
777
[ - 1, 11. This parameter remains fixed throughout the calculation. The transformed function fi(x, t) satisfies the LS equation:
The (unknown) solution of the LS equation is approximated by the linear combination
(14)
Substituting this approximation for the solution in the LS equation leads to a residual function:
where the dimensionless growth rate a(x) by:
is given
(22)
r(x, r) = m+l
a(x) = &*5@-‘(1 For the zero-collection ized so that:
+x)B(1 -X)*-8.
(15)
case rZ(x, t) is now normal-
C {Bj(x)fi(r) + [Bi(x)a’(x) + Bl(xMx)lJ(r)) j=O
i=U
I I --I
ii(x,t)dX=l.
(16)
I
ii@‘, T)ti(U, r)y(x’, U) --1
z)
I
l-u l_x
2
( >
dx’
1
--I
A(x’,
r)y(x’,
x) dx’,
(17)
where
r(n’,x)=
K(v’, 0)
08)
x
and u= u(x’, x) is the transformed ence , namely:
volume differ-
2(x-X’)-(l-X)(1-x’) (19)
U=2(x-x’)+(1-x)(l-x’)~
Discretization of equations (14) and (17) is carried out on a space of cubic splines. Cubic B-splines on x E [ - 1, l] are constructed as follows: the interval is partitioned by m nodal points {xi}K1 where - 1 = Xl-=ZX,< - . -
B!“(X)
=
o
Xj-lCXsXi,
(20)
otherwise.
7
The B-splines of order I (degree I- 1) are generated by the stable iterative method of Cox (1972) and deBoor (1972): *c,‘(x)
=
,
(x
--+,m~_,“w - (xi--wY-‘Yx)~
(21)
Xi-X;_/
Note that the index on the cubic B-splines runs from i=2,. . . , m + 3. A convenient notation is Bi (x) = B&(x) to denote the cubic B-spline which is nonzero over the interval (x~_~,_x~+~)so that the index now runs from i=O,. . . ,m+l.
2
j=(J
X
x
Bi(“)f,
(t)Bj(x’lfi(r)Y(x’,
u>
-I
m+l
2
l-_x
da?
(
x
-A@,
1
f-u x
The transformed collection term is
S(X, .,+;
xx-j
m+lm+l
-
>
+c
m+l c
idI
Ix
Bj(X')fi(t)y(X',
~icax~)
j=O
X)
dX’.
W)
--L
Note that the coefficients vi(Z)} are functions of the dimensionless time z. Where no confusion should arise this explicit dependence on time will be suppressed. Clearly the residual function r(x, r) should be small if the spline function (22) is to be a useful approximation to the solution of the original integro-differential equation. Let (q, r) denote the usual inner product: (+, r) =
1 Il(x)e
I -1
b,
w
where O(X) is some appropriately chosen weight factor. Furthermore, let {v*} be a chosen set of (m+2) test functions. The coefficients vi} satisfy the (m+ 1) system of ordinary differential equations: (7&,r)=O;
k=l,.
. . ,(m+l),
(25)
and the boundary condition (3). The approach usually taken when using spline methods for solving population balance equations is a collocation method. In the collocation method the test functions are chosen to be delta functions, ~~=6(x - u,), where uk E [ - 1, 11. An alternative approach, and the one taken in the present paper, is to choose the test functions to be cubic B-splines, t&k= Bk. This is the Galerkin method. Thus the residual is chosen to be zero only in the sense of an integral over the cubic B-splines. The problem now reduces to that of solving the system of ordinary differential equations: Lf(r)=-ivfr)+f(r)Cf(r),
(26)
L. D. ERASMUSet al.
778 7
2.0
6
AQ? 1.8
5
1.6
1.4
-0.6
-0.8
-0.4
0.4
0.2
0
-0.2
0.6
x Fig. 1.
Computed solution of the
LS
equation using the spline-Galerkin method. Since:
for the expansion coefficients f = (f&f,, . . . ,f,,,+t)T. The matrices in equation (26) are given by:
a’(x) Li.j=
'
-I
I
NjJ=
Bi(~)Bj(x)w(x)~9
=
%&_
B
[
(27)
(1-x)2-p (1 +x)*-B - (2 -/I)(1
1
I
-I
+
' I
+x)B(l
Bi(x)B;(n)a(x)w(x)dx
Bi(X)Bj(X)U'(X)CO(X)dX
(28)
-1
and
-x)‘-8
19
(30)
the integrand in the expression for Ni.j contains a weak singularity at x = - 1. This singularity can be removed by chasing an appropriate weight factor w(x) = (1 +x)‘-B. It should be remarked that this Table 1. Exact and computed moments MO(r) and M,(s) for the LS equation
' I
Bk(X)W(X)dX
-I
x I
Bi(U)Bj(X')y(X',
U)
MU
-1
5
0.2 0.4 0.6 0.8 1.0
l.orm 1.000 1.000 1.000 1.000 1.000
:.: 1:6 1.8 2.0
1.ooo l.OOU l.ooo l.oM) l.OOcl
0.0
X
l
--I
Bj(X’)y(X’,
X) dX’.
(29)
The integrals in equations (27-29) are evaluated numerically using a standard Gauss-Legendre quadrature formula.
(exact)
MU
(computed) l.OW l.ooo 1.000 l.uoO l.lXKl i.ooo l.olM l.ocKI l.ooO l.ooo l.tmO
M,
M,
(exact)
(computed)
2.ooo 3.224 4.656 6.265 8.036 9.951 12.002 14.180 16.478 18.887 21.406
2.000 3.067 4.412 5.979 7.733 9.651 11.714 13.912 16.240 18.674 21.219
Treatment of the population balance equation
779
0.8
0.6
-0.2
-0.4
-0.6
-0.8 __ -1
.I)
I
I
I
I
I
I
I
-0.8
I
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
x
Fig. 2. Error in the computed solution, ii(computed) - rl(exact).
2.5
2.0
1.5 G 5 IC 1.0
0.5
0 -“.t.
-‘.U
-“.4
-0.2
0
0.2
0.4
0.6
I Fig. 3. Computed
solution of the collection equation using the spline-Galerkin
method.
0.8
L. D. EWMUS et al.
780
Table2. Exact and computed moments M&) and M,(s) for the collection equation
T
MI (computed)
MO (exact)
MI (exact)
Ml (computed)
0.913 0.839
2.000 2.oiKl 2.000
2.000 2.aXl
0.723 0.675 0.634 0.597 0.564
?E 2:Oo0 2.Olm 2.wm 2.MKJ
0.508 0.535
f:E
where
1 Ai.j=
(33)
dx9
Bi(x)Bj(x) -I
0.0
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6
1.ooo
1.000 0.909 0.833 0.769 0.714 0.667 0.625 0.588 0.556 0.527 0.500
0.7n
f:E 2.000 2.ooo 2.aKl 2.000 2.aiw 2.ooo 2.oLul
approach is used to remove unwanted singularities in the integrand but that the integrand still remains a non-smooth function with possible singularities in the derivative. The system of equations (26) is subject to the initial condition:
gi=
Next we turn to the question of temporal discretization. We introduce the time step Ar. Time is now expressed in terms of an index r = r/At, while the coefficient fi, , denotes an approximation to J(t) , A classical one-step e-method gives the (m + 1) algebraic equations: m+l
m+l
NO)
-0.8
Fig. 4. Computed
-0.6
-0.4
-0.2
#PI+,
-I
L~.jfj.r+(l-8)Ar i=O
(31)
= g,
r
m+l
m+l
Any standard method, such as collocation, could be used to determine the coefficients (f;(O)}. However, in keeping with the Galerkin approach, the coefficients (f,(O)} are obtained by solving the matrix equation:
Bi (X)WO(X)dX*
-1
=c
~fit;-(")Bj(x)=~OSx)-
Il
X
=l i-0
Nk.i-
2 j=O
(35)
fi,r,
Ck,i.jf;..r 1
where k=l,. . . ,m+l. The value of B can be varied from 0= 0 (explicit scheme) to 0 = 1 (fully implicit scheme). The value of 8= l/2 gives the Crank-Nicolson method.
0
0.2
0.4
0.6
0.8
solution of the combined ripening plus collection equation using the spline-Galerkin method.
1.0
Treatment
of the population
'r
781
equation
chosen. A condition that satisfies the boundary conditions (2) and (3) is:
6-
/J ‘1
SF &
4-
I’ I :
1c
3-
:
2-
\
\
2v 4& v exp -, wo( u) = vo> ua (
\ \ \ \ \
where
0.2
0.3
0.4
0.5
0.6
0.7
0.1
I 0.9
Fig. 5. Computed solution after a time t = 2. The solid line shows the double distribution for the combined ripening plus collection equation. The broken line shows the result for ripening alone.
Supplementing equation boundary conditions:
(35) with the discrete
(36)
~fi,r+tsj(_1)=09
j=O
yields an algebraic system for (m +2) expansion coefficients fi, r+l. (Note that only four terms are summed in equation (36) because each subinterval is spanned by only four nonzero basis functions.) In the case of collection this is a nonlinear system which is solved using a standard quasi-Newton solver (ZSPOW from the IMSL package). The set fi_, is used as an initial estimate.
3. NUMERICAL
EXAMPLES
A number of simple test problems will be considered. For this purpose a simple initial condition is
Table 3. Computed moments M,,(r) and M,(r) for the combination of ripening plus collection equation
0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
is
the
average
particle
We begin with a numerical solution of the LS equation (ripening alone). The LS equation with growth law given by (7) has the solution:
x
0.0
M,,(O) = No and u.
(37)
volume.
1O0.1
balance
Ez
0:829 0.748 0.684 0.635 0.589 0.547 0.515 0.492 0.472
2.ooo 3.181 4.412 5.686 7.028 8.462 10.005 11.657 13.423 15.303 17.311
(1 +Bw -8)
n(u,t)=~
4Nol-(1-#3)$$ [ 1
For numerical calculations we have chosen the following parameters; the initial number of particles N,= 1 and the initial average volume of particles u0= 2. The initial condition now becomes: we(u) = u e-“.
(39)
We take /3= * and u,,~ = 1.547. Figure 1 shows the number of particles ii as a function of x for several values of the dimensionless time t. The curves represent a numerical solution of the LS equation using a basis of 19 cubic B-splines. The mapping parameter is 5 = 10 and the time step is At = 0.01. In this example there is a rapid growth of the small particles. In fact the left-hand-side of the particle spectrum is finite only for particle sizes which satisfy inequality (10). The larger particles grow more slowly as indicated by the small change on the right-hand-side of the particle spectrum. It should be remarked that only the positive part of the function is shown in the diagram. In practice a small negative part of the function appears on the lefthand-side of the particle spectrum. This negative value is only a few percent of the peak value and is the result of representing a nonsmooth function by a smooth cubic spline approximation. The negative values are set to zero for calculations at subsequent time intervals. The calculated moment using 19 spline basis functions and the exact moment MI are shown in Table 1. The calculated moment shows a small error (< 10%) which appears initially in the first step of the calculation and then remains constant throughout subsequent time steps. This error is reflected also in the difference between the numerical solution and the exact solution [given by equation (38)] as shown in Fig. 2. Again this error could be attributed to the smoothness properties of the cubic spline approximation. Next we turn to the collection equation. For the purpose of demonstrating the performance of the
L. D. ERASMUS er (II.
782
algorithm we consider the case of a constant kernel, K(u’, v) =$_ In this case the collection equation (Scott, 1968) has the solution:
I-_-
(W
where
MO=-
2No 2+gN,,t
and
M1=NovO.
(41)
Notice that in the case of collection alone the total number of particles is not conserved but that the total volume of particles is conserved. The solution of the collection equation should be well represented by a cubic spline approximation. We choose 5 = 1. Figure 2 shows the computed solution using 19 spline basis functions for several values of the time t. These curves are indistinguishable from the exact solution_ In this example there is little growth of the smaller particles, but rapid growth of the larger particles producing a long tail in the large particle end of the spectrum. It is clear that the effect of collection, while leading to growth in the average particle size, has a much different effect on the particle spectrum that does the growth due to ripening (LS equation). Table 2 shows the computed moments Mu and MI. These results agree extremely well with the exact values and illustrate the accuracy of the spline-Galerkin method for solving the collection equation. The numerical solution of the above two examples has been used to test the numerical procedures. An exact solution of the ripening plus collection equation for the growth law (6) and a constant (nonzero) collection kernel is not known. Thus one must rely on numerical approximations to obtain a solution of the integro-differential equation. A dimensionless measure of the collection/ripening rate is given by the ratio A = ~No~‘-~lu~. The relationship between g= 1 and cLj3 is such that A = 1 for a mapping parameter c= 10. The case A = 1 corresponds to particle growth from equal rates due to ripening and collection. Figure 4 shows the computed particle number spectrum. The rapid growth of the smaller particles can be attributed almost entirely to the effect of ripening. The rapid growth of the larger particles is an enhancement of that seen for the case of collection. An interesting feature of this spectrum is the formation of a double bump. Figure 5 shows the particle spectrum at time t = 2. The solution for both ripening and coltection (solid line) and the result for
ripening alone (broken line) are shown. The first peak in the solid line corresponds to the peak in the ripening spectrum. Table 3 shows the computed moments M,(t) and Ml(r). Although the exact solution for this last example is not known, for a combination of ripening and collection there is an expected decrease in the total number of particles MO and an increase in the total volume of particles Ml with increasing time. It can be remarked that calculations using more realistic kernels that describe sintering, such as those given by Ruckenstein and Pulvermacher (1973), also produce a double bump when combined with the growth law in equation (6). 4. CONCLUSIONS
Numerical solutions of the LS equation have been obtained in a case where the growth rate is a nonsmooth function of the particle volume. Numerical solutions for the combination of ripening plus collection have also been obtained from a numerical method that combines cubic splines and a Galerkin technique to determine the spline expansion coefficients. A novel feature of the ripening plus collection equation is the emergence of a double distribution in the particle size spectrum. NOMENCLATURE
Ai.i = Defined by equation (33) Bi = Cubic B-splines C *, r.i = Defined by equation (29) D = Particle diameter J(z) = Spline expansion coefficient f;. , = Expansion coefficient at discrete time step g, = Defined by equation (34) I( U) = Growth rate K(u’, U) = Collection kernel Li.i = Defined by equation (27) m = Number of nodal points MO = Number density defined by equation (5) MI = Volume density defined by equation (6) n(u, t) du = Average number of particles per unit volume NO = Initial average number of particles per unit volume Ni, j = Defined by equation (28) r = Residual function defined by equation (23) S = Source term defined by equation (6) t = Time u = Defined by equation (19) u, = Collocation points v = Particle volume B = Average particle volume u0= Initial average particle volume w,, = Initial volume distribution x = Dimensionless variable X, = Nodal points Greek letters a = Dimensionless growth rate /3= Exponent in growth term y = Dimensionless collection kernel 5 = Volume parameter
Treatment 8= A = 5= a, = T= x= vi = w =
of the population
Parameter in the e-method on rate Dimensionless ripening/collecti Constant collection kernel Growth parameter Dimensionless time Characteristic rate constant Test function Weight factor in equation (24) REFERENCES
Cox M. G.. The numerical evaluation of B-solines. J. Inst. _ Math. Abpf. 10, 134-149 (1972). Dadyburjor D. B. and E. Ruckenstein. Kinetics of Ostwald ripening. J. Cryst. Growth 40. 279-290 (1977). with B-splines. 1. Approx. deBoor C., On calculating Theory 6, 50-62 (1972). Evre D.. C. J. Wright and G. Reuter. Soline-collocation *with adaptive m&h grading for sob&g the stochastic collection equation. J. Comput. Phys. 78, 288-304 (1988). Gelbard F. and J. H. Seinfeld, Numerical solution of the dynamic equations for particulate systems. J. Comput. Phys. 23, 357-375 (1978).
balance
equation
783
Langmuir I., The evaporation of small drops. Phys. Reo. 12, 368 (1918). Lifshitz I. M. and V. V. Slyozov. The kinetics of precipitation from supersaturated solid solutions. J. Phys. Chem. Solids 19, 35-50 (l%l). Peterson T. W., F. Gelbard and J. H. Seinfeld, Dynamics of source-reinforced, coagulating, and condensing aerosols. J. Colloid Interface Sci. 63, 426-445 (1978). Ramabhadran T. E., T. W. Peterson and J. H. Seinfeld, Dynamics of aerosol coagulation and condensation. AIChE JI 22, 840-851 (1976). Ramkrishna D., The status of population balances. Reu. Chem. Engng 3,49-95 (1985). Ruckenstein E. and D. B. Dadyburjor, Mechanisms of aging of supported metal catalysts. J. Catalysis 48,73-80 (1977). Ruckenstein E. and B. Pulvermacher, Kinetics of crystallite sintering during heat treatment of supported metal catalysts. AZChE Jl 19, 356-364 (1973). Scott W. T., Analytic studies of cloud droplet coalescence I. J. Atmos. Sci. 35, 54-65 (1968). Viljoen H. J., D. Eyre and C.-J. W-right, Solving dynamic equations for the collection and evaporation of an aerosol. Can. J. Chem. Engng 68,938-942 (1990).