- Email: [email protected]
Fast method for solving nonequilibrium fixed-bed adsorption models with variable velocity and linear isotherm
0098-1354/92 $5.00+ 0.00 Copyright 0 1992 Pergamon Press Ltd
Computers them. Engng, Vol. 16, No. 6, pp. 535-544, 1992 Printed in Great Britain. All rights reserved
FAST METHOD FOR SOLVING NONEQUILIBRIUM FIXED-BED ADSORPTION MODELS WITH VARIABLE VELOCITY AND LINEAR ISOTHERM W. SUN Department (Received
and C. A. V. COSTAR
of Chemical Engineering, University of Porte, 4099 Porte Codex, Portugal
19 August 1991; final revision received 10 December 1991; received for publication 7 January 1992)
Abstract-The mass transfer equations for a fIxed-bed adsorber with linear isotherm and variable velocity are solved by combination of fast Fourier transform (FFT’) and method of lines techniques. The model considers axial dispersion and three diffusion steps in series: film, particle and microparticle. This is a three-region model where only the bulk phase mass transfer equations are nonlinear. The intraparticle diffusion equations are solved in the Laplace domain and this solution is the product of two functions: one is only dependent on model parameters and can be inverted by the FFT method and the other is just the time derivative of the axial bulk concentration. The mass transfer rate of the solute into the particle can then be calculated by the convolution integral of those two functions. Thus the three-region problem linked by boundary conditions is simplified to a single region problem. Finally the nonlinear one-region equations are solved by the method of lines using global orthogonal collocation for the axial coordinate. Several tests were done using a fairly large range of model parameters showing that: (a) mass balance relative errors are within 1%; (b) comparison between the solutions of the linear model (constant velocity) obtained with this method, by FFT and using a full orthogonal collocation solution shows that the relative error in the ordinates (outlet axial concentration) is always less than 0.5%; and (c) divergency or mass balance failure are observed when model parameters are such that a very steep axial profile develops. This failure is due to the use of global orthogonal collocation and thus a more suitable discretization method should be used in these conditions. Several plots display the histories of axial concentration and interstitial velocity showing the influence of model parameters and operating conditions (step, pulse and negative step inputs). The method of solution proposed in this paper is also suitable for solving other problems in various regions where the nonlinear equations only appear in some of them
three regions, rather than tridimensional, linked by boundary conditions. When the equilibrium isotherm is linear and the interstitial velocity is considered constant then the system transfer function can be obtained using Laplace transform methods and further inverted analytically (Rasmuson and Neretnieks, 1980; Rasmuson, 1982) or numerically by the FFT technique (Chen and Hsu, 1987). If the isotherm is nonlinear the usual way to solve the problem is the method of lines. The spatial variables are discretized by finite differences (Sun and Meunier, 1991), orthogonal collocation (Costa and Rodrigues, 1985; Raghavan and Ruthven, 1985) or finite elements (Sereno et al., 1991) and the resulting set of ordinary differential equations (ODES) solved by the Gear method (Hindmarsh, 1980). This was used to solve various adsorption models mostly in one region (axial coordinate) or two regions (axial coordinate + particle coordinate). The consideration of a complete equation for the third region (microparticle) will largely increase the number of ODES to be solved and thus the difficulty of the problem. That is why when by physical reasons it is important to consider microparticle diffusion,
INTRODUCTION
The important role of mathematical models to study, optimize and design fixed-bed adsorbers has stimulated the development of solution methods for the set of dynamic equations. Equilibrium models are represented by hyperbolic partial differential equations (PDEs); in the linear, semi-linear and quasi-linear cases an analytical or semi-analytical solution can be obtained by the method of characteristics (Rhee et al., 1986) whereas the nonlinear case requires a numerical solution (Loureiro and Rodrigues, 1991). The equilibrium models are useful for the qualitative understanding of the system behavior. The mathematical representation of a real fixedbed adsorption system requires that we take into account dispersive phenomena like axial dispersion and intraparticle mass transfer resistances. These models may have three spatial variables: axial coordinate of the fixed bed and particle and microparticle radial coordinates. The equations are in tTo whom all correspondence
should be addressed. 53s
W. SUN and C. A. V.
536
this mechanism is usually represented by a linear driving force approximation (Raghaven et al., 1985). In gas solid adsorption (solvent recovery, dehydration) the three-region model finds applications. In several cases the system shows two important characteristics: the equilibrium relationship is linear but the interstitial velocity is variable (LeVan et al., 1988) and thus the model is nonlinear. The objective of this wark is the development of a solution method for this case. We take advantage of the equation’s form, e.g. the nonlinear terms are only present in the bulk phase mass balance. The idea is to convert the three regions into one and thus dramatically reduce the problem dimensions and difficulty. This dimensionality reduction is accomplished by solving the linear equations in the Laplace domain and further obtaining the real-time solution by FFT. Then the problem is reduced to the region were the equations are nonlinear; these can now be solved by the method of lines using, for example, global orthogonal collocation for discretization of the spatial coordinate. This procedure saves computing time and problem solution difficulties when compared with the direct application of the method of lines to the three-region equations.
COSTA
regeneration
(negative step):
1 ac --
Pe
= C(O,@,
W
az Z_o
pulse: 1 ac -Pe az
= C(O,0) Z=o ac
r_l=o.
32 Macropore:
BC
s
(8)
= 0, Y-0
1, e)
qz,
-2% f
= c(z, e)
Y-l
.
(9)
Micropore:
MODEL Let us consider a fixed bed where solute adsorption from an inert is taking place. The adsorbent particles are spherical and made of an ensemble of spherical microparticles with pores amongst them. The diffusion path is as follows: film diffusion, macropore diffusion and micropore diffusion in series. The initial concentrations inside the particle and microparticle are constant and equal. The operation is isobaric and isothermal. The constitutive equations is dimensionless form are the following:
Bulk phase: cr-c-c,ac, c T
IC saturation
The dimensionless variables equations (1-l 3) are: C=(c
-co)/+,
Ci, = (Gincp
=
Cc,
-
N_D,r P- Rz’
(2)
C(Z, 0) = 0,
(3)
c, ay y_,’
(step input):
i ac -Pe az
I=
C(0, e) -
z_o
c,,
=
qz,
and
C,=c,Ic,,
cO)/cE
t
cO )lcE
3
Cpulsc
Ci
Z=z/L, z = LIuE,
=
(13)
parameters
in
c,=c,/c,,
=
Cci
Y, e).
(Cpu~u
-
-
cO )/CE
CO I/%
9
9
Y=RIRo, 9 = t/r.
(1)
ac,
344
(12)
y=l
-$=o. au -= dZ
=O,
c,tz,Y, i,e)
X= r/ro, aY
2 X-O
U=u/uE,
l?C ,+ug+344
BC
BC
(44
N-9 ‘-Ro’
where c, cE, co, cl,, cpubc. c,, c,, and ci are the concentrations of the solute in the bulk phase., reference, initial, inlet (step), inlet (pulse), total feed (constant), solute in the macropore and solute in the
Fast
method for solving adsorption models
micropore, respectively; II and U, are the interstitial velocity and a reference velocity, respectively; z is the fixed bed axial coordinate, R is the particle radial coordinate and r is the microparticle radial coordinate; L is the bed length, R, is the particle radius and r, is the microparticle radius; 6 is the dimensionless time and 0, is the dimensionless time length of the pulse; e is the external porosity and er is the macropore porosity; kr is the film mass transfer coefficient, D, is the macropore effective diffusivity, Q is the micropore effective diffusivity and D,, is the axial dispersion coefficient; K is the equilibrium constant. Equations (6) and (10) are linear parabolic equations with linear boundary conditions. The two nonlinear terms
and
537
Where:
s, =
NJ., N,+ N&'
S,=&coth(Jss)-1, s, =
s + 3drNiSr ’
NP S, = & Si =;(l
coth(JS,)
- 1,
+ K).
I
The function S, is independent of all spatial variables and depends only on mass transfer and equilibrium properties of the system. If the inverse function s5 of S, can be obtained then
can be calculated
by means of the convolution: 0
4(Z, 0) = only show up in equation (1). The first one is an implicit function of Z and 0 while the second is a function of the particle coordinates also. Thus if can be calculated as an explicit function acP/aYIr=l of Z and 0 only, the problem will be reduced to one region.
SOLUTION
METHOD
Let us take the Laplace (6-13)
transform
of equations
_ , (64
-3&,Niz
s
=
dY
@a)
0,
y=o
Pa)
e - A) dl.
(1%
Due to the complicated structure of S,, it is practically impossible to obtain an analytical inverse, s5. The direct use of numerical methods is also problematic because S, does not tend to zero when JsI increases unboundly. In order to meet this requirement S, has to bc transformed. Let G(s) = S,/s and F(s) = SC. The original function of F(s) is f(6) = aC/al?. The function G(s) tends zero when 1s 1 goes to infinity. Then the inverse of this function can be calculated using the Mellin’s formula (Bracewell, 1987): g(B) = 9-‘{G(s)}
r-1
>
s,(n)c(z,
s cl
= ;
a+rcOe”’ G(s) d.s, s n-k.2
(16)
where u is a real constant greater than the real parts of all function poles and that was set to zero here. The FFT technique (Hsu and Dranoff, 1987) can be used to approximate the Mellin’s integral and to obtain numerically g(e):
Y-l
(loa) dCi dX rco=O’ C,(Z, Y, 1, s) = Cr(Z, Y, s). This set of equations in the Laplace domain solved analytically to give:
dC, dY
=s,c. r-1
j = 1,2,. . . , N-l,
(17)
where
(124 (13a) can be
(14)
AT=%. T is the half period, 6 = j&T, s = ikla and N is the number of points used in FFT (should be a power of 2).
W. SUN and C. A. V. C~~TA
538
An important item to consider when inverting Laplace transforms by numerical methods is the rate at which the function in the Laplace domain tends to zero when IsI goes to infinity. This can be illustrated in the problem of reducing the Mellin’s integral to the Fourier integral. If we puts = iw in equation (16) we obtain the following complex Fourier integral: +a0 _ @‘G(h) s m
g(r)=&
do.
(18)
Most likely, the faster the function to be inverted decreases as 101 grows without bound, the more convenient the integral will be from a computational point of view and better accuracy will be achieved. We can attempt to accelerate the rate at which G(s) approaches zero when IsI goes to infinity, by decomposing this function into two summands: G(s)=G,(s)+G2(~).
+ G,(s),
(19a)
where D can be chosen as D-
limsG(s)=$. s-cc
P
N,N,S, - N;s N,(N,+ N,S,)s(l + s) ’
G,(s) tends to zero indeed faster than G, (8) by at least one order and can be efficiently inverted by FFT. The original function g(0) of G(s) is then:
g(B) = se-*
+ P-l{G2(s)].
(20)
c
The flux at the particle surface can then be obtained by solving:
q(Z, e) =
s@g(B 1
0
e
=
I
0
&de-1
The integral in equation the trapezoidal rule: 4(2,e)
= i yi’ J-0
After rearranging
(- . -.
-);
we get:
dC.
(21)
(21) can be approximated
{de -iae) + gre - (j +
n-2 +
by
l)ae]}
x {C[(j + l)AOl- WA@}.
(22)
C gj+l(Cm-j-Ca-j-d
j-0
9
1
(23)
where AB = O/n, gi = g(jA0) and C, = C(jAe) and n is the number of points considered. Let us analyze function g(B) in the vicinity of 0 = 0. The value of the function at that point can be calculated through g(O)=
Then:
G,(s) =
divident extrapolation (---); linear extrapolation 65,536 points (-).
(19)
These summands are chosen so that G,(s) tends to zero at the same rate as G(s) and G2(s) tends to zero faster than G(s). The function G(s) can be represented in the form: G(s) = &
Fig. 1. Calculation of g(0) in the vicinity of 0 = 0 (C,= 1, $ =O.S;r),=1;K= 100; N,= 1000;Np= 100;Ni= 10). 2048 points (0); g = 500/(1 + /XJ’R), c t .); Newton’s 3rd
limsG(s)=s. l-c0
P
Figure 1 shows the curve (full line) g(B) obtained with FF? (T = 50). The first point (0 = 0.001526) was obtained using N = 65,536. This calculation requires a lot of computational work and time and thus a less costly approximation is worthwhile. The curve g(0) can be calculated using FFT until a time that needs a reasonable N, say 2048 (0 = 0.048828 in the example shown) and then extrapolation to 0 = 0 or an asymptotic expression can be used. A third-order Newton divided difference formula (Kreyzig, 1988) can be used for extrapolation; in this case integration of equation (21) was done using the extrapolated g(B = 0) and g(B(N = 2048)) calculated by FFT; thus g(0) is approximated by a straight line (dashed line) between those two points. For small 0 an asymptotic expression can be derived (line * * .):
4
1
stel=Np 1+&/jj’ where B is a parameter calculated at 0(N = 2048). In this case integration between 0 = 0 and this point is analytical assuming that the concentration is piecewise linear in this interval.
Fast method for solving adsorption models To have an idea on the precision of these approximations the area under g(0) between 0 = 0 and B(N = 2048) can be calculated and compared. In the example of Fig. 1 those areas are 0.163 (solution by FFT), 0.198 (extrapolation) and 0.141 (asymptotic expression). The area for 6 > 8(N = 2048) is 3.59, and thus the differences between the results obtained by approximation and FFT are negligible. Note that if we take a straight line between 0(N = 2048) and g(0) = Nt/Nr, (line-.-.-) then the area under this line will be 12.3 meaning a gross error in the calculations with equation (21). Equation (23) can be substituted into equations (1) and (2) to yield:
c, - c - c, + 3bN, au az=
G -2
3bN c,
4(Z,6)=0
q(z, e).
(la)
(W
These equations are now written in one region. We can use orthogonal collocation (Villadsen and Michelsen, 1978) to discretize along the Z-direction transforming equation (la) into a set of ordinary differential equations. At a given 0 equation (2a) can be integrated numerically or analytically assuming that q(Z, 0) is piecewise linear in Z. RESULTS AND DISCUSSION
In order to validate the present method several tests were carried out with the following objectives: -tuning of the “numerical parameters”, e.g. the number of points in FFT, in order to obtain a convergent solution; -physical tests of this solution to ensure that this is the solution of the problem; -simulations to show the ability of the method to perform well in a fairly large range of model parameters and operating conditions. In the use of FFT two parameters must be set: the period 2T and the number of points N. The period is related to the time length of the solution; as an initial guess the period can be set as six times the stoichiometric time (8,, = 1 + $[l + &(l + K)]) and then adjusted to the minimum, The precision and computing time are related to N; large N means more precise
C(Z=l)=
539
results but more computing costs too. So a trade-off must be found and we advise 512 < N Q 2048. Of course the period should be as small as possible so that the same number of points give more precise results. In FFT calculations equation (17) is not used for the calculation of A(k) due to the presence of the infinite sum; instead a rounded-off approximation (Hsu and Dranoff, 1987) is preferred: A(+~$:+(k+;N)n],
k=O,l,...,N-1. (25)
The maximum number of points used in the integration of equation (21) by the trapezoidal rule [equation (23)] was determined by the value of N. All examples we ran show that using n = N gives enough precision. Probably n -c N can be used in some cases. The discretization of the axial coordinate can be done using one of the available tools. In this particular case we choose orthogonal collocation; this choice will be inappropriate, for example, in the case of high Pe numbers where plug flow is approached. Within the range of parameters we tested, 12-18 interior collocation points were enough. The integration of the resulting ODES was carried out by means of LSODE (Hindmarsh, 1980). The parameters of this package were set according to the instructions, i.e. stiff and numerical Jacobian options, relative error 5 x 10e3 and absolute error 5 x lo-‘. The integration proceeds step-by-step so that the convolution integral can be evaluated at each time step. Two “physical” tests were also done; firstly we must be sure that the numerical solution does not violate the global material balance and secondly that it does not add numerical dispersion relatively to the true solution. The step responses were used to carry out these tests. The relative error in the global mass balance was always inferior to 1% for the range of parameters we tested. The dispersion test needs a comparison with a known solution; for this purpose the related linear model (constant interstitial velocity) can be used. The dimensionless mass balance for the bulk phase is: %+&&$+34N,q(Z,8)=0.
(26)
The solution of equations (26) and (3-13) for a step input is, in the Laplace domain:
exp(Pe) (Y2- yi ) sI~~expQ,)(yrlPel)-y2cxp(72)(y1/Pe-
VI
(27)
W. SUN and C. A. V. COSTA
540
times corresponding to situations of high K or low intraparticle resistances, i.e. steeper profiles. The number of ODES to be solved in the of a full collocation discretization is case Narial 11+ N-r,, (1 + N,,,&] as compared with only Nti, for the combined solution. Raghavan and Ruthven (1985) proposed 2-6 collocation points for
where ?I,2 = f[Pe AZJPe’
+ 4Pe(s + 34N,S,)l.
It is difficult to find the original of equation (27) analytically and so FFT will be used. For this purpose the behavior of equation (27) can be improved by removing the singular point:
each intraparticle coordinate which means that theirsolution method generates 7-43 times more ODES than the method proposed in this work. This better performance of the combined method is thus likely to hold in the nonlinear case. Several simulations for the step input were carried out in the following range of parameters:
The original of the second term is 1 and FFT is applied to the first term. The comparison between the results obtained when the linear model is solved by FFT and by the method we proposed shows that the relative error in the ordinates (concentration) is always less than 0.5%. Raghaven and Ruthven (1985) solved this linear case using a full orthogonal collocation solution. In order to compare this method of solution with the one we are proposing, several simulations were ran for a fairly large range of model parameters and operating conditions. In those simulations we used 14 interior collocation points for the axial coordinate and four interior collocation points for the macroand microparticle radial coordinates in the case of the full collocation solution and 14 interior collocation points for the axial coordinate and 1024 points for the FFT in the case of the combined FFT-collocation solution. In both cases the ODES were solved using LSODE (Hindmarsh, 1980). There is practically no difference among the solutions obtained by direct use of FFT, full orthogonal collocation and the combined method. Figure 2 displays a typical case. The average CPU time for the combined solution is about 1 min in an IBM Risk 6000, whereas the solution using full orthogonal collocation takes about 3-15 times more in the same computer; higher CPU
c,=4;
c,=o;
C,=l;
9 =0.5;
&=
1,
K= l-200, Pe = 50-200, N, = 500-5000, Np = 0.1-100, Ni = 0.01-10, To give an idea of the physical meaning of these ranges we can say that on average for a molecular sieve, K > 10, Pe > 50, Nr = 5000, Np = lo-100 and Ni=O.l-1. Some of the simulations are shown in Figs 37 where the influence of K, Pe, Nr, Nr and Ni are depicted. These plots show that velocity changes are simultaneous with concentration changes, which agrees with equilibrium theory predictions (LeVan et al., 1988). Solute uptake during adsorption implies that the observed outlet flowrate is initially inferior to the feed flowrate; inlet and outlet flowrates only
CJ
::ii0
20
40
50
50
100
120
140
160
e
Fig. 2. Comparison among the solutions of the linear model obtained by full orthogonal collocation discretization (upper curve), direct inversion by FFT (middle curve) and by the combined method (lower curve). K = 80, Pe = 200, N,=~OOO,
N,=lO, N,=l.
0
50
100
150
200
250
e '
Fig. 3. InfIuence. of K in step responses (C,-; C,=l,Q
so.5;
&=
1; c,=4;
Co=Oh
U,---;
Fast method for solving adsorption
models
o.oV’ 0
541
I 40
80
120
160
200
0
Fig. 4. Influence of Pe in step responses (C, -; (I,---; c,, = I, f$ = 0.5; q$ = 1; CT = 4; c, = 0).
Fig. 6.
become the same when complete saturation of fixed bed occurs. In all these simulations an initially slight oscillation in the velocity profiles is observed due to the steepness of these profiles in early times. Figure 3 shows the influence of K; increasing K means increasing the equilibrium capacity of the fixed bed and thus the center of gravity of the concentration histories moves to longer times. Figures 4 and 5 show that in the range considered axial dispersion and film mass transfer have almost no effect on breakthrough curves. In Figs 6 and 7 the influence of intraparticle mass transfer resistances is depicted; higher resistances correspond to a smaller operating capacity of the fixed bed. Figures 8 and 9 show the combined influence of N,, Np and Ni and of Pe, Nt, and N,, respectively. These plots agree with the expected behavior of the system and add information on the way interstitial velocity changes at the fixed bed outlet. Whenever the axial concentration profile becomes very steep the method fails, i.e. divergency or material balance failure are observed. This happens when Pe 3 200 and the resistances to mass transfer are negligible or when adsorption capacity is small (K < 5). This is due to the use of global orthogonal collocation to discretize the axial coordinate. In these cases the use of a better suited discretization method, e.g. moving finite elements, will eliminate the prob-
lems. Other operating conditions were also tested. Figure 10 shows the results obtained for a pulse input 1) when the pulse time length (C,=4;C,=O;C,,,= is varied. In Fig. 11 the response to a negative step (Cr = 4; C, = 1) is plotted. The analysis of the simulations shows that the proposed method for reducing the number of spatial regions works well in a fairly large range of parameters.
of Np in step responses (C, --; U,- - -; C,” = I,4 = 0.5; & = 1; CT = 4; c, = 0).
Influence
CONCLUSIONS
A method is proposed for the solution of threeregion (microparticle, particle and bed) nonequilibrium fixed bed adsorption models with linear equilibrium and variable interstitial velocity. The method basis is the reduction of problem dimensions by using FFT inversion of Laplace domain functions in linear regions. The following algorithm is proposed for the solution of that problem: -solve the linear equations in the Laplace domain; -apply FFT to get the solution in the time domain; this implies: l
l
eventual transformation of the function to be inverted in order to meet FFT requirements, choice of a suitable half-period T and number of points N, 1.0
0.8
w
CJJ
0.6
0.4
0.2
0.0 0
40
80
120
180
200
e
Fig. 5. Influence of Nf in step responses (C, -; ci,= l,$ =0.5; 6,= I; c,=4; C,=O).
U,-- -;
Fig. 7. Influence of N, in step responses (C, -; U,-- -; c, = 1, f#l= 0.5; & = 1; CT = 4; c, = 0).
W. SUN and C. A. V. COSTA
542 1.0
KIMi P*=ZOO
0.8 s-o_-
Nf=5OO,Np=l,Nl=~
0.6 --+-
Nf=SOOO,
--C--
Nf=SOO,
Np=l, Nlml
c,u 0.4
0.2
---o--
Nf=SOOO.
Fig. 8. Combined influence of Nr, NP and N, (C, --;
-by means of a convolution calculate the desired real-time function to he used in the nonlinear equations; --solve the nonlinear one-region equations for example, the method of lines.
using,
Np’lOO.
U,---;
Ni=lO
NpolOO.Nl=lO
Ci, = I, r$ = 0.5; I&,= 1; C, = 4; C,, = 0).
This algorithm was tested using orthogonal collocation for the discretization of the 6xed-bed axial coordinate. These tests were performed over a large range of model parameters showing that: -mass
balance relative errors are within 1%;
1.0
0.0
0.6
Nf =5000
--“-
Pe=SO.Np=
--“-
P.130,
1. Nid
c,u 0.4
0.2
--&-
_
Np’lOO,
Pp=200,Np=100.
Ni=lO
Nl=10
0
Fig. 9. Combined influence of Fk, NP and Ni (C, -;
V,---;
Cb= 1, q5= 0.5; 4, = 1; CT = 4; C, = 0).
c.u
K=60 0.6
P6= 200 Nt = 6000 Np=l Nl-1
0
40
120
60
160
200
0
Fig. 10. Response to a pulse input. Influence of the pulse time length (C, -; &= 1; c,=4; C,=O). --comparison between the solutions of the linear model (constant velocity) obtained with this method and by FFT shows that the relative error in the ordinates (outlet solute axial concentration) is always less than 0.5%; -divergency or mass balance failure are observed when model parameters are such that a very steep axial profile develops. This is due to the use of global orthogonal collocation and thus a more suitable discretization method should be used in these conditions. It is worthwhile to mention that the use of this method saves a lot of computing time and numerical difficulties when compared to the direct use of the method of lines to solve this type of problem. 31 9.l2
This solution method is not limited to the class of problems treated. We think that it can be useful in the numerical solution of multiregion models where one or more regions are linear.
Acknowledgements-JNICT is gratefully acknowledged for financial support of this work. W. Sun also acknowledges Fundafio
c= e, = ca = ci = c, = cp = cpk = c, = C=
Cb = C, =
grant.
Solute concentration in bulk phase, gmol cm-) Initial concentration, gmol en-j Solute reference concentration, gmol cn+ Solute concentration in micropore, gmol cnd3 Inlet concentration (step input), gmol cnm3 Solute concentration in macropore, gmol crne3 Inlet concentration (Pulse), gmol cm-) Total concentration of solute and inert components in feed, gmol cmm3 Dimensionless solute concentration in bulk phase .Dimensionless initial concentration Dimensionless solute concentration in micropore Dimensionless inlet concentration (step input) Dimensionless solute concentration in macropore
CPvlr= Dimensionless inlet concentration (pulse) C, = Dimensionless total concentration of solute and inert components in bulk phase C = Dimensionless concentration variable of bulk phase cl
30
SO
60
120
160
El
Fig. 11. Response of a negative (C, -;
Oriente for a research
NOMENCLATURE
C, = %.__ *----..______ ------_-----____-________._________ Ci =
c.u
U,- - -; CPti = 1,4 = 0.5;
step of a presaturated bed (I.---; q5=OS; I$~= 1; CT = 4; C, = 0). K = 80, Pe=200. Nr=SODO, N,,= 1, Ni=l.
in Laplace domain Ci = Dimensionless concentration in Laplace domain
variable of micropore
C, = Dimensionless concentration variable of macropore in Laplace domain D, = Axial dispersion coeflicient, cm2s-l Di = Effective micropore diffusivity, cm2 S-I Dp = E&dive macropore diffusivity, cm2 s-l
W. SUN and C. A. V.
544 k, K L N Nr N,
= = = = = =
Np = Pe = r= R = r, = R, = g= t= T= U= u= us = X= Y= 2 =
External film mass transfer coefficient, cm s-’ Dimensionless adsorption equilibrium constant Bed length, cm Number of points used in FFT Dimensionless film mass transfer parameter Dimensionless mass transfer parameter im micropore Dimensionless mass transfer parameter in macropore Peclet number Radial coordinate for the microparticle, cm Radial coordinate for the pellet, cm Microparticle radius, cm Pellet radius, cm Dimensionless solute flux into the pellet Time, s Half operation period, s Dimensionless velocity Interstitial velocity, cm s-’ Reference interstitial velocity, cm s-’ Dimensionless microparticle radial coordinate Dimensionless pellet radial coordinate Dimensionless axial coordinate in the bed
Greek symbols a= /3 = L= ep = 4 = & = r = 0= w=
Constant of Mellin’s formula Parameter in equation (24) Bulk phase porosity, cm3 bed void/cm) bed Macropore porosity, cm’ macropores crnM3 pellet Dimensionless voidage parameter Dimensionless macropore voidage parameter Space time Dimensionless time variable Fourier integral variable
REFERENCES Bracewell R. N., The Fourier Transform and Its Applications, 3rd Edn. McGraw-Hill, New York (1987). Chen T. L. and J. T. Hsu, Prediction of breakthrough curves by the application of fast Fourier transform. AIChE J133, 1387-1390 (1987).
a,TA
Design of cyclic fixed bed Costa C. and A. Rodrigues, processes I. AIChE JI 31, 1655-1665 adsorption (1985). Hindmarsh A. C., LSODE and LSODI two new initial value ordinary differential equation solvers. ACM-SIGNUM News/err. 15, lo-11 (1980). Hsu J. T. and J. S. Dranoff. Numerical inversion of certain. Laplace transforms by ‘the direct application of fast Fourier transform (FFT) algorithm. Computers them. Engng II, 101-110 (1987). Kreyszig E., Advanced Engineering Mathematics, 6th Edn. Wiley, New York (1988). LeVan M. D. et al. Fixed bed adsorption of gases: effect of velocity variations on transition types. AIChE JI 34, 996-1005 (1988). Loureiro J. M. and A. Rodrigues, Two solution methods for hyperbolic systems of partial differential equations in chemical ,engineering. Chem. Engng Sci. 46, (1991). in press. Raghavan N. S. and D. M. Ruthven, Simulation of chromatographic response in columns packed with bidisperse structured particles. Chem. Engng Sci. 40,699-706 (1985). Raghavan N. S., M. M. Hassan and D. M. Ruthven, Numerical Simulation of a PSA system, Part I. AIChE Jl 31, 385-392 (1985). Rasmuson A., Time domain solution of a model for transport processes in bidisperse structured catalyst. Chem. /&pig Sci. 37, 787-788 (1982). Rasmuson A. and I. Neretnieks, Exact solution of a model for diffusion in particles and longitudinal dispersion in oacked beds. AIChE JI 26. 686-690 (1980). Rhee H. K., R. Aris and N. R. Amundson, First Order Partial Dtfirential Equations. Vol 1. Prentice-Hall, Englewood Cliffs, NJ (1986). Sereno C., A. Rodrigues and J. Villadsen, The moving finite element method with polynomial approximation of any degree. Comtmters them. Ennnn 15, 25-33 (1991). Sun L. M. and F. Meunier, An improved finite difference method for fixed bed multicomponent sorption. AIChE JI 37, 244-254 (1991). Villadsen J. and M. L. Michelson. Solution of Differential Eqaations Models by Polynomial Approximaiion. &enticeHall, Englewood Cliffs, NJ (1978).
Computers them. Engng, Vol. 16, No. 6, pp. 535-544, 1992 Printed in Great Britain. All rights reserved
FAST METHOD FOR SOLVING NONEQUILIBRIUM FIXED-BED ADSORPTION MODELS WITH VARIABLE VELOCITY AND LINEAR ISOTHERM W. SUN Department (Received
and C. A. V. COSTAR
of Chemical Engineering, University of Porte, 4099 Porte Codex, Portugal
19 August 1991; final revision received 10 December 1991; received for publication 7 January 1992)
Abstract-The mass transfer equations for a fIxed-bed adsorber with linear isotherm and variable velocity are solved by combination of fast Fourier transform (FFT’) and method of lines techniques. The model considers axial dispersion and three diffusion steps in series: film, particle and microparticle. This is a three-region model where only the bulk phase mass transfer equations are nonlinear. The intraparticle diffusion equations are solved in the Laplace domain and this solution is the product of two functions: one is only dependent on model parameters and can be inverted by the FFT method and the other is just the time derivative of the axial bulk concentration. The mass transfer rate of the solute into the particle can then be calculated by the convolution integral of those two functions. Thus the three-region problem linked by boundary conditions is simplified to a single region problem. Finally the nonlinear one-region equations are solved by the method of lines using global orthogonal collocation for the axial coordinate. Several tests were done using a fairly large range of model parameters showing that: (a) mass balance relative errors are within 1%; (b) comparison between the solutions of the linear model (constant velocity) obtained with this method, by FFT and using a full orthogonal collocation solution shows that the relative error in the ordinates (outlet axial concentration) is always less than 0.5%; and (c) divergency or mass balance failure are observed when model parameters are such that a very steep axial profile develops. This failure is due to the use of global orthogonal collocation and thus a more suitable discretization method should be used in these conditions. Several plots display the histories of axial concentration and interstitial velocity showing the influence of model parameters and operating conditions (step, pulse and negative step inputs). The method of solution proposed in this paper is also suitable for solving other problems in various regions where the nonlinear equations only appear in some of them
three regions, rather than tridimensional, linked by boundary conditions. When the equilibrium isotherm is linear and the interstitial velocity is considered constant then the system transfer function can be obtained using Laplace transform methods and further inverted analytically (Rasmuson and Neretnieks, 1980; Rasmuson, 1982) or numerically by the FFT technique (Chen and Hsu, 1987). If the isotherm is nonlinear the usual way to solve the problem is the method of lines. The spatial variables are discretized by finite differences (Sun and Meunier, 1991), orthogonal collocation (Costa and Rodrigues, 1985; Raghavan and Ruthven, 1985) or finite elements (Sereno et al., 1991) and the resulting set of ordinary differential equations (ODES) solved by the Gear method (Hindmarsh, 1980). This was used to solve various adsorption models mostly in one region (axial coordinate) or two regions (axial coordinate + particle coordinate). The consideration of a complete equation for the third region (microparticle) will largely increase the number of ODES to be solved and thus the difficulty of the problem. That is why when by physical reasons it is important to consider microparticle diffusion,
INTRODUCTION
The important role of mathematical models to study, optimize and design fixed-bed adsorbers has stimulated the development of solution methods for the set of dynamic equations. Equilibrium models are represented by hyperbolic partial differential equations (PDEs); in the linear, semi-linear and quasi-linear cases an analytical or semi-analytical solution can be obtained by the method of characteristics (Rhee et al., 1986) whereas the nonlinear case requires a numerical solution (Loureiro and Rodrigues, 1991). The equilibrium models are useful for the qualitative understanding of the system behavior. The mathematical representation of a real fixedbed adsorption system requires that we take into account dispersive phenomena like axial dispersion and intraparticle mass transfer resistances. These models may have three spatial variables: axial coordinate of the fixed bed and particle and microparticle radial coordinates. The equations are in tTo whom all correspondence
should be addressed. 53s
W. SUN and C. A. V.
536
this mechanism is usually represented by a linear driving force approximation (Raghaven et al., 1985). In gas solid adsorption (solvent recovery, dehydration) the three-region model finds applications. In several cases the system shows two important characteristics: the equilibrium relationship is linear but the interstitial velocity is variable (LeVan et al., 1988) and thus the model is nonlinear. The objective of this wark is the development of a solution method for this case. We take advantage of the equation’s form, e.g. the nonlinear terms are only present in the bulk phase mass balance. The idea is to convert the three regions into one and thus dramatically reduce the problem dimensions and difficulty. This dimensionality reduction is accomplished by solving the linear equations in the Laplace domain and further obtaining the real-time solution by FFT. Then the problem is reduced to the region were the equations are nonlinear; these can now be solved by the method of lines using, for example, global orthogonal collocation for discretization of the spatial coordinate. This procedure saves computing time and problem solution difficulties when compared with the direct application of the method of lines to the three-region equations.
COSTA
regeneration
(negative step):
1 ac --
Pe
= C(O,@,
W
az Z_o
pulse: 1 ac -Pe az
= C(O,0) Z=o ac
r_l=o.
32 Macropore:
BC
s
(8)
= 0, Y-0
1, e)
qz,
-2% f
= c(z, e)
Y-l
.
(9)
Micropore:
MODEL Let us consider a fixed bed where solute adsorption from an inert is taking place. The adsorbent particles are spherical and made of an ensemble of spherical microparticles with pores amongst them. The diffusion path is as follows: film diffusion, macropore diffusion and micropore diffusion in series. The initial concentrations inside the particle and microparticle are constant and equal. The operation is isobaric and isothermal. The constitutive equations is dimensionless form are the following:
Bulk phase: cr-c-c,ac, c T
IC saturation
The dimensionless variables equations (1-l 3) are: C=(c
-co)/+,
Ci, = (Gincp
=
Cc,
-
N_D,r P- Rz’
(2)
C(Z, 0) = 0,
(3)
c, ay y_,’
(step input):
i ac -Pe az
I=
C(0, e) -
z_o
c,,
=
qz,
and
C,=c,Ic,,
cO)/cE
t
cO )lcE
3
Cpulsc
Ci
Z=z/L, z = LIuE,
=
(13)
parameters
in
c,=c,/c,,
=
Cci
Y, e).
(Cpu~u
-
-
cO )/CE
CO I/%
9
9
Y=RIRo, 9 = t/r.
(1)
ac,
344
(12)
y=l
-$=o. au -= dZ
=O,
c,tz,Y, i,e)
X= r/ro, aY
2 X-O
U=u/uE,
l?C ,+ug+344
BC
BC
(44
N-9 ‘-Ro’
where c, cE, co, cl,, cpubc. c,, c,, and ci are the concentrations of the solute in the bulk phase., reference, initial, inlet (step), inlet (pulse), total feed (constant), solute in the macropore and solute in the
Fast
method for solving adsorption models
micropore, respectively; II and U, are the interstitial velocity and a reference velocity, respectively; z is the fixed bed axial coordinate, R is the particle radial coordinate and r is the microparticle radial coordinate; L is the bed length, R, is the particle radius and r, is the microparticle radius; 6 is the dimensionless time and 0, is the dimensionless time length of the pulse; e is the external porosity and er is the macropore porosity; kr is the film mass transfer coefficient, D, is the macropore effective diffusivity, Q is the micropore effective diffusivity and D,, is the axial dispersion coefficient; K is the equilibrium constant. Equations (6) and (10) are linear parabolic equations with linear boundary conditions. The two nonlinear terms
and
537
Where:
s, =
NJ., N,+ N&'
S,=&coth(Jss)-1, s, =
s + 3drNiSr ’
NP S, = & Si =;(l
coth(JS,)
- 1,
+ K).
I
The function S, is independent of all spatial variables and depends only on mass transfer and equilibrium properties of the system. If the inverse function s5 of S, can be obtained then
can be calculated
by means of the convolution: 0
4(Z, 0) = only show up in equation (1). The first one is an implicit function of Z and 0 while the second is a function of the particle coordinates also. Thus if can be calculated as an explicit function acP/aYIr=l of Z and 0 only, the problem will be reduced to one region.
SOLUTION
METHOD
Let us take the Laplace (6-13)
transform
of equations
_ , (64
-3&,Niz
s
=
dY
@a)
0,
y=o
Pa)
e - A) dl.
(1%
Due to the complicated structure of S,, it is practically impossible to obtain an analytical inverse, s5. The direct use of numerical methods is also problematic because S, does not tend to zero when JsI increases unboundly. In order to meet this requirement S, has to bc transformed. Let G(s) = S,/s and F(s) = SC. The original function of F(s) is f(6) = aC/al?. The function G(s) tends zero when 1s 1 goes to infinity. Then the inverse of this function can be calculated using the Mellin’s formula (Bracewell, 1987): g(B) = 9-‘{G(s)}
r-1
>
s,(n)c(z,
s cl
= ;
a+rcOe”’ G(s) d.s, s n-k.2
(16)
where u is a real constant greater than the real parts of all function poles and that was set to zero here. The FFT technique (Hsu and Dranoff, 1987) can be used to approximate the Mellin’s integral and to obtain numerically g(e):
Y-l
(loa) dCi dX rco=O’ C,(Z, Y, 1, s) = Cr(Z, Y, s). This set of equations in the Laplace domain solved analytically to give:
dC, dY
=s,c. r-1
j = 1,2,. . . , N-l,
(17)
where
(124 (13a) can be
(14)
AT=%. T is the half period, 6 = j&T, s = ikla and N is the number of points used in FFT (should be a power of 2).
W. SUN and C. A. V. C~~TA
538
An important item to consider when inverting Laplace transforms by numerical methods is the rate at which the function in the Laplace domain tends to zero when IsI goes to infinity. This can be illustrated in the problem of reducing the Mellin’s integral to the Fourier integral. If we puts = iw in equation (16) we obtain the following complex Fourier integral: +a0 _ @‘G(h) s m
g(r)=&
do.
(18)
Most likely, the faster the function to be inverted decreases as 101 grows without bound, the more convenient the integral will be from a computational point of view and better accuracy will be achieved. We can attempt to accelerate the rate at which G(s) approaches zero when IsI goes to infinity, by decomposing this function into two summands: G(s)=G,(s)+G2(~).
+ G,(s),
(19a)
where D can be chosen as D-
limsG(s)=$. s-cc
P
N,N,S, - N;s N,(N,+ N,S,)s(l + s) ’
G,(s) tends to zero indeed faster than G, (8) by at least one order and can be efficiently inverted by FFT. The original function g(0) of G(s) is then:
g(B) = se-*
+ P-l{G2(s)].
(20)
c
The flux at the particle surface can then be obtained by solving:
q(Z, e) =
s@g(B 1
0
e
=
I
0
&de-1
The integral in equation the trapezoidal rule: 4(2,e)
= i yi’ J-0
After rearranging
(- . -.
-);
we get:
dC.
(21)
(21) can be approximated
{de -iae) + gre - (j +
n-2 +
by
l)ae]}
x {C[(j + l)AOl- WA@}.
(22)
C gj+l(Cm-j-Ca-j-d
j-0
9
1
(23)
where AB = O/n, gi = g(jA0) and C, = C(jAe) and n is the number of points considered. Let us analyze function g(B) in the vicinity of 0 = 0. The value of the function at that point can be calculated through g(O)=
Then:
G,(s) =
divident extrapolation (---); linear extrapolation 65,536 points (-).
(19)
These summands are chosen so that G,(s) tends to zero at the same rate as G(s) and G2(s) tends to zero faster than G(s). The function G(s) can be represented in the form: G(s) = &
Fig. 1. Calculation of g(0) in the vicinity of 0 = 0 (C,= 1, $ =O.S;r),=1;K= 100; N,= 1000;Np= 100;Ni= 10). 2048 points (0); g = 500/(1 + /XJ’R), c t .); Newton’s 3rd
limsG(s)=s. l-c0
P
Figure 1 shows the curve (full line) g(B) obtained with FF? (T = 50). The first point (0 = 0.001526) was obtained using N = 65,536. This calculation requires a lot of computational work and time and thus a less costly approximation is worthwhile. The curve g(0) can be calculated using FFT until a time that needs a reasonable N, say 2048 (0 = 0.048828 in the example shown) and then extrapolation to 0 = 0 or an asymptotic expression can be used. A third-order Newton divided difference formula (Kreyzig, 1988) can be used for extrapolation; in this case integration of equation (21) was done using the extrapolated g(B = 0) and g(B(N = 2048)) calculated by FFT; thus g(0) is approximated by a straight line (dashed line) between those two points. For small 0 an asymptotic expression can be derived (line * * .):
4
1
stel=Np 1+&/jj’ where B is a parameter calculated at 0(N = 2048). In this case integration between 0 = 0 and this point is analytical assuming that the concentration is piecewise linear in this interval.
Fast method for solving adsorption models To have an idea on the precision of these approximations the area under g(0) between 0 = 0 and B(N = 2048) can be calculated and compared. In the example of Fig. 1 those areas are 0.163 (solution by FFT), 0.198 (extrapolation) and 0.141 (asymptotic expression). The area for 6 > 8(N = 2048) is 3.59, and thus the differences between the results obtained by approximation and FFT are negligible. Note that if we take a straight line between 0(N = 2048) and g(0) = Nt/Nr, (line-.-.-) then the area under this line will be 12.3 meaning a gross error in the calculations with equation (21). Equation (23) can be substituted into equations (1) and (2) to yield:
c, - c - c, + 3bN, au az=
G -2
3bN c,
4(Z,6)=0
q(z, e).
(la)
(W
These equations are now written in one region. We can use orthogonal collocation (Villadsen and Michelsen, 1978) to discretize along the Z-direction transforming equation (la) into a set of ordinary differential equations. At a given 0 equation (2a) can be integrated numerically or analytically assuming that q(Z, 0) is piecewise linear in Z. RESULTS AND DISCUSSION
In order to validate the present method several tests were carried out with the following objectives: -tuning of the “numerical parameters”, e.g. the number of points in FFT, in order to obtain a convergent solution; -physical tests of this solution to ensure that this is the solution of the problem; -simulations to show the ability of the method to perform well in a fairly large range of model parameters and operating conditions. In the use of FFT two parameters must be set: the period 2T and the number of points N. The period is related to the time length of the solution; as an initial guess the period can be set as six times the stoichiometric time (8,, = 1 + $[l + &(l + K)]) and then adjusted to the minimum, The precision and computing time are related to N; large N means more precise
C(Z=l)=
539
results but more computing costs too. So a trade-off must be found and we advise 512 < N Q 2048. Of course the period should be as small as possible so that the same number of points give more precise results. In FFT calculations equation (17) is not used for the calculation of A(k) due to the presence of the infinite sum; instead a rounded-off approximation (Hsu and Dranoff, 1987) is preferred: A(+~$:+(k+;N)n],
k=O,l,...,N-1. (25)
The maximum number of points used in the integration of equation (21) by the trapezoidal rule [equation (23)] was determined by the value of N. All examples we ran show that using n = N gives enough precision. Probably n -c N can be used in some cases. The discretization of the axial coordinate can be done using one of the available tools. In this particular case we choose orthogonal collocation; this choice will be inappropriate, for example, in the case of high Pe numbers where plug flow is approached. Within the range of parameters we tested, 12-18 interior collocation points were enough. The integration of the resulting ODES was carried out by means of LSODE (Hindmarsh, 1980). The parameters of this package were set according to the instructions, i.e. stiff and numerical Jacobian options, relative error 5 x 10e3 and absolute error 5 x lo-‘. The integration proceeds step-by-step so that the convolution integral can be evaluated at each time step. Two “physical” tests were also done; firstly we must be sure that the numerical solution does not violate the global material balance and secondly that it does not add numerical dispersion relatively to the true solution. The step responses were used to carry out these tests. The relative error in the global mass balance was always inferior to 1% for the range of parameters we tested. The dispersion test needs a comparison with a known solution; for this purpose the related linear model (constant interstitial velocity) can be used. The dimensionless mass balance for the bulk phase is: %+&&$+34N,q(Z,8)=0.
(26)
The solution of equations (26) and (3-13) for a step input is, in the Laplace domain:
exp(Pe) (Y2- yi ) sI~~expQ,)(yrlPel)-y2cxp(72)(y1/Pe-
VI
(27)
W. SUN and C. A. V. COSTA
540
times corresponding to situations of high K or low intraparticle resistances, i.e. steeper profiles. The number of ODES to be solved in the of a full collocation discretization is case Narial 11+ N-r,, (1 + N,,,&] as compared with only Nti, for the combined solution. Raghavan and Ruthven (1985) proposed 2-6 collocation points for
where ?I,2 = f[Pe AZJPe’
+ 4Pe(s + 34N,S,)l.
It is difficult to find the original of equation (27) analytically and so FFT will be used. For this purpose the behavior of equation (27) can be improved by removing the singular point:
each intraparticle coordinate which means that theirsolution method generates 7-43 times more ODES than the method proposed in this work. This better performance of the combined method is thus likely to hold in the nonlinear case. Several simulations for the step input were carried out in the following range of parameters:
The original of the second term is 1 and FFT is applied to the first term. The comparison between the results obtained when the linear model is solved by FFT and by the method we proposed shows that the relative error in the ordinates (concentration) is always less than 0.5%. Raghaven and Ruthven (1985) solved this linear case using a full orthogonal collocation solution. In order to compare this method of solution with the one we are proposing, several simulations were ran for a fairly large range of model parameters and operating conditions. In those simulations we used 14 interior collocation points for the axial coordinate and four interior collocation points for the macroand microparticle radial coordinates in the case of the full collocation solution and 14 interior collocation points for the axial coordinate and 1024 points for the FFT in the case of the combined FFT-collocation solution. In both cases the ODES were solved using LSODE (Hindmarsh, 1980). There is practically no difference among the solutions obtained by direct use of FFT, full orthogonal collocation and the combined method. Figure 2 displays a typical case. The average CPU time for the combined solution is about 1 min in an IBM Risk 6000, whereas the solution using full orthogonal collocation takes about 3-15 times more in the same computer; higher CPU
c,=4;
c,=o;
C,=l;
9 =0.5;
&=
1,
K= l-200, Pe = 50-200, N, = 500-5000, Np = 0.1-100, Ni = 0.01-10, To give an idea of the physical meaning of these ranges we can say that on average for a molecular sieve, K > 10, Pe > 50, Nr = 5000, Np = lo-100 and Ni=O.l-1. Some of the simulations are shown in Figs 37 where the influence of K, Pe, Nr, Nr and Ni are depicted. These plots show that velocity changes are simultaneous with concentration changes, which agrees with equilibrium theory predictions (LeVan et al., 1988). Solute uptake during adsorption implies that the observed outlet flowrate is initially inferior to the feed flowrate; inlet and outlet flowrates only
CJ
::ii0
20
40
50
50
100
120
140
160
e
Fig. 2. Comparison among the solutions of the linear model obtained by full orthogonal collocation discretization (upper curve), direct inversion by FFT (middle curve) and by the combined method (lower curve). K = 80, Pe = 200, N,=~OOO,
N,=lO, N,=l.
0
50
100
150
200
250
e '
Fig. 3. InfIuence. of K in step responses (C,-; C,=l,Q
so.5;
&=
1; c,=4;
Co=Oh
U,---;
Fast method for solving adsorption
models
o.oV’ 0
541
I 40
80
120
160
200
0
Fig. 4. Influence of Pe in step responses (C, -; (I,---; c,, = I, f$ = 0.5; q$ = 1; CT = 4; c, = 0).
Fig. 6.
become the same when complete saturation of fixed bed occurs. In all these simulations an initially slight oscillation in the velocity profiles is observed due to the steepness of these profiles in early times. Figure 3 shows the influence of K; increasing K means increasing the equilibrium capacity of the fixed bed and thus the center of gravity of the concentration histories moves to longer times. Figures 4 and 5 show that in the range considered axial dispersion and film mass transfer have almost no effect on breakthrough curves. In Figs 6 and 7 the influence of intraparticle mass transfer resistances is depicted; higher resistances correspond to a smaller operating capacity of the fixed bed. Figures 8 and 9 show the combined influence of N,, Np and Ni and of Pe, Nt, and N,, respectively. These plots agree with the expected behavior of the system and add information on the way interstitial velocity changes at the fixed bed outlet. Whenever the axial concentration profile becomes very steep the method fails, i.e. divergency or material balance failure are observed. This happens when Pe 3 200 and the resistances to mass transfer are negligible or when adsorption capacity is small (K < 5). This is due to the use of global orthogonal collocation to discretize the axial coordinate. In these cases the use of a better suited discretization method, e.g. moving finite elements, will eliminate the prob-
lems. Other operating conditions were also tested. Figure 10 shows the results obtained for a pulse input 1) when the pulse time length (C,=4;C,=O;C,,,= is varied. In Fig. 11 the response to a negative step (Cr = 4; C, = 1) is plotted. The analysis of the simulations shows that the proposed method for reducing the number of spatial regions works well in a fairly large range of parameters.
of Np in step responses (C, --; U,- - -; C,” = I,4 = 0.5; & = 1; CT = 4; c, = 0).
Influence
CONCLUSIONS
A method is proposed for the solution of threeregion (microparticle, particle and bed) nonequilibrium fixed bed adsorption models with linear equilibrium and variable interstitial velocity. The method basis is the reduction of problem dimensions by using FFT inversion of Laplace domain functions in linear regions. The following algorithm is proposed for the solution of that problem: -solve the linear equations in the Laplace domain; -apply FFT to get the solution in the time domain; this implies: l
l
eventual transformation of the function to be inverted in order to meet FFT requirements, choice of a suitable half-period T and number of points N, 1.0
0.8
w
CJJ
0.6
0.4
0.2
0.0 0
40
80
120
180
200
e
Fig. 5. Influence of Nf in step responses (C, -; ci,= l,$ =0.5; 6,= I; c,=4; C,=O).
U,-- -;
Fig. 7. Influence of N, in step responses (C, -; U,-- -; c, = 1, f#l= 0.5; & = 1; CT = 4; c, = 0).
W. SUN and C. A. V. COSTA
542 1.0
KIMi P*=ZOO
0.8 s-o_-
Nf=5OO,Np=l,Nl=~
0.6 --+-
Nf=SOOO,
--C--
Nf=SOO,
Np=l, Nlml
c,u 0.4
0.2
---o--
Nf=SOOO.
Fig. 8. Combined influence of Nr, NP and N, (C, --;
-by means of a convolution calculate the desired real-time function to he used in the nonlinear equations; --solve the nonlinear one-region equations for example, the method of lines.
using,
Np’lOO.
U,---;
Ni=lO
NpolOO.Nl=lO
Ci, = I, r$ = 0.5; I&,= 1; C, = 4; C,, = 0).
This algorithm was tested using orthogonal collocation for the discretization of the 6xed-bed axial coordinate. These tests were performed over a large range of model parameters showing that: -mass
balance relative errors are within 1%;
1.0
0.0
0.6
Nf =5000
--“-
Pe=SO.Np=
--“-
P.130,
1. Nid
c,u 0.4
0.2
--&-
_
Np’lOO,
Pp=200,Np=100.
Ni=lO
Nl=10
0
Fig. 9. Combined influence of Fk, NP and Ni (C, -;
V,---;
Cb= 1, q5= 0.5; 4, = 1; CT = 4; C, = 0).
c.u
K=60 0.6
P6= 200 Nt = 6000 Np=l Nl-1
0
40
120
60
160
200
0
Fig. 10. Response to a pulse input. Influence of the pulse time length (C, -; &= 1; c,=4; C,=O). --comparison between the solutions of the linear model (constant velocity) obtained with this method and by FFT shows that the relative error in the ordinates (outlet solute axial concentration) is always less than 0.5%; -divergency or mass balance failure are observed when model parameters are such that a very steep axial profile develops. This is due to the use of global orthogonal collocation and thus a more suitable discretization method should be used in these conditions. It is worthwhile to mention that the use of this method saves a lot of computing time and numerical difficulties when compared to the direct use of the method of lines to solve this type of problem. 31 9.l2
This solution method is not limited to the class of problems treated. We think that it can be useful in the numerical solution of multiregion models where one or more regions are linear.
Acknowledgements-JNICT is gratefully acknowledged for financial support of this work. W. Sun also acknowledges Fundafio
c= e, = ca = ci = c, = cp = cpk = c, = C=
Cb = C, =
grant.
Solute concentration in bulk phase, gmol cm-) Initial concentration, gmol en-j Solute reference concentration, gmol cn+ Solute concentration in micropore, gmol cnd3 Inlet concentration (step input), gmol cnm3 Solute concentration in macropore, gmol crne3 Inlet concentration (Pulse), gmol cm-) Total concentration of solute and inert components in feed, gmol cmm3 Dimensionless solute concentration in bulk phase .Dimensionless initial concentration Dimensionless solute concentration in micropore Dimensionless inlet concentration (step input) Dimensionless solute concentration in macropore
CPvlr= Dimensionless inlet concentration (pulse) C, = Dimensionless total concentration of solute and inert components in bulk phase C = Dimensionless concentration variable of bulk phase cl
30
SO
60
120
160
El
Fig. 11. Response of a negative (C, -;
Oriente for a research
NOMENCLATURE
C, = %.__ *----..______ ------_-----____-________._________ Ci =
c.u
U,- - -; CPti = 1,4 = 0.5;
step of a presaturated bed (I.---; q5=OS; I$~= 1; CT = 4; C, = 0). K = 80, Pe=200. Nr=SODO, N,,= 1, Ni=l.
in Laplace domain Ci = Dimensionless concentration in Laplace domain
variable of micropore
C, = Dimensionless concentration variable of macropore in Laplace domain D, = Axial dispersion coeflicient, cm2s-l Di = Effective micropore diffusivity, cm2 S-I Dp = E&dive macropore diffusivity, cm2 s-l
W. SUN and C. A. V.
544 k, K L N Nr N,
= = = = = =
Np = Pe = r= R = r, = R, = g= t= T= U= u= us = X= Y= 2 =
External film mass transfer coefficient, cm s-’ Dimensionless adsorption equilibrium constant Bed length, cm Number of points used in FFT Dimensionless film mass transfer parameter Dimensionless mass transfer parameter im micropore Dimensionless mass transfer parameter in macropore Peclet number Radial coordinate for the microparticle, cm Radial coordinate for the pellet, cm Microparticle radius, cm Pellet radius, cm Dimensionless solute flux into the pellet Time, s Half operation period, s Dimensionless velocity Interstitial velocity, cm s-’ Reference interstitial velocity, cm s-’ Dimensionless microparticle radial coordinate Dimensionless pellet radial coordinate Dimensionless axial coordinate in the bed
Greek symbols a= /3 = L= ep = 4 = & = r = 0= w=
Constant of Mellin’s formula Parameter in equation (24) Bulk phase porosity, cm3 bed void/cm) bed Macropore porosity, cm’ macropores crnM3 pellet Dimensionless voidage parameter Dimensionless macropore voidage parameter Space time Dimensionless time variable Fourier integral variable
REFERENCES Bracewell R. N., The Fourier Transform and Its Applications, 3rd Edn. McGraw-Hill, New York (1987). Chen T. L. and J. T. Hsu, Prediction of breakthrough curves by the application of fast Fourier transform. AIChE J133, 1387-1390 (1987).
a,TA
Design of cyclic fixed bed Costa C. and A. Rodrigues, processes I. AIChE JI 31, 1655-1665 adsorption (1985). Hindmarsh A. C., LSODE and LSODI two new initial value ordinary differential equation solvers. ACM-SIGNUM News/err. 15, lo-11 (1980). Hsu J. T. and J. S. Dranoff. Numerical inversion of certain. Laplace transforms by ‘the direct application of fast Fourier transform (FFT) algorithm. Computers them. Engng II, 101-110 (1987). Kreyszig E., Advanced Engineering Mathematics, 6th Edn. Wiley, New York (1988). LeVan M. D. et al. Fixed bed adsorption of gases: effect of velocity variations on transition types. AIChE JI 34, 996-1005 (1988). Loureiro J. M. and A. Rodrigues, Two solution methods for hyperbolic systems of partial differential equations in chemical ,engineering. Chem. Engng Sci. 46, (1991). in press. Raghavan N. S. and D. M. Ruthven, Simulation of chromatographic response in columns packed with bidisperse structured particles. Chem. Engng Sci. 40,699-706 (1985). Raghavan N. S., M. M. Hassan and D. M. Ruthven, Numerical Simulation of a PSA system, Part I. AIChE Jl 31, 385-392 (1985). Rasmuson A., Time domain solution of a model for transport processes in bidisperse structured catalyst. Chem. /&pig Sci. 37, 787-788 (1982). Rasmuson A. and I. Neretnieks, Exact solution of a model for diffusion in particles and longitudinal dispersion in oacked beds. AIChE JI 26. 686-690 (1980). Rhee H. K., R. Aris and N. R. Amundson, First Order Partial Dtfirential Equations. Vol 1. Prentice-Hall, Englewood Cliffs, NJ (1986). Sereno C., A. Rodrigues and J. Villadsen, The moving finite element method with polynomial approximation of any degree. Comtmters them. Ennnn 15, 25-33 (1991). Sun L. M. and F. Meunier, An improved finite difference method for fixed bed multicomponent sorption. AIChE JI 37, 244-254 (1991). Villadsen J. and M. L. Michelson. Solution of Differential Eqaations Models by Polynomial Approximaiion. &enticeHall, Englewood Cliffs, NJ (1978).