- Email: [email protected]
Analysis of transient models for catalytic tubular reactors by orthogonal collocation
Chemical Engineering Science, 197 l, Vol. 26, pp. 1555-1569.
Pergamon Press.
Printed in Great Britain
Analysis of transient models for catalytic tubular reactors by orthogonal collocation K N U D W,~DE H A N S E N Instituttet for Kemiteknik, Danmarks tekniske HCjskole, Lyngby, Denmark (Received 2 March 1971)
In the simulation of the dynamics of the catalytic fixed-bedreactor it is of great importance to use a correct model for the catalyst pellet. The various catalyst models are tested by a stationary as well as by a dynamic analysis. The results obtained with alternative models are then compared. Finally calculations of transients in a tubular reactor are Performed using a catalyst model, which includes both film and intraparticle resistances against heat and mass transport. The computation method is a combination of orthogonal collocation and the method of characteristics. This method allows the calculations to be performed under consumption of a reasonable amount of computer time. Abstract-
INTRODUCTION SIMULATION of the dynamics of distributed parameter systems has been subject to much attention during the last years. H o w e v e r , the study of complicated systems such as tubular reactors with catalytic reactions has been of a relatively moderate extent owing to the complicated mathematical description. But recently large scale c o m p u t e r and more efficient numerical solution procedure have made it possible to handle models for this very important, but complicated group of systems. T h e present paper deals with a tubular reactor packed with porous catalyst pellets o f spherical form. T h e reactant and products are gases. T h e reaction takes place when the reactant diffuses from the main stream into the pores o f the catalyst pellets. T h e conversion into products occurs when the reactant contacts the porewalls and it follows a first order expression. A ~ products r = k. c = kl. e -EmT. c where the rate constant k is a function of temperature according to an Arrhenius expression. After the reaction has taken place the products will diffuse through the pores back into the main stream. T h e reaction itself produces heat, which
is transported from the pellets to the main stream - m a i n l y through conduction in the solid material. Thus both concentration and temperature profiles develop in the pellets. T h e increase in temperature has the effect that the reaction rate grows exponentially, and it is particularly this non-linear coupling between mass and heat transfer, which causes difficulties in obtaining the solution. Earlier work has concentrated on the description of the gasflow in the external medium. Using a plug flow model Liu, Aris and A m u n d s o n [1-3] thus make computations comparing the results with calculations where an axial dispersed plug flow model has been used. Vanderveen, Luss and Amundson[4] and Deans and Lapidus[5] employ special flow models with networks of mixing elements in one and two dimensions respectively. In all these cases simple models have been used for the heat and mass transport in the catalyst pellets, as the transport resistance is thought mainly to be situated in the film around each single pellet, whereas internal gradients in concentration and temperature are neglected. M c G u i r e and Lapidus [6] include both types o f transport resistance and use Deans and Lapidus's[5] two dimensional flow model for the external medium. T h e solution method consists o f a discretisation technique followed by an
1555
K. W. HANSEN iterative solution of the matrix equations. This method permits computation of transients in very small reactors (5.15 catalyst diameters) to be carried out in about 3hr on an IBM 7090. Recently Feick and Quon [7] have improved the solution technique and reduced the computation time to 5-90min, depending on the chosen flowmodel (IBM 360•67). The reactor dimensions are in this case 10.50 particle diameters. The above mentioned papers all deal with cases in which the disturbances are large. Therefore it is necessary to use nonlinear models. However, when only small deviations from the steady state occur during the transient it is possible to linearize locally around the steady state and compute the frequency response of the model, Crider and Foss [8, 9]. Stangeland [ 10] uses such linear models to compute transfer functions for Film
Cat. pellet r=l
MI
Id2
M5
r=O
controllers which, based on measurements of the temperature in the inlet and at one point inside the reactor, manipulate the temperature or concentration in a secondary feed stream in order to eliminate temperature or concentration disturbancies in the outlet. Experimental verification of reactor models only has been carried out in a few studies. Lyche [11], and Hoiberg, Lyche and Foss [12] have compared measurements and theoretical computations. The models used are linear and do not include internal gradients in the pellets. The experiments were performed with the reaction between oxygen and hydrogen catalyzed by platinum on silica gel granules. Especially Hoiberg et al. find good agreement. Sinai and Foss[13] work with the liquid phase reaction Na~S203 + H20,, ~ products in a tubular reactor with an inert packing. In this case, too, a good agreement is found between the frequency reponse of the linearized model and the experimental results. However the conclusion must be that there is still a demand for both experimental and theoretical investigations in which the attention is drawn especially to the non-linear phenomena. This paper deals with valuation and comparison of alternative models for the single catalyst particle and finally a new computing method for the reactor transients will be presented.
ANALYSIS OF MODELS FOR THE CATALYTIC PELLET
M4
~s M5
M6
------
Temperature Concentration
Fig. 1. Catalyst models.
The catalyst pellets dominate the dynamic behavior of the fixed bed reactor when the reaction medium is a gas. It is therefore important to include all relevant transport resistances in the model for this element. The complicated models are, however, difficult to treat numerically and one is interested in using a model w h i c h - u n d e r observation of imposed conditions of a c c u r a c y is as simple as possible. A number of possible models, which include different combinations of transport resistances, will be listed in the following. From the results of the calculations it will appear in which parameter domain it is possible 1556
Analysisof transientmodels to substitute the complicated models with the more simple ones. Figure 1 shows the physical situation for 6 different models. Model No. 1 (M1) is the isothermal one, whereas the others are non-isothermal and non-linear. In model No. 2 (M2) it is assumed that only the transport resistances in the film are important. In M3 the film resistance against the heat transport and the intraparticle resistance against the mass transport are included. M4 is the reflection of M3 with internal heat diffusion and film resistance against mass transport. M5 only includes intraparticle resistances and in M6 both film and intraparticle resistances exist. M6 is the most complicated
model and agreement within a parameter domain between this model and one of the less complicated models indicates that the simpler model includes the relevant transport resistances. In Table 1 the mathematical formulation of each model is listed as differential equations. The dependent and independent variables and the parameters are all dimensionless quantities. In M6 5 parameters determine the steady state, while the 2 remaining parameters only appear in the dynamic description. The other models contain a smaller number of parameters which will appear through combination of the parameters in M6. Table 2 contains the definitions of the par-
Table 1. Mathematical representation o f the models for the catalyst pellet Model number M1
Differential equations ~Xk
e2--= Ot
Boundary conditions Xk(1) = I
V2Xk--tl~Xk
dxk
e2~-~- = 3o"( l - - x~ ) - - a x k e y(~-l/uk)
M2 d y k __ L e ~ - - - 3v( 1 - - y e ) + a f t er(l-llUk)Xk
OXk e . ) - - ~ ~2X k -- OfXk eY(1--11Yk) " Ot
M3
Xk(1 ) = 1 L e - d- ~y k -_- 3v(l --Yk) + 3 a f t e ytl-l/ut? f x k r 2 d r ~o
e2--~d x k = 3o" ( 1 -- xe) -- 3axk f o 1 e y~l-l~upr2 dr
M4
yk(l) = 1 L e -~ - = V2yk -}- otj~xk e't(l-ltU~ Ot OXk = V 2Xk - - OLXk e ~(1-lly k) ~2 O--t-
x k ( 1) = 1
L e v ~ = V2yk + aflXk C ( 1 - 1 t ~ Ot
yg ( 1 ) = 1
M5
M6
OXk ~2 _~ = V2Xk - - aXk eTtl-ltvk} Ot
VXk(1 ) = o ' ( l - - x k ( 1 ) )
L e - ~ = V 2yk + aflXk e yll-llul~ Ot
V y ~ ( 1 ) = v ( I -- y~ (1) )
1557
K. W. HANSEN Table 2. The parameters Parameter
Definition
Variationinterval
R 2 . kfe -y
a
0"01-100
De
/3
(-AH)c~oDe Tuok~
y
E R~T~
10-40
tr
kc.R De
20-20,000
v
h.R ke
1-200
0"0005-0"5
(pCp)k. De
Le
or is a modified Sherwood n u m b e r and v a modified Nusselt number. "Modified" signifies that the quantities contain parameters defined in different mediums. Large transfer coefficients are equivalent to small changes in potentials across the film. L e is a Lewis number for the catalyst, defined as the ratio between mass and thermal diffusivity. T h e values o f the parameters used in the following examples are summed up in Table 3.
0"5-100
ameters, and boundaries are stated within which one will normally expect the parameter values to be in physically realistic situations. R2/De is a characteristic diffusion time for the pellet and ks. e -~ the reaction rate at the gas temperature in the main stream. Consequently a is the ratio between the diffusion and the reaction time at the surface conditions. V~da= q~ is normally denominated the Thiele modulus. T h e parameter/3 is proportional to the reaction enthalpy, and for steady state condition/3 is the maximum dimensionless temperature rise inside a pellet when only internal resistances against heat and mass transport are considered. Y is proportional to the activation energy, which expresses the influence of the temperature on the reaction rate.
ANALYSIS OF THE STEADY STATE T h e treatment of stationary catalyst problems on a theoretical basis has dominated the literature during the last years. One has especially worked with the models M5 and M6 and tried to find the number of solutions as well as the stability of these and domains of existence and uniqueness as a function o f the parameters. In the early part o f the 1960's Tinkler and Metzner [14] and Weisz and Hicks [15] presented numerical computations which showed that 3 solutions could appear in M5 for sufficiently high values of the product y . / 3 within an interval of a. Hatfield and A d s [ 16] have used values for a,/3, 3' in the same range and calculated the steady states for M6, and they find that up to 5 solutions may appear when o- > v. Still with M6 Villadsen[17] and Cresswell [ 18] perform computations using a, /3, y values, which only give 1 solution for M5. F o r much lower y a n d / 3 values 3 steady states may now appear. Regarding the validity of the simpler models most of the published work only contains qualitative considerations dealing with the dominance of different transport resistances under various
Table 3. Values o f the parameters u s e d in preparation of Figs. 2 - 1 3 Figs.
c~
fl
Y
tr
v
2, 3 4, 5 6, 7 8, 9 11 12, 13
100 1 1 1 1 0-5
0"01 0" 1 0" 1 0"01 0"01 0"01
20 40 40 20 20 20
250 250 250 250 400 400
12-5 12-5 12"5 12"5 5 5
Le
1-100
~.z
U
W
A Y
AX
5 2"5
100 5
75 150
120 240
0.5 0.5
In the last Figs. (11-13) the stated values of a, fl and y refer to reactor inlet conditions and elsewhere to pellet surface conditions. 1558
Analysis of transient models
conditions. In the following two sections it is tried to obtain more exact information by computing the stationary profiles for all the formulated models using realistic parameter values inside the variation intervals shown in Table 2. COMPUTING
METHODS
The mathematical formulation of the stationary models appears from Table 1 by setting the accumulation terms on the left side of the equations equal to zero. M 1 is easily solved. The other models are nonlinear and more than one solution to the equations may occur. For M2 and M3 it is possible to construct algebraic algorithms, guaranteeing all solutions to be found, and then it is easy to determine the concentration and temperature profiles. The methods are described in Ref. [ 19]. For the solution of M4, M5, M6 orthogonal collocation is used in connection with a NewtonRaphson technique to solve the resulting algebraic equations. The method has been introduced by Villadsen [17, 20] and it has proved to be superior to a number of alternative methods. Orthogonal collocation is a solution method for differential equations and is one of several Methods of Weighted Residuals. The differential equation is satisfied in N predetermined points of the variation interval of the independent variable. These points are the zeros of an orthogonal polynomium of the order N. By choosing a suitable weighting function for these polynomials it is possible to fix the zeros at points where it is important to obtain an accurate solution. The catalyst problem considered is symmetrical, and an even set of Jacobi-polynomials can be used. The gradient Vy and the Laplace operator V~y are found as weighted averages of the value of the dependent variable at all collocation points and at the pellet surface:
The matrices A and B are only dependent on the chosen polynomials and can be found once for all. Details about the procedure and about the method itself can be found in Refs. [ 17] and [20]. In this case a weighting function (1 - F ) . r z is used. For the three models M4, M5 and M6 the number of solutions is determined by an approximate procedure consisting of orthogonal collocation with two internal points. Details are given in Ref. [ 19]. The two-point method produces fine results for M5 but should be used with caution for M6. It may be mentioned incidentally that Stewart and Villadsen [21] report to have obtained satisfactory results for model M5 with a method of lower o r d e r - a one-point orthogonal collocation. When the number of solutions and the corresponding approximative profiles are found, a high order method is used for the exact determination of the profiles. RESULTS
In Figs. 2-5 some examples are shown where the concentration and temperature profiles for
0"75 x x
0-5O a=lO0 o-=250 ~ + l
~-~1 u= 12.5
~=20
MI M2
(>25 l e l
MS
~×~
M6
N+I
Vy(l) = ~, A ( l , J ) y ( J ) J=l
I 0"~
N+I
0"50
0"75
r
V~y(1) = ~, B ( l , J ) y ( J ) .
Fig. 2. Concentrationprofile.Steady-state.
J=l
1559 C E S Vol. 26, No. 1 0 - - E
0"25
I.O
0-00
K. W. HANSEN concentration inside the pellet, whereas the temperature only shows a modest increase. It is observed that the concentration profiles for M 1, M3, M5 and M6 are close to each other, whereas M2 deviates strongly. F r o m this it is deduced: (1) the thermal effects have only little influence on the concentration profile, (2) the major mass transport resistance is distributed as diffusional resistance inside the pellet. As the mean temperature in M3 and M5 deviate from M6 to approximately the same degree it appears by observing the temperature profiles that the film resistance against the heat transport is of the same order o f magnitude as the intraparticle resistance. T h e s e results are true in many other cases, which means that M2 is a poor model in most practical cases. M 4 has been omitted from the figures because it will give the same result as M2. It is obvious that the concentration drop
I "020
I.OI5
× >, 1.010
IOO5
a-lO0 o-=250
M2 1.000' 0.00
~=0.01 v=12"5
~a,~
"),"• 20
\ M3
J 0"25
~l~
M5
~x~
I 050
M6 I 075
1.0
Fig. 3. Temperature profile.Steady-state.
,.or o.,,I-
+
MI 0"7
M2 M3 M5
0-6
05 0"0
~
~
M6
I
I
I
I
i
0"2
0"4
0-6
0"8
I'0
Fig. 4. Concentration profile. Steady-state. the various models are compared. F o r certain combinations of the parameters the solutions to all models are found and plotted in the same figure. Figures 2 and 3 are concentration and temperature profiles, respectively, computed for parameter values, which only generate one solution in all models, D u e to the large value of the Thiele modulus there is a pronounced decrease in the
across the film is much smaller than the decrease inside the pellet. T h e profiles in Figs. 4 and 5 are calculated for the same values of cr and v but for larger values o f the thermal parameters/3 and y and a smaller value of a. It would o c c u r for larger values o f the heat o f reaction and activation energy and for smaller catalyst dimensions. It appears that also in this case M2 deviates completely from the
1560
Analysis o f transient models
1.05 ~
x
~
a-I
x 1-04
:~ 1-03
~
x
M2
y-40
~A~
o'-250
~e~
M3 M5
v= 1 2 " 5
~x~
M6
°~o
1.02
x
1.01
I'00 0.0
•
I 0"2
I •
I
0"4
06
I 0"8
~'~,Q I~)
t
Fig. 5. Temperature profile. Steady-state.
The conclusion of these remarks and a series other models. The increase of temperature is still modest, but owing to the strong influence of of computations presented in Ref. [19] is that temperature on the reaction rate it gives rise to for mass transfer the internal diffusion resistance a considerable decrease in concentration. Here will dominate the external film resistance, whereas the internal and external resistances are often M5 gives a better result than M3. While only one solution was obtained in all of equal magnitude for the heat transport. Thus models for the parameter values used in Figs. 2 in certain cases it will be possible to substitute and 3 the parameter values in Figs. 4 and 5 will M6 with M5 and in other cases with M3. generate three solutions in M2 and M3, but only ANALYSIS OF THE DYNAMIC STATE one in M5 and M6. It is the "lowest" steady state Very few computational studies of transients in M2 and M3 which is shown in the figures. The in a single catalyst pellet have previously been "higher" steady states give very high conversion reported. and temperature. It is now obvious that one has Hlavacek and Marek[22] present computato be very careful when using simpler models tions in which M5 is used. They investigate the for the catalyst pellets in a fixed-bed reactor, influence of the value of the Lewis number when because the conditions and parameters (a,/3, y) different initial conditions are used in both the vary along the reactor. To be valid within a multiand single-solution domains. parameter range the simpler models must provide Ferguson and Finlayson[23] use almost the the same number of solutions and approximately the same spatial mean value of temperature and same computational technique as in this paper to calculate transients for M5 and M6. They apply concentration in each steady state as M6. As mentioned above, one might very well find orthogonal collocation combined with an exmore than one solution to M5 and M6. The film plicit integration method in the time domain resistances--especially against heat transport-- using polynomials for each model. Orthogonal collocation is also used in this are much more important in the "higher" steady states. For these solutions there will exist an study. The polynomials are even and belong to agreement between M6 and M3 rather than the Jacobi family as in the stationary case. The resulting ordinary differential equations are between M6 and M5. 1561
K. W. HANSEN integrated by a fourth order R u n g e - K u t t a method with a variable step size. T h e adjustment o f the step size takes place through calculation of the truncation error according to principles, which are recommended by Zonneveld [24]. T h e computations which are presented below only show the results for models M1, 3, 5, 6, since as mentioned in the previous sections M2 and M4 are not applicable in the relevant parameter domain.
and concentration for the models M 1, 3, 5, 6 are shown. As initial conditions profiles deviating from the steady state have been inserted. 8 collocation points have been used and this proved to give an accuracy o f at least 3 - 4 figures. Figures 6 and 7 show the transient response o f the spatial mean concentration and temperature after x = 0.8 and y = 1.15 have been inserted as initial values at all the collocation points. T h e stationary parameter values correspond to those of Figs. 4 and 5, and only one stationary state occurred in M5 and M6. T h e diffusion resistance inside the pellet is dominating
RESULTS In Figs. 6 - 9 some transients o f temperature I'O--
MI
0"8
A~
~
01=I B-0.1 o"=250 v=12.5 Le= I
g 0-6 x
o
~
"
× 0.2
0.¢
-00
I
0"05
I
I
O"I0
01,5
I 0,20
t
Fig. 6. Mean concentration. Transient response. Initial values: xk = 0.8, y, = 1.15. 1.20 -
;:0.' )'=40 ~" I't0
1.00--
000
V=12"5 Le=l C2=0-5
I
0"05
I
0-10 t
-A"
I
015
¢
020
Fig. 7. Mean temperature. Transient response. Initial values: xk = 0.8, Yk= 1.15. 1562
Analysisof transientmodels
oa
× Le=l
~
x
~
×_
x
/
B-o~, /
/ @4
Le-5
i
e
,'~,*'zs° )t=20
~
1,'=12.5
+
x
I
o.2
(>0
t
~
+
Le=lO0 I 0.1
I 0.2
I 0-3
I 04
•
(>5
t
Fig. 8. Effect of the Lewis number. M e a n concentration. T r a n s i e n t response of M6. Initial values: xk - Stationary profile, yk = 1.3.
1.3
La=lO0 Le=25
~
+
i
°-,
.8=001
~"
e
~
• I" I --
I .0 043
X ~ L e =I
I 0.1
I 0.2
I 0"3
7-20 o'=250
~ e
-0.4
vii2'5
E 2 =0.5
0.5
t
Fig. 9. Effect of the Lewis number. M e a n temperature. T r a n s i e n t response of M6. lnitial values: xk ~ Stationary profile,yk = l '3.
for both heat and mass transport. In the stationary state M5 is thus preferable to M3. The temperature perturbation causes the mean concentration in Fig. 6 to decrease heavily at the beginning of the transient for all the models shown, except of course for the isothermal model M 1, which goes straight back to the stationary state. M3 has a smaller minimum value than M6, but agrees better with M6 than M5 at the beginning of the transient. Later on, M3 goes quickly back to its inaccurate stationary state so that M5 is the best model at the latter part of
the transient. The temperature plot shows a maximum due to the strong heat generation. During this period it appears that also with regard to temperature M3 agrees better with M6 than M5 does. Consequently it must be assumed that the dominating resistance against the heat transport is to be found in the external film at the beginning of the transient, whereas the internal heat diffusion resistance later on gradually becomes dominant just as in the stationary state. The conclusion of these and other results
1563
K. W. H A N S E N
presented in Ref. [19] must be that obviously none of the simple models are sufficiently correct to substitute M6 throughout the parameter space investigated. For the computations of transients in a reactor it will therefore be necessary to use M6 as catalyst model. Actual reactor computations will show whether considering a large number of coupled catalyst pellets might after all make the use of simpler models feasible. Finally Figs. 8 and 9 show the influence of variations in the Lewis number on the transient response of M6 after a disturbance in the stationary temperature profile. For large Le-values the temperature dynamics is very slow causing a heavy decrease in the concentration. In that case it is possible to use a simpler method. Presumably a quasistationary solution, in which the accumulation term in the mass balance has been eliminated will be acceptable apart from the very first period.
Catalyst U OxA.
(III)
W Oyk = V 2 y k q.. otflXk e y(1-1/uk) Or
(IV)
0~" = V2Xk-- axk e ytl-1/uk)
where U=
E2. R 2 • Vs De. L
W = (pCp)k. R z • vs ke. L
The boundary conditions are Oxk(1 Or "-) = o'(xo--xk(1))
(V)
Oyk -~-r (1) = u(yg--yk(l))
(VI)
together with the symmetry conditions ANALYSIS OF THE REACTOR DYNAMICS
In this investigation the interest is mainly concentrated on the formulation of models for the catalyst. Therefore a number of simplifying assumptions will be introduced for the other elements in the reactor: (a) Adiabatic operation of the reactor. (b) No heat capacity of the reactor wall. (c) The density of the gas is constant and independent of temperature. (d) Other physical quantities are constant. (e) No axial and radial mixing in the external medium. With M6 as catalyst model the describing equations will be the following: External gas medium OYuo~.t-VO0-~ = A Y ( y k ( 1 ) --Y~')
(I)
Ox~ t _ v O x g = A X ( x k ( 1 ) - - x o ) Or Oz
(II)
where AY=
a.h.L ~
• (oC,,)o
AX=a.ke.L . v,
E~ . v ,
oxk -~r (0) = 0
(0) = O.
(VII, VIII)
The dimensionless parameters U and W are the ratio of the diffusion time for the mass transport to the reactor residence time and the diffusion time for the heat transport to the reactor residence time, respectively. High values of U and W cause the catalyst dynamics to be much slower than the transportation lag in the external medium. The ratio W~ U is equal to Le/e., and is thus of importance for the dynamics of the single pellets. Three independent variables occur in the m o d e l - t w o in space and one in time. The catalyst concentration xk and temperature Yk are functions of space variables r and z whereas the gas concentration xo and gas temperature yg only depend on the axial position z in the reactor. With the assumptions previously introduced the coefficients A X , A Y, U, W, o- and u are constants. This is a rather crude assumption, which, however, may be justified here since the basic transport processes are the central investigation object.
1564
Analysis of transient models SOLUTION METHOD
The numerical technique used in this paper has been specifically developed to try to meet the requirement of a low consumption of computer time. The equations for the catalyst (III)-(VIII) are transformed to a set of 2N ordinary differential equations, as in the dynamic analysis of the single pellet, by using orthogonal collocation with N internal points. The collocation equations in point No. (I) are: dXk (1)
dr
1 N+I
=-U "~ B(I,J).Xk(J) J=l
--Uxk(I)
(IX)
e y(1-11u~/))
d ~ = A Y ( y k ( N + 1) - y g ) .
The characteristics of the equations are plotted in the (r, z) plane in Fig. 10. Each Cline here represents the time variation of the concentration and temperature profiles at a certain point in the reactor. Similarly the variation of concentration and temperature in a flow element in the external medium is recorded by following a g-line from z = 0 until z = 1. The surface values Xk(N+ 1) and y k ( N + 1) can be eliminated by the boundary conditions (XI) and (XII). (IX), (X) and (XIII), (XIV) can then be rewritten. (IX) and (XIII) thus become:
d x k ( 1 ) - - ~ BB(l,J) xk(J)+CC(l) xo dr
dYk (1)
d'r
"
J=l
I~ ~ B(I,J) . Yk(J) + wXk(1)
O/
- - "-~Xk ( 1 ) e ytl-vuktm
eY(1-11Yk(l)).
N+I
(J) =
l/)
(XI)
,]=1 N+I
~, A (N+ 1,J)y~(J) = ~'(yg--yk(N+ 1)). (XlI) ./=1
The solution of these equations and the coupled equations for the external medium is performed by the method of characteristics[26]. The characteristics corresponding to the catalyst equations are straight lines parallel to the r-axis in the (r, z) plane, whereas the characteristics corresponding to the equations for the external medium are curves, which for constant flow are straight lines with a slope equal to V. With the relation dz/dr = V the equations for the external medium can be rewritten as (XIII) and (XIV).
~ HH(J). xk(J) + H H ( N + 1)xo. dr = J=l (XVI)
In an adequate number of points on the z-axis the Eqs. (XV) and its temperature analogue are set up. These equations must be integrated along the Ccharacteristics whereas (XVI) and its temperature analogue must be integrated along the 8-characteristics. The calculations are started along the 8characteristic, which is nearest to (r, z) = (0, 1). The integration from one intersection of this line ~-Iines ~0
g
- lines
o
0-0 Title,
dxu = AX(xk(N + 1) --xg)
(XV)
dx v
(X)
The boundary conditions for r = 1 are transformed to
dr
"
1 N+I J=l
Z A ( N + 1,J)x
(XIV)
(XlII) 1565
l"
Fig. 10. Characteristic curves for the reactor.
K. W. H A N S E N
with a ~-characteristic to another intersection is performed by a second order predictor-corrector scheme except for the first step, where it is necessary to use a first order predictor. After having finished the calculations along this line the integration is continued correspondingly along the nearest 8-line. The developed computer program is able to calculate the transients after step changes and sinus disturbances in the inlet concentration and temperature as well as step changes in the flow. A considerable part of the program is devoted to administration of the enormous data material produced in just one single transient. During the computation about ¼ million points are thus passed in the (r, z) plane. The computation of each transient is made in several stages, between which step length in time and space may be adjusted if required.
proximately reached the stationary profile of the isothermal model. Then it changes slowly following the increase in temperature which gradually causes the concentration profile to change curvature. The dotted profile in Fig. 11 is the stationary profile when using M5 as catalyst model instead of M6. The profiles are quite close to each other and consequently the film resistance against the heat transport is of only minor importance for the conversion. If the catalyst profiles are studied during the transient period it will appear that this result is not due to the fact that the external film resistance against the heat transport is insignificant compared with the internal resistance, but the increase of temperature in one pellet is rather small compared with the increase of temperature along the reactor. For this computation a step size of 0.02 in the axial direction was used together with time increments of 0.05 at the beginning and later on 0.1.
RESULTS
In the following some results are shown which illustrate the most important dynamic characteristics of the reactor. Furthermore it will be possible to draw conclusions about the applicability of the computation method. The parameter values for each of the examples are summed up in Table 3. By computations, described in Ref. [19], it has been ascertained in all cases that only one stationary state exists for the catalyst pellets in the inlet. Thus only one possible stationary state exists for all pellets in the reactor, as the local 3/and/3 values will both be smaller than 3/and/3 calculated on inlet conditions. The stationary state for the reactor was found by introducing initial values equal to 1 for concentration and temperature in both the catalyst and the gas medium and observing the transients until no change took place. Figure 11 shows an example, in which the concentration and temperature profiles in the reactor are plotted. The ratio between W and U is 20. The concentration profile in each single pellet will therefore at the beginning of the transient change much faster than the temperature profile. After 5-10 residence times the concentration profile has ap-
1"08
~®_
/:////
I
r~4
1"06 >'~ 1'04 1"02
I
1.00 .
r.O.4
r=O-8
------ M5 0"92
Q=I r-2o 0"=400 0.6
0 0"820-00
T.,J
N~
AX-,20
X\
"
8 0"25
4 0"50
0"75
I'00
Z
Fig. 11. T r a n s i e n t response. Reactor profiles. Initial conditions: x~ = x, = Yk= Y u = 1.
1566
Analysis of transient models
The catalyst profiles were computed with two collocation points. On an IBM 7094 the computations could be performed in about 80 min. Figure 12 shows the result of a perturbation of the inlet temperature from 1 to 1.1 during 6 residence times. The stationary profiles have been obtained by an integration procedure similar to the one mentioned above. The parameter values are somewhat different from those on Fig. 11. It is especially noted that the heat capacity of the catalyst is a factor 10 smaller. This causes a much faster dynamic behavior. (Le/e2 = 2). After the 6 residence times the stationary state a is changed to the b-profiles. A temperature maximum has developed in the inlet while the concentration drops throughout the reactor. Hereafter, it will be seen that the temperature wave is moving slowly towards the reactor outlet with a velocity approximately equal to El (pCp)o/ [ (1-- ~1) (pCp)k]. Gradually the wave changes its form becoming wider and getting a still higher maximum value. The temperature perturbation is immediately followed by a drop in the outlet concentration, whereas the outlet temperature only changes slowly. At the beginning of the transient the rate of conversion will thus be lower in the last half of the reactor,
and consequently the heat generation will be smaller in this part too. During this period the outlet temperature will therefore drop and reach the minimum value after 22 residence times. The deviations from the stationary profiles are pronounced for the case shown in Fig. 12. In order to investigate whether the temperature difference between surface and centre of the single catalyst pellet still remains small relative to the overall temperature changes in the reactor, the time variation of the catalyst profiles is shown at the reactor outlet in Fig. 13. The pulse disturbance in the inlet temperature causes the concentration and temperature values to decrease until r = 27, whereafter the temperature starts increasing rapidly. After 48 residence times the temperature approaches the maximum value, and at the same time the concentration profile is very steep. During the part of the transient shown in Fig. 13 the temperature of the external medium is higher than the temperature of the pellets, and
r=48
14 45 1'3 ,~
i.z
41 36
115f
g/~h Profile
>,~ 1,4
Time
o b c d e f 9 h -6 0 6 22 32414851
1,3
d
e
Iq i l0
I
I
~
/
~
075 Stationary state
/
o
-
x 0"8
a=0.5 b ~=0.01 AY=ISO 7"=20 AX=240 0"61-- o",400 U=2"5 U=5 W=5
~
j
j
~ 0,25
@4 0"00
I
"r=O
f ~
F2
I'0
\27
Stationary state and
I 0"25
0"50
0"75
I'00 ~
Fig. 12. Transient responsel Reactor profiles. Initial conditions: Stationary state. Inlet temperature equal to 1.1 for 6 residence times.
1567
:i i .J. 025
0.50 r
0-75
Fig. 13. Catalyst profiles at z = I for the case in Fig. 12.
[
K. W. H A N S E N
thus the temperature decreases towards the centre of the pellet, but the profile is, however, still almost fiat. Integration step sizes between 0.01 and 0.005 were used in the axial direction and from 0.01 to 0.04 in the time direction. The computer time until r = 51 amounted to 50 min.
NOTATION
A A (I, J)
AX AY a B (1, J)
CONCLUSION
The reactor dynamics for gaseous reaction media are dominated by the heat transport processes. In the single catalyst pellet the concentration profile will often follow the temperature profile in a quasistationary state. A temperature disturbance moves like a slow wave through the reactor due to the small heat capacity of the gas medium compared with that of the catalyst. The temperature variations in the individual pellets are rather small compared with the temperature changes throughout the reactor, at least when only one stationary state of the catalyst is possible for all pellets in the reactor. These observations make it possible to use simpler models, and the validity of these models will be reported elsewhere (Ref. [25]). The numerical technique presented has proved efficient for the computation of transients in reactors. The method seems to allow a considerable reduction of the computation time as compared to previous works. Furthermore the computations have given important information about the best way to build up the routines just as it has become obvious that certain modifications of the method will reduce the computation time even more. This is due to the fact that stability problems arise when very steep concentration profiles occur in the catalyst pellets. This may be remedied either by changing the weight function or by a segmentation of the variation interval of the radial pellet cordinate using a spline-collocation approach. Acknowledgements--The author wants to express his gratitude to Messrs. J. Villadsen, M. Kiimmel and H. Livbierg for helpful discussions during the work and for their valuable comments on the presentation of the results in this paper.
BB(I,J) Cp
CC(1,J) c De E --AH
HH(1)
reactant weight coefficient in the expansion of first derivative
akcL/(elvs) ahL/(e,(pC,)ov,) catalyst surface/reactor volume = 3(1--Ea)/R weight coefficient in the expansion of the Laplace operator
(I/U)(B(I, J ) - - A ( N + I , J ) . B ( I , N + 1)](o'+A(N+ 1 , N + 1))) specific heat capacity
B(I,N+I).o-/(U(o-+A(N+I, N + 1))) concentration of A effective diffusion coefficient activation energy heat of reaction
- - A X . A ( N + I,I)/(o-+ A ( N + I,
N+I)) h heat-transfer coefficient kc mass-transfer coefficient ke effective heat conductivity kt frequency factor L length of reactor Le Lewis number (pCp)kDdke I position in reactor N number of internal collocation points R radius of catalyst pellet Ro gas constant r dimensionless radius r'/R or reaction rate r' radial position in catalyst T temperature t dimensionless time t'De/R 2 t' time
U V
~zR2vs/(De. L) v/v,
v average velocity in reactor
W x y z
(pCp)kR2vs/(ke . L) dimensionless concentration c/coo dimensionless temperature T/Too dimensionless position in reactor 1/L
Greek symbols
1568
Ot
R 2 .
ks.
e--E/~Roroo)]De
Analysis of transient models
fl (--AH)cgoDd(Tooke) Y E/(RoTuo) 8 def. in Fig. 10 def. in Fig. 10 E1 void fraction of reactor E2 void fraction of catalyst p density Thiele modulus V ~ v hR/ke [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21 ] [22] [23] [24] [25] [26]
or
kcR/De
1"
dimensionless
time t'vs/L
Subscripts gas catalyst 0 s u r f a c e c o n d i t i o n s ( F i g s . 2 - 9 ) o r inlet conditions (Figs. 11-13) s stationary
g k
REFERENCES L1U S. L. and A M U N D S O N N. R., lndl Engng Chem. Fundls 1962 1 200. L I U S. L., A M U N D S O N N. R. and ARIS R., ibid. 1963 2 12. LIU S. L. and A M U N D S O N N. R.,ibid. 1963 2 183. V A N D E R V E E N J. W., LUSS D. and A M U N D S O N N. R.,A.I.Ch.E. J11968 14 637. D E A N S H. A. and L A P I D U S L., ibid. 1960 6 656. M c G U I R E M. L. and L A P I D U S L., ibid. 1965 1185. F E I C K J. and Q U O N D., Can. J. chem. Engng 1970 48 205. C R I D E R J. E. and FOSS A. S.,A.I.Ch.E.J11966 12 514. C R I D E R J. E. and FOSS A. S., ibid. 1968 14 77. S T A N G E L A N D B. E., Ph.D. Thesis, University of California, Berkeley 1967. L Y C H E B. J., Ph.D. Thesis, University of California, Berkeley 1968. H O I B E R G J. A., L Y C H E B. J. and FOSS A. S., Experimental evaluation of dynamic models for a fixed-bed catalytic reactor, Paper presented at The 1st Int. Symp. on Chem. Reaction Engng, Washington, June 1970. SINAI J. and FOSS A. S.,A.I.Ch.E. J11970 16 658. T I N K L E R J. D. and M E T Z N E R A. B., Ind. Engng Chem. 1961 53 663. W E I S Z P. B. and HICKS J. S., Chem. Engng Sci. 1962 17 265. H A T F I E L D B. and ARIS R., ibid. 1969 24 1213. V I L L A D S E N J. V., Selected Approximation Methods for Chemical Engng Problems. Danmarks tekniske HCjskole 1970. CRESSWELL D. L., Chem. Engng Sci. 1970 25 267. W , ~ D E H A N S E N K., Ph.D. Thesis (in Danish), Danmarks tekniske HCjskole 1970. V I L L A D S E N J. V. and S T E W A R T W. E., Chem. Engng Sci. 1967 22 1483. STEWART W. E. and V I L L A D S E N J. V., A .I.Ch.E. J11969 15 28. H L A V A C E K V. and M A R E K M., paper presented at The 4th European Syrup. on Chemical Reaction Engng, paper 2.5, September 1968. F E R G U S O N N. B. and F I N L A Y S O N B. A., Chem. EngngJ. 1970 1 327. Z O N N E V E L D J. A.,Automatic Numerical Integration. Mathematical Centrum, Amsterdam, 1964. W/EDE H A N S E N , K. To be published. AMES W. F., Nonlinear Partial Differential Equations in Engineering. Academic Press, New York 1965. R 6 s u m 6 - D a n s la simulation de la dynamique du r6acteur catalytique b, lit fix6, il est important d'utiliser un mod61e correct pour le comprim6 de catalyseur. Les diff6rents mod61es de catalyseur sont test6s d'apr~s une analyse statique ainsi que dynamique. Les r6sultats obtenus avec les deux modEles sont ensuite compar6s. Enfin, des calculs de transition sont faits dans un r6acteur tubulaire en utilisant un modEle de catalyseur qui tient compte des r6sistances de film et de particules contre le transport de masse et de chaleur. La m6thode de calcul est une combinaison de la collocation orthogonale et de la m6thode des caract6ristiques. Cette m6thode permet d'effectuer les calculs, 5. l'alde d'un ordinateur, dans un temps raisonnable. Zusammenfassung-Bei der Nachahmung der Dynamik des katalytischen Festbettreaktors ist es besonders wichtig ein richtiges Modell fiir die Katalysatorkugel zu verwenden. Die verschiedenen Katalysatormodelle werden durch eine station~e sowie durch eine dynamische Analyse gepriift. Die mit Modellen verschiedener Art erhaltenen Ergebnisse werden dann verglichen. Schliesslich werden Berechnungen von f0bergangsvorgSngen in einem Rohrreaktor unter Verwendung eines Katalysatormodells ausgefiihrt, das sowohl F i l m - a l s anch Zwischenteilchenwiderst~nde gegen W ~ m e - u n d Stofftransport einschliesst. Die Berechnungsmethode ist eine Kombination orthogonaler Zusammenstellung und der Methode der charakteristischen Merkmale. Mittels dieser Methode kiSnnen die Berechnungen unter Aufwendung eines verniinftigen Ausmasses an Computerzeit durchgefiihrt werden. 1569
Pergamon Press.
Printed in Great Britain
Analysis of transient models for catalytic tubular reactors by orthogonal collocation K N U D W,~DE H A N S E N Instituttet for Kemiteknik, Danmarks tekniske HCjskole, Lyngby, Denmark (Received 2 March 1971)
In the simulation of the dynamics of the catalytic fixed-bedreactor it is of great importance to use a correct model for the catalyst pellet. The various catalyst models are tested by a stationary as well as by a dynamic analysis. The results obtained with alternative models are then compared. Finally calculations of transients in a tubular reactor are Performed using a catalyst model, which includes both film and intraparticle resistances against heat and mass transport. The computation method is a combination of orthogonal collocation and the method of characteristics. This method allows the calculations to be performed under consumption of a reasonable amount of computer time. Abstract-
INTRODUCTION SIMULATION of the dynamics of distributed parameter systems has been subject to much attention during the last years. H o w e v e r , the study of complicated systems such as tubular reactors with catalytic reactions has been of a relatively moderate extent owing to the complicated mathematical description. But recently large scale c o m p u t e r and more efficient numerical solution procedure have made it possible to handle models for this very important, but complicated group of systems. T h e present paper deals with a tubular reactor packed with porous catalyst pellets o f spherical form. T h e reactant and products are gases. T h e reaction takes place when the reactant diffuses from the main stream into the pores o f the catalyst pellets. T h e conversion into products occurs when the reactant contacts the porewalls and it follows a first order expression. A ~ products r = k. c = kl. e -EmT. c where the rate constant k is a function of temperature according to an Arrhenius expression. After the reaction has taken place the products will diffuse through the pores back into the main stream. T h e reaction itself produces heat, which
is transported from the pellets to the main stream - m a i n l y through conduction in the solid material. Thus both concentration and temperature profiles develop in the pellets. T h e increase in temperature has the effect that the reaction rate grows exponentially, and it is particularly this non-linear coupling between mass and heat transfer, which causes difficulties in obtaining the solution. Earlier work has concentrated on the description of the gasflow in the external medium. Using a plug flow model Liu, Aris and A m u n d s o n [1-3] thus make computations comparing the results with calculations where an axial dispersed plug flow model has been used. Vanderveen, Luss and Amundson[4] and Deans and Lapidus[5] employ special flow models with networks of mixing elements in one and two dimensions respectively. In all these cases simple models have been used for the heat and mass transport in the catalyst pellets, as the transport resistance is thought mainly to be situated in the film around each single pellet, whereas internal gradients in concentration and temperature are neglected. M c G u i r e and Lapidus [6] include both types o f transport resistance and use Deans and Lapidus's[5] two dimensional flow model for the external medium. T h e solution method consists o f a discretisation technique followed by an
1555
K. W. HANSEN iterative solution of the matrix equations. This method permits computation of transients in very small reactors (5.15 catalyst diameters) to be carried out in about 3hr on an IBM 7090. Recently Feick and Quon [7] have improved the solution technique and reduced the computation time to 5-90min, depending on the chosen flowmodel (IBM 360•67). The reactor dimensions are in this case 10.50 particle diameters. The above mentioned papers all deal with cases in which the disturbances are large. Therefore it is necessary to use nonlinear models. However, when only small deviations from the steady state occur during the transient it is possible to linearize locally around the steady state and compute the frequency response of the model, Crider and Foss [8, 9]. Stangeland [ 10] uses such linear models to compute transfer functions for Film
Cat. pellet r=l
MI
Id2
M5
r=O
controllers which, based on measurements of the temperature in the inlet and at one point inside the reactor, manipulate the temperature or concentration in a secondary feed stream in order to eliminate temperature or concentration disturbancies in the outlet. Experimental verification of reactor models only has been carried out in a few studies. Lyche [11], and Hoiberg, Lyche and Foss [12] have compared measurements and theoretical computations. The models used are linear and do not include internal gradients in the pellets. The experiments were performed with the reaction between oxygen and hydrogen catalyzed by platinum on silica gel granules. Especially Hoiberg et al. find good agreement. Sinai and Foss[13] work with the liquid phase reaction Na~S203 + H20,, ~ products in a tubular reactor with an inert packing. In this case, too, a good agreement is found between the frequency reponse of the linearized model and the experimental results. However the conclusion must be that there is still a demand for both experimental and theoretical investigations in which the attention is drawn especially to the non-linear phenomena. This paper deals with valuation and comparison of alternative models for the single catalyst particle and finally a new computing method for the reactor transients will be presented.
ANALYSIS OF MODELS FOR THE CATALYTIC PELLET
M4
~s M5
M6
------
Temperature Concentration
Fig. 1. Catalyst models.
The catalyst pellets dominate the dynamic behavior of the fixed bed reactor when the reaction medium is a gas. It is therefore important to include all relevant transport resistances in the model for this element. The complicated models are, however, difficult to treat numerically and one is interested in using a model w h i c h - u n d e r observation of imposed conditions of a c c u r a c y is as simple as possible. A number of possible models, which include different combinations of transport resistances, will be listed in the following. From the results of the calculations it will appear in which parameter domain it is possible 1556
Analysisof transientmodels to substitute the complicated models with the more simple ones. Figure 1 shows the physical situation for 6 different models. Model No. 1 (M1) is the isothermal one, whereas the others are non-isothermal and non-linear. In model No. 2 (M2) it is assumed that only the transport resistances in the film are important. In M3 the film resistance against the heat transport and the intraparticle resistance against the mass transport are included. M4 is the reflection of M3 with internal heat diffusion and film resistance against mass transport. M5 only includes intraparticle resistances and in M6 both film and intraparticle resistances exist. M6 is the most complicated
model and agreement within a parameter domain between this model and one of the less complicated models indicates that the simpler model includes the relevant transport resistances. In Table 1 the mathematical formulation of each model is listed as differential equations. The dependent and independent variables and the parameters are all dimensionless quantities. In M6 5 parameters determine the steady state, while the 2 remaining parameters only appear in the dynamic description. The other models contain a smaller number of parameters which will appear through combination of the parameters in M6. Table 2 contains the definitions of the par-
Table 1. Mathematical representation o f the models for the catalyst pellet Model number M1
Differential equations ~Xk
e2--= Ot
Boundary conditions Xk(1) = I
V2Xk--tl~Xk
dxk
e2~-~- = 3o"( l - - x~ ) - - a x k e y(~-l/uk)
M2 d y k __ L e ~ - - - 3v( 1 - - y e ) + a f t er(l-llUk)Xk
OXk e . ) - - ~ ~2X k -- OfXk eY(1--11Yk) " Ot
M3
Xk(1 ) = 1 L e - d- ~y k -_- 3v(l --Yk) + 3 a f t e ytl-l/ut? f x k r 2 d r ~o
e2--~d x k = 3o" ( 1 -- xe) -- 3axk f o 1 e y~l-l~upr2 dr
M4
yk(l) = 1 L e -~ - = V2yk -}- otj~xk e't(l-ltU~ Ot OXk = V 2Xk - - OLXk e ~(1-lly k) ~2 O--t-
x k ( 1) = 1
L e v ~ = V2yk + aflXk C ( 1 - 1 t ~ Ot
yg ( 1 ) = 1
M5
M6
OXk ~2 _~ = V2Xk - - aXk eTtl-ltvk} Ot
VXk(1 ) = o ' ( l - - x k ( 1 ) )
L e - ~ = V 2yk + aflXk e yll-llul~ Ot
V y ~ ( 1 ) = v ( I -- y~ (1) )
1557
K. W. HANSEN Table 2. The parameters Parameter
Definition
Variationinterval
R 2 . kfe -y
a
0"01-100
De
/3
(-AH)c~oDe Tuok~
y
E R~T~
10-40
tr
kc.R De
20-20,000
v
h.R ke
1-200
0"0005-0"5
(pCp)k. De
Le
or is a modified Sherwood n u m b e r and v a modified Nusselt number. "Modified" signifies that the quantities contain parameters defined in different mediums. Large transfer coefficients are equivalent to small changes in potentials across the film. L e is a Lewis number for the catalyst, defined as the ratio between mass and thermal diffusivity. T h e values o f the parameters used in the following examples are summed up in Table 3.
0"5-100
ameters, and boundaries are stated within which one will normally expect the parameter values to be in physically realistic situations. R2/De is a characteristic diffusion time for the pellet and ks. e -~ the reaction rate at the gas temperature in the main stream. Consequently a is the ratio between the diffusion and the reaction time at the surface conditions. V~da= q~ is normally denominated the Thiele modulus. T h e parameter/3 is proportional to the reaction enthalpy, and for steady state condition/3 is the maximum dimensionless temperature rise inside a pellet when only internal resistances against heat and mass transport are considered. Y is proportional to the activation energy, which expresses the influence of the temperature on the reaction rate.
ANALYSIS OF THE STEADY STATE T h e treatment of stationary catalyst problems on a theoretical basis has dominated the literature during the last years. One has especially worked with the models M5 and M6 and tried to find the number of solutions as well as the stability of these and domains of existence and uniqueness as a function o f the parameters. In the early part o f the 1960's Tinkler and Metzner [14] and Weisz and Hicks [15] presented numerical computations which showed that 3 solutions could appear in M5 for sufficiently high values of the product y . / 3 within an interval of a. Hatfield and A d s [ 16] have used values for a,/3, 3' in the same range and calculated the steady states for M6, and they find that up to 5 solutions may appear when o- > v. Still with M6 Villadsen[17] and Cresswell [ 18] perform computations using a, /3, y values, which only give 1 solution for M5. F o r much lower y a n d / 3 values 3 steady states may now appear. Regarding the validity of the simpler models most of the published work only contains qualitative considerations dealing with the dominance of different transport resistances under various
Table 3. Values o f the parameters u s e d in preparation of Figs. 2 - 1 3 Figs.
c~
fl
Y
tr
v
2, 3 4, 5 6, 7 8, 9 11 12, 13
100 1 1 1 1 0-5
0"01 0" 1 0" 1 0"01 0"01 0"01
20 40 40 20 20 20
250 250 250 250 400 400
12-5 12-5 12"5 12"5 5 5
Le
1-100
~.z
U
W
A Y
AX
5 2"5
100 5
75 150
120 240
0.5 0.5
In the last Figs. (11-13) the stated values of a, fl and y refer to reactor inlet conditions and elsewhere to pellet surface conditions. 1558
Analysis of transient models
conditions. In the following two sections it is tried to obtain more exact information by computing the stationary profiles for all the formulated models using realistic parameter values inside the variation intervals shown in Table 2. COMPUTING
METHODS
The mathematical formulation of the stationary models appears from Table 1 by setting the accumulation terms on the left side of the equations equal to zero. M 1 is easily solved. The other models are nonlinear and more than one solution to the equations may occur. For M2 and M3 it is possible to construct algebraic algorithms, guaranteeing all solutions to be found, and then it is easy to determine the concentration and temperature profiles. The methods are described in Ref. [ 19]. For the solution of M4, M5, M6 orthogonal collocation is used in connection with a NewtonRaphson technique to solve the resulting algebraic equations. The method has been introduced by Villadsen [17, 20] and it has proved to be superior to a number of alternative methods. Orthogonal collocation is a solution method for differential equations and is one of several Methods of Weighted Residuals. The differential equation is satisfied in N predetermined points of the variation interval of the independent variable. These points are the zeros of an orthogonal polynomium of the order N. By choosing a suitable weighting function for these polynomials it is possible to fix the zeros at points where it is important to obtain an accurate solution. The catalyst problem considered is symmetrical, and an even set of Jacobi-polynomials can be used. The gradient Vy and the Laplace operator V~y are found as weighted averages of the value of the dependent variable at all collocation points and at the pellet surface:
The matrices A and B are only dependent on the chosen polynomials and can be found once for all. Details about the procedure and about the method itself can be found in Refs. [ 17] and [20]. In this case a weighting function (1 - F ) . r z is used. For the three models M4, M5 and M6 the number of solutions is determined by an approximate procedure consisting of orthogonal collocation with two internal points. Details are given in Ref. [ 19]. The two-point method produces fine results for M5 but should be used with caution for M6. It may be mentioned incidentally that Stewart and Villadsen [21] report to have obtained satisfactory results for model M5 with a method of lower o r d e r - a one-point orthogonal collocation. When the number of solutions and the corresponding approximative profiles are found, a high order method is used for the exact determination of the profiles. RESULTS
In Figs. 2-5 some examples are shown where the concentration and temperature profiles for
0"75 x x
0-5O a=lO0 o-=250 ~ + l
~-~1 u= 12.5
~=20
MI M2
(>25 l e l
MS
~×~
M6
N+I
Vy(l) = ~, A ( l , J ) y ( J ) J=l
I 0"~
N+I
0"50
0"75
r
V~y(1) = ~, B ( l , J ) y ( J ) .
Fig. 2. Concentrationprofile.Steady-state.
J=l
1559 C E S Vol. 26, No. 1 0 - - E
0"25
I.O
0-00
K. W. HANSEN concentration inside the pellet, whereas the temperature only shows a modest increase. It is observed that the concentration profiles for M 1, M3, M5 and M6 are close to each other, whereas M2 deviates strongly. F r o m this it is deduced: (1) the thermal effects have only little influence on the concentration profile, (2) the major mass transport resistance is distributed as diffusional resistance inside the pellet. As the mean temperature in M3 and M5 deviate from M6 to approximately the same degree it appears by observing the temperature profiles that the film resistance against the heat transport is of the same order o f magnitude as the intraparticle resistance. T h e s e results are true in many other cases, which means that M2 is a poor model in most practical cases. M 4 has been omitted from the figures because it will give the same result as M2. It is obvious that the concentration drop
I "020
I.OI5
× >, 1.010
IOO5
a-lO0 o-=250
M2 1.000' 0.00
~=0.01 v=12"5
~a,~
"),"• 20
\ M3
J 0"25
~l~
M5
~x~
I 050
M6 I 075
1.0
Fig. 3. Temperature profile.Steady-state.
,.or o.,,I-
+
MI 0"7
M2 M3 M5
0-6
05 0"0
~
~
M6
I
I
I
I
i
0"2
0"4
0-6
0"8
I'0
Fig. 4. Concentration profile. Steady-state. the various models are compared. F o r certain combinations of the parameters the solutions to all models are found and plotted in the same figure. Figures 2 and 3 are concentration and temperature profiles, respectively, computed for parameter values, which only generate one solution in all models, D u e to the large value of the Thiele modulus there is a pronounced decrease in the
across the film is much smaller than the decrease inside the pellet. T h e profiles in Figs. 4 and 5 are calculated for the same values of cr and v but for larger values o f the thermal parameters/3 and y and a smaller value of a. It would o c c u r for larger values o f the heat o f reaction and activation energy and for smaller catalyst dimensions. It appears that also in this case M2 deviates completely from the
1560
Analysis o f transient models
1.05 ~
x
~
a-I
x 1-04
:~ 1-03
~
x
M2
y-40
~A~
o'-250
~e~
M3 M5
v= 1 2 " 5
~x~
M6
°~o
1.02
x
1.01
I'00 0.0
•
I 0"2
I •
I
0"4
06
I 0"8
~'~,Q I~)
t
Fig. 5. Temperature profile. Steady-state.
The conclusion of these remarks and a series other models. The increase of temperature is still modest, but owing to the strong influence of of computations presented in Ref. [19] is that temperature on the reaction rate it gives rise to for mass transfer the internal diffusion resistance a considerable decrease in concentration. Here will dominate the external film resistance, whereas the internal and external resistances are often M5 gives a better result than M3. While only one solution was obtained in all of equal magnitude for the heat transport. Thus models for the parameter values used in Figs. 2 in certain cases it will be possible to substitute and 3 the parameter values in Figs. 4 and 5 will M6 with M5 and in other cases with M3. generate three solutions in M2 and M3, but only ANALYSIS OF THE DYNAMIC STATE one in M5 and M6. It is the "lowest" steady state Very few computational studies of transients in M2 and M3 which is shown in the figures. The in a single catalyst pellet have previously been "higher" steady states give very high conversion reported. and temperature. It is now obvious that one has Hlavacek and Marek[22] present computato be very careful when using simpler models tions in which M5 is used. They investigate the for the catalyst pellets in a fixed-bed reactor, influence of the value of the Lewis number when because the conditions and parameters (a,/3, y) different initial conditions are used in both the vary along the reactor. To be valid within a multiand single-solution domains. parameter range the simpler models must provide Ferguson and Finlayson[23] use almost the the same number of solutions and approximately the same spatial mean value of temperature and same computational technique as in this paper to calculate transients for M5 and M6. They apply concentration in each steady state as M6. As mentioned above, one might very well find orthogonal collocation combined with an exmore than one solution to M5 and M6. The film plicit integration method in the time domain resistances--especially against heat transport-- using polynomials for each model. Orthogonal collocation is also used in this are much more important in the "higher" steady states. For these solutions there will exist an study. The polynomials are even and belong to agreement between M6 and M3 rather than the Jacobi family as in the stationary case. The resulting ordinary differential equations are between M6 and M5. 1561
K. W. HANSEN integrated by a fourth order R u n g e - K u t t a method with a variable step size. T h e adjustment o f the step size takes place through calculation of the truncation error according to principles, which are recommended by Zonneveld [24]. T h e computations which are presented below only show the results for models M1, 3, 5, 6, since as mentioned in the previous sections M2 and M4 are not applicable in the relevant parameter domain.
and concentration for the models M 1, 3, 5, 6 are shown. As initial conditions profiles deviating from the steady state have been inserted. 8 collocation points have been used and this proved to give an accuracy o f at least 3 - 4 figures. Figures 6 and 7 show the transient response o f the spatial mean concentration and temperature after x = 0.8 and y = 1.15 have been inserted as initial values at all the collocation points. T h e stationary parameter values correspond to those of Figs. 4 and 5, and only one stationary state occurred in M5 and M6. T h e diffusion resistance inside the pellet is dominating
RESULTS In Figs. 6 - 9 some transients o f temperature I'O--
MI
0"8
A~
~
01=I B-0.1 o"=250 v=12.5 Le= I
g 0-6 x
o
~
"
× 0.2
0.¢
-00
I
0"05
I
I
O"I0
01,5
I 0,20
t
Fig. 6. Mean concentration. Transient response. Initial values: xk = 0.8, y, = 1.15. 1.20 -
;:0.' )'=40 ~" I't0
1.00--
000
V=12"5 Le=l C2=0-5
I
0"05
I
0-10 t
-A"
I
015
¢
020
Fig. 7. Mean temperature. Transient response. Initial values: xk = 0.8, Yk= 1.15. 1562
Analysisof transientmodels
oa
× Le=l
~
x
~
×_
x
/
B-o~, /
/ @4
Le-5
i
e
,'~,*'zs° )t=20
~
1,'=12.5
+
x
I
o.2
(>0
t
~
+
Le=lO0 I 0.1
I 0.2
I 0-3
I 04
•
(>5
t
Fig. 8. Effect of the Lewis number. M e a n concentration. T r a n s i e n t response of M6. Initial values: xk - Stationary profile, yk = 1.3.
1.3
La=lO0 Le=25
~
+
i
°-,
.8=001
~"
e
~
• I" I --
I .0 043
X ~ L e =I
I 0.1
I 0.2
I 0"3
7-20 o'=250
~ e
-0.4
vii2'5
E 2 =0.5
0.5
t
Fig. 9. Effect of the Lewis number. M e a n temperature. T r a n s i e n t response of M6. lnitial values: xk ~ Stationary profile,yk = l '3.
for both heat and mass transport. In the stationary state M5 is thus preferable to M3. The temperature perturbation causes the mean concentration in Fig. 6 to decrease heavily at the beginning of the transient for all the models shown, except of course for the isothermal model M 1, which goes straight back to the stationary state. M3 has a smaller minimum value than M6, but agrees better with M6 than M5 at the beginning of the transient. Later on, M3 goes quickly back to its inaccurate stationary state so that M5 is the best model at the latter part of
the transient. The temperature plot shows a maximum due to the strong heat generation. During this period it appears that also with regard to temperature M3 agrees better with M6 than M5 does. Consequently it must be assumed that the dominating resistance against the heat transport is to be found in the external film at the beginning of the transient, whereas the internal heat diffusion resistance later on gradually becomes dominant just as in the stationary state. The conclusion of these and other results
1563
K. W. H A N S E N
presented in Ref. [19] must be that obviously none of the simple models are sufficiently correct to substitute M6 throughout the parameter space investigated. For the computations of transients in a reactor it will therefore be necessary to use M6 as catalyst model. Actual reactor computations will show whether considering a large number of coupled catalyst pellets might after all make the use of simpler models feasible. Finally Figs. 8 and 9 show the influence of variations in the Lewis number on the transient response of M6 after a disturbance in the stationary temperature profile. For large Le-values the temperature dynamics is very slow causing a heavy decrease in the concentration. In that case it is possible to use a simpler method. Presumably a quasistationary solution, in which the accumulation term in the mass balance has been eliminated will be acceptable apart from the very first period.
Catalyst U OxA.
(III)
W Oyk = V 2 y k q.. otflXk e y(1-1/uk) Or
(IV)
0~" = V2Xk-- axk e ytl-1/uk)
where U=
E2. R 2 • Vs De. L
W = (pCp)k. R z • vs ke. L
The boundary conditions are Oxk(1 Or "-) = o'(xo--xk(1))
(V)
Oyk -~-r (1) = u(yg--yk(l))
(VI)
together with the symmetry conditions ANALYSIS OF THE REACTOR DYNAMICS
In this investigation the interest is mainly concentrated on the formulation of models for the catalyst. Therefore a number of simplifying assumptions will be introduced for the other elements in the reactor: (a) Adiabatic operation of the reactor. (b) No heat capacity of the reactor wall. (c) The density of the gas is constant and independent of temperature. (d) Other physical quantities are constant. (e) No axial and radial mixing in the external medium. With M6 as catalyst model the describing equations will be the following: External gas medium OYuo~.t-VO0-~ = A Y ( y k ( 1 ) --Y~')
(I)
Ox~ t _ v O x g = A X ( x k ( 1 ) - - x o ) Or Oz
(II)
where AY=
a.h.L ~
• (oC,,)o
AX=a.ke.L . v,
E~ . v ,
oxk -~r (0) = 0
(0) = O.
(VII, VIII)
The dimensionless parameters U and W are the ratio of the diffusion time for the mass transport to the reactor residence time and the diffusion time for the heat transport to the reactor residence time, respectively. High values of U and W cause the catalyst dynamics to be much slower than the transportation lag in the external medium. The ratio W~ U is equal to Le/e., and is thus of importance for the dynamics of the single pellets. Three independent variables occur in the m o d e l - t w o in space and one in time. The catalyst concentration xk and temperature Yk are functions of space variables r and z whereas the gas concentration xo and gas temperature yg only depend on the axial position z in the reactor. With the assumptions previously introduced the coefficients A X , A Y, U, W, o- and u are constants. This is a rather crude assumption, which, however, may be justified here since the basic transport processes are the central investigation object.
1564
Analysis of transient models SOLUTION METHOD
The numerical technique used in this paper has been specifically developed to try to meet the requirement of a low consumption of computer time. The equations for the catalyst (III)-(VIII) are transformed to a set of 2N ordinary differential equations, as in the dynamic analysis of the single pellet, by using orthogonal collocation with N internal points. The collocation equations in point No. (I) are: dXk (1)
dr
1 N+I
=-U "~ B(I,J).Xk(J) J=l
--Uxk(I)
(IX)
e y(1-11u~/))
d ~ = A Y ( y k ( N + 1) - y g ) .
The characteristics of the equations are plotted in the (r, z) plane in Fig. 10. Each Cline here represents the time variation of the concentration and temperature profiles at a certain point in the reactor. Similarly the variation of concentration and temperature in a flow element in the external medium is recorded by following a g-line from z = 0 until z = 1. The surface values Xk(N+ 1) and y k ( N + 1) can be eliminated by the boundary conditions (XI) and (XII). (IX), (X) and (XIII), (XIV) can then be rewritten. (IX) and (XIII) thus become:
d x k ( 1 ) - - ~ BB(l,J) xk(J)+CC(l) xo dr
dYk (1)
d'r
"
J=l
I~ ~ B(I,J) . Yk(J) + wXk(1)
O/
- - "-~Xk ( 1 ) e ytl-vuktm
eY(1-11Yk(l)).
N+I
(J) =
l/)
(XI)
,]=1 N+I
~, A (N+ 1,J)y~(J) = ~'(yg--yk(N+ 1)). (XlI) ./=1
The solution of these equations and the coupled equations for the external medium is performed by the method of characteristics[26]. The characteristics corresponding to the catalyst equations are straight lines parallel to the r-axis in the (r, z) plane, whereas the characteristics corresponding to the equations for the external medium are curves, which for constant flow are straight lines with a slope equal to V. With the relation dz/dr = V the equations for the external medium can be rewritten as (XIII) and (XIV).
~ HH(J). xk(J) + H H ( N + 1)xo. dr = J=l (XVI)
In an adequate number of points on the z-axis the Eqs. (XV) and its temperature analogue are set up. These equations must be integrated along the Ccharacteristics whereas (XVI) and its temperature analogue must be integrated along the 8-characteristics. The calculations are started along the 8characteristic, which is nearest to (r, z) = (0, 1). The integration from one intersection of this line ~-Iines ~0
g
- lines
o
0-0 Title,
dxu = AX(xk(N + 1) --xg)
(XV)
dx v
(X)
The boundary conditions for r = 1 are transformed to
dr
"
1 N+I J=l
Z A ( N + 1,J)x
(XIV)
(XlII) 1565
l"
Fig. 10. Characteristic curves for the reactor.
K. W. H A N S E N
with a ~-characteristic to another intersection is performed by a second order predictor-corrector scheme except for the first step, where it is necessary to use a first order predictor. After having finished the calculations along this line the integration is continued correspondingly along the nearest 8-line. The developed computer program is able to calculate the transients after step changes and sinus disturbances in the inlet concentration and temperature as well as step changes in the flow. A considerable part of the program is devoted to administration of the enormous data material produced in just one single transient. During the computation about ¼ million points are thus passed in the (r, z) plane. The computation of each transient is made in several stages, between which step length in time and space may be adjusted if required.
proximately reached the stationary profile of the isothermal model. Then it changes slowly following the increase in temperature which gradually causes the concentration profile to change curvature. The dotted profile in Fig. 11 is the stationary profile when using M5 as catalyst model instead of M6. The profiles are quite close to each other and consequently the film resistance against the heat transport is of only minor importance for the conversion. If the catalyst profiles are studied during the transient period it will appear that this result is not due to the fact that the external film resistance against the heat transport is insignificant compared with the internal resistance, but the increase of temperature in one pellet is rather small compared with the increase of temperature along the reactor. For this computation a step size of 0.02 in the axial direction was used together with time increments of 0.05 at the beginning and later on 0.1.
RESULTS
In the following some results are shown which illustrate the most important dynamic characteristics of the reactor. Furthermore it will be possible to draw conclusions about the applicability of the computation method. The parameter values for each of the examples are summed up in Table 3. By computations, described in Ref. [19], it has been ascertained in all cases that only one stationary state exists for the catalyst pellets in the inlet. Thus only one possible stationary state exists for all pellets in the reactor, as the local 3/and/3 values will both be smaller than 3/and/3 calculated on inlet conditions. The stationary state for the reactor was found by introducing initial values equal to 1 for concentration and temperature in both the catalyst and the gas medium and observing the transients until no change took place. Figure 11 shows an example, in which the concentration and temperature profiles in the reactor are plotted. The ratio between W and U is 20. The concentration profile in each single pellet will therefore at the beginning of the transient change much faster than the temperature profile. After 5-10 residence times the concentration profile has ap-
1"08
~®_
/:////
I
r~4
1"06 >'~ 1'04 1"02
I
1.00 .
r.O.4
r=O-8
------ M5 0"92
Q=I r-2o 0"=400 0.6
0 0"820-00
T.,J
N~
AX-,20
X\
"
8 0"25
4 0"50
0"75
I'00
Z
Fig. 11. T r a n s i e n t response. Reactor profiles. Initial conditions: x~ = x, = Yk= Y u = 1.
1566
Analysis of transient models
The catalyst profiles were computed with two collocation points. On an IBM 7094 the computations could be performed in about 80 min. Figure 12 shows the result of a perturbation of the inlet temperature from 1 to 1.1 during 6 residence times. The stationary profiles have been obtained by an integration procedure similar to the one mentioned above. The parameter values are somewhat different from those on Fig. 11. It is especially noted that the heat capacity of the catalyst is a factor 10 smaller. This causes a much faster dynamic behavior. (Le/e2 = 2). After the 6 residence times the stationary state a is changed to the b-profiles. A temperature maximum has developed in the inlet while the concentration drops throughout the reactor. Hereafter, it will be seen that the temperature wave is moving slowly towards the reactor outlet with a velocity approximately equal to El (pCp)o/ [ (1-- ~1) (pCp)k]. Gradually the wave changes its form becoming wider and getting a still higher maximum value. The temperature perturbation is immediately followed by a drop in the outlet concentration, whereas the outlet temperature only changes slowly. At the beginning of the transient the rate of conversion will thus be lower in the last half of the reactor,
and consequently the heat generation will be smaller in this part too. During this period the outlet temperature will therefore drop and reach the minimum value after 22 residence times. The deviations from the stationary profiles are pronounced for the case shown in Fig. 12. In order to investigate whether the temperature difference between surface and centre of the single catalyst pellet still remains small relative to the overall temperature changes in the reactor, the time variation of the catalyst profiles is shown at the reactor outlet in Fig. 13. The pulse disturbance in the inlet temperature causes the concentration and temperature values to decrease until r = 27, whereafter the temperature starts increasing rapidly. After 48 residence times the temperature approaches the maximum value, and at the same time the concentration profile is very steep. During the part of the transient shown in Fig. 13 the temperature of the external medium is higher than the temperature of the pellets, and
r=48
14 45 1'3 ,~
i.z
41 36
115f
g/~h Profile
>,~ 1,4
Time
o b c d e f 9 h -6 0 6 22 32414851
1,3
d
e
Iq i l0
I
I
~
/
~
075 Stationary state
/
o
-
x 0"8
a=0.5 b ~=0.01 AY=ISO 7"=20 AX=240 0"61-- o",400 U=2"5 U=5 W=5
~
j
j
~ 0,25
@4 0"00
I
"r=O
f ~
F2
I'0
\27
Stationary state and
I 0"25
0"50
0"75
I'00 ~
Fig. 12. Transient responsel Reactor profiles. Initial conditions: Stationary state. Inlet temperature equal to 1.1 for 6 residence times.
1567
:i i .J. 025
0.50 r
0-75
Fig. 13. Catalyst profiles at z = I for the case in Fig. 12.
[
K. W. H A N S E N
thus the temperature decreases towards the centre of the pellet, but the profile is, however, still almost fiat. Integration step sizes between 0.01 and 0.005 were used in the axial direction and from 0.01 to 0.04 in the time direction. The computer time until r = 51 amounted to 50 min.
NOTATION
A A (I, J)
AX AY a B (1, J)
CONCLUSION
The reactor dynamics for gaseous reaction media are dominated by the heat transport processes. In the single catalyst pellet the concentration profile will often follow the temperature profile in a quasistationary state. A temperature disturbance moves like a slow wave through the reactor due to the small heat capacity of the gas medium compared with that of the catalyst. The temperature variations in the individual pellets are rather small compared with the temperature changes throughout the reactor, at least when only one stationary state of the catalyst is possible for all pellets in the reactor. These observations make it possible to use simpler models, and the validity of these models will be reported elsewhere (Ref. [25]). The numerical technique presented has proved efficient for the computation of transients in reactors. The method seems to allow a considerable reduction of the computation time as compared to previous works. Furthermore the computations have given important information about the best way to build up the routines just as it has become obvious that certain modifications of the method will reduce the computation time even more. This is due to the fact that stability problems arise when very steep concentration profiles occur in the catalyst pellets. This may be remedied either by changing the weight function or by a segmentation of the variation interval of the radial pellet cordinate using a spline-collocation approach. Acknowledgements--The author wants to express his gratitude to Messrs. J. Villadsen, M. Kiimmel and H. Livbierg for helpful discussions during the work and for their valuable comments on the presentation of the results in this paper.
BB(I,J) Cp
CC(1,J) c De E --AH
HH(1)
reactant weight coefficient in the expansion of first derivative
akcL/(elvs) ahL/(e,(pC,)ov,) catalyst surface/reactor volume = 3(1--Ea)/R weight coefficient in the expansion of the Laplace operator
(I/U)(B(I, J ) - - A ( N + I , J ) . B ( I , N + 1)](o'+A(N+ 1 , N + 1))) specific heat capacity
B(I,N+I).o-/(U(o-+A(N+I, N + 1))) concentration of A effective diffusion coefficient activation energy heat of reaction
- - A X . A ( N + I,I)/(o-+ A ( N + I,
N+I)) h heat-transfer coefficient kc mass-transfer coefficient ke effective heat conductivity kt frequency factor L length of reactor Le Lewis number (pCp)kDdke I position in reactor N number of internal collocation points R radius of catalyst pellet Ro gas constant r dimensionless radius r'/R or reaction rate r' radial position in catalyst T temperature t dimensionless time t'De/R 2 t' time
U V
~zR2vs/(De. L) v/v,
v average velocity in reactor
W x y z
(pCp)kR2vs/(ke . L) dimensionless concentration c/coo dimensionless temperature T/Too dimensionless position in reactor 1/L
Greek symbols
1568
Ot
R 2 .
ks.
e--E/~Roroo)]De
Analysis of transient models
fl (--AH)cgoDd(Tooke) Y E/(RoTuo) 8 def. in Fig. 10 def. in Fig. 10 E1 void fraction of reactor E2 void fraction of catalyst p density Thiele modulus V ~ v hR/ke [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21 ] [22] [23] [24] [25] [26]
or
kcR/De
1"
dimensionless
time t'vs/L
Subscripts gas catalyst 0 s u r f a c e c o n d i t i o n s ( F i g s . 2 - 9 ) o r inlet conditions (Figs. 11-13) s stationary
g k
REFERENCES L1U S. L. and A M U N D S O N N. R., lndl Engng Chem. Fundls 1962 1 200. L I U S. L., A M U N D S O N N. R. and ARIS R., ibid. 1963 2 12. LIU S. L. and A M U N D S O N N. R.,ibid. 1963 2 183. V A N D E R V E E N J. W., LUSS D. and A M U N D S O N N. R.,A.I.Ch.E. J11968 14 637. D E A N S H. A. and L A P I D U S L., ibid. 1960 6 656. M c G U I R E M. L. and L A P I D U S L., ibid. 1965 1185. F E I C K J. and Q U O N D., Can. J. chem. Engng 1970 48 205. C R I D E R J. E. and FOSS A. S.,A.I.Ch.E.J11966 12 514. C R I D E R J. E. and FOSS A. S., ibid. 1968 14 77. S T A N G E L A N D B. E., Ph.D. Thesis, University of California, Berkeley 1967. L Y C H E B. J., Ph.D. Thesis, University of California, Berkeley 1968. H O I B E R G J. A., L Y C H E B. J. and FOSS A. S., Experimental evaluation of dynamic models for a fixed-bed catalytic reactor, Paper presented at The 1st Int. Symp. on Chem. Reaction Engng, Washington, June 1970. SINAI J. and FOSS A. S.,A.I.Ch.E. J11970 16 658. T I N K L E R J. D. and M E T Z N E R A. B., Ind. Engng Chem. 1961 53 663. W E I S Z P. B. and HICKS J. S., Chem. Engng Sci. 1962 17 265. H A T F I E L D B. and ARIS R., ibid. 1969 24 1213. V I L L A D S E N J. V., Selected Approximation Methods for Chemical Engng Problems. Danmarks tekniske HCjskole 1970. CRESSWELL D. L., Chem. Engng Sci. 1970 25 267. W , ~ D E H A N S E N K., Ph.D. Thesis (in Danish), Danmarks tekniske HCjskole 1970. V I L L A D S E N J. V. and S T E W A R T W. E., Chem. Engng Sci. 1967 22 1483. STEWART W. E. and V I L L A D S E N J. V., A .I.Ch.E. J11969 15 28. H L A V A C E K V. and M A R E K M., paper presented at The 4th European Syrup. on Chemical Reaction Engng, paper 2.5, September 1968. F E R G U S O N N. B. and F I N L A Y S O N B. A., Chem. EngngJ. 1970 1 327. Z O N N E V E L D J. A.,Automatic Numerical Integration. Mathematical Centrum, Amsterdam, 1964. W/EDE H A N S E N , K. To be published. AMES W. F., Nonlinear Partial Differential Equations in Engineering. Academic Press, New York 1965. R 6 s u m 6 - D a n s la simulation de la dynamique du r6acteur catalytique b, lit fix6, il est important d'utiliser un mod61e correct pour le comprim6 de catalyseur. Les diff6rents mod61es de catalyseur sont test6s d'apr~s une analyse statique ainsi que dynamique. Les r6sultats obtenus avec les deux modEles sont ensuite compar6s. Enfin, des calculs de transition sont faits dans un r6acteur tubulaire en utilisant un modEle de catalyseur qui tient compte des r6sistances de film et de particules contre le transport de masse et de chaleur. La m6thode de calcul est une combinaison de la collocation orthogonale et de la m6thode des caract6ristiques. Cette m6thode permet d'effectuer les calculs, 5. l'alde d'un ordinateur, dans un temps raisonnable. Zusammenfassung-Bei der Nachahmung der Dynamik des katalytischen Festbettreaktors ist es besonders wichtig ein richtiges Modell fiir die Katalysatorkugel zu verwenden. Die verschiedenen Katalysatormodelle werden durch eine station~e sowie durch eine dynamische Analyse gepriift. Die mit Modellen verschiedener Art erhaltenen Ergebnisse werden dann verglichen. Schliesslich werden Berechnungen von f0bergangsvorgSngen in einem Rohrreaktor unter Verwendung eines Katalysatormodells ausgefiihrt, das sowohl F i l m - a l s anch Zwischenteilchenwiderst~nde gegen W ~ m e - u n d Stofftransport einschliesst. Die Berechnungsmethode ist eine Kombination orthogonaler Zusammenstellung und der Methode der charakteristischen Merkmale. Mittels dieser Methode kiSnnen die Berechnungen unter Aufwendung eines verniinftigen Ausmasses an Computerzeit durchgefiihrt werden. 1569