solution of population balance equations by mwr - Semantic Scholar
Computers and Chemical Engineering, Vol. I, pp. 23-31.
Pergamon Press, 1977.
Printed in Great Britain
SOLUTION OF POPULATION BALANCE EQUATIONS BY MWR P. N. SINGH* and D. RAMKRIsHNAt Department of Chemical Engineering, Indian Institute of Technology, Kanpur, U.P. 208016,India (Received21March 1975) Abstract-The population balance equation for crystallization in a continuous mixed suspension and mixed product removal crystallizer accounting for the effects of arbitrary crystal breakage (an outstanding problem) has been solved by the method of weighted residuals. The trial functions used were problem-specific polynomials generated by the Gram-Schmidt orthogonalization process with a suitable weight function in the definition of the inner product. The weight function emerged from the analytical solution of multilated versions of the original population balance equation. The methods used also included a Vorobyev’s method of moments which is essentially a variation of the method of weighted residuals. Accurate solutions were obtained with only four problem-specific polynomials thus rating these trial functions above the more standard (and traditional) choices like Laguerre polynomials with which no satisfactory solutions could be obtained even when a large number of functions were employed. Scope-Population balances[l], which are important to the description of the behavior of dispersed phase systems, have found limited applications in chemical engineering because of the difficulty in obtaining solutions to the integral or integro-differential equations which arise therefrom. Attempts in the past have concentrated in estimating the moments of the number density function appearing in the equation through the solution of moment equations[l, 2,9]. This approach has powerful limitations[3]. In this paper, efficient methods of solution have been presented to solve selected population balance equations arising in cystallization. The examplesincludeone that has defied solution in the past. The method of solution presented in this paper is an application of the general method of weighted residuals, of which the method of moment equations referred to earlier is a special case (see, e.g. [3]),with trial functions that were constructed specifically[6].These specific trial functions are obtained by Gram-Schmidt orthogonalization of the set (~‘7using suitably weighted inner products. The nature of the weight functions in the inner product has been discussed before [6]. 1. INTRODUCTION
Population balance equations (PBE) which arise in the description of particulate systems, have been of great interest to many diverse areas of chemical engineering. The general structure of these equations was discussed by Hulburt & Katz[l]. The solution of PBE has however been a chief deterrent to their application in practice. Hulburt & Katz[l] and Hulburt & Akiyama[2] suggested the generation of moment equations from the PBE and to construct the number density function, i.e. the solution by a suitable Laguerre function expansion in terms of coefficients that were known linear combinations of the moments. Where such moment equations are not obtainable either because they are unclosed or because they contain terms involving the solution not expressible as integral moments, they proposed the substitution of the Laguerre function expansion of the solution in such terms. Ramkrishna[3] has shown that this procedure is equivalent to a particular application of the method of weighted residuals (MWR) to the solution of PBE. Subramanian & Ramkrishna[4] further demonstrated the applicability of MWR to the solution of PBE by analyzing the behavior of a biological population. MWR which has found increasing applications of late is the subject matter of an entire book by Finlayson[S]. Essentially, the method consists in expanding the unknown solution in terms of a set of known trial functions, and evaluating the coefficients of expansion in some best way. This is often *Present address: Department of Chemical Engineering, SurathkalEngineeringCollege, Surathkal, Karnataka, India. tTo whom correspondence should be addressed: Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, MN 55455,U.S.A.
accomplished by orthogonalizing the residual (obtained by substituting the trial solution into the equation) with a set of weighting functions that may or may not be identical to the set of trial functions. The success of MWR depends much on the choice of the trial functions. It is desirable that they originate from a complete set although this by itself does not assure rapid convergence to the solution (if convergence is there at all). In this context, Ramkrishna[6] has pointed out the possibility of using problem-specific polynomials (PSP) as trial functions. PSP are generated by Gram-Schmidt orthogonalizing the sequence {x”} with a suitable weight function which presumably has the same general shape and trend as that of the solution sought. In this paper, we consider the applicability of this suggestion to solving a population balance equation for a particle .breakage process. Thus we consider the problem of obtaining the crystal size distribution in a continuous, mixed suspension, mixed product removal crystallizer (CMSMPRC) accounting for the effects of crystal breakage[7]. A method due to Vorobyev[8], called the method of moments, which represents a variation of MWR, appears as yet inadequately tapped in the solution of engineering problems. Since Vorobyev’s method of moments (VMM) consists in generating a set of trial functions from an initial function, it fits neatly into the scheme of this investigation. We present a sketchy description of VMM below leaving the details to Vorobyev[8]. Vorobyeo’s method of moments Consider a bounded linear operator A defined on a real Hilbert space X and the equation w=oAw+f, 23
(I)
P. N.
24
SINGH
and D.
RAMKRISHNA
concentration and free of solids while simultaneously the crystal suspension is withdrawn at an identical rate. Supersaturation by, say cooling of the solution in crystallizer provides for the required nucleation. Crystal breakage invariably occurs in such crystallizers and Randolph[9] provides the population balance equation for a CMSMPRC for the case of binary breakage. The steady (2) state population balance is given by
where w, f E X*; u is a known real parameter. Equation (1) is to be solved for the unknown w given f. MWR consists in expanding the trial solution w, in a finite linear combination of n trial vectors (pa,pl,. . . , +I (normally originating from a set complete in X) given by n-l
W” =
c
QjQj
i=”
and substituting into Eq. (1). The different terms appearing in Eq. (1) after substitution of (2) are collated into a residual term which is to be as near zero as possible by appropriate choice of the coefficients of expansion of the trial solution. This is normally accomplished by orthogonalizing the residual with a set of functions {@i}/-’ generally derived from a complete set in X A particular application of MWR is the Galerkin technique of orthogonalizing the residual from Eqs. (1 and 2) with the vectors +Q,,qpl,. . . , 4pn-,.It is easy to show that this procedure obtains the exact solution of the equation wn =P,AP,w,
+P,f,
(3)
where P, is the orthogonal projection of X onto the linear manifold spanned by po, (pl,. . . , pn-I, denoted by H,,.? The approximation w, is said to converge to the actual solution w of Eq. (1) if lim ((w,- w((= 0, “‘l
(4)
@WI=-;w+(B-D), where D is the “death” function relating to disappearance of crystals of size 1 by breakage assumed to be D(I) = k’w(l)l”,
in which k’ and M are constants; B represents the “birth” function which for binary breakage is found to be B(I)=~‘P(+(;)+D(&)]dr.
qj+l = Ap,
j =O, 1,2 ,..,n-2.
+f.
(6)
Vorobyev[8] establishes the convergence of w, to w when A is a completely continuous operator. Since the operators which are encountered in this paper are completely continuous the convergence of VMM may be assured. Besides Galerkin’s technique, other methods of weighted residuals have also been employed in which the orthogonalization of the residual is performed with a set of vectors {I&,}~“-’ not necessarily the same as the trial set {Qj)o’-‘. 2. PBE FOR A CONTINUOUS
x=E
w(O y=*
k, = k’f3(Ih)“,
(10)
where it has been tacitly assumed that the McCabe “AL law” holds which implies that i is a constant. Equation (7) then transforms to the dimensionless version $+(l+k,x”)y=k,
I”
‘P(e)
(5)
When q0 is taken to be f then Eq. (3) may be replaced by w,, =P,AP,w,
(9)
The function P(E) is a probability density such that P(E) de represents the probability that a crystal of size EI to (e t de)l is formed by breakage of a crystal of size 1. For binary breakage P(r) is symmetric about l = l/2. Following Randolph[9], we switch to the dimensionless variables and constants I
where /].I(represents the vector norm generated by the inner product on X We disregard here the question of when such convergence may ensue; it will be presumed that the method will converge leaving the test to actual computation. The method due to Vorobyev is essentially Galerkin’s technique in which the trial vectors (pO,(p,, . . . , B-, are generated by the recursive relation
(8)
CRYSTALLIZER
A continuous, mixed suspension, mixed product removal crystallizer (CMSMPRC) consists of a well-stirred vessel to which is fed the solution of constant solute
Equation (11) is subject to the initial condition y(0) = 1. Before Eq. (11) can be solved the function P(E) must be specified, which depends on the breakage process. Thus for the attrition model described by Randolph [9] Eq. (11) degenerates into g+y
(12)
so that y(x) = emX.For the binary equal breakage model one has P(E) = 8(~ -i) which produces the steady state PBE 2+(1-t
klx”)y =2”+‘k,xmy(2x).
(13)
The binary uniform breakage model assumes P(E) = 1.0 so that g+(l+k,x”‘)y=k,
I’ K:Tyt)
*Boldface letters such as w, f denote vectors in the abstract Hilbert space X When K is a function space w - w(x), f = j(x). tFor any vector 9 E X the orthogonal projection P,(p is given by
=o,
+
($r y(e)] de.
(14)
This paper undertakes to solve Eqs. (13 and 14) by MWR using problem-specific polynomials.
25
Solutionof populationbalanceequationsby MWR 3.SOLlJTION BY MWR
We consider the solution of Eqs. (13 and 14) by MWR using as trial functions “appropriately” generated problem-specific polynomials. What is appropriate is indeed a matter of intuitive judgement. Thus for the present we may assume that the birth term in Eq. (11) which constitutes the right hand side of the equation is non-existent yielding the differential equation =O,
$+(l+k,x”)o
e-lr+W”“)‘(m+l)l g
=
(16)
We use w(x), clearly positive-valued, as the weight function for generating polynomials by GSOP on the set {x”}. Thus the nth degree polynomial p.(x) is given by
(20)
so that the trial solution w,(x) is given by wn(x)
e-lr+(k,x~+')/lm+l)l~
=
cpi(x)= e-lx+lP,~~+‘)/(m+~)lp,(X)
(15)
where we have replaced the dependent variable y by w. Equation (15), when solved subject to o(O) = 1, yields w(x)
Increased breakage is characterized by higher values of Q. Computations have been made for Q = 0.2 and 0.5, which are the values used by Randolph [9] in the solution of Eq. (13) obtained by the Adams-Moulton-Shell technique. However, no solution had been possible for Eq. (14) in the past. For the application of MWR the trial functions {cpi(x)}a”..’ based on the polynomials pi(x) are expressed as
(21)
aipj(x),
Alternatively, the exponential solution e-” which emerges from the attrition model may also serve as a weight function for generating polynomials. The polynomials so
m
IJI
PO(X)= 1
e-15+(t,5
0
x” -
m+‘v(m+l)l d5
gopr(x) [
{“Pi
e-‘*“k~‘“““‘m+“‘d~
(17)
P.(X) =
n-1 ,,,” _-&(D)
,-oz[“pk(() e-,“+‘k,6”‘+‘)ib”+,), dt]’ e-d;;l+(k,rim*‘)/(m+l)l
Obviously the polynomial coefficients are to be obtained by numerical computation. Table 1 presents the coefficients for these problem-specific polynomials p.(x) represented by P.(X)
ankXk.
= i
k=O
Table 1 Constant Coefficients coeff.
Coeff. ofx
Coeff.
CO&.
Coeff.
0fX2
ofx)
ofx”
0.0
0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.810 0.0 0.0 0.0 0.0 A
generated are none other than the Laguerre polynomials denoted L,(x). Thus the trial functions for this case may be called Laguerre functions which are obtained by multiplying the Laguerre polynomials for e-X. Calculations were made for both Eqs. (13 and 14) using MWR with trial functions given by (20) and the Laguerre functions. The task that remains is of course the determination of the n constants (YO,ol,. . . , (Y”-~,which calls for n equations. The initial condition y(0) = 1 instantly identifies one of these equations. Thus n-l
PSP with weight function (16) Q = 0.5
2.022 - 1.151 1.297 - 1.419 1.524 1.0 Laguerre polynomials 1.0 with weight 1.0 1.0 function 1.0 eC”
1.314 - 3.429 6.312 - 10.161 0.0 - 1.0 -2.0 - 3.0 -4.0
1.244 - 5.218 13.808 0.0 0.0 0.5 1.5 3.0
1.044 - 6.052 0.0 0.0 0.0 -; -f
Evidently the coefficients {a,t}o’ would depend on the parameters kr and m. The value of m is assun+ed to be 4, following Randolph[9]; k,is determined by the extent of breakage and has been characterized by Randolph[9] by a quantity Q, which represents the fraction of particles of the “dominant” size broken to those that are washed out of the crystallizer.* Larson & Randolph[lO] determined the dominant size as given by x = 3. Thus
2 wm =1.
(22)
The residual for either of the Eqs. (13 and 14) is obtained by substituting the trial solution into the equation and algebraically summing all the terms together. We denote the residual function by R(x; ,(Y~, oI,. . . , an-J. Both Eqs. (13 and 14) imply that as x+0 y’(O)+-1, which is a useful condition to impose on the solution. This may be seen to be equivalent to equating to zero the residual at x = 0. Thus R(0; CXO, (YI,.. . , a.-,) = 0.
(23)
It remains to identify (n - 2) more equations in the coefficients ao, aI,. . . , CZ-~.We therefore set
I
oi R(x;ao,al,...,
(Y.-,)&(X)dx = 0
0
(24)
j=O,1,2 ,...,n-3, Q
=
kJ”‘y(3) Y(3)
Or
k, = Q/3”.
(19)
*The dominant crystal size is that which has the maximumweight fraction in the slurry.
where {I,$(x)}~“-’ represent a set of functions not necessarily identical to the trial functions. In the present work, two choices were made for $i(x). The first assumed +(x)= pi(x) while the second had Jll(x) = pi(x). The
26
P. N. SINGH and D. RAMKRISHNA
results of calculation are presented in the next section. In order to apply the Vorobyev method of moments we transform Eqs. (13 and 14) into their respective integral forms. Thus for Eq. (13) we have
while for Eq. (14) one obtains
xexp
[
5 -x +k’(S”nT’,IXm+‘)]~m-2Y(~). (26)
In Eqs. (25 and 26) the inhomogeneous term f(x) is given by f(x) = e-Ir+(~,xm+‘)/(m+l)l, (27) which is the same as (16). The integral representations (25 and 26), free from the unbounded operators present in their counterparts (13 and 14), may now be written in operator form as y = f + k,Ay,
(28)
where A is a linear operator defined on _ZJO,m) and given by c5”
formula. To obtain the coefficients of expansion of the trial solution, (~0, al, . . . anml,Eqs. (24), which can be condensed into a single matrix equation, were solved by a standard matrix-inversion subroutine. The successive application of Eq. (5) using the integral operators in Eq. (29) involves the evaluation of iterated integrals. These integrals were evaluated by using Simpson’s rule and checked to four decimal place accuracy with the evaluation of the same integrals by Romberg’s method[l3]. The computation of & from &, requires knowledge of the latter at values of its arguments different from those for which &I has been evaluated previously. This was done by using a cubic spline interpolation procedure [ 141. The computation times presented herein include the compilation time and the time taken for the solution with the number of trial functions varying from three to the maximum until convergence is obtained. Approximate solutions obtained by MWR using PSP as trial functions took ca. 4 min for the uniform breakage model and ca. 3 min for the binary equal breakage model. Vorobyev’s method required shorter times. Thus only 2 min were sufficient to obtain very accurate solutions with n varying from 2 to 6. 5. RESULTS AND DISCUSSION Whether the approximate solution obtained is indeed an acceptable solution is not entirely a trivial matter. It is common to assume that if the solution does not materially change with increasing number of trial functions, the solution on hand is a satisfactory solution. In this work, this criterion was supplemented by an additional check
exp [(f-x)+~(~~“-xm”]Y(5)dS
[forEq.(25)]
Ay =
(2%
2~d~~d~Sexp[5-x+k’(5~‘~~m”)]$m-2Y(~)
The operator A is completely continuous which grants the convergence of VMM [8]. The trial functions are now generated as suggested by Eq. (5). When {+J~}~“-’ have been obtained the method is identical to Galerkin’s technique. The orthogonalization of the residual with {pj}o”-’is then to be performed using the regular inner product of &[O, m). Notice that the inhomogenous term f in Eq. (28) is the same as that used as the weight function for PSP in the previous methods.
[forEq.(26)].
which consists in partitioning the equation into a right hand side and a left hand side (such that either one on substitution of the approximate solution is small in relation to unity) and demanding that the ratio of the two be sufficiently close to unity. Thus the Eqs. (13 and 14) as also (25 and 26), written as they are, were found to be adequate partitions. The results for the binary equal breakage and the binary uniform breakage models are presented and discussed separately.
4. COMPUTATIONAL PROCEDURES
The computations, which were performed on an IBM-7044 computer, essentially consisted of (i) determination of the coefficients of the problem-specific polynomials, (ii) orthogonalization of the residual with weighting functions, and (iii) determination of the coefficients of expansion of the trial solution in terms of the trial functions. These pertain to the method of weighted residuals, while calculations using Vorobyev’s method of moments involved successive application of Eq. (5) using the integral operators defined by (29). For the method of weighted residuals, the polynomial coefficients are obtained by performing the integrations indicated in Eqs. (17). Also, the orthogonalization procedure in Eqs. (24), amounts to a set of integrations over the semi-infinite interval. These integrals were evaluated by a 32 point Gauss-Laguerre quadrature
Binary equal breakage model
Either of Eqs. (13 and 25) represents the binary equal breakage model. The solutions by MWR using problem-specific polynomials (17) for Q = 0.5 have been presented in Figs. 1 and 2. It was found that when the orthogonalizing function +j(x) is the same as pi(x), the solution settles down to within five decimal places with only three problemspecific polynomials, i.e. n = 3. However, n = 5 displays a strange aberration (Fig. 1). Check of the numerical results for possible errors confirmed the result shown. With increasing II the solution again is restored to that at n = 3; further Table 2 presents the ratio of left hand side to right hand side for this solution. It is clear from Table 2 that this ratio is indeed close to unity. For this reason the solution obtained is referred to as the “accurate” solution. Figures
21
Solutionofpopulationbalanceequationsby MWR
Accurate
solution
solution
-Accurate
Bmory
A
equal breakage
Fig. 1. Effect of crystal breakage on population density in a CMSMPRC.
Fig. 2. Effect of crystal breakage on population density in a CMSMPRC,comparisonof weighting functions.
1 and 2 represent the accurate solution by a continuous curve. Discretized representation has been used to denote solutions in general for visual convenience. Figure 2 shows the results of computation using PSP and $i = & The resulting approximate solution scatters about the accurate solution without settling down. Thus it appears that the choice of I+$= pi is more adequate. Figures 3 and 4 show the calculations with Laguerre polynomials for two different choices of I+%. It is seen that as many as fourteen Laguerre polynomials are required to obtain a reasonably acceptable solution (Fig. 3). In Fig. 4 where $i = +, the solution is not obtained even with fourteen or more Laguerre functions. Figure 5 shows the results of computation using VMM. It is found that excellent solutions are obtained for n 2 2. Similar calculations have been made by Singh[lO] for
Q = 0.2, and the results possess features identical to those at Q = 0.5. Figure 6 summarizes the calculations for Q = 0.5 and 0.2 obtained by MWR using PSP and VMM. Calculations by Randolph have also been shown in the figure. The figure quantitatively brings out the obvious effect of the value of Q on the solution. Binary
uniform
breakage
model
The model is represented by either of Eqs. (14 and 26). No solution has been obtained in the literature for this situation. Also the equation is of general interest to particulate systems in which the particles grow and undergo splitting. The calculations made are identical to those for the binary equal breakage model in the sense the same set of
Table 2. Comparison of accurate solutions to the PBE for binary equal breakage model Q = 0.5 MWR(Laguerre polynomials) ” = 14 VMMn=2
MWR (PSP) n = 6 x
Y
Ratio
0.0 0.2 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0
1@300 0.8178 0.6062 0.3772 0.2413 0.1528 09036 x 10-l 0.4667x 10-l 0.1950x 10-l 06468 x lo-* 0.1399x lo-’ 0.1748x lo-’ 0.1008x lo-’ 0.1864x lo-”
1.000 0.990 1.020 1.050 1.010 10Nl l+tO 1.050 I.001 1.055 1.040 1.035 1.045 1.050
Y
Ratio
1NKtO 1GlO 0.8163 1.100 0.6085 1.080 0.3798 1.090 0.2438 1.030 1.1506 1.120 0.8823x IO-’ 1.200 0.4570x 10-l 0.980 0.1995x 10-l 1.200 0.8098x lo-* 0.890 0.1802x 1O-2 1.500 0.2200x lo-’ 1.450 0.1135x 1o-4 1.380 0.2836x lo-” 1.530
Y
Ratio
l.OtlOO 0.8187 0.6070 0.3750 0.2419 0.1552 0.8961x 10-l 0.4550x 10-l 0.1950x 10-l 0.6383x lo-* 0.1404x 1om2 0.1755x lo-’ 0.1010x lo-’ 0.2074x 10m6
1.000 1.000 1.000 0998 0999 l+nxl 1.001 I@31 1.001 1.001 1.001
1ml 1.010 1.010
P. N.
28
SINCH
and D.
RAMKRISHNA
IO Qz05 -Accurate
solution
n=3
*
Binary
equal
breakage
-Accurate
Blnory
0
J
IO
solution
equal
breakage
2 3 DlmenslonlQss
I
4 5 SIZP. x
6
7
7.5
Fig. 5. Effect of crystal breakage on population density in a CMSMPRCVor Obyev’s method of moments.
Fig. 3. Effect of crystal breakage on population density in a CMSMPRC,comparisonof weight functions.
0105 0.02
CA=05
x MWR . -VMM -Accurate .
.
solution
n:3
-VMM Binary
0 n:4
MWR
A Numerical solution (Randolph)
equal
breakage
x n:5 d n=6 0
n =I4
A
O,(X)=P-XL,w
4 -I
P,(x)=e-xL,(x)
0
Binary
equal
I
2
breakage 3
4
5
6
7
75
Fig. 4. Effect of crystal breakage on population density in a
Fig. 6. Comparison of results obtained by VMM, MWR and numerical integration (Randolph).
CMSMPRC,comparison of weighting functions.
trial functions and orthogonalizing functions have been used. Figures 7 and 8 show the resuits of MWR using PSP for Q = 0.5. The solution again settles down to a final value within n = 4. Table 3 shows the ratio of left hand side to right hand side which is indeed close to unity. Again the accurate solution is represented by a continuous curve. On the other hand, the use of Laguerre polynomials produced no satisfactory results even when fourteen
functions are used as is established by Figs. 9 and 10. Figure 11 represents the results for VMM which again show that very accurate solutions are obtained for n 12. Figure 12 presents the calculations for Q = 0.2 and 0.5 obtained by MWR and VMM. 6. CONCLUSIONS The results in this paper have established that appropriately generated problem-specific polynomials are distinctly superior to the standard polynomials which
29
Solution of population balance equations by MWR
Q=O5 -Accurate .
solution
n.3
en:4 x n:5 A n-6 0 ” :14
E
0
mi(x)=p,(x)
,B
@i(X)=ewX
& n.6
i
6(‘*02klx5)
Vi (X):P,(X) Bi;
L,(X)
$Ji(X)=Ll(X)
.
&nary
uniform
I 1
I 2
breakage
t un,ttcrm
, \
;kag;
,
,
0
,
IO4 0
I
I 4
I
3
Dlmenslonless 0
I
2
5
4
3
Dlmonslonless
6
7
75
51ze.x
SIZ’2,
I
I
I
5
6
7
75
x
Fia. 9. Effect of crystal breakage on population density in a CMSMPRC,comparison of weight functions.
Fig. 7. Effect of crystal breakage on population density in a CMSMPRC.
.
o n.6 @,(x)_p,(x) +i(x)=pi(x) B~nory
,-(x*02 ,-(x+02 “nltorm
klX5) klX5) braakag@
l
I
.
. %
l \
104.
0
I
I
I
1
2
3
Dlmonslonloss
\. 75
Fig. 8. Effect of crystal breakage on population density in a CMSMPRC,comparison of weightingfunctions.
have been used in the past. Singh[lO] has also used PSP to solve nonlinear population balance equations arising in the description of agglomerative particulate systems. That the appropriateness of the weight function depends on its similarity in shape and trend to the solution in question[6] is borne out by Fig. 13 which shows the weight functions e-’ and (16) plotted alongside the solutions for the two models.
II 4
5 SIZQ.
6
75
x
Fig. 10. Effect of crystal breakage on population density in a CMSMPRC,comparison of weightingfunctions. Also Vorobyev’s method of moments has proved particularly effective in the solution of the problems
considered herein. The extension of VMM to nonlinear problems considered herein. The extension of VMM to nonlinear problems however does not appear to have been attempted. In recent times, the orthogonal collocation procedure[ll, 121 has proved a valuable tool in the solution of mathematical equations of physical phenomena. In a sense, the effort here has incorporated the collocation procedure with x = 0 as a collocation point [see Eq. (23)]. Singh[lO] has used the collocation technique in its entirety in combination with problem-
30
P. N.
SINGH
and D.
RAMKRISHNA
Table 3. Comparison of accurate solutions to the PBE for binary uniform breakage
model
Q =O.S MWR(Laguerrepolynomials) n = 14
MWR(PSP) II = 6 x 0.0
0.2 0.5 1.0 1.5
2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.4 6.0
Ratio
Y
14000 0.8191 0.6085 0.3743 0.2332 0.1461 0.8972 x 0.5141 x 0.2545 x 0.9701 x 0.2529 x 0.3700 x 0.4565 x 0.5746 x
IO-’ 10-l 10-l lo-* 10m2 lo-’ 10m4 lo-”
l+OO 0997 0999 1.005 1.001 0997 0.997 1.001 1.005 1.003 0998 0997 0.995 0.990
Ratio
Y
1@000 0.8187 06076 0.3751 0.2335 0.1463 0.8730 x 0.4925 x 0.2623 x 0.1122 x 0.4083 x 0.7938 x 0.1829 x -
10-l 10-l IO-’ 10-l lo-’
IO-’ lo-’
0.4504x lo-’
n=4
VMM
Ratio
Y
1.000
1.0000
I@0
0.950 0,965
0.8191 06086 0.3743 0.2332
0997
1.100 1.004 1.060 0.850 0.800 1.200 1.884 3.199 15.930 10.930 7.500
0
0.1461 0.8979 0.5143 0.2542 0.9697 0.2494 0.3611 0.4532 0.5689
1
2
l@O
I@05
x x x x x
10-l 10-l 10-l lo-’ IO-’ x lo-’ x IO-’ x lo-”
3
1.OOl l+NXl 1.000 l.OlM 1,002 1.003 1.011 1.022 1,032 1.035
4
5
Dimensionless
see,
6
Fig. 12. Comparison of results obtained by VMM and
Fig.
11.Effect of crystal breakage on population density in a CMSMPRC, Vorobyev’smethodof moments.
specific polynomials successfully in the solution of nonlinear population balance equations arising in the theory of agglomerative processes. CONCLUSIONS
AND SIGNIFICANCE
Solutions were obtained for integro-differential equations of population balance models of a continuous crystallizer. These are presented in Tables l-3 and in Figs. 1, 5-8 and 11. Accurate solutions were obtained by using
7 -;15
x MWR.
the specially constructed trial functions referred to in the preceding section. only a small number of which were necessary. The weight function for the orthogonalization procedure, which generates the trial functions, was chosen to be the solution of a simplified version of the population balance model. Figures 3,4,9 and 10 show that standard Laguerre function expansion of the solution suggested by earlier workers [l, 21 has serious limitations in striking contrast with the methods presented herein. Figures 5, 11 and 12 expose an additional method called Vorobyev’s method of moments[8] as a very attractive prospect for the solution of linear population balance equations of dispersed phase systems in which particle breakage is predominant.
31
Solution of population balance equations by MWR
0=05 -~-c-x ___ ,-(x+02
Q R w(l) w, x y
kl x5)
constant defined by Eq. (19) residual number density of crystals nth approximate of w dimensionless crystal size dimensionless number density
Greek symbols a. coefficient of nth approximating function in trial solution
cp. tin 0 6
nth approximating function in trial solution nth function used to orthogonalize residual average holding time Dirac delta function constant parameter in Eq. (1) E fraction of daughter crystal size referred to parent size o weight function in GSOP.
Abbreviations used CMSMPR continuous, mixed suspension, mixed product removal GSOP Gram-Schmidt orthogonalization process MWR method of weighted residuals PBE population balance equation VMM Vorobyev’s method of moments. REFERENCES
0
1
2
3
4
Dlmonslonloss
5 sze,
6
7
15
x
Fig. 13. Shapes and trends of weight functions and accurate solution.
NOMENCLATURE
A a.* B(I) D(I) X k’ k.l L L”
1 P(e) P. P”
bounded linear operator coefficient of xk in the nth degree polynomial number birth rate of crystals of size 1 number death rate of crystals of size I Hilbert space constant in Eq. (8) dimensionless constant in Eq. (10) linear growth rate of crystal nth degree Laguerre polynomials length of crystal orthogonal projection operator orthogonal projection operator nth degree problem-specific polynomial
1. H. M. Hulburt and S. L. Katz, Chem. Engng Sci. 19, 555 (1964). 2. H. M. Hulburt and T. Akiyama, lnd. Engng Chem. Fund/s. 8, 319 (1%9). 3. D. Ramkrishna, Chem. Engng Sci. 26, 1134(1971). 4. G. Subramanian and D. Ramkrishna, Math. Biosci. 10, 1 (1971). 5. B. A. Finlayson, The Method of Weighted Residuals and Variational Principles. Academic Press, New York (1972). 6. D. Ramkrishna, Chem. Engng Sci. 28, 1362(1973). 7. A. D. Randolph and M. A. Larson, Theory of Particulare Processes. Academic Press, New York (1971). 8. Yu. V. Vorobyev, Method of Moments in Applied Mathematics (Translated from the Russian by B. Seckler). Gordon & Breach, New York (1965). 9. A. D. Randolph, lnd. Engng Chem. Fundls. 8, 53 (1%9). 10. P. N. Singh, Ph.D. Dissertation, Indian Institute of Technology, Kanpur (1974). 11. B. A. Finlayson, Chem. Engng Sci. 26, 1081(1971). 12. J. V. Villadsen and W. E. Stewart, Chem. Engng Sci. 22, 1483 (1967). 13. B. Carnahan, H. A. Luther and J. 0. Wilkes, Applied Numerical Methods. Wiley, New York (1969). 14. .I. H. Ahlberg, E. N. Nilson and .I. L. Walsh, The Theory of Splines and Their Applications. Academic Press, New York (1967).
Pergamon Press, 1977.
Printed in Great Britain
SOLUTION OF POPULATION BALANCE EQUATIONS BY MWR P. N. SINGH* and D. RAMKRIsHNAt Department of Chemical Engineering, Indian Institute of Technology, Kanpur, U.P. 208016,India (Received21March 1975) Abstract-The population balance equation for crystallization in a continuous mixed suspension and mixed product removal crystallizer accounting for the effects of arbitrary crystal breakage (an outstanding problem) has been solved by the method of weighted residuals. The trial functions used were problem-specific polynomials generated by the Gram-Schmidt orthogonalization process with a suitable weight function in the definition of the inner product. The weight function emerged from the analytical solution of multilated versions of the original population balance equation. The methods used also included a Vorobyev’s method of moments which is essentially a variation of the method of weighted residuals. Accurate solutions were obtained with only four problem-specific polynomials thus rating these trial functions above the more standard (and traditional) choices like Laguerre polynomials with which no satisfactory solutions could be obtained even when a large number of functions were employed. Scope-Population balances[l], which are important to the description of the behavior of dispersed phase systems, have found limited applications in chemical engineering because of the difficulty in obtaining solutions to the integral or integro-differential equations which arise therefrom. Attempts in the past have concentrated in estimating the moments of the number density function appearing in the equation through the solution of moment equations[l, 2,9]. This approach has powerful limitations[3]. In this paper, efficient methods of solution have been presented to solve selected population balance equations arising in cystallization. The examplesincludeone that has defied solution in the past. The method of solution presented in this paper is an application of the general method of weighted residuals, of which the method of moment equations referred to earlier is a special case (see, e.g. [3]),with trial functions that were constructed specifically[6].These specific trial functions are obtained by Gram-Schmidt orthogonalization of the set (~‘7using suitably weighted inner products. The nature of the weight functions in the inner product has been discussed before [6]. 1. INTRODUCTION
Population balance equations (PBE) which arise in the description of particulate systems, have been of great interest to many diverse areas of chemical engineering. The general structure of these equations was discussed by Hulburt & Katz[l]. The solution of PBE has however been a chief deterrent to their application in practice. Hulburt & Katz[l] and Hulburt & Akiyama[2] suggested the generation of moment equations from the PBE and to construct the number density function, i.e. the solution by a suitable Laguerre function expansion in terms of coefficients that were known linear combinations of the moments. Where such moment equations are not obtainable either because they are unclosed or because they contain terms involving the solution not expressible as integral moments, they proposed the substitution of the Laguerre function expansion of the solution in such terms. Ramkrishna[3] has shown that this procedure is equivalent to a particular application of the method of weighted residuals (MWR) to the solution of PBE. Subramanian & Ramkrishna[4] further demonstrated the applicability of MWR to the solution of PBE by analyzing the behavior of a biological population. MWR which has found increasing applications of late is the subject matter of an entire book by Finlayson[S]. Essentially, the method consists in expanding the unknown solution in terms of a set of known trial functions, and evaluating the coefficients of expansion in some best way. This is often *Present address: Department of Chemical Engineering, SurathkalEngineeringCollege, Surathkal, Karnataka, India. tTo whom correspondence should be addressed: Department of Chemical Engineering and Materials Science, University of Minnesota, Minneapolis, MN 55455,U.S.A.
accomplished by orthogonalizing the residual (obtained by substituting the trial solution into the equation) with a set of weighting functions that may or may not be identical to the set of trial functions. The success of MWR depends much on the choice of the trial functions. It is desirable that they originate from a complete set although this by itself does not assure rapid convergence to the solution (if convergence is there at all). In this context, Ramkrishna[6] has pointed out the possibility of using problem-specific polynomials (PSP) as trial functions. PSP are generated by Gram-Schmidt orthogonalizing the sequence {x”} with a suitable weight function which presumably has the same general shape and trend as that of the solution sought. In this paper, we consider the applicability of this suggestion to solving a population balance equation for a particle .breakage process. Thus we consider the problem of obtaining the crystal size distribution in a continuous, mixed suspension, mixed product removal crystallizer (CMSMPRC) accounting for the effects of crystal breakage[7]. A method due to Vorobyev[8], called the method of moments, which represents a variation of MWR, appears as yet inadequately tapped in the solution of engineering problems. Since Vorobyev’s method of moments (VMM) consists in generating a set of trial functions from an initial function, it fits neatly into the scheme of this investigation. We present a sketchy description of VMM below leaving the details to Vorobyev[8]. Vorobyeo’s method of moments Consider a bounded linear operator A defined on a real Hilbert space X and the equation w=oAw+f, 23
(I)
P. N.
24
SINGH
and D.
RAMKRISHNA
concentration and free of solids while simultaneously the crystal suspension is withdrawn at an identical rate. Supersaturation by, say cooling of the solution in crystallizer provides for the required nucleation. Crystal breakage invariably occurs in such crystallizers and Randolph[9] provides the population balance equation for a CMSMPRC for the case of binary breakage. The steady (2) state population balance is given by
where w, f E X*; u is a known real parameter. Equation (1) is to be solved for the unknown w given f. MWR consists in expanding the trial solution w, in a finite linear combination of n trial vectors (pa,pl,. . . , +I (normally originating from a set complete in X) given by n-l
W” =
c
QjQj
i=”
and substituting into Eq. (1). The different terms appearing in Eq. (1) after substitution of (2) are collated into a residual term which is to be as near zero as possible by appropriate choice of the coefficients of expansion of the trial solution. This is normally accomplished by orthogonalizing the residual with a set of functions {@i}/-’ generally derived from a complete set in X A particular application of MWR is the Galerkin technique of orthogonalizing the residual from Eqs. (1 and 2) with the vectors +Q,,qpl,. . . , 4pn-,.It is easy to show that this procedure obtains the exact solution of the equation wn =P,AP,w,
+P,f,
(3)
where P, is the orthogonal projection of X onto the linear manifold spanned by po, (pl,. . . , pn-I, denoted by H,,.? The approximation w, is said to converge to the actual solution w of Eq. (1) if lim ((w,- w((= 0, “‘l
(4)
@WI=-;w+(B-D), where D is the “death” function relating to disappearance of crystals of size 1 by breakage assumed to be D(I) = k’w(l)l”,
in which k’ and M are constants; B represents the “birth” function which for binary breakage is found to be B(I)=~‘P(+(;)+D(&)]dr.
qj+l = Ap,
j =O, 1,2 ,..,n-2.
+f.
(6)
Vorobyev[8] establishes the convergence of w, to w when A is a completely continuous operator. Since the operators which are encountered in this paper are completely continuous the convergence of VMM may be assured. Besides Galerkin’s technique, other methods of weighted residuals have also been employed in which the orthogonalization of the residual is performed with a set of vectors {I&,}~“-’ not necessarily the same as the trial set {Qj)o’-‘. 2. PBE FOR A CONTINUOUS
x=E
w(O y=*
k, = k’f3(Ih)“,
(10)
where it has been tacitly assumed that the McCabe “AL law” holds which implies that i is a constant. Equation (7) then transforms to the dimensionless version $+(l+k,x”)y=k,
I”
‘P(e)
(5)
When q0 is taken to be f then Eq. (3) may be replaced by w,, =P,AP,w,
(9)
The function P(E) is a probability density such that P(E) de represents the probability that a crystal of size EI to (e t de)l is formed by breakage of a crystal of size 1. For binary breakage P(r) is symmetric about l = l/2. Following Randolph[9], we switch to the dimensionless variables and constants I
where /].I(represents the vector norm generated by the inner product on X We disregard here the question of when such convergence may ensue; it will be presumed that the method will converge leaving the test to actual computation. The method due to Vorobyev is essentially Galerkin’s technique in which the trial vectors (pO,(p,, . . . , B-, are generated by the recursive relation
(8)
CRYSTALLIZER
A continuous, mixed suspension, mixed product removal crystallizer (CMSMPRC) consists of a well-stirred vessel to which is fed the solution of constant solute
Equation (11) is subject to the initial condition y(0) = 1. Before Eq. (11) can be solved the function P(E) must be specified, which depends on the breakage process. Thus for the attrition model described by Randolph [9] Eq. (11) degenerates into g+y
(12)
so that y(x) = emX.For the binary equal breakage model one has P(E) = 8(~ -i) which produces the steady state PBE 2+(1-t
klx”)y =2”+‘k,xmy(2x).
(13)
The binary uniform breakage model assumes P(E) = 1.0 so that g+(l+k,x”‘)y=k,
I’ K:Tyt)
*Boldface letters such as w, f denote vectors in the abstract Hilbert space X When K is a function space w - w(x), f = j(x). tFor any vector 9 E X the orthogonal projection P,(p is given by
=o,
+
($r y(e)] de.
(14)
This paper undertakes to solve Eqs. (13 and 14) by MWR using problem-specific polynomials.
25
Solutionof populationbalanceequationsby MWR 3.SOLlJTION BY MWR
We consider the solution of Eqs. (13 and 14) by MWR using as trial functions “appropriately” generated problem-specific polynomials. What is appropriate is indeed a matter of intuitive judgement. Thus for the present we may assume that the birth term in Eq. (11) which constitutes the right hand side of the equation is non-existent yielding the differential equation =O,
$+(l+k,x”)o
e-lr+W”“)‘(m+l)l g
=
(16)
We use w(x), clearly positive-valued, as the weight function for generating polynomials by GSOP on the set {x”}. Thus the nth degree polynomial p.(x) is given by
(20)
so that the trial solution w,(x) is given by wn(x)
e-lr+(k,x~+')/lm+l)l~
=
cpi(x)= e-lx+lP,~~+‘)/(m+~)lp,(X)
(15)
where we have replaced the dependent variable y by w. Equation (15), when solved subject to o(O) = 1, yields w(x)
Increased breakage is characterized by higher values of Q. Computations have been made for Q = 0.2 and 0.5, which are the values used by Randolph [9] in the solution of Eq. (13) obtained by the Adams-Moulton-Shell technique. However, no solution had been possible for Eq. (14) in the past. For the application of MWR the trial functions {cpi(x)}a”..’ based on the polynomials pi(x) are expressed as
(21)
aipj(x),
Alternatively, the exponential solution e-” which emerges from the attrition model may also serve as a weight function for generating polynomials. The polynomials so
m
IJI
PO(X)= 1
e-15+(t,5
0
x” -
m+‘v(m+l)l d5
gopr(x) [
{“Pi
e-‘*“k~‘“““‘m+“‘d~
(17)
P.(X) =
n-1 ,,,” _-&(D)
,-oz[“pk(() e-,“+‘k,6”‘+‘)ib”+,), dt]’ e-d;;l+(k,rim*‘)/(m+l)l
Obviously the polynomial coefficients are to be obtained by numerical computation. Table 1 presents the coefficients for these problem-specific polynomials p.(x) represented by P.(X)
ankXk.
= i
k=O
Table 1 Constant Coefficients coeff.
Coeff. ofx
Coeff.
CO&.
Coeff.
0fX2
ofx)
ofx”
0.0
0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.810 0.0 0.0 0.0 0.0 A
generated are none other than the Laguerre polynomials denoted L,(x). Thus the trial functions for this case may be called Laguerre functions which are obtained by multiplying the Laguerre polynomials for e-X. Calculations were made for both Eqs. (13 and 14) using MWR with trial functions given by (20) and the Laguerre functions. The task that remains is of course the determination of the n constants (YO,ol,. . . , (Y”-~,which calls for n equations. The initial condition y(0) = 1 instantly identifies one of these equations. Thus n-l
PSP with weight function (16) Q = 0.5
2.022 - 1.151 1.297 - 1.419 1.524 1.0 Laguerre polynomials 1.0 with weight 1.0 1.0 function 1.0 eC”
1.314 - 3.429 6.312 - 10.161 0.0 - 1.0 -2.0 - 3.0 -4.0
1.244 - 5.218 13.808 0.0 0.0 0.5 1.5 3.0
1.044 - 6.052 0.0 0.0 0.0 -; -f
Evidently the coefficients {a,t}o’ would depend on the parameters kr and m. The value of m is assun+ed to be 4, following Randolph[9]; k,is determined by the extent of breakage and has been characterized by Randolph[9] by a quantity Q, which represents the fraction of particles of the “dominant” size broken to those that are washed out of the crystallizer.* Larson & Randolph[lO] determined the dominant size as given by x = 3. Thus
2 wm =1.
(22)
The residual for either of the Eqs. (13 and 14) is obtained by substituting the trial solution into the equation and algebraically summing all the terms together. We denote the residual function by R(x; ,(Y~, oI,. . . , an-J. Both Eqs. (13 and 14) imply that as x+0 y’(O)+-1, which is a useful condition to impose on the solution. This may be seen to be equivalent to equating to zero the residual at x = 0. Thus R(0; CXO, (YI,.. . , a.-,) = 0.
(23)
It remains to identify (n - 2) more equations in the coefficients ao, aI,. . . , CZ-~.We therefore set
I
oi R(x;ao,al,...,
(Y.-,)&(X)dx = 0
0
(24)
j=O,1,2 ,...,n-3, Q
=
kJ”‘y(3) Y(3)
Or
k, = Q/3”.
(19)
*The dominant crystal size is that which has the maximumweight fraction in the slurry.
where {I,$(x)}~“-’ represent a set of functions not necessarily identical to the trial functions. In the present work, two choices were made for $i(x). The first assumed +(x)= pi(x) while the second had Jll(x) = pi(x). The
26
P. N. SINGH and D. RAMKRISHNA
results of calculation are presented in the next section. In order to apply the Vorobyev method of moments we transform Eqs. (13 and 14) into their respective integral forms. Thus for Eq. (13) we have
while for Eq. (14) one obtains
xexp
[
5 -x +k’(S”nT’,IXm+‘)]~m-2Y(~). (26)
In Eqs. (25 and 26) the inhomogeneous term f(x) is given by f(x) = e-Ir+(~,xm+‘)/(m+l)l, (27) which is the same as (16). The integral representations (25 and 26), free from the unbounded operators present in their counterparts (13 and 14), may now be written in operator form as y = f + k,Ay,
(28)
where A is a linear operator defined on _ZJO,m) and given by c5”
formula. To obtain the coefficients of expansion of the trial solution, (~0, al, . . . anml,Eqs. (24), which can be condensed into a single matrix equation, were solved by a standard matrix-inversion subroutine. The successive application of Eq. (5) using the integral operators in Eq. (29) involves the evaluation of iterated integrals. These integrals were evaluated by using Simpson’s rule and checked to four decimal place accuracy with the evaluation of the same integrals by Romberg’s method[l3]. The computation of & from &, requires knowledge of the latter at values of its arguments different from those for which &I has been evaluated previously. This was done by using a cubic spline interpolation procedure [ 141. The computation times presented herein include the compilation time and the time taken for the solution with the number of trial functions varying from three to the maximum until convergence is obtained. Approximate solutions obtained by MWR using PSP as trial functions took ca. 4 min for the uniform breakage model and ca. 3 min for the binary equal breakage model. Vorobyev’s method required shorter times. Thus only 2 min were sufficient to obtain very accurate solutions with n varying from 2 to 6. 5. RESULTS AND DISCUSSION Whether the approximate solution obtained is indeed an acceptable solution is not entirely a trivial matter. It is common to assume that if the solution does not materially change with increasing number of trial functions, the solution on hand is a satisfactory solution. In this work, this criterion was supplemented by an additional check
exp [(f-x)+~(~~“-xm”]Y(5)dS
[forEq.(25)]
Ay =
(2%
2~d~~d~Sexp[5-x+k’(5~‘~~m”)]$m-2Y(~)
The operator A is completely continuous which grants the convergence of VMM [8]. The trial functions are now generated as suggested by Eq. (5). When {+J~}~“-’ have been obtained the method is identical to Galerkin’s technique. The orthogonalization of the residual with {pj}o”-’is then to be performed using the regular inner product of &[O, m). Notice that the inhomogenous term f in Eq. (28) is the same as that used as the weight function for PSP in the previous methods.
[forEq.(26)].
which consists in partitioning the equation into a right hand side and a left hand side (such that either one on substitution of the approximate solution is small in relation to unity) and demanding that the ratio of the two be sufficiently close to unity. Thus the Eqs. (13 and 14) as also (25 and 26), written as they are, were found to be adequate partitions. The results for the binary equal breakage and the binary uniform breakage models are presented and discussed separately.
4. COMPUTATIONAL PROCEDURES
The computations, which were performed on an IBM-7044 computer, essentially consisted of (i) determination of the coefficients of the problem-specific polynomials, (ii) orthogonalization of the residual with weighting functions, and (iii) determination of the coefficients of expansion of the trial solution in terms of the trial functions. These pertain to the method of weighted residuals, while calculations using Vorobyev’s method of moments involved successive application of Eq. (5) using the integral operators defined by (29). For the method of weighted residuals, the polynomial coefficients are obtained by performing the integrations indicated in Eqs. (17). Also, the orthogonalization procedure in Eqs. (24), amounts to a set of integrations over the semi-infinite interval. These integrals were evaluated by a 32 point Gauss-Laguerre quadrature
Binary equal breakage model
Either of Eqs. (13 and 25) represents the binary equal breakage model. The solutions by MWR using problem-specific polynomials (17) for Q = 0.5 have been presented in Figs. 1 and 2. It was found that when the orthogonalizing function +j(x) is the same as pi(x), the solution settles down to within five decimal places with only three problemspecific polynomials, i.e. n = 3. However, n = 5 displays a strange aberration (Fig. 1). Check of the numerical results for possible errors confirmed the result shown. With increasing II the solution again is restored to that at n = 3; further Table 2 presents the ratio of left hand side to right hand side for this solution. It is clear from Table 2 that this ratio is indeed close to unity. For this reason the solution obtained is referred to as the “accurate” solution. Figures
21
Solutionofpopulationbalanceequationsby MWR
Accurate
solution
solution
-Accurate
Bmory
A
equal breakage
Fig. 1. Effect of crystal breakage on population density in a CMSMPRC.
Fig. 2. Effect of crystal breakage on population density in a CMSMPRC,comparisonof weighting functions.
1 and 2 represent the accurate solution by a continuous curve. Discretized representation has been used to denote solutions in general for visual convenience. Figure 2 shows the results of computation using PSP and $i = & The resulting approximate solution scatters about the accurate solution without settling down. Thus it appears that the choice of I+$= pi is more adequate. Figures 3 and 4 show the calculations with Laguerre polynomials for two different choices of I+%. It is seen that as many as fourteen Laguerre polynomials are required to obtain a reasonably acceptable solution (Fig. 3). In Fig. 4 where $i = +, the solution is not obtained even with fourteen or more Laguerre functions. Figure 5 shows the results of computation using VMM. It is found that excellent solutions are obtained for n 2 2. Similar calculations have been made by Singh[lO] for
Q = 0.2, and the results possess features identical to those at Q = 0.5. Figure 6 summarizes the calculations for Q = 0.5 and 0.2 obtained by MWR using PSP and VMM. Calculations by Randolph have also been shown in the figure. The figure quantitatively brings out the obvious effect of the value of Q on the solution. Binary
uniform
breakage
model
The model is represented by either of Eqs. (14 and 26). No solution has been obtained in the literature for this situation. Also the equation is of general interest to particulate systems in which the particles grow and undergo splitting. The calculations made are identical to those for the binary equal breakage model in the sense the same set of
Table 2. Comparison of accurate solutions to the PBE for binary equal breakage model Q = 0.5 MWR(Laguerre polynomials) ” = 14 VMMn=2
MWR (PSP) n = 6 x
Y
Ratio
0.0 0.2 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0
1@300 0.8178 0.6062 0.3772 0.2413 0.1528 09036 x 10-l 0.4667x 10-l 0.1950x 10-l 06468 x lo-* 0.1399x lo-’ 0.1748x lo-’ 0.1008x lo-’ 0.1864x lo-”
1.000 0.990 1.020 1.050 1.010 10Nl l+tO 1.050 I.001 1.055 1.040 1.035 1.045 1.050
Y
Ratio
1NKtO 1GlO 0.8163 1.100 0.6085 1.080 0.3798 1.090 0.2438 1.030 1.1506 1.120 0.8823x IO-’ 1.200 0.4570x 10-l 0.980 0.1995x 10-l 1.200 0.8098x lo-* 0.890 0.1802x 1O-2 1.500 0.2200x lo-’ 1.450 0.1135x 1o-4 1.380 0.2836x lo-” 1.530
Y
Ratio
l.OtlOO 0.8187 0.6070 0.3750 0.2419 0.1552 0.8961x 10-l 0.4550x 10-l 0.1950x 10-l 0.6383x lo-* 0.1404x 1om2 0.1755x lo-’ 0.1010x lo-’ 0.2074x 10m6
1.000 1.000 1.000 0998 0999 l+nxl 1.001 I@31 1.001 1.001 1.001
1ml 1.010 1.010
P. N.
28
SINCH
and D.
RAMKRISHNA
IO Qz05 -Accurate
solution
n=3
*
Binary
equal
breakage
-Accurate
Blnory
0
J
IO
solution
equal
breakage
2 3 DlmenslonlQss
I
4 5 SIZP. x
6
7
7.5
Fig. 5. Effect of crystal breakage on population density in a CMSMPRCVor Obyev’s method of moments.
Fig. 3. Effect of crystal breakage on population density in a CMSMPRC,comparisonof weight functions.
0105 0.02
CA=05
x MWR . -VMM -Accurate .
.
solution
n:3
-VMM Binary
0 n:4
MWR
A Numerical solution (Randolph)
equal
breakage
x n:5 d n=6 0
n =I4
A
O,(X)=P-XL,w
4 -I
P,(x)=e-xL,(x)
0
Binary
equal
I
2
breakage 3
4
5
6
7
75
Fig. 4. Effect of crystal breakage on population density in a
Fig. 6. Comparison of results obtained by VMM, MWR and numerical integration (Randolph).
CMSMPRC,comparison of weighting functions.
trial functions and orthogonalizing functions have been used. Figures 7 and 8 show the resuits of MWR using PSP for Q = 0.5. The solution again settles down to a final value within n = 4. Table 3 shows the ratio of left hand side to right hand side which is indeed close to unity. Again the accurate solution is represented by a continuous curve. On the other hand, the use of Laguerre polynomials produced no satisfactory results even when fourteen
functions are used as is established by Figs. 9 and 10. Figure 11 represents the results for VMM which again show that very accurate solutions are obtained for n 12. Figure 12 presents the calculations for Q = 0.2 and 0.5 obtained by MWR and VMM. 6. CONCLUSIONS The results in this paper have established that appropriately generated problem-specific polynomials are distinctly superior to the standard polynomials which
29
Solution of population balance equations by MWR
Q=O5 -Accurate .
solution
n.3
en:4 x n:5 A n-6 0 ” :14
E
0
mi(x)=p,(x)
,B
@i(X)=ewX
& n.6
i
6(‘*02klx5)
Vi (X):P,(X) Bi;
L,(X)
$Ji(X)=Ll(X)
.
&nary
uniform
I 1
I 2
breakage
t un,ttcrm
, \
;kag;
,
,
0
,
IO4 0
I
I 4
I
3
Dlmenslonless 0
I
2
5
4
3
Dlmonslonless
6
7
75
51ze.x
SIZ’2,
I
I
I
5
6
7
75
x
Fia. 9. Effect of crystal breakage on population density in a CMSMPRC,comparison of weight functions.
Fig. 7. Effect of crystal breakage on population density in a CMSMPRC.
.
o n.6 @,(x)_p,(x) +i(x)=pi(x) B~nory
,-(x*02 ,-(x+02 “nltorm
klX5) klX5) braakag@
l
I
.
. %
l \
104.
0
I
I
I
1
2
3
Dlmonslonloss
\. 75
Fig. 8. Effect of crystal breakage on population density in a CMSMPRC,comparison of weightingfunctions.
have been used in the past. Singh[lO] has also used PSP to solve nonlinear population balance equations arising in the description of agglomerative particulate systems. That the appropriateness of the weight function depends on its similarity in shape and trend to the solution in question[6] is borne out by Fig. 13 which shows the weight functions e-’ and (16) plotted alongside the solutions for the two models.
II 4
5 SIZQ.
6
75
x
Fig. 10. Effect of crystal breakage on population density in a CMSMPRC,comparison of weightingfunctions. Also Vorobyev’s method of moments has proved particularly effective in the solution of the problems
considered herein. The extension of VMM to nonlinear problems considered herein. The extension of VMM to nonlinear problems however does not appear to have been attempted. In recent times, the orthogonal collocation procedure[ll, 121 has proved a valuable tool in the solution of mathematical equations of physical phenomena. In a sense, the effort here has incorporated the collocation procedure with x = 0 as a collocation point [see Eq. (23)]. Singh[lO] has used the collocation technique in its entirety in combination with problem-
30
P. N.
SINGH
and D.
RAMKRISHNA
Table 3. Comparison of accurate solutions to the PBE for binary uniform breakage
model
Q =O.S MWR(Laguerrepolynomials) n = 14
MWR(PSP) II = 6 x 0.0
0.2 0.5 1.0 1.5
2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.4 6.0
Ratio
Y
14000 0.8191 0.6085 0.3743 0.2332 0.1461 0.8972 x 0.5141 x 0.2545 x 0.9701 x 0.2529 x 0.3700 x 0.4565 x 0.5746 x
IO-’ 10-l 10-l lo-* 10m2 lo-’ 10m4 lo-”
l+OO 0997 0999 1.005 1.001 0997 0.997 1.001 1.005 1.003 0998 0997 0.995 0.990
Ratio
Y
1@000 0.8187 06076 0.3751 0.2335 0.1463 0.8730 x 0.4925 x 0.2623 x 0.1122 x 0.4083 x 0.7938 x 0.1829 x -
10-l 10-l IO-’ 10-l lo-’
IO-’ lo-’
0.4504x lo-’
n=4
VMM
Ratio
Y
1.000
1.0000
I@0
0.950 0,965
0.8191 06086 0.3743 0.2332
0997
1.100 1.004 1.060 0.850 0.800 1.200 1.884 3.199 15.930 10.930 7.500
0
0.1461 0.8979 0.5143 0.2542 0.9697 0.2494 0.3611 0.4532 0.5689
1
2
l@O
I@05
x x x x x
10-l 10-l 10-l lo-’ IO-’ x lo-’ x IO-’ x lo-”
3
1.OOl l+NXl 1.000 l.OlM 1,002 1.003 1.011 1.022 1,032 1.035
4
5
Dimensionless
see,
6
Fig. 12. Comparison of results obtained by VMM and
Fig.
11.Effect of crystal breakage on population density in a CMSMPRC, Vorobyev’smethodof moments.
specific polynomials successfully in the solution of nonlinear population balance equations arising in the theory of agglomerative processes. CONCLUSIONS
AND SIGNIFICANCE
Solutions were obtained for integro-differential equations of population balance models of a continuous crystallizer. These are presented in Tables l-3 and in Figs. 1, 5-8 and 11. Accurate solutions were obtained by using
7 -;15
x MWR.
the specially constructed trial functions referred to in the preceding section. only a small number of which were necessary. The weight function for the orthogonalization procedure, which generates the trial functions, was chosen to be the solution of a simplified version of the population balance model. Figures 3,4,9 and 10 show that standard Laguerre function expansion of the solution suggested by earlier workers [l, 21 has serious limitations in striking contrast with the methods presented herein. Figures 5, 11 and 12 expose an additional method called Vorobyev’s method of moments[8] as a very attractive prospect for the solution of linear population balance equations of dispersed phase systems in which particle breakage is predominant.
31
Solution of population balance equations by MWR
0=05 -~-c-x ___ ,-(x+02
Q R w(l) w, x y
kl x5)
constant defined by Eq. (19) residual number density of crystals nth approximate of w dimensionless crystal size dimensionless number density
Greek symbols a. coefficient of nth approximating function in trial solution
cp. tin 0 6
nth approximating function in trial solution nth function used to orthogonalize residual average holding time Dirac delta function constant parameter in Eq. (1) E fraction of daughter crystal size referred to parent size o weight function in GSOP.
Abbreviations used CMSMPR continuous, mixed suspension, mixed product removal GSOP Gram-Schmidt orthogonalization process MWR method of weighted residuals PBE population balance equation VMM Vorobyev’s method of moments. REFERENCES
0
1
2
3
4
Dlmonslonloss
5 sze,
6
7
15
x
Fig. 13. Shapes and trends of weight functions and accurate solution.
NOMENCLATURE
A a.* B(I) D(I) X k’ k.l L L”
1 P(e) P. P”
bounded linear operator coefficient of xk in the nth degree polynomial number birth rate of crystals of size 1 number death rate of crystals of size I Hilbert space constant in Eq. (8) dimensionless constant in Eq. (10) linear growth rate of crystal nth degree Laguerre polynomials length of crystal orthogonal projection operator orthogonal projection operator nth degree problem-specific polynomial
1. H. M. Hulburt and S. L. Katz, Chem. Engng Sci. 19, 555 (1964). 2. H. M. Hulburt and T. Akiyama, lnd. Engng Chem. Fund/s. 8, 319 (1%9). 3. D. Ramkrishna, Chem. Engng Sci. 26, 1134(1971). 4. G. Subramanian and D. Ramkrishna, Math. Biosci. 10, 1 (1971). 5. B. A. Finlayson, The Method of Weighted Residuals and Variational Principles. Academic Press, New York (1972). 6. D. Ramkrishna, Chem. Engng Sci. 28, 1362(1973). 7. A. D. Randolph and M. A. Larson, Theory of Particulare Processes. Academic Press, New York (1971). 8. Yu. V. Vorobyev, Method of Moments in Applied Mathematics (Translated from the Russian by B. Seckler). Gordon & Breach, New York (1965). 9. A. D. Randolph, lnd. Engng Chem. Fundls. 8, 53 (1%9). 10. P. N. Singh, Ph.D. Dissertation, Indian Institute of Technology, Kanpur (1974). 11. B. A. Finlayson, Chem. Engng Sci. 26, 1081(1971). 12. J. V. Villadsen and W. E. Stewart, Chem. Engng Sci. 22, 1483 (1967). 13. B. Carnahan, H. A. Luther and J. 0. Wilkes, Applied Numerical Methods. Wiley, New York (1969). 14. .I. H. Ahlberg, E. N. Nilson and .I. L. Walsh, The Theory of Splines and Their Applications. Academic Press, New York (1967).
