- Email: [email protected]
Non-isothermal multicomponent adsorption processes and their numerical treatment by means of integro-differential equations
Compners
&em.
Engng,
Vol. 13, No. 9, PP. 1025-1030, 1989
0098-I354/89 $3.00 f 0.00 Copyright0 1989 PergamonPressplc
Printed in Great Britain. All tights reserved
NON-ISOTHERMAL MULTICOMPONENT PROCESSES AND THEIR NUMERICAL MEANS OF INTEGRO-DIFFERENTIAL
ADSORPTION TREATMENT EQUATIONS
BY
G. NAGEL’ and G. KLUGE* ‘Sektion Mathematik and %ektion Chemie, Karl-Marx-UniversitPt, Leipzig 7010, G.D.R.. (Received
28 October
1988; final revision received 7 March
Karl-Marx-Platz,
1989; receiued for publicarion 3 1 March
1989)
Abstract--Considering a solid phase diffusion model for a multicomponent adsorption process, the mass and the energy balances are described by linear partial differential equations, but the boundary conditions of the mass balances are non-linear because of the non-linearity of the adsorption isotherms. In this case an analytical solution of the differential equations for the Dirichlet problem can be given. Applying this analytical solution, the prescribed boundary conditions change to integro-differential equations containing only the unknown values on the outer particle surface. However, these non-linear integro-differential equations have to be treated numerically. To solve this problem a numerical method is presented. Such a combined analytical and numerical treatment of the partial differential equations for the adsorbent particle is more effective than a merely numerical one.
INTRODUCTION
MODEL
One of the basic problems of modelling multicomponent adsorption processes in the fixed bed is the numerical treatment of the mass and the energy balance equations in the particle. The collocation method has proved to be an effective and generally applicable solution procedure (Villadsen and Stewart, 1967; Nagel and Adler, 1971; Nagel et al., 1987). A rougher approximation method, the so-called “linear driving force model” can be understood as the most simplified case of the collocation method. It is represented by Harwell et al. (1980), for example. In our recent paper (Nagel et al., 1987) a parallel diffusion model was used for the mass transfer in the particle. That means two diffusion fluxes progress in parallel in the direction of the particle centre, one of them (D,ac,,/~?r) in the macropores and the other one (D,aq,/ar) in the adsorbed phase. In many examples which have been investigated, the diffusion flux in the macropores was found to be negligible compared to the flux in the adsorbed phase. For this reason it should be possible to apply a simpler diffusion model which considers only the diffusion flux in the adsorbed phase, and one can take advantage of the special shape of the model equations in the mathematical treatment. On this assumption, the balance equations within the particle can be solved analytically. We have to employ only a numerical method to fulfil the boundary conditions on the external surface of the particle. This paper deals only with the balance equations for the particle. The coupling with the balance equations for the fluid phase in the fixed bed can be performed by wellknown methods (Nagel ei al., 1987).
EQUATIONS
In the development of the model equations, the temperature and concentration fields are assumed to be symmetrical within the spherical particle. The transport mechanism of the adsorbate is simplified by assuming a diffusion flux only in the adsorbed phase. The film mass transfer and the film heat transfer are involved by the boundary conditions on the external surface of the particle. All transport coefficients are assumed to be independent of concentration and temperature, respectively. As the process of adsorption onto the solid is rarely rate limiting, we can assume the existence of a local equilibrium between the concentrations in the macropores cPi and in the adsorbed qi phase throughout the particle. In the multicomponent case this equilibrium is described by the non-linear relations: cpi = c,;(4,,
i=l,...,n.
. . 74n. n,
Based on these suppositions, the balance equations (1) and (2) with the boundary conditions (3) and (4) and the initial conditions (5) and (6) are:
1025
dqi = at a=, at
’( -a
O
Tp
a24i I
D.
(
tTr*
;aq, rdr
>’
a=T, -+‘ ;$ dr*
+ >
i=l,...,n,
(1)
5 i=
(2)
tzo, bounded for r-0, i=l,...,n,
aT ar
2
&,
I
=
Bi,(T,-
T,,)
(3) (4)
1026
G.
r = 1,
for
i=l,...,n,
The
the relations:
(5)
T, = T,, independent O
KLUGE
Furthermore,
t > 0,
cpi = qi = 0,
for
and G.
NAGEL
of r
4Jt)
(6) =
S’0
are defined
dr
s
r=o.
coefficients
o’qi(r, t)3rz
=
by:
3D,&7i(l, -= dr
I)
O,$
(D,r)q:(t
-
’Diti (D;r)q:‘(t f0
t) dr,
7)
-
(14)
dr (15)
= q;(r),
or,R, /5 , Q;=y.
Bi,=
-AN-
ANALYTICAL
3a 8T,(l,
(7)
s
e
t)
ar
=. s
SOLUTION
’ 4
(e-27) T;(t
- 2 Qiw j)
The differential equations (1) and (2) as well as the boundary condition (4) are linear. Only the boundary condition (3) is non-linear. Due to the non-linear boundary condition (3) an analytical solution for the total problem (l-6) cannot be given. Assuming the values ~~(1, t) and T,(l, f) on the particle surface are given, an analytical solution of problems (1) and (2) can be stated satisfying the initial conditions (5) and (6) and the condition for the boundedness for r--+0. Applying the Laplace transformation one obtains the solution:
(0;~) - $ (a711
;= ,
x
are obtained,
t)dr
-
q;‘(r - r) d7.
(16)
with
q(p)=
Vde-P’$(t)dz i 0 =- ; (J;;coth,/$
-
=1-$+%-
1)
&+
-...
(17)
and q,(r, t) =
D,v
S’0
(r, QT)qT(t
-
7)
dr,
i=l,...,n.
(8) = 2
I T,(r,
r) = r,, +
s0
B,$
(B,z)q:‘(~
with
t,=-.
1’0
D, a -
and
the values
= 4i(l,
t>,
T;(r)
F (k=l
ljA _
transform
= r,c1,
(rl P) =
s 0
(18)
(as)
conditions
r;‘(i
- c;,(t)], -
(19)
i = I, . . , n,
T) dz
t>,
(11) respect
to
(12)
is:
e-p’ q (r, r) dt =
I’o [+ (oil>
- J/ (as)l~f’(r
-
7)
d7
particle
’2kx sin (her) eck2,+ I
7.
+
&
- ,$, QiW
on the external
=
its Laplace
e-@:‘.
(10)
Di
and the prime denotes the derivation with the time d/dt. The function 9 (r, t) is defined by: q (r, t) =
,g,
- r)dr
=3pi[c,(l)
&(t)
3+ 5
Using relations (15) and (16) the boundary (3) and (4) become:
J’0
T; and q: denote surface:
-
- T) - Tin] dr
arp (r, ar)[T;(t
sinh(r +‘$ r sinh Jd
)
’
(13)
3n Bi,[r,(t)
-
r:(r)],
t > 0.
(20)
The original problem (l-6) reduces to the treatment of the integro-differential equations (IDES) (I 9) and (20) with the initial conditions (5) and (6) (for r = 1). In the IDES, only the values on the particle surface q), ct, TS, or their derivatives with respect to time appear. Contrary to the “linear driving force model” relations (19) and (20) are exact. The IDES (19) for the mass balances have already been stated by Solotarev and Dubinin (1973). As far as we know, such IDES have not yet been employed for the calculation of breakthrough curves or kinetic curves from single-particle experiments. The reason is that
Non-isothermal multicomponent adsorption processes equations of this type have rarely been treated in numerical mathematics, and standard methods for their treatment have not yet been published. In the next section an effective method for their numerical treatment will be presented. To apply this method the energy balance (20) has to be transformed such that only one convolution integral occurs as in the mass balances (19). After multiplying equations (19) with Q,(l +ci), i = 1,. . . , n, and adding the results to equation (20) and considering definition (IO), we obtain:
f s0
a$ (m)yC(t
- r) dr = 3a Bi,[T,(t)
NUMERICAL METHODS Let us consider the IDE (problem
t >O,
c;,tt)l,
I > 0,
(21)
(28)
i b,f(t k=O
-kh,y(t
2 Qiwl(r)
vw(t) = T”,(t) - T, +
(22)
,=I
and the initial condition Y,(O) = 0.
(23)
Now the IDES (19) and (21) are to be treated numerically. For the moment we concentrate on a relation between the mean loading qi and the loading on the particle surface q: , which enables us to estimate when the diffusion resistance can be neglected. From relation (17) it follows: m e-@D& (Dir) dt s0
--2P2 15D, +3lSD;’
=1-A!-
+‘*.
(24)
On the other hand, for the Laplace transform of the Dirac distribution b (t - to) yields: ic e-p’s (t - I,) dr I0 2 =e
2
-loY&foP++
+...
(25)
The Dirac distribution S (r - 1,) with 1
for the function represents an approximation D& (Dir) in relation (14) as far as their 0th and 1st moment agree and their Laplace transforms differ only to a term 0 (0;‘) for Di--t co. Replacing the function D,JI (Dir) in relation (14) by the Dirac distribution 6 (r - T,,).one obtains the approximation: IjJt) =
I 6 (r - to)q;(t
-
7)
= 0,
f(t,y(t))
dr = q;(t - to).
s0
For large values of D, the mean loading 4, lags behind the loading on the particle surface q: by the time distance 1”.
-kh)), bo
=!=0,
(29)
= 0
for
t < 0.
(30)
The function y(t) is to be found out in both problems, where the kernel k (t) and the right-handside f(r, Y (t)) are given as well as the constant coefficients akr b, and the step size h. For simplicity’s sake we suppose that the functions k (t), y(t), y’(f), f(t, y (t)) are Laplace transformable, that means: k(p)=
coe -P’k (I) dt,
s0
m
Y(P)
e+y
=
(r) dt,
P_F (P) =
e-P'y'(f) dt,
s 0
f(p)
=
sa
eFP’f(t,y
(t)) dt,
0
may exist for Re a(p)). Obviously This dependence, Furthermore, we
p > c with a constant c (c < 0 for the functiony(p) depends on y(t). however, is not especially denoted. suppose the normalization: X k(t)dt = 1, J0 that means E(O) = 1. Let us ask: are there such functions k(t) that problems 1 and 2 are equivalent? This question can simply be answered by applying the Laplace transformation to both of the problems. Problem 1: ff(P)PY
(26)
*o = m
difference
with the initial condition: y(t)
with
(27)
Y (0) = 0 and, on the other hand, the implicit equation of order I (problem 2)
t =- 0, + 5 Qi(I + c,) 38Jcfi(r) ,= I
1):
‘k(r)y’(r --r)dr =f(t,y(t)), I0 with the initial condition:
i a,y(t-kh)=h k=O
- T;(I)]
1027
Problem
(P) -J(P)-
(31)
2:
(a0 + (I, emhp+ . . + a,e-‘hp)P (p) =-h [b, + b,e-“p +.
. . + b,e-‘*p]f(p).
(32)
From equation (32) it follows:
Thus, if E(p) with K(P)
=
problems
has the particular shape E(p)
= l?,(p)
aoz’+a,z’-‘+...+a, z = ehp, (34) hp (b,z’ + b, z’- I + . . . + 6,) ’ 1 and 2 are equivalent.
G. NAGEL~~~G.KLUGE
1028
Here the polynomials in the numerator and the denominator of term (34) are supposed to be relatively prime. For the existence of an original function K,(t) to K,(p) the coefficient b, has to be #O. The original function K,(t) is a step function with a step size h. We can dispense with a detailed discussion of the function K,(r) since these facts are known. However, i?,(p) may be supposed to be analytical in zero p = 0 and, furthermore, z,(O) = 1, which is analogous to c(p). Hence, R,(p) must satisfy the consistence conditions:
which are well-known from the ordinary differential equations. For a kernel k(t) with a Laplace transform f;(p) not of the shape (34) we can approximate G(p) by a function K,(p) of the shape (34) whose coefficients ak, bk, k = 0, . ,1, are disposed so that z,(p) agrees with g(p) ‘*as well as possible”. By “as well as possible” we will understand that the power series expansions of both functions z,(p) and k”(p) in zero p = 0 agree up to a powerp” as high as possible. That means the derivatives &J’(O) and /%‘(O), j = 1, . ,)r”, agree as well as the moments (if they exist) original the functions K,(I) and k(r): L s0
K,(r) t’dt =
‘* i 0
k(l)Pdt,
j=l,...,
tn.
Asforthe21+2parametersak,b,,k=0,1,...,I,of function (34), one of them is redundant and two of them are disposed by the consistence conditions, but 21 - 1 parameters are available. In the applications we will restrict to 2 < 2 because the expansion of functions E(p) and I?,(p) into the power series is rather troublesome. More important than a good approximation of the power series of E(p) and K,(p) is the stability of the arising difference equation (29). It can prove to be necessary to agree not all powers to YPZ = 21 - 1, and to dispose of the free parameters so that the difference equation (29) gets good stability properties. Remark-If kernel k(t) in IDE (27) is the Dirac distribution S(t), then (27) represents an ordinary differential equation; J”(f)
=.f(G r(t)),
with f;(p) = 1. In this case the approximation method described above represents the well-known implicit difference method for the numerical treatment of ordinary differential equations. The difference method (29) for the treatment of an IDE (27) behaves in a much more good-natured way with respect to the stability and to a good approximability, compared to this method for ordinary differential equations.
Obviously, one can generalize the difference method described above to systems of IDES:
r‘ki(~)y;(r-
. .l r,(t)>,
r)dT =.f;(t,y,(l).
Jo
i=i.....n.
t
ro,
y;(O) = 0.
In the case of IDES (19) and (21) it yields [cf. (17)]: p(p)=1
-&P*+&JJ*~--&~*‘+
-“.
(36)
with p* =g
for IDES
(19)
p* = $
for IDE
(21).
Thus, for the one-step coefficients: a,=
1.
a,=
i = I,. . , , n, (37)
method 1
-1,
b,=2+-
(1 = 1, 171= 1) the
1 lSh+’
b, = 1 - b,,
(3X)
result in: h* = lzDi
for IDES (19)
h*=ha
for IDE
i = 1,
. , n, (39)
(21).
However, the one-step method is a somewhat rough approximation (approximately as exact as the linear driving force model). Therefore, we have applied a two-step method (I =2). For that one the coefficients: a, = 1 -2u,.
az=a,-
1,
6, = 1 -bb,--bhz,
(40)
result from the consistence conditions. where the coefficients a,,. b,, h, have to satisfy the conditions for the moment agreements:
a,=-+
1 2
2
(3rd moment),
23625wh*
(43)
with w
=&+%.
*’
(W
For little values h* (z lop2 or smaller) one can agree all three moments (m = 3). Then, the coefficient a, results from condition (43). If h* increases, the coefficient a, decreases. Keeping the monotonic stability of the method, the coefficients a2 and b2 may not be too small. Therefore, for larger values h* (x IO--‘) it is better to agree only two moments (m = 2), and to choose a larger value for n, than would result from condition (43). The largeness chosen for a, depends on the right-hand-side of the IDE. Moreover, attention is to be paid to a particularity of the functionf(t. _J~(If) at zero I = 0. For a bounded
Non-isothermal multicomponent adsorption processes kernel k(l) not changing its sign for 0 < t < to the left-hand-side of the IDE (27) becomes:
s
1029
Table 1. Comparison of relative CPU times for example calculations by various methods (m = number of coHocation points) Collocation method This method
m =2
??I=3
m=4
Linear driving force model
I
a3
z5.5
29.5
El
I
lim r-+0 Thus, functionf(t, tion
0
- 7) dt = 0.
k(s)y’(r
v(t))
(45)
has also to satisfy the condi-
lim f(r, r(t)) I-+0
= 0.
For a singular kernel k(t) as in the cases of IDES (19) and (21) discussed above it is possible that:
f lim k(r)y’(t t-.+0 s 0
- 7) dz + 0
for the graphical representation of the breakthrough curves, it is referred to the previous papers mentioned above. Unfortunately, for this example, the numerical differences by various methods listed in Table 1 were too small to illustrate these very well in the graphical representation of the breakthrough curves.
and therefore:
CONCLUSIONS lim
I-+0
fCi.
y (t 1) + 0.
(46)
For example, from: k(t) N C,t-“2, and condition
y’(t)
_ C*l-“2,
t*
+o
(28) follows:
, s0
WT)Y’(~
-7)d7 y(t)
=f(&y@))-
C,C,z,
2C2r1’2
for r -b+O. As the approximative kernel K,(r) for r - 4-O is always bounded and does not change the sign for 0 < r c h, it is always necessary to set: f(r, y(r))
= 0
for r = 0,
The combined analytical and numerical treatment of the balance equations described above results in a considerable reduction of the calculation expense compared to the collocation method. Just the collocation method is actually an effective and accurate method. Certainly, the new method is only qualified for the solid phase diffusion model, but the collocation method is applicable to a larger extent. The numerical method for the treatment of integrodifferential equations proposed in this work would also be of interest independently of adsorpiion problems. Furthermore, the integro-differential equations have proved to be easier to treat by implicit difference methods than ordinary differential equations of first order.
in difference equation (29), also in cases where condition (46) is valid. This is to be considered for the start steps. Comparison
of numerical
methods
To test the method described here we have chosen the adsorption of a two-component system in a fixed bed as an example [the adsorption of the solutes of ethane and carbon dioxide in a bed of 5A molecular sieves (Basmadjian and Wright, 198111.This example has already been used by applying the collocation method in our recent papers (Nagel er al., 1987; Nagel and Kluge, 1987). Table 1 contains the relative CPU times for the calculation of the concentrations and temperature breakthrough curves of this example. Compared to the collocation method, the new method reduced the computing time considerably. As for the numerical accuracy, the calculated concentrations and temperature values of both methods agreed practically. On the other hand, the computing expense with the new method was found to be comparable with that applying the simpler “linear driving force model”, but the accuracy of the numerical data was better. Moreover, it should be noticed that one has to neglect the resistance for the mass and heat transfer on the external particle surface (pi, Bi, - co) when the “linear driving force model” is applied. As
AND SIGNIFICANCE
NOMENCLATURE
of heat transfer within the adsorbent particle, =I,l(p,C,R;), s-’ Biot number for the heat transfer, particle fluid = uhR,//l, Heat capacity of the adsorbent particle, kJ kg-‘K-’ Concentration of i in the flowing medium, kg of i me3 Concentration of i in the fluid phase within the particle, kg of i mm3 Coefficient of the mass transfer of i in the SK’ adsorbed phase =D,,lRi, Effective diffusion coefficient of i in the adsorbed
a = Coefficient
Bi, = C, = cl,.= cPi= Dj = D,, = -Ali,
phase, m* s-’
= Adsorption heat of component i, kJ kg-’ h = Step size k, = Film mass transfer coefficient of i, m S-’ n = Number of adsorbable components p = Laplace transform variable, s-’ e, = -AH,,lC,, K (I, _.= Concentration of i in the adsorbed phase (loading), kg of i kg--’ solid q, = Integral mean value of q,, kg of i kg-’ solid R, = Radius of the spherical adsorbent particle, m r = Normalized radial coordinate of the particle T,= Temperature of the flowing medium, K T,, = Initial temperature of the particle, K TP= Temperature of the particle, K t = Time, s t, = cf. (26) y, = cf. (22)
G. NAGEL and G. KLUGE
1030 Greek
ferrers ah = Heat transfer
coeficient,
particle-fluid,
kW m-r
P, = :,!;P R ) c,=cf. (16) p E., = Effective thermal conductivity of the particle, kW m-’ K-t pp = Density of the adsorbent, kg m-’ particle SUperSCriptS -
s = Particle = Laplace ’= d/dr
surface transform
REFERENCES
Basmadjian D. and D. W. Wright, Non-isothermal sorption of ethanecarbon dioxide mixture in beds of SA molecular sieves. Chem. Engng Sri. 34, 937-940 (1981).
Harwell J. H., A. I. Liapis, R. Litchfield and D. T. Hanson, A non-equilibrium model for fixed-bed multi-component adiabatic adsorption. Chem. Engng Sci. 35, 2287-2296 (1980). Nagel G. and R. Adler, Anwedung der Kollokationsmethode auf die Berechnung polytroper Rohrreaktoren. Chem. Technik. 23, 335-341 (1971). Nagel G., G. Kluge and W. Flock, Modelling of nonisothermal multi-component adsorption in adiabatic fixed beds--I. The numerical solution of the parallel diffusion model. Chem. Engng Sci. 42, 143- 153 (1987). Nagel G. and G. Kluge, Modelling of non-isothermal multi-component adsorption in non-adiabatic fixed beds. Nrtn. J. Ind. Chem. 15, 63--71 (1987). Solo&rev P. P. and M. M. Dubinin, On equations describing the internal diffusivity within granular adsorbents (in (1973). Russian). DokI. Akad. Nauk. SSSR 210, 136139 Villadsen J. V. and W. E. Stewart. Solution of boundarvvalue problems by orthogonal collocation. Chem. Enngng Sci. 22, 148331501 (1967).
&em.
Engng,
Vol. 13, No. 9, PP. 1025-1030, 1989
0098-I354/89 $3.00 f 0.00 Copyright0 1989 PergamonPressplc
Printed in Great Britain. All tights reserved
NON-ISOTHERMAL MULTICOMPONENT PROCESSES AND THEIR NUMERICAL MEANS OF INTEGRO-DIFFERENTIAL
ADSORPTION TREATMENT EQUATIONS
BY
G. NAGEL’ and G. KLUGE* ‘Sektion Mathematik and %ektion Chemie, Karl-Marx-UniversitPt, Leipzig 7010, G.D.R.. (Received
28 October
1988; final revision received 7 March
Karl-Marx-Platz,
1989; receiued for publicarion 3 1 March
1989)
Abstract--Considering a solid phase diffusion model for a multicomponent adsorption process, the mass and the energy balances are described by linear partial differential equations, but the boundary conditions of the mass balances are non-linear because of the non-linearity of the adsorption isotherms. In this case an analytical solution of the differential equations for the Dirichlet problem can be given. Applying this analytical solution, the prescribed boundary conditions change to integro-differential equations containing only the unknown values on the outer particle surface. However, these non-linear integro-differential equations have to be treated numerically. To solve this problem a numerical method is presented. Such a combined analytical and numerical treatment of the partial differential equations for the adsorbent particle is more effective than a merely numerical one.
INTRODUCTION
MODEL
One of the basic problems of modelling multicomponent adsorption processes in the fixed bed is the numerical treatment of the mass and the energy balance equations in the particle. The collocation method has proved to be an effective and generally applicable solution procedure (Villadsen and Stewart, 1967; Nagel and Adler, 1971; Nagel et al., 1987). A rougher approximation method, the so-called “linear driving force model” can be understood as the most simplified case of the collocation method. It is represented by Harwell et al. (1980), for example. In our recent paper (Nagel et al., 1987) a parallel diffusion model was used for the mass transfer in the particle. That means two diffusion fluxes progress in parallel in the direction of the particle centre, one of them (D,ac,,/~?r) in the macropores and the other one (D,aq,/ar) in the adsorbed phase. In many examples which have been investigated, the diffusion flux in the macropores was found to be negligible compared to the flux in the adsorbed phase. For this reason it should be possible to apply a simpler diffusion model which considers only the diffusion flux in the adsorbed phase, and one can take advantage of the special shape of the model equations in the mathematical treatment. On this assumption, the balance equations within the particle can be solved analytically. We have to employ only a numerical method to fulfil the boundary conditions on the external surface of the particle. This paper deals only with the balance equations for the particle. The coupling with the balance equations for the fluid phase in the fixed bed can be performed by wellknown methods (Nagel ei al., 1987).
EQUATIONS
In the development of the model equations, the temperature and concentration fields are assumed to be symmetrical within the spherical particle. The transport mechanism of the adsorbate is simplified by assuming a diffusion flux only in the adsorbed phase. The film mass transfer and the film heat transfer are involved by the boundary conditions on the external surface of the particle. All transport coefficients are assumed to be independent of concentration and temperature, respectively. As the process of adsorption onto the solid is rarely rate limiting, we can assume the existence of a local equilibrium between the concentrations in the macropores cPi and in the adsorbed qi phase throughout the particle. In the multicomponent case this equilibrium is described by the non-linear relations: cpi = c,;(4,,
i=l,...,n.
. . 74n. n,
Based on these suppositions, the balance equations (1) and (2) with the boundary conditions (3) and (4) and the initial conditions (5) and (6) are:
1025
dqi = at a=, at
’( -a
O
Tp
a24i I
D.
(
tTr*
;aq, rdr
>’
a=T, -+‘ ;$ dr*
+ >
i=l,...,n,
(1)
5 i=
(2)
tzo, bounded for r-0, i=l,...,n,
aT ar
2
&,
I
=
Bi,(T,-
T,,)
(3) (4)
1026
G.
r = 1,
for
i=l,...,n,
The
the relations:
(5)
T, = T,, independent O
KLUGE
Furthermore,
t > 0,
cpi = qi = 0,
for
and G.
NAGEL
of r
4Jt)
(6) =
S’0
are defined
dr
s
r=o.
coefficients
o’qi(r, t)3rz
=
by:
3D,&7i(l, -= dr
I)
O,$
(D,r)q:(t
-
’Diti (D;r)q:‘(t f0
t) dr,
7)
-
(14)
dr (15)
= q;(r),
or,R, /5 , Q;=y.
Bi,=
-AN-
ANALYTICAL
3a 8T,(l,
(7)
s
e
t)
ar
=. s
SOLUTION
’ 4
(e-27) T;(t
- 2 Qiw j)
The differential equations (1) and (2) as well as the boundary condition (4) are linear. Only the boundary condition (3) is non-linear. Due to the non-linear boundary condition (3) an analytical solution for the total problem (l-6) cannot be given. Assuming the values ~~(1, t) and T,(l, f) on the particle surface are given, an analytical solution of problems (1) and (2) can be stated satisfying the initial conditions (5) and (6) and the condition for the boundedness for r--+0. Applying the Laplace transformation one obtains the solution:
(0;~) - $ (a711
;= ,
x
are obtained,
t)dr
-
q;‘(r - r) d7.
(16)
with
q(p)=
Vde-P’$(t)dz i 0 =- ; (J;;coth,/$
-
=1-$+%-
1)
&+
-...
(17)
and q,(r, t) =
D,v
S’0
(r, QT)qT(t
-
7)
dr,
i=l,...,n.
(8) = 2
I T,(r,
r) = r,, +
s0
B,$
(B,z)q:‘(~
with
t,=-.
1’0
D, a -
and
the values
= 4i(l,
t>,
T;(r)
F (k=l
ljA _
transform
= r,c1,
(rl P) =
s 0
(18)
(as)
conditions
r;‘(i
- c;,(t)], -
(19)
i = I, . . , n,
T) dz
t>,
(11) respect
to
(12)
is:
e-p’ q (r, r) dt =
I’o [+ (oil>
- J/ (as)l~f’(r
-
7)
d7
particle
’2kx sin (her) eck2,+ I
7.
+
&
- ,$, QiW
on the external
=
its Laplace
e-@:‘.
(10)
Di
and the prime denotes the derivation with the time d/dt. The function 9 (r, t) is defined by: q (r, t) =
,g,
- r)dr
=3pi[c,(l)
&(t)
3+ 5
Using relations (15) and (16) the boundary (3) and (4) become:
J’0
T; and q: denote surface:
-
- T) - Tin] dr
arp (r, ar)[T;(t
sinh(r +‘$ r sinh Jd
)
’
(13)
3n Bi,[r,(t)
-
r:(r)],
t > 0.
(20)
The original problem (l-6) reduces to the treatment of the integro-differential equations (IDES) (I 9) and (20) with the initial conditions (5) and (6) (for r = 1). In the IDES, only the values on the particle surface q), ct, TS, or their derivatives with respect to time appear. Contrary to the “linear driving force model” relations (19) and (20) are exact. The IDES (19) for the mass balances have already been stated by Solotarev and Dubinin (1973). As far as we know, such IDES have not yet been employed for the calculation of breakthrough curves or kinetic curves from single-particle experiments. The reason is that
Non-isothermal multicomponent adsorption processes equations of this type have rarely been treated in numerical mathematics, and standard methods for their treatment have not yet been published. In the next section an effective method for their numerical treatment will be presented. To apply this method the energy balance (20) has to be transformed such that only one convolution integral occurs as in the mass balances (19). After multiplying equations (19) with Q,(l +ci), i = 1,. . . , n, and adding the results to equation (20) and considering definition (IO), we obtain:
f s0
a$ (m)yC(t
- r) dr = 3a Bi,[T,(t)
NUMERICAL METHODS Let us consider the IDE (problem
t >O,
c;,tt)l,
I > 0,
(21)
(28)
i b,f(t k=O
-kh,y(t
2 Qiwl(r)
vw(t) = T”,(t) - T, +
(22)
,=I
and the initial condition Y,(O) = 0.
(23)
Now the IDES (19) and (21) are to be treated numerically. For the moment we concentrate on a relation between the mean loading qi and the loading on the particle surface q: , which enables us to estimate when the diffusion resistance can be neglected. From relation (17) it follows: m e-@D& (Dir) dt s0
--2P2 15D, +3lSD;’
=1-A!-
+‘*.
(24)
On the other hand, for the Laplace transform of the Dirac distribution b (t - to) yields: ic e-p’s (t - I,) dr I0 2 =e
2
-loY&foP++
+...
(25)
The Dirac distribution S (r - 1,) with 1
for the function represents an approximation D& (Dir) in relation (14) as far as their 0th and 1st moment agree and their Laplace transforms differ only to a term 0 (0;‘) for Di--t co. Replacing the function D,JI (Dir) in relation (14) by the Dirac distribution 6 (r - T,,).one obtains the approximation: IjJt) =
I 6 (r - to)q;(t
-
7)
= 0,
f(t,y(t))
dr = q;(t - to).
s0
For large values of D, the mean loading 4, lags behind the loading on the particle surface q: by the time distance 1”.
-kh)), bo
=!=0,
(29)
= 0
for
t < 0.
(30)
The function y(t) is to be found out in both problems, where the kernel k (t) and the right-handside f(r, Y (t)) are given as well as the constant coefficients akr b, and the step size h. For simplicity’s sake we suppose that the functions k (t), y(t), y’(f), f(t, y (t)) are Laplace transformable, that means: k(p)=
coe -P’k (I) dt,
s0
m
Y(P)
e+y
=
(r) dt,
P_F (P) =
e-P'y'(f) dt,
s 0
f(p)
=
sa
eFP’f(t,y
(t)) dt,
0
may exist for Re a(p)). Obviously This dependence, Furthermore, we
p > c with a constant c (c < 0 for the functiony(p) depends on y(t). however, is not especially denoted. suppose the normalization: X k(t)dt = 1, J0 that means E(O) = 1. Let us ask: are there such functions k(t) that problems 1 and 2 are equivalent? This question can simply be answered by applying the Laplace transformation to both of the problems. Problem 1: ff(P)PY
(26)
*o = m
difference
with the initial condition: y(t)
with
(27)
Y (0) = 0 and, on the other hand, the implicit equation of order I (problem 2)
t =- 0, + 5 Qi(I + c,) 38Jcfi(r) ,= I
1):
‘k(r)y’(r --r)dr =f(t,y(t)), I0 with the initial condition:
i a,y(t-kh)=h k=O
- T;(I)]
1027
Problem
(P) -J(P)-
(31)
2:
(a0 + (I, emhp+ . . + a,e-‘hp)P (p) =-h [b, + b,e-“p +.
. . + b,e-‘*p]f(p).
(32)
From equation (32) it follows:
Thus, if E(p) with K(P)
=
problems
has the particular shape E(p)
= l?,(p)
aoz’+a,z’-‘+...+a, z = ehp, (34) hp (b,z’ + b, z’- I + . . . + 6,) ’ 1 and 2 are equivalent.
G. NAGEL~~~G.KLUGE
1028
Here the polynomials in the numerator and the denominator of term (34) are supposed to be relatively prime. For the existence of an original function K,(t) to K,(p) the coefficient b, has to be #O. The original function K,(t) is a step function with a step size h. We can dispense with a detailed discussion of the function K,(r) since these facts are known. However, i?,(p) may be supposed to be analytical in zero p = 0 and, furthermore, z,(O) = 1, which is analogous to c(p). Hence, R,(p) must satisfy the consistence conditions:
which are well-known from the ordinary differential equations. For a kernel k(t) with a Laplace transform f;(p) not of the shape (34) we can approximate G(p) by a function K,(p) of the shape (34) whose coefficients ak, bk, k = 0, . ,1, are disposed so that z,(p) agrees with g(p) ‘*as well as possible”. By “as well as possible” we will understand that the power series expansions of both functions z,(p) and k”(p) in zero p = 0 agree up to a powerp” as high as possible. That means the derivatives &J’(O) and /%‘(O), j = 1, . ,)r”, agree as well as the moments (if they exist) original the functions K,(I) and k(r): L s0
K,(r) t’dt =
‘* i 0
k(l)Pdt,
j=l,...,
tn.
Asforthe21+2parametersak,b,,k=0,1,...,I,of function (34), one of them is redundant and two of them are disposed by the consistence conditions, but 21 - 1 parameters are available. In the applications we will restrict to 2 < 2 because the expansion of functions E(p) and I?,(p) into the power series is rather troublesome. More important than a good approximation of the power series of E(p) and K,(p) is the stability of the arising difference equation (29). It can prove to be necessary to agree not all powers to YPZ = 21 - 1, and to dispose of the free parameters so that the difference equation (29) gets good stability properties. Remark-If kernel k(t) in IDE (27) is the Dirac distribution S(t), then (27) represents an ordinary differential equation; J”(f)
=.f(G r(t)),
with f;(p) = 1. In this case the approximation method described above represents the well-known implicit difference method for the numerical treatment of ordinary differential equations. The difference method (29) for the treatment of an IDE (27) behaves in a much more good-natured way with respect to the stability and to a good approximability, compared to this method for ordinary differential equations.
Obviously, one can generalize the difference method described above to systems of IDES:
r‘ki(~)y;(r-
. .l r,(t)>,
r)dT =.f;(t,y,(l).
Jo
i=i.....n.
t
ro,
y;(O) = 0.
In the case of IDES (19) and (21) it yields [cf. (17)]: p(p)=1
-&P*+&JJ*~--&~*‘+
-“.
(36)
with p* =g
for IDES
(19)
p* = $
for IDE
(21).
Thus, for the one-step coefficients: a,=
1.
a,=
i = I,. . , , n, (37)
method 1
-1,
b,=2+-
(1 = 1, 171= 1) the
1 lSh+’
b, = 1 - b,,
(3X)
result in: h* = lzDi
for IDES (19)
h*=ha
for IDE
i = 1,
. , n, (39)
(21).
However, the one-step method is a somewhat rough approximation (approximately as exact as the linear driving force model). Therefore, we have applied a two-step method (I =2). For that one the coefficients: a, = 1 -2u,.
az=a,-
1,
6, = 1 -bb,--bhz,
(40)
result from the consistence conditions. where the coefficients a,,. b,, h, have to satisfy the conditions for the moment agreements:
a,=-+
1 2
2
(3rd moment),
23625wh*
(43)
with w
=&+%.
*’
(W
For little values h* (z lop2 or smaller) one can agree all three moments (m = 3). Then, the coefficient a, results from condition (43). If h* increases, the coefficient a, decreases. Keeping the monotonic stability of the method, the coefficients a2 and b2 may not be too small. Therefore, for larger values h* (x IO--‘) it is better to agree only two moments (m = 2), and to choose a larger value for n, than would result from condition (43). The largeness chosen for a, depends on the right-hand-side of the IDE. Moreover, attention is to be paid to a particularity of the functionf(t. _J~(If) at zero I = 0. For a bounded
Non-isothermal multicomponent adsorption processes kernel k(l) not changing its sign for 0 < t < to the left-hand-side of the IDE (27) becomes:
s
1029
Table 1. Comparison of relative CPU times for example calculations by various methods (m = number of coHocation points) Collocation method This method
m =2
??I=3
m=4
Linear driving force model
I
a3
z5.5
29.5
El
I
lim r-+0 Thus, functionf(t, tion
0
- 7) dt = 0.
k(s)y’(r
v(t))
(45)
has also to satisfy the condi-
lim f(r, r(t)) I-+0
= 0.
For a singular kernel k(t) as in the cases of IDES (19) and (21) discussed above it is possible that:
f lim k(r)y’(t t-.+0 s 0
- 7) dz + 0
for the graphical representation of the breakthrough curves, it is referred to the previous papers mentioned above. Unfortunately, for this example, the numerical differences by various methods listed in Table 1 were too small to illustrate these very well in the graphical representation of the breakthrough curves.
and therefore:
CONCLUSIONS lim
I-+0
fCi.
y (t 1) + 0.
(46)
For example, from: k(t) N C,t-“2, and condition
y’(t)
_ C*l-“2,
t*
+o
(28) follows:
, s0
WT)Y’(~
-7)d7 y(t)
=f(&y@))-
C,C,z,
2C2r1’2
for r -b+O. As the approximative kernel K,(r) for r - 4-O is always bounded and does not change the sign for 0 < r c h, it is always necessary to set: f(r, y(r))
= 0
for r = 0,
The combined analytical and numerical treatment of the balance equations described above results in a considerable reduction of the calculation expense compared to the collocation method. Just the collocation method is actually an effective and accurate method. Certainly, the new method is only qualified for the solid phase diffusion model, but the collocation method is applicable to a larger extent. The numerical method for the treatment of integrodifferential equations proposed in this work would also be of interest independently of adsorpiion problems. Furthermore, the integro-differential equations have proved to be easier to treat by implicit difference methods than ordinary differential equations of first order.
in difference equation (29), also in cases where condition (46) is valid. This is to be considered for the start steps. Comparison
of numerical
methods
To test the method described here we have chosen the adsorption of a two-component system in a fixed bed as an example [the adsorption of the solutes of ethane and carbon dioxide in a bed of 5A molecular sieves (Basmadjian and Wright, 198111.This example has already been used by applying the collocation method in our recent papers (Nagel er al., 1987; Nagel and Kluge, 1987). Table 1 contains the relative CPU times for the calculation of the concentrations and temperature breakthrough curves of this example. Compared to the collocation method, the new method reduced the computing time considerably. As for the numerical accuracy, the calculated concentrations and temperature values of both methods agreed practically. On the other hand, the computing expense with the new method was found to be comparable with that applying the simpler “linear driving force model”, but the accuracy of the numerical data was better. Moreover, it should be noticed that one has to neglect the resistance for the mass and heat transfer on the external particle surface (pi, Bi, - co) when the “linear driving force model” is applied. As
AND SIGNIFICANCE
NOMENCLATURE
of heat transfer within the adsorbent particle, =I,l(p,C,R;), s-’ Biot number for the heat transfer, particle fluid = uhR,//l, Heat capacity of the adsorbent particle, kJ kg-‘K-’ Concentration of i in the flowing medium, kg of i me3 Concentration of i in the fluid phase within the particle, kg of i mm3 Coefficient of the mass transfer of i in the SK’ adsorbed phase =D,,lRi, Effective diffusion coefficient of i in the adsorbed
a = Coefficient
Bi, = C, = cl,.= cPi= Dj = D,, = -Ali,
phase, m* s-’
= Adsorption heat of component i, kJ kg-’ h = Step size k, = Film mass transfer coefficient of i, m S-’ n = Number of adsorbable components p = Laplace transform variable, s-’ e, = -AH,,lC,, K (I, _.= Concentration of i in the adsorbed phase (loading), kg of i kg--’ solid q, = Integral mean value of q,, kg of i kg-’ solid R, = Radius of the spherical adsorbent particle, m r = Normalized radial coordinate of the particle T,= Temperature of the flowing medium, K T,, = Initial temperature of the particle, K TP= Temperature of the particle, K t = Time, s t, = cf. (26) y, = cf. (22)
G. NAGEL and G. KLUGE
1030 Greek
ferrers ah = Heat transfer
coeficient,
particle-fluid,
kW m-r
P, = :,!;P R ) c,=cf. (16) p E., = Effective thermal conductivity of the particle, kW m-’ K-t pp = Density of the adsorbent, kg m-’ particle SUperSCriptS -
s = Particle = Laplace ’= d/dr
surface transform
REFERENCES
Basmadjian D. and D. W. Wright, Non-isothermal sorption of ethanecarbon dioxide mixture in beds of SA molecular sieves. Chem. Engng Sri. 34, 937-940 (1981).
Harwell J. H., A. I. Liapis, R. Litchfield and D. T. Hanson, A non-equilibrium model for fixed-bed multi-component adiabatic adsorption. Chem. Engng Sci. 35, 2287-2296 (1980). Nagel G. and R. Adler, Anwedung der Kollokationsmethode auf die Berechnung polytroper Rohrreaktoren. Chem. Technik. 23, 335-341 (1971). Nagel G., G. Kluge and W. Flock, Modelling of nonisothermal multi-component adsorption in adiabatic fixed beds--I. The numerical solution of the parallel diffusion model. Chem. Engng Sci. 42, 143- 153 (1987). Nagel G. and G. Kluge, Modelling of non-isothermal multi-component adsorption in non-adiabatic fixed beds. Nrtn. J. Ind. Chem. 15, 63--71 (1987). Solo&rev P. P. and M. M. Dubinin, On equations describing the internal diffusivity within granular adsorbents (in (1973). Russian). DokI. Akad. Nauk. SSSR 210, 136139 Villadsen J. V. and W. E. Stewart. Solution of boundarvvalue problems by orthogonal collocation. Chem. Enngng Sci. 22, 148331501 (1967).