Machine Computation of Transients in Fixed Beds with Intraparticle Diffusion and Nonlinear Kinetics James M. Wheeler Eastman Kodak Co., Rochester, N . Y . 14650
Stanley Middleman' University of Massachusetts, Amherst, Mass. OiOOd
The dynamics of fixed-bed performance was studied numerically. The model ignores axial diffusion in the fluid phase, but accounts for intraparticle diffusion and particle to fluid convective mass transfer. A homogeneous chemical reaction i s assumed to occur within the particles, the rate being governed b y a nonlinear Michaelis-Menten expression. The algorithm proposed for numerical solution lends itself readily to use with the IBM 360 continuous systems modeling program. A stability analysis gives good estimates of a maximum time step which ensures stability while minimizing computing time.
111 the course of developing models for describing the dynamics of transport phenomena in metabolically active tissue, a model was studied which can be cast in a form identical to one which describes the dynamics of ion exchange beds, chromatographic columns, and packed-bed reactors. This paper presents these results in a context suitable to the display of various aspects of column performance. Of particular interest is the inclusion of nonlinear kinetics, aiid an evaluation of the interaction of kinetics and intraparticle diffusion in determining breakthrough curves. Of equai importance is the development of a numerical scheme which has certain advantages over others previously used for this class of problems.
Mathematical Formulation
The system of interest is a packed bed of uiiiform spherical particles about and among which a liquid flows with steady uniform velocity. d material balance on some chemical species dissolved in the liquid phase leads to Equation 1,
dc dc 3 ~ 7 - + ~ - + - J ~ ( l - e ) dz at R
= O
is the flux of material (whose concentration is labeled c in the fluid) from the fluid into the solid. Diffusion in the fluid phase has been neglected relative to convection. W t h i n each particle the same species (whose concentration is labeled q within the particle) satisfies the material balance given in Equation 2 ,
This may be used to rewrite Equation 1 in the form
Equations 2 aiid 4 constitute a determinate set of equations once initial and boundary conditions on c and q are specified, and after a model for @(q) is assumed. At the interface an appropriate condition gives the flux in terms of a mass transfer coefficient,
(5) Spherical symmetry requires dq/br = 0 a t r = 0. Concentration c * is assumed to be in equilibrium with concentration q in the solid, a t the interface, so that c * = e* [q(R)]. For the case7 to be demonstrated a linear equilibrium is assumed, so that c * = b q(R) (6)
If the bed is initially a t solute concentration c o and q o , and is then subjected to an input e t , the initial and boundary conditions are c(z,O) =
J R
c0(2)
q(z,r,O) = qo(z,r) for all t
20
(7)
c(0,t) = c,(t) a t z = 0, for all t
>0
(8)
If the initial state is a t equilibrium, c o
=
b q"iR).
The problem is now cast in terms of dimensionless independent variables, giving Equations 9 and 10. The intraparticle diffusion coefficient is assumed to be a constant. A homogeneous chemical reaction is assumed to occur within each particle, and is characterized by a rate Gt which may depend on q. At the particle-fluid interface the flux, J R ,must balance the diffusion on the solid side of the interface, or
J R = D br 2 i 7=R
(3)
bc
bc
-
ax
+
G
+
= O
a
(10) Parameters LY aiid aP are Peclet numbers, defined by a = 3D(1
- s)L/R2U E
1
To whom correspondence should be sent.
624 Ind. Eng. Chem. Fundam., Vol. 9, No. 4, 1970
(9)
ffP
=
3(1 -
e)
ff
(12)
using the Duhamel theorem to convert Equation 1 to an integrodifferential equation. The present work differs from these earlier studies in two respects. The conversion of Equations 1 and 2 to a single integrodifferential equation seems to give rise to severe computational problems, even when a formal analytical solution of the resulting equation may be effected. Hence, the philosophy of this type of approach is rejected a priori. The secoud major difference 1s the inclusion of intraparticle chemical reaction. Since the problem will be treated numerically, there is little additional complication in considering nonlinear kinetics.
(13)
Numerical Method
(14)
One derivative is removed from Equation 9 by defining a dimensionless contact time
The independent variables are normalized to the bed length, L , sphere radius, R, and residence time, t o = E L / C ,and are given by
x
=
z/L
r/R
p =
T
=
t/t,
The reaction rate becomes
P(q) =
t
o
w
The initial and boundary conditions now take the form
c(x,O)
dX,P,O)
gO(x,p)
=
n(0,t) =
-
(11)
= C0(X)
Qi
(bq/b/)),=1=
y(bqIp-1
- c)
drlldPI,=O = 0
(15)
T = T - x
I n Equation 14 a Sherwood number appears, y =
hR/D
and changing variables. Then the system becomes
which is based on the inti aparticle diffusivity. Inspection shows that c and q depend on parameters E , cy, y , and b, and any parameters that appear in the reaction rate term, P ( q ) . T o complete the model, some statement regarding kinetics must be made. Becauqe of the biological context which motivated this work, the assumption was made that Michaelis-llenten enzyme kinetics were obeyed, so that
Kq
+ K,
= ___
q
or
Thus the assumption of Xichaelia-llenten kinetics introduces two parameters, K = ~t~ and K,. If K,, is either very large or very small, with reqpect to q , P(q) gives first-order or zeroorder kinetics, respectively.
K i t h the simple finite-difference approxiniations 2 =
m(Az), m
=
O,l,. . .,JI
Related Studies
Rosen (1952, 1954) studied the problem formulated here, but without the complication of chemical reaction. The problem he formulated, then, was linear and could be solved analytically. Basically, his method of approach was to use Duhamel's principle to write q as a function of g@). With Equations 5 aud 6 one could write y(R)in terms of c. As a result, q was expressed as an infinite series of integrals of time and length derivatives of c (z,t). Equation 1, for c, then takes the form of a n integrodifferential equation, where the integral involves an infinite series. Computationally, Rosen's forniulation requires the use of numerical methods for evaluation of the solution. The complexity of this problem is such that some of the advantage of having a formal analytical solution is washed out in the wake of generatiug numerical results. However, several useful asymptotic expressions are given which take fairly simple analytical forms, and from which considerable information can be extracted. More recently, Colwell and Dranoff (1969) examined the effects of axial dispersion (as characterized by the addition of a term - DL b 2c / b 9to the left-hand side of Equation 1 ) and introduced a nonlinear equilibrium at the particle-fluid boundary (by taking l / b to be a linear function of c* in Equation 6). They solved the differential equation for c after first
the above differential equations are represented by
for
m
=
1 , Ji'
j = l , J and where
Ind. Eng. Chem. Fundam., Vol. 9, No. 4, 1970
625
Figure 2.
T Effect of kinetic parameters on washout curves K, = 0.01, E = 0.1, b = 1 , y = 2 , a = 45
T
Figure 1 . Comparison of Rosen's analytical solutions with numerical solutions of this work t
= 0.1, b = 1 3 . 5 , ~=~ 13.5, y = 2.96(A) and 0.74(8)
Approximation of t'he initial and boundary conditions gives
C,(O) q,,j(O)
C",,
=
=
712 =
1,JI
(27)
q a m , j , m = 1,Jf
j = 1,J
CdT)
-
qm ,J+1
-
AP
qmJ
=
C i ( T ) ,T
b qm'Jc' =y[
(28)
>0
$- q n ' J -
2
(29)
.;I,
m T
=
1,JI
>0
(30)
upon the mesh size. A rather coarse mesh, J I = J = 10, was judged adequate for illustrating the solution technique and interpret'ing solutions in subsequent work. However, te&s with other mesh sizes indicated that solution accuracy could be greatly increased a t the expense of computer time, A sizable reduction in computer time call be realized by appropriately choosing AT, the incremental step in T . Ideally, AT should be near the maxiinurn value, (AT).vf,for which the solution is stable. (AT), may be estimat'ed by standard matrix techniques (Smith, 1965) if p,, j depends linearly upon qm, j , an assumption which was made in t,he stability analysis. For simplicity, our estimate of (AT),,{ assumed the use of Euler's integration method. The result of the stability analysis is the restriction (where ki =
K/Km)
qm,O = qm,l, m = 1,X
T>O
(31)
+
The iiodes located a t j = 0 and j = J 1 lie outside the sphere and are eliminated from Equations 25 and 26 by using Equations 30 and 31. Equation 26 represents -11 X J nonlinear first-order differential equations which must be integrated while simultaneously solving the M algebraic equations of Equation 25. Digital simulation languages developed in recent gears permit solving such a system with relatively little programming effort. I n this work t'he IBN 360 continuous systems modeling program (CSYLP) was used. The computation method was checked by setting P,,j = 0 and comparing washout curves with the corresponding exact solutions obtained by Rosen for two sets of bed parameters (Figure 1). An Euler integration routine with a fixed increment in T of 0.001 was used. Several other integration techniques %ere tested, and in each case the solution accuracy was only slightly improved while requiring substantially more computer time. For example, a variable-step, fourth-order Runge-Kutta method increased accuracy by less than 0.5'% but required 4 minutes us. 1 minute for Euler's method, Solution accuracy in this case is limited not by the integration technique but by the truncation error, which depends 626 Ind. Eng. Chem. Fundam., Vol. 9, No. 4, 1970
I3ecau.e the finite difference scheme use3 an average value of q a t the sphere surface (Equation 30), the Sherwood number, y , does not appear in the bounds on AT, and thus has no effect on stability. The conclusion of Equation 32 was ernpirically confirmed by integrating the equations with successively larger values of AT until the solution became unstable. This test was per~ kl. I n each trial informed for several values of ~ u p / ( A p )aiid stability occurred a t a value of AT 10 to 2070 larger than predicted by Equation 32. Such results are explained by the fact that the stability analysis considers all round-off errors to be of like sign, whereas in actuality they are not. Results
Figure 1 gives "washout" curves which correspond to the case co(s) 1 and C, ( T ) = 0, in the absence of chemical reaction. Because of the large number of combinations of parameters (such as a,y,K , aiid K,) possible when chemical reaction is considered, only an incomplete illustration of other iiumerical results can be presented here. Figure 2 shows the effect of reaction rate on washout. T h e
a=45 0010 w’w
‘0
5
10
-
15
20
20
15
3
T
25
Figure 3. Effect of particle Peclet number on breakthrough curves, in absence of chemical reaction e = 0.1,b = 1 , y =
I
IO
5
Figure 4. Effect of particle Peclet number on breakthrough curves, with chemical reaction K,
= 0.01,
K
= 0.05,
e = 0.1, b = 1, y =
2
2
initial state of the bed is taken to be the equilibrium distribution corresponding to a steady inlet coneenbration of unity. For no reaction, c o ( x ) 1, but for K # O theinitial distribution, e o @ ) , is not uniform, and depends on t’he kinetics. The initial outlet concentration, c(l@), depends on K , as can be seen in Figure 2. .Isone might expect, increasing the rate of reaction increases the rate a t which the bed is cleared, all other parameters being held constant. Figure 3 shows the effect of the particlel’eclet number on the approach to equilibrium, in the absence of chemical reaction. For this case the bed was taken to be initially enipty, co(r)5 0, and “breakthrough” curves are shown for the response to a steady unit input, C t ( T ) = 1. Figure 4 s h o w breakthrough curves for the same paranieters as in Figure 3, except that chemical reaction is included. In the absence of reaction (Figure 3), the hreakthrough curves approach equilibrium faster for the larger Peclet numbers, as espected. Under the conditions illustrated, however, chemical reaction can reverse this sit’uation,as Figure 4 sho~vs. Conclusions
The proposed algorithm for solving Equations 1 and 2 offers an esplicit technique for dealing with nonlinear source terms as well as a variety of boundary and initial conditions. With little difficulty, the method could be extended to account for other effects such as interparticle diffusion and temperature variations. Programming the algorithm is simple when a digital simulation program is available. Otherwise, a routine for integrating the large number of simultaneous ordinary differential equations would need to be written. The degree of solution accuracy depends upon the choices of A x , A p , and AT. It was found that solutions sufficiently accurate for many applications are obtained using Ax = A p = 0.1 and with AT slightly less than its upper bound of Equation 32.
Acknowledgment
The authors ale grateful to the Eastmaii Kodak Co. for generously contributing the computer time used in this qtudy. Nomenclature
b
=
C
= =
equilibrium coefficient, Equation 6 concentration i n fluid phase, nioles/cm3 C* concentration in eauilibriuni with solid surface, Equation 6 D = intraiiarticle diffu.ivitv. “ , em2, aec h = fluid-particle mass transfer coefficient, cm;/.?ec JR = surface flux, nioles,icm*-sec K , Knt = kinetic coefficients, Equation 16, sec-I, nioles,/cm3 L = bed length, cm P = dimensionless reaction rate = concent’ration in particle, moles,!em3 P r = radial coordinate in particle, cm R = radius of particle, em @ = reaction rate, Equation 16, moles,’cin3-sec t = time, sec t” = residence time, sec T = dimensionless contact time = superficial velocity, cm’sec X = dimensionless axial coordinate z = axial coordinate, cm CY, CYP = Peclet numbers Y = Sherm-ood number E = void fraction K = Kt, = dimensionless kinetic constant P = dimensionless radial coordinate T = dimensionless time I
u
literature Cited
COlWell, c. J., Dranoff, J. (1969).
s.,I K D . ENG.CHEhI. FUNDAM. 8, 193
\ - - - - I
Rosen, J. B., J . Chem. Phys. 20,387 (1952). Rosen, J. B., Ind. Eng. Chem. 46, 1590 (1964). Smith, G. D., “Sumerical Solution of Partial Differential Equations,” Oxford University Press, Cambridge, 1965. RECEIVED for review November 4, 1969 ACCEPTEDAugust 3, 1970 Ind. Eng. Chem. Fundam., Vol. 9, No. 4, 1970
627
