- Email: [email protected]
Integral method of analysis for chemical reaction in a nonisothermal finite cylindrical catalyst pellet-I. Dirichlet problem
Chemical Printed
Engineering Science, in Great Elritam.
Vol.
42,
No.
I, pp.
27-33,
1987.
0009-2509/87 Pergamon
$3.00 + 0.00 Journals Ltd.
INTEGRAL METHOD OF ANALYSIS FOR CHEMICAL REACTION IN A NONISOTHERMAL FINITE CYLINDRICAL CATALYST PELLET-I. DIRICHLET PROBLEM SURYANARYANA
MUKKAVILLI,t LAWRENCE L. TAVLARIDESts and CHARLES V. WITTMANNI Department of Chemical Engineering, Illinois Institute of Technology Chicago, IL60616, U.S.A. and tDepartment of Chemical Engineering and Materials Science, Syracuse University, Syracuse, NY 13210, U.S.A. (Received
22 October
1984)
Abstract-An integral equation method is developed to solve the problem of diffusion and reaction in a porous nonisothermal finite cylindrical pellet in the absence of external transport resistances. Green’s function method is applied to transform the partial differential equation for concentration into a Fredholm integral equation. A modified Green’s function method is developed to accelerate the convergence of the partial eigen series by an order of two while exploiting the symmetry properties of the classical Green’s function formulation. The resulting integral equation is solved by a Newton-Kantorovich iteration scheme to obtain effectiveness factors for various nonlinear reaction rate forms.
It appears from the above review that only a limited effort has been made in the modeling of finite cylindrical pellets. Further, no attempts at solving the corresponding problem with external transport resistances (Robin problem) have been reported. In this paper we develop integral equation formulations using the Green’s function approach (Carslaw and Jaeger, 1973). Amundson and Schilson (1961) obtained the Green’s function for reaction in a sphere and solved the resulting linear Fredholm integral equation using a successive approximation technique. Kesten (1969) applied the Green’s function analysis to obtain concentration profiles for ammonia decomposition in a spherical catalyst pellet. Dente ef al. (1969a, b) also used the successive approximation
INTRODUCTION
Even though the finite cylinder is a common geometry for industrial catalysts only a few researchers attempted to model diffusion and reacfion in porous finite cylindrical catalyst pellets. The more involved problems of uniqueness and stability of the steady state in finite cylinders have attracted even less attention. The earliest work on finite cylindrical pellets is that of Aris (1975), in which he obtained analytical expressions for reaction in an isothermal finite cylinder with no external transport resistances (Dirichlet boundary conditions). Gunn (1967) and Amundson and Luss (1967) obtained similar analytical expressions. Gunn (1967) solved the corresponding problem for hollow cylindrical pellets and gives effectiveness factor vs Thiele modulus curves for several values of length to diameter ratio of the cylinder. Ho and Hsiao (1977) obtained approximate effectiveness factors for a finite cylindrical pellet based on singular perturbation on an infinite cylinder. Their work is restricted to first order reaction in an isothermal pellet as are the results of
method integral Dixit
to obtain multiplicity limits using a Votterra equation formulation. and Tavlarides (1982) were the first to use a
cylindrical catalysts is that of Sorenson et al. (1973). They used the orthogonal collocation technique to solve the problem of first order reaction in a nonisothermal pellet with Dirichlet boundary conditions. A predictor
Newton-Kantorovich iteration scheme to solve nonlinear Fredholm integral equations for the analysis of reaction in a sphere and an infinite cylinder. They demonstrated the utility of this technique for the case of Fischer-Tropsch synthesis. In this paper we present an integral equation formulation for reaction in a finite cylinder subjected to Dirichlet boundary conditions. We begin with the mathematical description of the Dirichlet Boundary Value Problem (DBVP) and the corresponding classical Green’s function formulation. This classical Green’s function method will be modified to accelerate the convergence of the partial eigen series using a modified Green’s function formulation. This modified
+Present address: Department of Chemical Engineering and Materials Science, Syracuse University, Syracuse, NY 13210, U.S.A. 5Author to whom correspondence should he addressed. “Present address: Hercules Research Center, Wilmington, DE 19894, U.S.A.
Green’s function formulation illustrates the convergence of the eigen series by TWO orders of magnitude while preserving the numerically advantageous symmetry properties of the classical Green’s function solution. The Newton-Kantorovich iteration scheme and some of the numerical results will be presented in
previous investigations (Aris, 1975; Amundson Luss, 1967; Gunn, 1967). The most useful work in the. modeling of
and finite
27
SURYANARYANA MUKKAWLLI
28 terms
of
effectiveness
factor
plots.
In a succeeding
paper (Mukkavilii et al., 1987) we present the corresponding results for the case with external transport limitations (Robin problem). The advantages of this integral solution are: (i) the accelerated convergence of the eigen series for the modified Green’s function formulation, and (ii) the quadratic rate of convergence 1971) of the iteration scheme.
(Baker
et al.,
et al.
and (6) where y = C/C, and 0 = T/T, are the dimensionless concentration and temperature, C, and TB are the bulk concentration and temperature, R = R’(T,C)/ R’(T,, C,) is the dimensionless reaction rate, &=
MATHEMATICAL
FORMULATION
Consider a chemical reaction taking place inside a finite cylindrical catalyst pellet, shown in Fig. 1. Transport of heat and mass through the gas transport zone to the external surface of the catalyst pellet, denoted by XI, is characterized by the respective Nusselt and Sherwood numbers. Pore diffusion and catalytic reaction take place in the pellet interior, denoted by 0. Following
the classical
formulation
described
by
Aris (1975), which assumes a homogeneous porous pellet and using effective transport coefficients, the steady state concentration and temperature profiles for a single chemical reaction are given by, D,V’C
= pbSgR’(T,C)
in D u X2
(1)
is the well known
is the Prater
in Q uXI
(2)
where C and Tare respectively the concentration and temperature, D, and K, are the effective diffusivity and thermal conductivity, pb is the pellet density, S, is the catalytic area per unit mass, ( -AH,) is the heat of reaction and R’ is the rate of reaction per unit area. For a finite cylindrical pellet the Laplacian operator
where r’ and Z’ are the radial and axial coordinates respectively, as shown in Fig. 1. Non-dimensionalizing
the eqs (l), (2) and (3) results,
- (AH,)D,C,
number,
and
V*O=
-/?c$2R(l?,y)
in R uX! in fi uaQ
(4) (5)
is the aspect ratio, with R. pellet and L its length.
(9)
the radius
of the catalyst
For the DBVP, negligible interphase transport resistances would mean that the concentration and temperature on the external surface (i30) of the catalyst pellet are the same as the bulk values. The appropriate boundary conditions then are for the DBVP on aR
(IOa)
and 8 = Note and/or
1
that dR is the external
on XI.
(lob)
surface denoted
by r =
1
z = 0, 1.
Prater analysis Prater (1958) simplified the system of two differential equations describing temperature and concentration profiles into a single differential equation and an algebraic equation relating temperature and concentration in a linear fashion. The Prater analysis shows that for the DBVP, inQuX!.
(11)
Hence the DBVP for diffusion and reaction in a finite cylinder is defined by eq. (4), (lOa) and (11).
CLASSICAL
Fig. 1. Finite cylindrical catalyst pellet
Ro L
Q=l+p(l-y) = @R(&y)
(8)
J%Ts
EC-
in V’y
(7)
modulus
y=l -(-AH,)pbS,R’(T,C)
%C,)
D,C,
Thiele B=
and, K,V’T=
R$p,S,R’( J
GREEN’S
FUNCTION
FORMULATION
The classical Green’s function is defined as the solution for a linear differential equation, with a unit source, accompanied by linear homogeneous boundary conditions (Baker et al., 1971; Carslaw and Jaeger, 1959). The source is usually taken as a dirac delta function. Numerous examples of the use of Green’s functions with applications in diverse areas of physical sciences may be found in Carslaw and Jaeger (1973). The problem of diffusion and reaction inside a catalyst pellet may be viewed as one with an infinite number of “sinks”, corresponding to the consumption of the
Integral method of analysis-I reactant,
and,
as such is an analogue
of the Green’s
function with a unit point “source”. In the context of the DBVP defined earlier, the classical Green’s function, G,(r, z; r,,, q,), with a unit source at (ro,zo) is defined by the following set of
29
It can be seen symmetric about
from
eqs (15)
and (18)
that
Cd is
(r,, z,), i.e. the source point Gd (r, z; ra, zo) = Cd (re, zO; r, z). Further Gd has discontinuous first derivatives dG,/az and aG,,/ar at z = z0 and r = r. respectively [see eqs (12) and (16)].
equations V=G,, =
+(r
with homogeneous
MODIFIED
--r,,)d(z-z,)
boundary Cd =0
(12)
condition on &2.
(13)
In order to solve the system of eq. (12) and (13) for Cd, we use a partial eigen series expansion in the axial direction, z. Define the Transform f. =-
1
r
6 nm s 0
j’sinnnzdz
(14a)
and, the Inverse f=
g “==
f,sinnnz;
(14b)
= ;.
(14c)
I
where a,,
Applying the transform (14a) to eq. (12) and using the boundary conditions Cd = 0 at z = 1 and at z = 0 from eq. (13), we obtain, Gd =
f n=
1
gn(r; ro) sin nxz sin rz7~z~
(15)
where g. (r; ro) satisfies,
1d dgn --(r+-c’ n2n’y. rdr dr along
with the boundary g.(r;re)
=
(16)
(17a)
and g. (r; rO) is finite in r E (0,l).
(17b)
In addition gn(r; r-a), the one dimensional Green’s function in the radial direction satisfies, at the source point r = ro, s”lr=r+
-~ dgn dr
=
r=rO+
converges
=
--,
2 i-0
(17d)
1[
Z,(&nnr)Z,(&n~r,)K,(en~)
modified
I[ -2
Z,(en~r)Z,(ennr,)K,(Enlc) Z,(~nxr,)& Z&in)
of coefficientsf, as l/n’
+
in the
1 as n -
a.
Green’s
function =
-;
1
for the DBVP,
Gd,. as
6(r-rr,)~d(z;.70)
(19)
along with the boundary conditions, eq. (13). In order to accelerate the convergence of the partial eigen series for G,, we require that the source functions &, (z; zO) be such that SdECk with preferably as high development that follows
a “k” as possible. we choose
In
the
SdEQ i.e. [,, has a piecewise continuous The resulting Gd is such that, GdECJ
z,(&nnr)K&nxr,)
Io(enn)
FORMULATION
Hence the coefficients in the eigen series for Cd converge only as l/n2 as n - co. Clearly this slow rate of convergence must be improved upon for numerical purposes. In order to accomplish this we define a
1
J
1.
(En7rr)
r < r0
(18) r 1 r0
derivative
(with respect
of eqs (16) and (17) is given by
-2
9Jr;rd
(x) E Ck, then the sequence
(17c)
= Snlr-r,,
and
The solution
Iff
eigen series for f(x)
V*G, 1
FUNCTION
Theorem
conditions = 0 at r =
GREEN’S
At this point it is instructive to examine the partial eigen series obtained for Cd. We note that Cd EC’, C“ being the class of functions with piecewise continuous derivatives up to and including the kth order. This discontinuous nature of the derivatives of G,, has an important bearing on the convergence of the eigen series obtained above. Since only a partial eigen series expansion in the axial direction, z, is used in the above development, the discontinuity in the derivative, dG,/ik, need not be considered. However, if, a full eigen series expansion, i.e. in both axial (z) and radial (r) directions is used, the following observations regarding the convergence of the partial eigen series are appropriately extended. The following theorem (Tolstov, 1962) on the convergence of the coefficients of the eigen series is of importance.
to z).
of order
1.
30
SURYANARYANA
i%fUKKAVILLI
This would mean that the eigen series for G, will converge at a rate two orders of magnitude faster than that for Gd. It should be noted that the inherent symmetry in the Green’s function, Cd, is of great numerical advantage. This symmetry property, resulting from the symmetry of the dirac-delta source function [S (z - ze) = 6 (zO -z)], must bc preserved in Cd. In the Appendix, we show that the source id must also satisfy the homogeneous boundary conditions in order to accelerate the convergence of the partial eigen these limitations we choose
L(z; 20) =
series
z(l -zo) { %(l-2)
for
6,.
z <
2,;
z >
ZfJ,
then
given
where
Green’s
function
for
DBVP,
cm
gd,
is again
given
by eq. (18).
is
Comparing
eqs (15) and (21) we note that the partial eigen series for G,, converges at a rate two orders faster than that for Gd. Further, C?, (r, 1; r,,, zO) = 6, ( rO, zO; r, z).
INTEGRAL
METHOD
OF ANALYSIS
The solution for the DBVP for reaction and diffusion inside the finite cylindrical catalyst pellet may now be expressed in terms of the modified Green’s function, Cd. The inhomogeneous boundary conditions, posed by eq. (lOa), may be reduced into homogeneous boundary conditions by using the transformation, 8=1-y.
ss 1
1
?= ‘
-1
NUMERICAL
by,
g,,(r; re)
Once eq. (24) is solved for the concentration profile, y, inside the catalyst pellet, for a given Thiele modulus, 4, and reaction rate expression, R(y), the effectiveness factor, q, for the pellet may be calculated. The effectiveness factor, q. is given by
Within
i.e. the source in the axial direction for G, is the onedimensional Green’s function in the axial direction for the DBVP. The modified
e?t Cd.
0
dr dz.
R(y(r,z))r
PROCEDURE
AND
(25)
RESULTS
The nonlinear Fredholm integral equation, eq. (24), is solved using a Prime 400 computer at the Illinois Institute of Technology. The iteration scheme and some of the typical results are presented below.
Numerical
procedure
The modified Green’s function, cd, given by eqs (21) and (18) is calculated for a given value of the aspect ratio, c, using the well known Euler-McLauren summation technique (Dahlquist and Bjorck, 1974), which is essentially a repeated averaging method. The integrals in eqs (24) and (25) are obtained by using Gauss-Radau quadrature with an 8 x 8 point interior grid. Accuracy of the integration scheme is checked by using a 16 x 16 grid. The nonlinear integration equation, eq. (24), is solved by using the Newton-Kantorovich iteration scheme, described below. In operator form, eq. (24) may be written as F(y) The
= 0.
Newton-Kantorovich
(26)
iteration
scheme
is given
F(Yi)
by (27)
(22)
This transformation known as the “Superposition Method”, overcomes the difficulties involved in differentiating Gd on the boundary aR, as would be the case for the inhomogeneous boundary conditions. Substituting eq. (22) in eq. (4) and multiplying the
where, 6F, is the variation the present case is given
result with integrating
The Newton-Kantorovich iteration technique is quadratically convergent and is a superior iteration scheme to the linearly convergent successive substitution tnethod used by Dente et al. (1969a, b). When the iteration scheme satisfies a root mean
& and subtracting from eq. (19) x J, and with respect to r and z gives,
1 s 0
Y(r,,z)id(z;zo)dz
= 4’
1
1
ss0
Cl
x &(r,
R(1
-Y’(r,z))
z; rO, z,,)r dr dz.
Since &, and Gd are symmetric we may transpose the (re,ze) surface to the (r,z) surface. Thus
(23) from
1 s 0
Y(r,z,)i,(z;z,)dz, 1 X
ss0
= fz(l
-z)-@
I 0
R(Y
(re, ze))C% (r.z;
re,ze)re
dr0 dze.
(24)
Equation (24) gives the solution for the concentration profils inside the finite cylindrical pellet in a closed form.
of F with respect by
to JJand for
6F,,=$$
square
criteria
J[k
4 C (yi + 1
yij2] G (error)
-
m
the corresponding
effectiveness
Results For the DBVP
the rate
dimensionless
concentration, R(Y)
= Y” exP
factor,
expression y, would
YB(1 -Y) l-fB(l-Y)
q, is calculated.
in terms be,
of the
Integralmethod of analysis-1 where n is the reaction order and
E Y=RT, is the Arrhenius number. In Fig. 2, effectiveness factor vs Thiele modulus curves are given for a first order reaction (n = 1) in an isothermal (/I = 0) pellet. The parameter for the curves is the aspect ratio, E, in the range of 0.25-2.0. The analytical results of Gunn (1967) and the perturbation approximation of Ho and Hsiao (1977) are shown for comparison. It is noted that the numerical results are in good agreement with those of Gunn (1967) for Thiele modulus values less than 10.0. As the Thiele modulus increased above 10.0, the numerical calculations become less accurate, with a relative error of about 3 % at I$ = 40. This is due to the fact that for large values of c#, corresponding to severe internal diffusional limitations, the reaction zone is confined to a narrow shell
31
near the outer surface of the catalyst pellet. Such a problem with steep concentration gradients can be solved by using two regions of integration, with the interior of the pellet containing a relatively sparse grid. It must be noted that the range of Thiele modulus for which stiffness problems arise is within the asymptotic range where the effectiveness factor q a l/4_ It is shown (Aris, 197.5; Amundson and Luss, 1967) that, by modifying the definition of the Thiele modulus to include as the scaling parameter for length, the ratio volume/surface area of the pellet, the effectiveness factor curves for various E values can be brought together. This leads to the modified Thiele modulus, rh
+m =(2;E) In Fig. 3, the effectiveness factor versus modified Thiele modulus curves are given for various nonlinear reaction rate models. The results of the present study
2% l
0.25
f 4
Y:? 2.0
0 Gunn (1967) R&.(41 l HO ond Hsiao (1977) Ref.61
Fig. 2. Effectiveness factor vs Thiele modulus for first-order reaction-Dirichlet
_
0.01
*
SorenPen
problem.
g & (1973) Ref. (7)
0. I MODIFIED
THIEtJC MODULUS.
Fig. 3. Effectiveness factor vs modified Thiele modulus problem.
1.0 #‘m
ICI0
for nonlinear reaction rate forms-Dirichlet
32
SLJRYANARYANA
are compared with those from orthogonal collocation calculations (Sorensen et al., 1973). It is noted that our numerical calculations agree within 1%. Sorensen et al. (1973) calculated the bounds on the multiplicity analysis by using an initial value problem approach similar to that due to Dente et al. (1969a). The present study can be readily extended to obtain the corresponding Volterra integral equation formulation.
T
temperature
Ts Y
temperature of the bulk, K concentration of reactant in the pellet, dimensionless axial coordinate, dimensionless axial source coordinate, dimensionless
z 20 Greek P
CONCLUSIONS
An integral equation method is developed to model diffusion and reaction in a finite cylindrical catalyst pellet with no external transport limitations. The Green’s function kernel is modified to accelerate the convergence of the partial eigen series by two orders of magnitude. A Newton-Kantorovich iteration scheme is
used to solve the resulting nonlinear integral equation. Our numerical results for various reaction rate expressions compare well with previous results (Gunn, 1967; Ho and Hsiao 1977). It is noted that the modified Green’s function method exploits the symmetry properties of the classical Green’s function formulation while the convergence of the partial eigen series is accelerated by two orders of magnitude. Further, the convergence of the rate of quadratic Newton-Kantorovich technique is faster than the successive substitution. Acknowledgments-Thr general support and the assistantship (S.M.) provided by the Department of Chemical Engineering of Illinois Institute of Technology through Dr. D. T. Wasan and the Department of ChemicaI Engineering and Materials Science of Syracuse University are gratefully acknowledged. NOTATION
concentration of reactant in the pellet, g mol/cm3 concentration of reactant in the bulk, g mol/cm’ effective diffusivity of reactant in the porous catalyst pellet, cm’/s Green’s function for the finite cylindrical pellet for DBVP modified Green’s function for the finite cylindrical pellet for DBVP modified Besel function modified Bessel function effective thermal conductivity in the porous catalyst pellet, Cal/cm/K/s number of grid points in the radial direction length of the finite cylindrical catalyst pellet, cm number of grid points in the axial direction radial coordinate, dimensionless radial source coordinate, dimensionless dimensionless reaction rate reaction rate, g mol/cm’/s radius of the finite cylindrical catalyst pellet, cm catalytic surface area. cm’/n
al.
MUKKAVILLI~~
Y s -AH, E L ;: :, Pb
R an
of the pellet, K
letters
Prater temperature, dimensionless Arrhenius number, dimensionless dirac delta function heat of reaction, Cal/g mol aspect ratio of the finite cylindrical pellet, ROIL one dimensional Green’s function in the axial direction effectiveness factor, dimensionless temperature of the pellet, dimensionless Thiele modulus, dimensionless modified Thieie modulus, dimensionless bulk density of the pellet, g/cm3 interior of the pellet external surface of the pellet
REFERENCES
Amundson, N. R. and Schilson, R. E., 1961,
Intraparticle diffusion and conduction in porous catalysts-I. Single reactions. Chem. Engng Sci. 13, 226-236. Amundson, N. R. and Luss, D., 1967, On a conjecture of Aris: proof and remarks. A.I.Ch.E. J. 13, 759-761. Aris, R., 1975, The Mathrmaticul Theory of Diffusion and Reaction in Permeable Cutulysts, Vol. I. Calendron Press, Oxford. Baker, B. S., Gidaspow, D. and Wasan, D. T., 1971, Thermal phenomena in fuel cells and batteries. Advances in Electrochemistry and Electrochemical Engineering, Vol. 8, pp. 63-156. interscience Publishers, New York. Carslaw, H. S. and Jaeger, J. C., 1959, Conducrion ofHear in Solids. Calendron Press, Oxford. Dahlquist, G. and Bjiirck, A., 1974, Numerical Methods (translated by Anderson, C.). Prentice Hall, Englewood Cliffs, New Jersey. Derite, M., Biardi, G. and Ranzi, E., 1969a, Un nuovo methods per il calcolo dei profili di concentrazione a de1 fattore di utilivazione per reazioni complesse evolvenlise i granuli catalitici non isothermi. Quad. Ing. Chem. Ital. 5, 65-72.
Dente, M., Biardi, G. and Ranzi, G., 1969b, Problemi riguardantila determinazionediretta dei limiti statici di
‘ignizione’ e ‘spegnimento’ per reazioni cineticamente semplici evolventici in granuli catalitici non isothermi. Quad. Ing. Chem. Ital. 5, 122-129. Dixit, R. S. and Tavlarides, L. L., 1982, Integral method of analysis of Fischer-Tropsch synthesis reactions in a catalyst pellet. Chem. Engng Sci. 37, 539-544. Gunn, D. J., 1967, Diffusion and chemical reaction in catalysis and absorption. Chem. Engng Sci. 22, 1439-1455. Ho, T. C. and Hsiao, G. C., 1977, Estimation of the effectiveness factor for a cylindrical catalyst support: a singular perturbation approach. Chem. Engng Sci. 32, 63-66. Kesten, A. S., 1969, An integral equation method for evaluating the effects of film and pore diffusion of heat and mass on reaction rates in porous catalyst particles. A.I.Ch.E. J 15, 128-131. Mukkavilli, S., Tavlarides, L. L. and Wittmann, C. V., 1987, _. . . . ^ ..^ -.Integral method ot analysis tor chemtcal reaction in a
Integral method of analysis-1
33
nonisothermal finite cylindrical catalyst pellet-II. Robin problem. Chem. Engng Sci. 42, 3540. Sorensen, J. Q., Guertin, E. W. and Stewart, W. E., 1973, Computational models for cylindrical catalyst particles. A.I.Ch.E.
J. 19, 969-975.
Tolstov, G. P., 1962, Fourier Series (translated by Silverman, R. A.). Dover Publ., New York. Consider
APPENDIX
Consider f(x)
??f$?&-
such that,
f;(x)
7
dx.
(A-3)
f, and fz such that,
fi(O)=JX)=O
0 < x < a;
f(x)= {h(x) a=sx<
1,
(A-1)
and define the transform-inverse pair J =2JiJ(x)sinnnxdx, J(x) =
s
and
fi (a) =X(a)-
(A-4)
Then, the first two terms on the right-hand side of eq. (A-3) cancel each other and As converge as I/I?. If in addition. (A-2a)
F f. sin nrtx. (A-2b) n=1 Using the product integration and assuming thatfr andh are continuously differentiable in the respective intervals, we have from (A-Za),
f;(o)
=fi(l)
= 0
and
f;(a)
=&(a)
(A-5)
then, f. s converge as 1/n3. Thus we see that in addition to continuity of the function (and its derivatives) we require that the function (and its derivatives) also satisfy the homogeneous boundary conditions in order to accelerate the convergence of the eigen series.
Engineering Science, in Great Elritam.
Vol.
42,
No.
I, pp.
27-33,
1987.
0009-2509/87 Pergamon
$3.00 + 0.00 Journals Ltd.
INTEGRAL METHOD OF ANALYSIS FOR CHEMICAL REACTION IN A NONISOTHERMAL FINITE CYLINDRICAL CATALYST PELLET-I. DIRICHLET PROBLEM SURYANARYANA
MUKKAVILLI,t LAWRENCE L. TAVLARIDESts and CHARLES V. WITTMANNI Department of Chemical Engineering, Illinois Institute of Technology Chicago, IL60616, U.S.A. and tDepartment of Chemical Engineering and Materials Science, Syracuse University, Syracuse, NY 13210, U.S.A. (Received
22 October
1984)
Abstract-An integral equation method is developed to solve the problem of diffusion and reaction in a porous nonisothermal finite cylindrical pellet in the absence of external transport resistances. Green’s function method is applied to transform the partial differential equation for concentration into a Fredholm integral equation. A modified Green’s function method is developed to accelerate the convergence of the partial eigen series by an order of two while exploiting the symmetry properties of the classical Green’s function formulation. The resulting integral equation is solved by a Newton-Kantorovich iteration scheme to obtain effectiveness factors for various nonlinear reaction rate forms.
It appears from the above review that only a limited effort has been made in the modeling of finite cylindrical pellets. Further, no attempts at solving the corresponding problem with external transport resistances (Robin problem) have been reported. In this paper we develop integral equation formulations using the Green’s function approach (Carslaw and Jaeger, 1973). Amundson and Schilson (1961) obtained the Green’s function for reaction in a sphere and solved the resulting linear Fredholm integral equation using a successive approximation technique. Kesten (1969) applied the Green’s function analysis to obtain concentration profiles for ammonia decomposition in a spherical catalyst pellet. Dente ef al. (1969a, b) also used the successive approximation
INTRODUCTION
Even though the finite cylinder is a common geometry for industrial catalysts only a few researchers attempted to model diffusion and reacfion in porous finite cylindrical catalyst pellets. The more involved problems of uniqueness and stability of the steady state in finite cylinders have attracted even less attention. The earliest work on finite cylindrical pellets is that of Aris (1975), in which he obtained analytical expressions for reaction in an isothermal finite cylinder with no external transport resistances (Dirichlet boundary conditions). Gunn (1967) and Amundson and Luss (1967) obtained similar analytical expressions. Gunn (1967) solved the corresponding problem for hollow cylindrical pellets and gives effectiveness factor vs Thiele modulus curves for several values of length to diameter ratio of the cylinder. Ho and Hsiao (1977) obtained approximate effectiveness factors for a finite cylindrical pellet based on singular perturbation on an infinite cylinder. Their work is restricted to first order reaction in an isothermal pellet as are the results of
method integral Dixit
to obtain multiplicity limits using a Votterra equation formulation. and Tavlarides (1982) were the first to use a
cylindrical catalysts is that of Sorenson et al. (1973). They used the orthogonal collocation technique to solve the problem of first order reaction in a nonisothermal pellet with Dirichlet boundary conditions. A predictor
Newton-Kantorovich iteration scheme to solve nonlinear Fredholm integral equations for the analysis of reaction in a sphere and an infinite cylinder. They demonstrated the utility of this technique for the case of Fischer-Tropsch synthesis. In this paper we present an integral equation formulation for reaction in a finite cylinder subjected to Dirichlet boundary conditions. We begin with the mathematical description of the Dirichlet Boundary Value Problem (DBVP) and the corresponding classical Green’s function formulation. This classical Green’s function method will be modified to accelerate the convergence of the partial eigen series using a modified Green’s function formulation. This modified
+Present address: Department of Chemical Engineering and Materials Science, Syracuse University, Syracuse, NY 13210, U.S.A. 5Author to whom correspondence should he addressed. “Present address: Hercules Research Center, Wilmington, DE 19894, U.S.A.
Green’s function formulation illustrates the convergence of the eigen series by TWO orders of magnitude while preserving the numerically advantageous symmetry properties of the classical Green’s function solution. The Newton-Kantorovich iteration scheme and some of the numerical results will be presented in
previous investigations (Aris, 1975; Amundson Luss, 1967; Gunn, 1967). The most useful work in the. modeling of
and finite
27
SURYANARYANA MUKKAWLLI
28 terms
of
effectiveness
factor
plots.
In a succeeding
paper (Mukkavilii et al., 1987) we present the corresponding results for the case with external transport limitations (Robin problem). The advantages of this integral solution are: (i) the accelerated convergence of the eigen series for the modified Green’s function formulation, and (ii) the quadratic rate of convergence 1971) of the iteration scheme.
(Baker
et al.,
et al.
and (6) where y = C/C, and 0 = T/T, are the dimensionless concentration and temperature, C, and TB are the bulk concentration and temperature, R = R’(T,C)/ R’(T,, C,) is the dimensionless reaction rate, &=
MATHEMATICAL
FORMULATION
Consider a chemical reaction taking place inside a finite cylindrical catalyst pellet, shown in Fig. 1. Transport of heat and mass through the gas transport zone to the external surface of the catalyst pellet, denoted by XI, is characterized by the respective Nusselt and Sherwood numbers. Pore diffusion and catalytic reaction take place in the pellet interior, denoted by 0. Following
the classical
formulation
described
by
Aris (1975), which assumes a homogeneous porous pellet and using effective transport coefficients, the steady state concentration and temperature profiles for a single chemical reaction are given by, D,V’C
= pbSgR’(T,C)
in D u X2
(1)
is the well known
is the Prater
in Q uXI
(2)
where C and Tare respectively the concentration and temperature, D, and K, are the effective diffusivity and thermal conductivity, pb is the pellet density, S, is the catalytic area per unit mass, ( -AH,) is the heat of reaction and R’ is the rate of reaction per unit area. For a finite cylindrical pellet the Laplacian operator
where r’ and Z’ are the radial and axial coordinates respectively, as shown in Fig. 1. Non-dimensionalizing
the eqs (l), (2) and (3) results,
- (AH,)D,C,
number,
and
V*O=
-/?c$2R(l?,y)
in R uX! in fi uaQ
(4) (5)
is the aspect ratio, with R. pellet and L its length.
(9)
the radius
of the catalyst
For the DBVP, negligible interphase transport resistances would mean that the concentration and temperature on the external surface (i30) of the catalyst pellet are the same as the bulk values. The appropriate boundary conditions then are for the DBVP on aR
(IOa)
and 8 = Note and/or
1
that dR is the external
on XI.
(lob)
surface denoted
by r =
1
z = 0, 1.
Prater analysis Prater (1958) simplified the system of two differential equations describing temperature and concentration profiles into a single differential equation and an algebraic equation relating temperature and concentration in a linear fashion. The Prater analysis shows that for the DBVP, inQuX!.
(11)
Hence the DBVP for diffusion and reaction in a finite cylinder is defined by eq. (4), (lOa) and (11).
CLASSICAL
Fig. 1. Finite cylindrical catalyst pellet
Ro L
Q=l+p(l-y) = @R(&y)
(8)
J%Ts
EC-
in V’y
(7)
modulus
y=l -(-AH,)pbS,R’(T,C)
%C,)
D,C,
Thiele B=
and, K,V’T=
R$p,S,R’( J
GREEN’S
FUNCTION
FORMULATION
The classical Green’s function is defined as the solution for a linear differential equation, with a unit source, accompanied by linear homogeneous boundary conditions (Baker et al., 1971; Carslaw and Jaeger, 1959). The source is usually taken as a dirac delta function. Numerous examples of the use of Green’s functions with applications in diverse areas of physical sciences may be found in Carslaw and Jaeger (1973). The problem of diffusion and reaction inside a catalyst pellet may be viewed as one with an infinite number of “sinks”, corresponding to the consumption of the
Integral method of analysis-I reactant,
and,
as such is an analogue
of the Green’s
function with a unit point “source”. In the context of the DBVP defined earlier, the classical Green’s function, G,(r, z; r,,, q,), with a unit source at (ro,zo) is defined by the following set of
29
It can be seen symmetric about
from
eqs (15)
and (18)
that
Cd is
(r,, z,), i.e. the source point Gd (r, z; ra, zo) = Cd (re, zO; r, z). Further Gd has discontinuous first derivatives dG,/az and aG,,/ar at z = z0 and r = r. respectively [see eqs (12) and (16)].
equations V=G,, =
+(r
with homogeneous
MODIFIED
--r,,)d(z-z,)
boundary Cd =0
(12)
condition on &2.
(13)
In order to solve the system of eq. (12) and (13) for Cd, we use a partial eigen series expansion in the axial direction, z. Define the Transform f. =-
1
r
6 nm s 0
j’sinnnzdz
(14a)
and, the Inverse f=
g “==
f,sinnnz;
(14b)
= ;.
(14c)
I
where a,,
Applying the transform (14a) to eq. (12) and using the boundary conditions Cd = 0 at z = 1 and at z = 0 from eq. (13), we obtain, Gd =
f n=
1
gn(r; ro) sin nxz sin rz7~z~
(15)
where g. (r; ro) satisfies,
1d dgn --(r+-c’ n2n’y. rdr dr along
with the boundary g.(r;re)
=
(16)
(17a)
and g. (r; rO) is finite in r E (0,l).
(17b)
In addition gn(r; r-a), the one dimensional Green’s function in the radial direction satisfies, at the source point r = ro, s”lr=r+
-~ dgn dr
=
r=rO+
converges
=
--,
2 i-0
(17d)
1[
Z,(&nnr)Z,(&n~r,)K,(en~)
modified
I[ -2
Z,(en~r)Z,(ennr,)K,(Enlc) Z,(~nxr,)& Z&in)
of coefficientsf, as l/n’
+
in the
1 as n -
a.
Green’s
function =
-;
1
for the DBVP,
Gd,. as
6(r-rr,)~d(z;.70)
(19)
along with the boundary conditions, eq. (13). In order to accelerate the convergence of the partial eigen series for G,, we require that the source functions &, (z; zO) be such that SdECk with preferably as high development that follows
a “k” as possible. we choose
In
the
SdEQ i.e. [,, has a piecewise continuous The resulting Gd is such that, GdECJ
z,(&nnr)K&nxr,)
Io(enn)
FORMULATION
Hence the coefficients in the eigen series for Cd converge only as l/n2 as n - co. Clearly this slow rate of convergence must be improved upon for numerical purposes. In order to accomplish this we define a
1
J
1.
(En7rr)
r < r0
(18) r 1 r0
derivative
(with respect
of eqs (16) and (17) is given by
-2
9Jr;rd
(x) E Ck, then the sequence
(17c)
= Snlr-r,,
and
The solution
Iff
eigen series for f(x)
V*G, 1
FUNCTION
Theorem
conditions = 0 at r =
GREEN’S
At this point it is instructive to examine the partial eigen series obtained for Cd. We note that Cd EC’, C“ being the class of functions with piecewise continuous derivatives up to and including the kth order. This discontinuous nature of the derivatives of G,, has an important bearing on the convergence of the eigen series obtained above. Since only a partial eigen series expansion in the axial direction, z, is used in the above development, the discontinuity in the derivative, dG,/ik, need not be considered. However, if, a full eigen series expansion, i.e. in both axial (z) and radial (r) directions is used, the following observations regarding the convergence of the partial eigen series are appropriately extended. The following theorem (Tolstov, 1962) on the convergence of the coefficients of the eigen series is of importance.
to z).
of order
1.
30
SURYANARYANA
i%fUKKAVILLI
This would mean that the eigen series for G, will converge at a rate two orders of magnitude faster than that for Gd. It should be noted that the inherent symmetry in the Green’s function, Cd, is of great numerical advantage. This symmetry property, resulting from the symmetry of the dirac-delta source function [S (z - ze) = 6 (zO -z)], must bc preserved in Cd. In the Appendix, we show that the source id must also satisfy the homogeneous boundary conditions in order to accelerate the convergence of the partial eigen these limitations we choose
L(z; 20) =
series
z(l -zo) { %(l-2)
for
6,.
z <
2,;
z >
ZfJ,
then
given
where
Green’s
function
for
DBVP,
cm
gd,
is again
given
by eq. (18).
is
Comparing
eqs (15) and (21) we note that the partial eigen series for G,, converges at a rate two orders faster than that for Gd. Further, C?, (r, 1; r,,, zO) = 6, ( rO, zO; r, z).
INTEGRAL
METHOD
OF ANALYSIS
The solution for the DBVP for reaction and diffusion inside the finite cylindrical catalyst pellet may now be expressed in terms of the modified Green’s function, Cd. The inhomogeneous boundary conditions, posed by eq. (lOa), may be reduced into homogeneous boundary conditions by using the transformation, 8=1-y.
ss 1
1
?= ‘
-1
NUMERICAL
by,
g,,(r; re)
Once eq. (24) is solved for the concentration profile, y, inside the catalyst pellet, for a given Thiele modulus, 4, and reaction rate expression, R(y), the effectiveness factor, q, for the pellet may be calculated. The effectiveness factor, q. is given by
Within
i.e. the source in the axial direction for G, is the onedimensional Green’s function in the axial direction for the DBVP. The modified
e?t Cd.
0
dr dz.
R(y(r,z))r
PROCEDURE
AND
(25)
RESULTS
The nonlinear Fredholm integral equation, eq. (24), is solved using a Prime 400 computer at the Illinois Institute of Technology. The iteration scheme and some of the typical results are presented below.
Numerical
procedure
The modified Green’s function, cd, given by eqs (21) and (18) is calculated for a given value of the aspect ratio, c, using the well known Euler-McLauren summation technique (Dahlquist and Bjorck, 1974), which is essentially a repeated averaging method. The integrals in eqs (24) and (25) are obtained by using Gauss-Radau quadrature with an 8 x 8 point interior grid. Accuracy of the integration scheme is checked by using a 16 x 16 grid. The nonlinear integration equation, eq. (24), is solved by using the Newton-Kantorovich iteration scheme, described below. In operator form, eq. (24) may be written as F(y) The
= 0.
Newton-Kantorovich
(26)
iteration
scheme
is given
F(Yi)
by (27)
(22)
This transformation known as the “Superposition Method”, overcomes the difficulties involved in differentiating Gd on the boundary aR, as would be the case for the inhomogeneous boundary conditions. Substituting eq. (22) in eq. (4) and multiplying the
where, 6F, is the variation the present case is given
result with integrating
The Newton-Kantorovich iteration technique is quadratically convergent and is a superior iteration scheme to the linearly convergent successive substitution tnethod used by Dente et al. (1969a, b). When the iteration scheme satisfies a root mean
& and subtracting from eq. (19) x J, and with respect to r and z gives,
1 s 0
Y(r,,z)id(z;zo)dz
= 4’
1
1
ss0
Cl
x &(r,
R(1
-Y’(r,z))
z; rO, z,,)r dr dz.
Since &, and Gd are symmetric we may transpose the (re,ze) surface to the (r,z) surface. Thus
(23) from
1 s 0
Y(r,z,)i,(z;z,)dz, 1 X
ss0
= fz(l
-z)-@
I 0
R(Y
(re, ze))C% (r.z;
re,ze)re
dr0 dze.
(24)
Equation (24) gives the solution for the concentration profils inside the finite cylindrical pellet in a closed form.
of F with respect by
to JJand for
6F,,=$$
square
criteria
J[k
4 C (yi + 1
yij2] G (error)
-
m
the corresponding
effectiveness
Results For the DBVP
the rate
dimensionless
concentration, R(Y)
= Y” exP
factor,
expression y, would
YB(1 -Y) l-fB(l-Y)
q, is calculated.
in terms be,
of the
Integralmethod of analysis-1 where n is the reaction order and
E Y=RT, is the Arrhenius number. In Fig. 2, effectiveness factor vs Thiele modulus curves are given for a first order reaction (n = 1) in an isothermal (/I = 0) pellet. The parameter for the curves is the aspect ratio, E, in the range of 0.25-2.0. The analytical results of Gunn (1967) and the perturbation approximation of Ho and Hsiao (1977) are shown for comparison. It is noted that the numerical results are in good agreement with those of Gunn (1967) for Thiele modulus values less than 10.0. As the Thiele modulus increased above 10.0, the numerical calculations become less accurate, with a relative error of about 3 % at I$ = 40. This is due to the fact that for large values of c#, corresponding to severe internal diffusional limitations, the reaction zone is confined to a narrow shell
31
near the outer surface of the catalyst pellet. Such a problem with steep concentration gradients can be solved by using two regions of integration, with the interior of the pellet containing a relatively sparse grid. It must be noted that the range of Thiele modulus for which stiffness problems arise is within the asymptotic range where the effectiveness factor q a l/4_ It is shown (Aris, 197.5; Amundson and Luss, 1967) that, by modifying the definition of the Thiele modulus to include as the scaling parameter for length, the ratio volume/surface area of the pellet, the effectiveness factor curves for various E values can be brought together. This leads to the modified Thiele modulus, rh
+m =(2;E) In Fig. 3, the effectiveness factor versus modified Thiele modulus curves are given for various nonlinear reaction rate models. The results of the present study
2% l
0.25
f 4
Y:? 2.0
0 Gunn (1967) R&.(41 l HO ond Hsiao (1977) Ref.61
Fig. 2. Effectiveness factor vs Thiele modulus for first-order reaction-Dirichlet
_
0.01
*
SorenPen
problem.
g & (1973) Ref. (7)
0. I MODIFIED
THIEtJC MODULUS.
Fig. 3. Effectiveness factor vs modified Thiele modulus problem.
1.0 #‘m
ICI0
for nonlinear reaction rate forms-Dirichlet
32
SLJRYANARYANA
are compared with those from orthogonal collocation calculations (Sorensen et al., 1973). It is noted that our numerical calculations agree within 1%. Sorensen et al. (1973) calculated the bounds on the multiplicity analysis by using an initial value problem approach similar to that due to Dente et al. (1969a). The present study can be readily extended to obtain the corresponding Volterra integral equation formulation.
T
temperature
Ts Y
temperature of the bulk, K concentration of reactant in the pellet, dimensionless axial coordinate, dimensionless axial source coordinate, dimensionless
z 20 Greek P
CONCLUSIONS
An integral equation method is developed to model diffusion and reaction in a finite cylindrical catalyst pellet with no external transport limitations. The Green’s function kernel is modified to accelerate the convergence of the partial eigen series by two orders of magnitude. A Newton-Kantorovich iteration scheme is
used to solve the resulting nonlinear integral equation. Our numerical results for various reaction rate expressions compare well with previous results (Gunn, 1967; Ho and Hsiao 1977). It is noted that the modified Green’s function method exploits the symmetry properties of the classical Green’s function formulation while the convergence of the partial eigen series is accelerated by two orders of magnitude. Further, the convergence of the rate of quadratic Newton-Kantorovich technique is faster than the successive substitution. Acknowledgments-Thr general support and the assistantship (S.M.) provided by the Department of Chemical Engineering of Illinois Institute of Technology through Dr. D. T. Wasan and the Department of ChemicaI Engineering and Materials Science of Syracuse University are gratefully acknowledged. NOTATION
concentration of reactant in the pellet, g mol/cm3 concentration of reactant in the bulk, g mol/cm’ effective diffusivity of reactant in the porous catalyst pellet, cm’/s Green’s function for the finite cylindrical pellet for DBVP modified Green’s function for the finite cylindrical pellet for DBVP modified Besel function modified Bessel function effective thermal conductivity in the porous catalyst pellet, Cal/cm/K/s number of grid points in the radial direction length of the finite cylindrical catalyst pellet, cm number of grid points in the axial direction radial coordinate, dimensionless radial source coordinate, dimensionless dimensionless reaction rate reaction rate, g mol/cm’/s radius of the finite cylindrical catalyst pellet, cm catalytic surface area. cm’/n
al.
MUKKAVILLI~~
Y s -AH, E L ;: :, Pb
R an
of the pellet, K
letters
Prater temperature, dimensionless Arrhenius number, dimensionless dirac delta function heat of reaction, Cal/g mol aspect ratio of the finite cylindrical pellet, ROIL one dimensional Green’s function in the axial direction effectiveness factor, dimensionless temperature of the pellet, dimensionless Thiele modulus, dimensionless modified Thieie modulus, dimensionless bulk density of the pellet, g/cm3 interior of the pellet external surface of the pellet
REFERENCES
Amundson, N. R. and Schilson, R. E., 1961,
Intraparticle diffusion and conduction in porous catalysts-I. Single reactions. Chem. Engng Sci. 13, 226-236. Amundson, N. R. and Luss, D., 1967, On a conjecture of Aris: proof and remarks. A.I.Ch.E. J. 13, 759-761. Aris, R., 1975, The Mathrmaticul Theory of Diffusion and Reaction in Permeable Cutulysts, Vol. I. Calendron Press, Oxford. Baker, B. S., Gidaspow, D. and Wasan, D. T., 1971, Thermal phenomena in fuel cells and batteries. Advances in Electrochemistry and Electrochemical Engineering, Vol. 8, pp. 63-156. interscience Publishers, New York. Carslaw, H. S. and Jaeger, J. C., 1959, Conducrion ofHear in Solids. Calendron Press, Oxford. Dahlquist, G. and Bjiirck, A., 1974, Numerical Methods (translated by Anderson, C.). Prentice Hall, Englewood Cliffs, New Jersey. Derite, M., Biardi, G. and Ranzi, E., 1969a, Un nuovo methods per il calcolo dei profili di concentrazione a de1 fattore di utilivazione per reazioni complesse evolvenlise i granuli catalitici non isothermi. Quad. Ing. Chem. Ital. 5, 65-72.
Dente, M., Biardi, G. and Ranzi, G., 1969b, Problemi riguardantila determinazionediretta dei limiti statici di
‘ignizione’ e ‘spegnimento’ per reazioni cineticamente semplici evolventici in granuli catalitici non isothermi. Quad. Ing. Chem. Ital. 5, 122-129. Dixit, R. S. and Tavlarides, L. L., 1982, Integral method of analysis of Fischer-Tropsch synthesis reactions in a catalyst pellet. Chem. Engng Sci. 37, 539-544. Gunn, D. J., 1967, Diffusion and chemical reaction in catalysis and absorption. Chem. Engng Sci. 22, 1439-1455. Ho, T. C. and Hsiao, G. C., 1977, Estimation of the effectiveness factor for a cylindrical catalyst support: a singular perturbation approach. Chem. Engng Sci. 32, 63-66. Kesten, A. S., 1969, An integral equation method for evaluating the effects of film and pore diffusion of heat and mass on reaction rates in porous catalyst particles. A.I.Ch.E. J 15, 128-131. Mukkavilli, S., Tavlarides, L. L. and Wittmann, C. V., 1987, _. . . . ^ ..^ -.Integral method ot analysis tor chemtcal reaction in a
Integral method of analysis-1
33
nonisothermal finite cylindrical catalyst pellet-II. Robin problem. Chem. Engng Sci. 42, 3540. Sorensen, J. Q., Guertin, E. W. and Stewart, W. E., 1973, Computational models for cylindrical catalyst particles. A.I.Ch.E.
J. 19, 969-975.
Tolstov, G. P., 1962, Fourier Series (translated by Silverman, R. A.). Dover Publ., New York. Consider
APPENDIX
Consider f(x)
??f$?&-
such that,
f;(x)
7
dx.
(A-3)
f, and fz such that,
fi(O)=JX)=O
0 < x < a;
f(x)= {h(x) a=sx<
1,
(A-1)
and define the transform-inverse pair J =2JiJ(x)sinnnxdx, J(x) =
s
and
fi (a) =X(a)-
(A-4)
Then, the first two terms on the right-hand side of eq. (A-3) cancel each other and As converge as I/I?. If in addition. (A-2a)
F f. sin nrtx. (A-2b) n=1 Using the product integration and assuming thatfr andh are continuously differentiable in the respective intervals, we have from (A-Za),
f;(o)
=fi(l)
= 0
and
f;(a)
=&(a)
(A-5)
then, f. s converge as 1/n3. Thus we see that in addition to continuity of the function (and its derivatives) we require that the function (and its derivatives) also satisfy the homogeneous boundary conditions in order to accelerate the convergence of the eigen series.