64
Ind. Eng. Chem. Fundam.
0 = permeation time lag, s 7
= D t / L 2 ,dimensionless
Literature Cited Ash, R.; Barrer, R. M.; Nicholson, D. 2.Phys. Chem. 1083, 37, 257-288. Astarita. G.; Savage, D. W.; Blsio, A. "Gas Treating with Chemical Solvents": Wiley-Interscience: New York, 1983. Danckwerts, P. V.; McNeil, K. M. Trans. Insf. Chem. Eng. 1887, 4 5 , T32T49. Folkner, C. A.; Noble, R. D. J . Membr. Sci. 1083. 12. 303-312 Frisch, H. L. J . Phys. Chem. 1957, 61, 93-95. Frisch, H. L. Po&". Eng. Sci. 1080, 20, 2-13. Goddard, J. D.; Schultz, J. S.:Suchdeo, S. R. AIChEJ. 1074, 2 0 , 625-645. Goodknight, R. C.;Fan, I. J . Phys. Chem. 1081, 65,1709-1712. &os, G.:Moll, W.; Hoppe. H.; &os, H. J . Gen. Physiol. 1076, 6 7 , 773-790. Henls, J. M. S.;Tripodi, M. K. Science 1083, 220, 11-17. Lonsdale, H. K. J . Membr. Sci. 1082, 10, 81-181. Matson, S. L.; Herrick, C. S.;Ward, W. J. Ind. Eng. Chem. Process Res. D e v . 1077. 16. 370-374. Matson, S. L.;Lopez, J.; Quinn, J. A. Chem. Eng. Scl. 1983, 3 8 , 503-524. Meldon, J. H.; DeKoning, J.; Stroeve. P. Bioelect. Bioenerg. 1078, 5 , 77-87. Meldon, J. H.; Kang, Y. S. AIChE Symp. Ser. 1083, 79, No. 227, 38-42.
1985,2 4 , 64-77
Meldon, J. H.: Smith, K. A.; Conon, C. K. Chem. Eng. Sci. 1977, 3 2 , 939-950. Michaels, A. S.; Bixler, H. J.: b i n . H. L. J . ADD/. Phvs. 1084.. 3 5 .. ,. 3185-3178. Paul, D. R. J . PO/ym. Sci. A2 1080, 7 , 1811-1818. Petropoulos, J. H.; Roussls, P. P. In "Organic Solid State Chemistry"; Adler. G., Ed.; Gordon and Breach: New York, 1969; pp 343-357. Sasidhar, V.; Ruckenstein, E. J . Membr. Sci. 1083, 13, 67-84. Tajar, J. G.; Miller, I. F. A1Ch.E J . 1072, 18, 78-83. Tshudy, J. A.; Von Frankenberg, C. J . Po/ym. Sci., Polym. Phys. 1973, 1 1 , 2027-2037. Way, J. D.; Noble, R. D.; Flynn, T. M.; Sloan, E. D. J . Membr. Sci. 1082, 12, 239-259. Weber, 0.W.; Miller, 1. F.; Gregor, H. P. AIChE J . 1070, 16, 609-614. Wolfe, R. G.; Miller, I. F.: Gregor, H. P. J . Biomed. Mater. Res. 1970, 4 , 295-3 12.
Received for review December 5 , 1983 Accepted August 8, 1984
This work was supported by research grant No. DAAG29-83-K0137 from the U.S. Army Research Office.
Solution Methods for Problems with Discontinuous Boundary Conditions in Heat Conduction and Diffusion with Reaction Patrick L. Mills Chemical Reaction Systems Group, Corporate Research Laboratory, Monsanto Company, St. Louis, Missouri 63 167
Steven S. La1 and Mllorad P. DudukovlC" Chemical Reaction Engineering Laboratory, Department of Chemical Engineering, Washington Universiv, St. Louis, Missouri 63130
Boundary value problems that arise in heat conduction and diffusion with reaction involving linear differential operators and discontinuous boundary conditions along a smooth boundary are considered. I t Is shown that problems of this type can be effectively solved by a number of techniques such as those based upon the method of weighted residuals, integral equation methods, and the finite element method. Comparisons among results for a number of test cases show that the integral equation method is superior to others since it requires solution of a much smaller system of linear equations whlle still achieving more accurate results. The example problems treated yield the effectiveness factor for a discretely distributed catalyst on a nonadsorbing support, the correction factor needed for mass transfer Coefficients for transport to nonuniform surfaces, and the corrected effectiveness factor for spherical pores such as that encountered in zeolites.
Introduction A number of problems in transport phenomena are associated with discontinuous boundary conditions, i.e., with the situation when there is an abrupt change in the type of boundary condition a t a point or line along a smooth boundary curve or surface. These can be called mixed or split boundary value problems (Mills and DudukoviE, 1980, 1982) which are to be distinguished from problems with mixed or Robin (third kind) boundary conditions. The variety of chemical engineering situations that lead to problems with discontinuous boundary conditions is well represented by a number that have been examined recently and which include: (i) diffusion and reaction in porous catalyst pellets which have partial external and/or internal wetted surfaces (Goto et al., 1981; Herskowitz et al., 1979; Martinez et al., 1981; Mills and DudukoviE, 1979, 1980b; Ramachandran and Smith, 1979; Tan and Smith, 1980); (ii) diffusion of reactants to a surface with nonuniform
(patchy) catalyst distribution and surface reaction (Loffler and Schmidt, 1975); (iii) determination of the correction factor for diffusion coefficients measured for spherical catalyst pellets imbedded in a tube by monitoring the concentration response at steady state (Meyer et al. 1976) or at transient conditions (DudukoviE et al., 1981); (iv) catalyst effectiveness of irregularly shaped pores (Chu and Chon, 1970); (v) analysis of conduction heat transfer between solids which have partial contact along a common boundary (Dundurs and Panek, 1976; Sadhal, 1980,1981; Schneider et al., 1977);and (vi) fluid mechanical modeling of a particular type of Stokes flow phenomena (Takagi, 1974). Problems of this type arise in other branches of applied science (acoustics, crack propagation) and various methods of solution have been investigated. The few analytical solutions that are available have been reported in a monograph by Sneddon (1966) and classical papers by Tranter 0 1985 American
Chemical Society
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985
(1959, 1964) and Collins (19621, to name a few. The difficulties in solving nonstandard problems of this type and the need for accurate and reliable solution algorithms has been examined by Kelman (1967a,b, 1970, 1979, 1980), Kelman and Koper (1973), Kelman and Simpson (1975), and Kelman et al. (1977), who have made significant contributions in this field. In contrast, the chemical engineering literature has ignored, until recently, the computational aspects of this problem. For example, when discontinuous boundary conditions arose in the mathematical description of a chemical engineering problem, they were sometimes handled by standard techniques such as finite differences, as illustrated by Loffler and Schmidt (1975), or the Galerkin method, applied by Chu and Chon (1970). In these cases, no significant investigation regarding the accuracy of the method was reported. On other occasions, the discontinuous conditions which naturally arose during modeling of the system were eliminated by modifying the system geometry, e.g., from a sphere to a cube to yield a simpler problem (Herskowitz et al., 1979). In a related case, Herskowitz et al. (1979), in order to reconcile their data with theory, estimated that the mass-transfer coefficient for exchange between pore mouths and the bulk fluid for pellets of porosity e is equal to the mass transfer coefficient for exchange from the whole external surface (obtained by a dissolution method) divided by the pellet's porosity. Questions can be raised regarding both the accuracy of the numerical methods used and the validity of the analytical or heuristic approximations. We have reported previously (Mills and DudukoviE, 1980) that finite difference schemes (Schneideret al., 1977; Goto et al., 1981; Loffler and Schmidt, 1975; Martinez et al., 1981), the artifical interface method (DudukoviE and Mills, 1978; Kelman, 1967a; 1967b; 1970),and the method of weighted residuals (Kelman and Koper, 1973; Mills and DudukoviE, 1980,1982) have been used for solution of split boundary value problems. It has been established that finite difference schemes often require an extremely fine mesch size, are frequently inaccurate in the vicinity of the boundary singularity (Whiteman, 1968) and are, in general, much more inefficient than other techniques. The artificial interface method is readily applied only for certain problems posed in Cartesian or polar coordinates, which limits its usefulness. In contrast, the method of weighted residuals (least squares, Galerkin, and collocation) provides a systematic approach to solution of split boundary value problems in any geometry, and appropriate general formulas and computer software packages have been developed (Mills and DudukoviE, 1980, 1982). Convergence theorems for the least-squares approach exist (Feinerman and Kelman, 1974; Feinerman et al., 1974) and extensive numerical testing has been done with the other MWRs to suggest that they are the most reliable technique available today. Nevertheless, for certain values of the boundary condition discontinuity, MWRs still require the solution of a large set (N I 150) of linear equations in order to achieve desired accuracy and, hence, can only be used on large computers at the expense of substantial computational time. Having demonstrated that problems with discontinuous boundary conditions are ubiquitous and important in chemical engineering, it is our goal to determine accurate yet readily useable methods of solution which preferably do not require extensive computations or large computers to obtain reliable results. We surmise that if such methods are available, problems of this type would occur more frequently rather than being handled by engineering es-
65
timates and approximations of dubious accuracy. We submit that anyone with sufficient stamina, a large computer at his disposal, and some rudimentary knowledge of algebra and calculus can employ a brute-force approach and obtain a finite difference solution of split boundary value problems with linear operators. This, however, is not necessarily good engineering practice and we would like to explore and present the best methods for solution of problems of this type. As might be anticipated, we will show that there is a definite trade-off. To achieve computational efficiency, more effort involving analytical methods is necessary before using the computer. In the integral equation method, for example, simpler computational schemes are obtained, which results in less computer time and storage requirements. We will present, on the basis of selected examples, each of which has a physical significance, a comparison between the method of weighted residuals, the integral equation approach, and the finite element method for solution of split boundary value problems. While the development of a reliable method for solution of split boundary value problems is one of our objectives, the results obtained in solving various example problems are of equal importance. These example problems provide us with information on how to estimate the correction factor for mass transfer coefficients at partially adsorbing surfaces and how to calculate the effectiveness fador for heterogeneous catalytic surfaces, both for metal dispersions on supports and for spherically shaped pores such as that found in zeolite type structures. Problem Formulation and Solution Methods The example problems to be treated later in this paper are based generally on the following two-dimensional transport equation V 2 u- c2u = 0 in D
(1)
where D is a region enclosed by a sectionally smooth boundary nd
a 0 = CdDi i=I
Generally, discontinuous boundary conditions occur along one of the boundaries corresponding to i = s. Normally, there are one or two discontinuities so that this particular boundary dD, can be further divided into boundary segments dDSi m
dD, = CdD,,; i=l
m =2 or 3
(2)
The boundary condition is of the same type along each dosi and can be given by au Xiu + vi- = wi on dosi;i = 1, ..., m (3) an Equations of type 3 naturally apply to all other smooth boundaries, dDi, for i = 1, ..., n d , and i # s. Equation 1is reduced to Laplace's equation for c = 0, which describes steady-state diffusion of heat or mass in an isotropic medium. The case c > 0 leads to the Emden-Fowler equation with f ( ~ =) u (Mehta and Ark, 1971) and is encountered in the theory of diffusion with reaction where c is the Thiele modulus or its equivalent. The quantity u is either dimensionless concentration or temperature. Parameters Xi,vi, and aiare all related to appropriate physical quantities such as conductivity, heat transfer coefficients or diffusivity, mass transfer coefficients, and rate constants and are defined on (0, a). Appropriate rearrangement of these parameters leads to key dimensionless groups such as the Biot number and the
66
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985
surface Damkoehler number. Their origin and physical significance in heat conduction and diffusion with reaction are summarized well by Ozisik (1968) and Aris (1975), respectively. The Helmholtz equation is obtained for imaginary values of c where Icl denotes the wave number. This case will not be treated here since it is discussed in detail elsewhere (Collins, 1962). A. Method of Weighted Residuals. The starting point as shown earlier (Mills and Dudukovit, 1980, 1982) has been to utilize the separation-of-variables solution for eq 1 using all the boundary conditions other than those on dD, where the split conditions occur. This allows representation of the unknown function u by the following infinite series
B. Integral Equation Method. Our success in obtaining efficient solutions to a particular split boundary value problem for the effectiveness factor for a partially externally wetted catalyst pellet (Mills and DudukoviE, 1982b)by an integral equation formulation led us to generalize this approach to all problems of the type given by eq 5. Since the general case is treated elsewhere (Mills and Dudukovit, 1984), we restrict our attention here only to dual series problems ( m = 2) which correspond to the examples to be treated in this paper. Following Sneddon (19661, we extend one of the series in eq 5, for i = j = 1or 2 and defied for .$-l i t 5 t j , into the remaining region k by introducing an unknown function &(E) in tj i 5 I tj+l m
m
C a n ~ n j r n ( t )=
(4)
where X n ( g )and I’n([) are the eigenfunctions associated with the separation-of-variables solution in the g and 5 coordinates, respectively, while a, are the unknown coefficients to be determined from the split boundary conditions on aD,. Application of these conditions yields the following dual ( m = 2) or triple ( m = 3) series equations m
C a , ~ ~ i r ~= (fi(t); t ) t i - 1 < t < 5; for i = 1,..., m
n=O
n=O
= &!k(t); t k - 1 4 k=l,2; k # j
...,
t
5
tj
(7)
tk
Formation of the inner product ( w ( t ) r n ( [ ) r,([)) , over (0, 1)followed by interchanging the order of summation and integration leads to the following expression for the series coefficients a,
k#j
In this equation, the weighting function w(4) assures orthogonality of the r,(f)’sover (0, 1). Using the orthogonality property and substituting eq 8 for an’s into eq 7 results in the following Fredholm integral equation of the first kind dt
&>kj(t,r)&!k(t)
(6)
General expressions for the A j ’ s and Fj’s of each MWR method were developed by Mills and Dudukovit (1980, 1982). A computer software package which allows the user to select a particular method of weighted residuals for a specified kernel and series modifiers is also available (Mills and DudukoviE, 1980; 1982). In order to solve a particular linear operator problem with discontinuous boundary conditions by the MWR technique, the user must first reduce the problem to the standard form given by eq 5. Next, the user must specify the series modifiers pni, the series kernel [sin A,,& cos A,(, or P,(cos e), where A, = n, n + l f 2 , etc.], the boundary forcing functions fi(t),and the size of the matrix N. In all cases, the infinite set of eq 6 is reduced to a finite set for j = 0, 1, ...,N - 1, which can be solved by Gaussian elimination or LU decomposition to yield the first N series coefficients aj’s. Our computer package is used for solution of some of the example problems presented here. Although the above described algorithm is a powerful software package for solution of split boundary value problems, there are several drawbacks to the method. General error estimates are not available, convergence can be quite slow in some problems, and solution to a large number of equations may be required.
Fk(r);
t’k-1
< 5’ < 5 / k , k
# j
(9) where Fk(5’)
m
= Fi; i = 0, 1, 2,
5
(5)
In these equations, the pni are the so-called series modifiers (Kelman and Koper, 1973), r,(t)is the series kernel, fi(t) are functions prescribed by boundary conditions, and a, are the unknown series coefficients. Appropriate dimensionalization of [ ensures that to= 0 and E, = 1. The series kernel r,([)is orthogonal over (0,l) with respect to some weighting function but not over ti). The working equations for the method of weighted residuals have been developed (Mills and DudukoviE, 1980; 1982) by formation of the residual from eq 5 followed by application of an appropriate inner product of the residual and selected trial functions. In all cases (least-squares, Galerkin, collocation, and the method of moments) an infinite system of algebraic linear equations is obtained for the a j ’ s CajAij j=O
fj(t); t j - 1
(10)
= f k ( r ) - hkj(t’)
The integral equation kernel & j ( [ , r ) and the forcing function h k j ( r ) are given by the following infinite series which contain the series modifiers pni and the series kernel
rn(t) Kkj(t,P)
=
2[
n=O
s,arn(t)r,(E’)/PnjSE*
U(f)r,%)
dt]
EO
(11)
Closed form inverses to the integral equation given by eq 9 are seldom possible for the kernel and forcing functions usually encountered in chemical engineering problems so that numerical inversion is necessary. Use of stretched coordinates x = ( E - tk-l)/(tk - t k - 1 ) allows standard quadrature methods over the interval (0, 1) to be used. Trapezoidal, Simpson’s, and Gaussian quadrature rules have been tested and found to give comparable results (Mills and Dudukovic, 1984). In chemical engineering transport and reaction problems, the kernal frequently contains a logarithmic singularity at x = y and a modified quadrature method (Baker, 1977; Kantorovich and Krylov, 1964) must be used. This leads to the following set of N linear equations for the values of the unknown function g k ( t m ) evaluated at the quadrature rule abscissas, tm
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985 67
where within this element, u ( ~is) ,given as a linear combination of the nodal values, ul k
u ( ~=) CNiul = NTu
(16)
l=i
where N T = ( N i N j N k )and uT = (uiujuk).The elements of the vector N assume the following form for a triangular element
Ni = (ai + bit
+ e17)/2Ae
(17)
where (184 - .5+251+1;~ = i, j , k (18b) bi = 1i+1 - 711+2; 1 = i, j , k (18c) ci = b+z - &+I; 1 = i, j , k with the convention that j = i + 1 and k = i + 2. The area a1 = E1+1711+2
Cases where the kernel does not contain a pole for f = may be solved by the usual quadrature method (Baker, 1977) so that eq 14a and 15 are omitted, with eq 14b now being valid for all 1 and m. Solution of the linear system defined by eq 13 may be performed by standard techniques such as Gaussian elimination or LU decomposition. Various mathematical software packages such as the HSL or IMSL library contain excellent subroutines for this purpose. As shown above, the integral equation method allows the problem to be reduced to the solution of a set of linear algebraic equations. We will show here that this set is much smaller than the one required by the method of weighted residuals. The solution of eq 13 seems simple, and yet it should be pointed out that it assumes one is able to sum the kernel accurately, thereby providing the elements of matrix B with 6 to 7 digit accuracy. Various studies have shown (Baker, 1977; SIAM Conference, 1981) that the condition numbers for matrix B (Le., the ratio of the largest to smallest eigenvalue) are adversely affected when the matrix elements Blm’s contain errors or are known with low accuracy. In many chemical engineering problems, the integral equation kernel is normally given as a slowly convergent infinite series. Early numerical experiments (Mills and DudukoviE, 1979; 1982b) attempted to use direct numerical summation of the kernel and confirmed that the required accuracy was not achieved. However, in each problem encountered so far, an equivalent closed form expression for the kernel can be developed which provides sufficient accuracy for the Blm’sand makes the integral equation method the most powerful. Such closed form expressions for the kernel series are given for the example problems in this paper. The detailed transformations leading from the infinite series to closed form representation of the kernel are given elsewhere (Mills and DudukoviE, 1982b, 1984). C. The Finite Element Method. The advantage over other computational techniques of the finite element method for the solution of boundary value problems with irregular boundaries or discontinuous boundary conditions has been exploited for quite some time in civil engineering, particularly in structural analysis. The method has recently been introduced for solution of chemical engineering problems by Finlayson (1980). An abbreviated outline of the method as applied to mixed boundary value problems is given here. The reader interested in further details is referred to the excellent and thorough user-oriented monograph of Rao (1982). For simplicity and illustration purposes, the methodology for a two-dimensional Cartesian coordinate system with discontinuous boundary conditions is given. The domain of interest, D, is divided into NE triangular elements. The three nodes of our general mth element are , (tk, q k ) , relabeled i, j , k at positions (&, ai), (.&, ~ j ) and spectively. The value of the dependent variable every-
of the element is given by 1 A(e) = - E j ) ( T j - 71k) - ( E j - E k ) ( V i (19) 2 Upon using the Galerkin method, one requires that the residual of eq 1be orthogonal to each trial function Nifor the whole area of the element which, with the use of Gauss’s theorem, results in eq 20. Substitution of eq 16
-I(&
?*))
and its derivatives allows eq 20 to be reduced further to a set of linear equations of the form K(e) = p ( e ) (21) where the matrix K(e)is given by the following sum K(e) = Kl(d + &(e) + &(e) (22) The matrix Kl(e)results from term I in eq 21 and is given by Kl(e)= JBTDB dA
1
where the matrix B is given by
The diagonal matrix D is the identify matrix I when V2 is the conventional Laplace operator, but if one deals with scaled problems or anisotropic media and interprets
then the matrix D becomes D=(:
:)
(25)
Term I1 in eq 20 is zero if the element is completely within region D. This term is also assigned a zero value for those elements on the boundary of D at which either au/an = 0 (oi= Xi = 0) or u is specified. In the latter case, i.e., for Dirichlet conditions, the appropriate nodal values of u are specified. If the general boundary condition given by eq 3 applies for the portion of the element on the
68
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985
Table I. Elements of Matrices and Vectors for Eq 22 bEbj2 f a+i2 atbibj
edge
condition
+ a+jcj
aEbjbk t C U V C ~ C J
p(e) :
ij
Nh =
T = Ts
Sij = length of ij
jk
I C
0
0
c+
E
X*
Figure 1. Heat conduction in a semiinfinite plate. The lower edge is maintained at T = T , for 0 5 z* < c* and is insulated for c* < x*
L-2
5 1.
hri = 0
( y s)j k2[ p 11
3
0 2 1 -(-)-
Si, = length of j k
boundary Ci of region D with nonzero values of Xi, vi, wi, then term I1 of eq 20 generates the matrix K2(e)and the vector P(e)
Finally, matrix K3(e)originates from term I11 of eq 20 and has the form
K 3 ( e ) = c2 s g N N T dA For the trial function given by eq 17, the elements of K and P can be readily calculated following the procedure given by Rao (1982). These are summarized in Table I for ease of reference. The above procedure yields a local set of equations for the values of the unknown function at the nodes of a particular element. A n appropriate bookkeeping procedure for numbering of nodes (Rao, 1982) is handled in a subroutine (Lai, 1983) and allows the system of equations to be augmented for all the N N nodes in domain D to give the following system of linear equations KNNXNNUNN = PNN (29) The computer algorithm further reduces the system of equations by the number of specified nodal values, ND, due to specified Dirichlet conditions at some boundaries and solves the resulting set of linear equations by the Choleski decomposition method for symmetric matrices. Details for the example problems treated here are given by Lai (1983). The above development illustrates that a standard application of the finite element method reduces the discontinuous or split boundary value problem to a set of linear equations. Once a general node numbering routine, a subroutine for evaluation of elements of K and P, and a linear equation solver are available (Lai, 1983), then the user need only supply the number of elements desired, NE, and the appropriate boundary conditions and their locations. The finite element method can then be used as readily as the pre-programmed version of MWRs treated earlier. The question arises as to how well the results
obtained compare to those calculated by either the MWR's and by the integral equation method and which seems to be the simplest and most reliable method to use. This question is addressed by introduction of various example problems given below. Numerical Examples Several examples involving heat conduction and diffusion reaction systems with discontinuous boundary conditions are given below to illustrate application of the above proposed methods for solution of split boundary value problems. Example 1. Heat Conduction in a Partially Insulated, Semiinfinite Plate. This example is selected to compare the integral equation and MWR solutions to the available analytical solution. The physical problem is that of two-dimensional, steady-state heat conduction in a semi-infinite plane whose lower edge y = 0 is maintained at a prescribed temperature T = T, for 0 I x* < c* but is insulated for c* < x* I 1 as shown in Figure 1. The problem is described by eq 1in Cartesian coordinates with c = 0, i.e.
with boundary conditions Io , t = 0; aula( =
o
7=O,u=l; O I ( < C ?I = 0, au/ag = 0; c < 5 I
(31) (32) (334 (33b)
lim(u)= 0
(34)
ql0,[=7r; u = o
r)--
The dimensionless temperature is denoted here by u ( T - T,)/(T, - T J , while the dimensionless coordinates are = 7 r x * / l , q = ry*/l. The boundary discontinuity is c = ac*/l. The remaining variables are defined in the Nomenclature section. A separation-of-variables solution which satisfies eq 30 and the boundary conditions given by eq 31, 32, and 34 is m
u(l,q) =
n=0
an(n+ 72)-1e-(n+1/2)~cos [ ( n+
Y&1
(35)
Application of the discontinuous boundary conditions given by eq 33a and 33b leads to the following dual-series equations m
c a,(n +
y2)-l cos
n=O
[ ( n+
y2)[] = 1;
0 IE
1
(52b)
72
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985
w(S) = 1
(55)
The infinite series representation of the kernel and the forcing function obtained by substituting eq 51-55 into eq 11 and 12 are
KIZ(E,I')=
-
"( + +-) + W
1
tanh ( n d / w ) cos (n7rE) cos (n7r.$') n7r n=1
2c
Da
(i, 2 cos (n7rE)cos w
+ N
(n7rQ) (56a)
= 20 QUADRATURE POINTS
10-4-
I
10-1
I
100
1 ' 1
,
101
,
I
l ' l l ' l l
102
'
I
' I '
If
103
DAMKOEHLER NUMBER, Da
Figure 4. Effectiveness factor as a function of Damkoehler number at various plw and 6 f w values.
-
All the above expressions given by eq 48-56b can be obtained when the limit 4 0 is applied to the formulas developed for the general case of 4 # 0 by Mills and DudukoviE (1984). The form of the integral equation kernel and the forcing function given by eq 56a and 56b, respectively, cannot be directly summed with sufficient accuracy on any existing digital computer. Transformations of the kernel and the forcing function that are necessary to yield a readily summable form are shown in detail for the general case of 4 # 0 by Mills and DudukoviE (1984). Here, we summarize the fiial representations for our case of 4 = 0 for the kernel K&,F), forcing function F 1 ( r ) ,and the integrated kernel A ( f ) ,all of which can be readily evaluated with high accuracy. 6 1 1 Klz(E,f')= ; -; In lcos (TE) - cos (7rr)l- - In 2 + -
tanh
($)
7r
-1
cos (n7rS)cos (n7rp)
2 2
n7r
n=l
+
The function denoted by $ in the above expressions refers to Clausen's integral (1832) defined by sin (ne 0 IO I27r $(e) = -; n=l nz
c m
=
-Lo
In (2 sin
i)
d.2
(59)
This function was evaluated in this work by use of the series expansion in terms of Bernoulli numbers B2k given by Ashour and Sabri (1956). Calculated values for $(e) agreed with those given by Clausen to 16 decimals on our IBM 370165 using the first 17 terms with double precision arithmetic. For the case where mass transfer through the boundary layer is controlling, the rate can be obtained by letting Da m in eq 50-58. This corresponds to mass transfer to a surface that has a fraction plw that is permeable to the diffusing species, while the remaining fraction 1 - ( p / w ) is impermeable. The quantity of interest is the ratio of the observed mass-transfer coefficient k , to the conventional mass transfer coefficient based on film theory k , = D / 6 , which can be obtained by applying the following limiting process to eq 49 N = lim {q.Da) = k , / k , = 1 - a o ( 6 / w ) (60)
-
Da--.;
W
. 1J J
(n7r)"
n=l
(57th
n=l
'I$( f
(nx)" 7r2
+ I ' ) + $[ T (
f - I ' ) ] } (58)
The results for example problem 2 are computed by three methods: (i) using the MWR's where N linear equations of the type given by eq 6 are solved; (ii) using the integral equation method and solving N linear eq 13 at N Gaussian quadrature points with KI2(E,F)and F l ( p ) , and A(,$') given by eq 57a, 57b, and 58; and (iii) by the finite element method solving linear eq 29 for NN nodal values with cq = ( 6 / w ) z and CY? = 1 in the expression for Kl(e),( - X / u ) = Da in the expression for Kz(e)and P = 0, K3(e)= 0 for 4 = 0 all in Table I. In Figure 4 the catalyst effectiveness factor is plotted against the Damkoehler number at different plw and 6/w values. The results presented in Figure 4 agree in a semiquantitative sense with those of Loffler and Schmidt (1975). An absolute comparison is not readily performed because these authors did not present their values in tabular form. Based on our numerical experience with all other methods and other related problems, a very fine mesh size and a large number of equations would be necessary to obtain finite-difference results which are
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985 73 10= 1 - --pw=oe -pw
OB--
-p'w
= 02
1
- - - p / w = 0002 PARAMETER
06
_ e - -
b
w
F
DAMKOEHLER NUMBER,
Da
Figure 5. The observable TDa as a function of the Damkoehler number, Da. I N T E G R A L E Q U A T l O N ( N = 20) --
a
-INTEGRAL
EQUATION (N = 20)
-
= 0.1
8iW
10-4
...FINITE ELEMENT (NN = 1681)
0.0
0.2
0.6
0.4
0.8
0
'
PMI
Figure 7. Effect of the finite element step size on the calculated effectiveness factor results.
I
/
......... flm ELEUEWT ph = 0.75
......... I
0.0 I
0.0
0.2
0.6
0.4
0.0
1.0
PIW
......
b 0.010
p/w = ..................................... ..... .......................... ................. ... ..................... 0.05
0.1 FINITE ELEMENT (NN = 1681) - LEAST-SQUARE (N = 150) 8IW = 5
00
02
04
00
--
-08
0
10
PIW Figure 6. Comparison between effectiveness factor values calculated by various methods as a function of plw. (a) Da = 1, 10; (b) Da =
2 4 6 8 10 12 14 16 18 20 NUMBER OF GAUSSIAN QUADRATURE POINTS, N
I
l
0
15
I
0
l
,
l
l
l
l
,
,
,
l
30 45 60 75 90 105 120 135 150 165 LEAST SQUARE METHOD MATRIX SIZE, N
' do '
1280
1920
2580
3200
NUMBER OF ELEMENTS, NE
100.
Figure 8. Comparison between the integral equation, least-squares, and finite element methods for evaluation of 01 as a function of the order-of-approximation.
comparable to those from the integral equation. Figure 4 quantifies the intuitively expected results which show that at fixed Damkoehler number (fixed ratio of diffusion and surface reaction resistance) and fixed boundary layer thickness (6/w),the larger fraction of the surface is catalytically active (Le., the larger p/w) the higher the effectiveness factor. As expected, at Da < 1 and p / w < 0.5, qo plw. At fixed Da and plw, decreasing values for 6/w result in a subsequent decrease in the effectiveness factor. Figure 5 presents the above findings in a slightly different form. Here, the Damkoehler number is plotted against the product qDa which represents the ratio of the observed ration rate and maximum rate for mass transfer with p / w and 6/w as parameters. Figure 6 shows a comparison between the results for the effectiveness factor vs. plw obtained by the integral
equation, least-squares, and finite element methods. The agreement among all methods is good except at very small values of plw, which was discussed in more detail elsewhere (Mills and DudukoviE, 1984). However, a large system of equations ( N = 150) for least-squares and NN = 1681nodes have to be used for the finite element method to obtain accuracy comparable to that of the integral equation method with only 20 quadrature points. The approach of the finite element results to the integral equation results is shown in Figure 7. Figure 8 gives a comparison between mass transfer coefficient correction factors calculated by the integral equation, least-squares, and finite element methods as a function of the number of Gaussian quadrature points N , least-squares matrix sizes N , and number of nodes N N , respectively. The number of equations needed for the
F;:
74
Ind. Eng. Chem. Fundam., Vol. 24, No. 1, 1985 i o
r"
/ o w = i 0
09
08 07
eu
a 06
a
2
05
/ l W
I
= 01
9
Ea
E
Figure 10. Illustration of a sperical ink-bottle type pore.
o4 03
INTEGRAL EQUATION (N = 20)LEAST SQUARE (N = 150)
02
FINITE ELEMENT("
= 1681)
___._
01
00 00
02
0.6
04
08
10
PW
Figure 9. Mass transfer correction factor as a function of the fractional active area for various numerical methods.
finite element method is NN - int(p/A) where A is the spacing selected. The aspect ratio 6 / w is fixed at 0.1 in Figure 8, while p l w is varied from 0.05 to 0.75. An analogous behavior is observed a t other 6 / w values. Correction factors a are calculated from eq 60. When the coefficient a. is evaluated by the integral equation method, the correction factors are essentially constant at each p / w value when four or more quadrature points are used ( N > 4). Those calculated by the least-squares method keep increasing, even for matrix sizes of 100 or more. Similar behavior is observed for the finite element method. Figure 9 shows that the only significant deviation between the integral equation method (N = 20), least-squares method (N = 150), and the finite element method (NN= 1681) occurs a t very low p / w (plw < 0.2) and high 6Jw ( 6 / w 1 10). Based on our experience with example problem 1,this is to be expected. By analogy, we expect the integral equation results to yield the best approximation to the true answer, although an absolute proof of that is lacking. The interpretation of Figure 9 is of particular importance. The parameter a represents the ratio of the actual mass transfer coefficient k , for a surface having an active fraction p / w that participates in the mass transfer process to the mass transfer coefficient for a uniformly accessible surface k,". Based on two-film theory, k," = D / 6 . Referring to Figure 3, we see that if the boundary y = 6 is maintained at c = c b while the fraction of the surface for 0 < x < p / w at y = 0 is at c = 0 and the rest of the surface for p / w < x I 1impenetrable to the diffusing species, then the flux per unit area of total surface is kmcb while the flux per unit area of active surface is k,,Cb(p/W). The flux to the uniform surface is kmoCb. The ratio of the two fluxes is a = k , / k m o . Figure 9 demonstrates that the two-dimensional model predicts k , = k," at 6 / w > 50 and p / u > 0.2. Thus, when 6 f p = 0 (lo3),i.e., the boundary layer thickness is much larger than the size of the discontinuities at the surface, the mass flux per unit area of total surface remains the same as for a uniform surface, k , = k," and the flux per unit area of active surface is enhanced by an amount k,a = k," (wlp). At the other extreme, when the boundary layer thickness is small compared to discontinuities at the surface and 6 / w = 0.1, the flux on the total surface is diminished almost in proportion to the fraction
of active surface so that k , = ( p / w ) k , " and the flux on the active area is not augmented compared to the uniform surface so that k , = k,". Only under this condition does the assertion of Herskowitz et al. (1979) seem to be valid. Example 3. Effectiveness Factor in Spherical Ink-Bottle Type Catalyst Pores. Here we consider a spherical ink-bottle type pore whose surface is maintained at c = c , for 0 5 8 < a but undergoes first-order reaction rs = kscs for a < 8 I T at r = R on the walls of the pore. A diagram of the pore is given in Figure 10 which can be used to represent an idealization of a cavity frequently found in zeolite or related catalysts. This problem was first treated by Chu and Chon (1970), who evaluated the pore effectiveness factor by the Galerkin method. The objective here is to compare their results with those obtained by the various methods presented here. The following differential equations describe the system
(
p2$)
+i a (sin
8
$) =O
(61)
sin 8 88
aP
1 au
p = 1; --
D~ ap
+ u =O;
au
8=o;T-=o; a0
LY
