- Email: [email protected]
Moving reaction zones in fixed bed reactors under the influence of various parameters
ChemicalEngineering Science, 1972. Vol. 27, pp. 1485- 1496.
Pergamon Press.
Printed in Great Britain
Moving reaction zones in fixed bed reactors under the influence of various parameters D. VORTMEYER and W. JAHNELt Technische Universit;it Miinchen, Germany (Received 28 September 197 1) Abstract-Moving reaction zones are obtained by the solution of the interdependent and unsteady state set of the energy and material balance differential equations for a fixed bed catalytic reactor with a fast exothermic reaction. Parametric studies concerning the influence of gas velocity and reactant concentration, of heat and mass diffusivities, of heat radiation, of the kinetic parameters and of the particle diameter were performed. A comparison with experimental work turned out to be very satisfactory.
IT IS established that in the case of fast exothermic reactions in fixed bed reactors the total conversion takes place over a few particle diameters range known as the reaction zone. Experimental investigations by Wicke and Vortmeyer [l] have shown that the reaction zone may move through the turbular reactor in either axial direction (Fig. 1). It was further demonstrated that direction and velocity w of the zone depend heavily on the gas velocity u. and on the reactant concentration cGo. The experimental relationship is:
where CYand p depend on other parameters of the system. Positive values of w denote movement against the gas flow while negative values indicate that the reaction zone is being blown out of the reactor. The exponents 0.77 of u,, and 0.5 of cc,, were later confirmed by further measurements by Wicke et al. [9, lo]. The movement against the flow is the most interesting feature and can only be explained by axial diffusion processes of heat and mass. Since it may safely be assumed that the transport of energy is the dominating process, a qualitative treatment of the energy equation in [l]
and an extended one in [2] by analogy to laminar flames gave very satisfactory results for the exponents of u. and cco. In addition, it led to interpretation of the coefficients OLand p in formula (1). ff is a coefficient of great complexity while /? turns out as the heat capacity ratio of fluid to solid. However many questions still remained open concerning the influence of various other parameters such as radiation, mass diffusion, kinetics, particle diameter etc. Since it is very difficult to change some of these parameters experimentally it seemed worthwhile to solve the complete set of the unsteady state differential equations numerically and thus simulate the dynamics of the system on the computer. Results by Vortmeyer and Jahnel[3] for the energy equation alone turned out to be very encouraging, since the exponents 0.77 of u. and 0.5 of cGo from the Computer solutions were in perfect agreement with formula (1). 1. THE EQUATIONS OF THE HOMOGENEOUS MODEL FOR A FIXED BED ADIABATIC REACTOR
Reactors with fast exothermic reactions exhibit substantial differences in temperature between solid and gas, particularly in the range of the reaction zone where the reaction rate is
tPresent address: Motoren-Turbinen-Union, Miinchen, Germany.
1485 CESVol27Na8-A
D. VORTMEYER
and W. JAHNEL
= a(g,aT/ax)/ax-
a(upGcpGT)/ax
.
- a(Z - K)/ax - AH . i.
T,=473
(2)
Although not necessary from the computational point of view, the average values A$, and c were introduced in the energy equation and the heat capacity of the gas was neglected in Eq. (2).
[K]
Energy balance: P,C~aTiat = ~&a2Tlax2 --uopco~pc aTlax - a( 1 - K)/ax - AH?. -5
-
-6
-
-7
1 5
0
1 IO
15
10%
0.
I 20
I 25
30
[m/secl
Fig. 1. Reaction zone velocity as a function of u0 and cGO,
very high. Nevertheless, the one dimensional homogeneous model was applied. Since the movement of the reaction zone depends on heat flux from the reactant to the non-reactant part of the bed, the application of the method appears to be appropriate in this range of relatively low temperatures with only minor temperature differences between solid and gas. The following additional assumptions were made: (i) A quasi homogeneous reaction is assumed to take place with no change in the number of moles. (ii) The axial conduction and diffusion processes under unsteady conditions may be evaluated from the Fourier and Fick equations. (iii) The radiative flux is approximated by the two-flux model as developed in [4] for fixed beds. The equations obtained were: Energy balance: (P,c,+
Although confirmed only for isothermal conditions, the effective axial diffusion coefficient was evaluated from Pe’ = udlD& = const. and was applied to the non isothermal reaction zone. Then the following material balance is obtained. E,aylat = D& a2ylax2 - u,,ay/ax - i
(4)
with (3 The transport equations for radiation [4] are: dZ/d_r+ al = bC,( TI 100)4 + gK
(6)
dZC/d_r- UK = bC,(T/ 1OO)4- gl.
(7)
The coefficients a, b, g depend on the emissivity and on the radiative transmission number of the bed. a=(1+B)2+(1+B)(1-B)(1-+J1-B)2 (1 +B)2-(1 -B)2(1-E)2
(l+B)d (8)
b= (1+B)2-(1+B)(1-B)(1-e)~E,(1-B)2 (1+B)2-(1-B)2(1-e)2
(l+B)d /a\ \7/
g=
lphd
(3)
am
(1+BY+(1+B)(1-B) (1+B)2(l-B)*(l-_E)*
(I_E)
(0;;;;. (10)
1486
Moving reactionzones
E is the emissivity of the particle surface and B is a dimensionless radiation transmission number for a fixed bed [5]. 2. THE
BOUNDARY
CONDITIONS
In order to calculate the dynamic behaviour of the tubular reactor and particularly the transit velocity of the reaction zone, boundary conditions have to be fixed and must be chosen in such a way that they are in agreement with the experimental system. The experiments [ l] were performed in such a way that the transit velocity of the reaction zone was measured at a great distance from the inlet and outlet of the bed. Mathematically this means an infinitely long reactor with the reactions taking place somewhere in the middle. However, only the finite reactor was treated. Consequently boundary conditions had to be introduced at the outlet and the second derivatives were chosen for this purpose. The conditions for the radiative transport equation at inlet and outlet are that there are no net fluxes of radiation. The following system of boundary conditions arises: t=O
--m
T,s
T(x)<
TE
YGO5 y(x) 2 0 (11) tz
0
x=-m
y=yGO T= TG I-K=0
t 2 0
x=+m
(12)
a2yla9=o a2TlaX2 = 0 I-K=O.
(13)
For ignition to take place, temperature and concentration profiles T(x) and y(x) have to be present somewhere in the reactor between plus and minus infinity at time level zero. At all other times the only gas stream conditions we know of are the inlet concentration and inlet temperature when the reactor zone is far away. Since the
reactor system in the infinitely long reactor is an initial value problem, the conditions at x = + m would come out by computation as aTlax = aylax = 0 or aTlaX
= a2yla2 = 0
for the adiabatic system. However for computational est condition
reasons the weak-
(14)
a2Tla.e = a2ylax2 = 0 was used at the now finite end of the reactor. Preliminary calculations for various reactor lengths have shown that in the case of particle diameters d = 0404 m, a length of O-1 m was sufficient to give results which were almost in exact agreement with those obtained when using very long reactors. Therefore all calculations were performed with this shorter reactor length. The influence of the boundary conditions at x = +m was also tested, with no change in the results. 3. THE
NUMERICAL
METHOD
The physical problem is described by a nonlinear system of two partial differential equations of the parabolic type (3), (4) and two ordinary first order differential Eqs. (6), (7). These four interdependent differential equations require two different types of solutions. Several methods for solving non-linear equations that result in linear finite difference equations, which can be solved by existing algorithms and which do not involve excessive iteration, have been developed. In this paper only forward Taylor Series projections are used. These methods do not need any iterations to treat the nonlinearity. Only in the case of ordinary differential equations does an iteration have to be made in order to satisfy suitable boundary conditions. The previously mentioned system of differential equations causes a combined initial-boundary-value problem. The well-known Runge-Kutta method was used to solve the ordinary differential equations
1487
D. VORTMEYER
which describe the transport mechanism of radiation. For this purpose it is convenient to write the equations in the following form: d+/dx = q($, x) + PVJ.
(1%
In this case I,IJrepresents the dependent variable while x is the independent variable. The nonlinear term q($,x) contains the coupling expression and p is a constant (see Eqs. 16, 17). If the radiation flux from x = - 00 to x = + CQis considered, I/J= Z(x) and p = --a
q=
bC,(T(~)/10O)~+gK(x).
Otherwise $ = K(x), the radiation opposite direction, and
and W. JAHNEL
After solving the differential equations (6) and (7) by the Runge-Kutta method we obtain a new distribution for Z(x) and K(x). The boundary conditions Z--K = 0 are now checked. If the condition is fulfilled and two successive distributions differ by less than a predicted error, the iteration is stopped and the correct Z(x) and K(x) distribution values are found. If the condition is not fulfilled the obtained distribution values are used to continue the iteration process. The radiation flux appears in the energy balance equation in the form
(16)
e(Z(x) - K(x))lax = 2X,( T(x)/ 100)4 - (g - a)(Z(x) - K(x)) (18)
flux in the
which can easily be determined by subtraction of Eq. (7) from Eq. (6). After having found the p = a q = - bC,(T(x)l loo)4 - gZ(x). ( 17) radiational heat flux distribution and the reaction rate, the parabolic equations for mass and The initial conditions for x = -CO and x = + 03 energy balance can be solved by an explicit are Z-K = 0 in both cases. forward method for the next time step using a Let us consider an arbitrary time step t at Taylor Series projection. This method has the which the temperature and concentration distriadvantage that no iteration is necessary because bution over the reactor length are known. Now the non-linearity only corresponds to the points the energy flux and the reaction rate which in the previous time step. The possibility exists depend on the distributions can be calculated that numerical instability arises when using this for this time step. method and therefore a stability analysis is The calculation of reaction rate Eq. (5) is easy necessary. to perform whereas the computation of radiaFor solving the partial differential equations it tional heat fluxes is much more difficult. In this is more convenient to write the equation in a case two interdependent ordinary differential general form such as the following: equations are to be solved. Because of the expression bC, ( T(x) I 100) 4, which contains the W) = a$ + P4t + Y& + 84, = A$& -&0. ( 19) interdependent variables, an iteration becomes necessary. The coefficients (Y,/3, y, 6 and f are determined In order to begin the iteration, Z(x) and K(x) by the special differential equations. In case of have to be chosen. For this purpose we use a energy balance, the dependent variable 4 is linear distribution of the form equal to the temperature T and the coefficients Z(x) = K(x) are in which the boundary values are set as Z(x = 0) = C,( T(x = O)/ 100)4 and
a=0 f=
ZW = x~) = C,(T(x = xiv)/ 1OO)4.
p=15,&.
y=l.lOpG&~
- S(Z(x) - K(x))bx
The mass balance 1488
- AHi.
equation
6=--h&
(20) is treated in the
Moving reaction zones
same way. In this case we set C#J = y and the coefficients are a=0
p=1
y=uo
a,,=-2&+6
(25)
a=-p&.o
(26)
(21)
f=-;.
For the energy equation is:
The boundary conditions can be written as
For
balance the finite difference
2
and
lsj
T IJtl=-H~*J(~)+Ti~(l-2~~)
It will be remembered that the reaction rate, the nonlinear term of the equations, contains the coupling elements as well. Because of the properties of the forward difference method, iteration is unnecessary. Both of the partial differential equations are to be solved independently from each other. After transforming the partial differential equation into a finite difference equation by using a forward Taylor Series projection over three discrete points shown in Fig. 2 we obtain the following equation:
4 lJ+l
=
+
--1 (2
Ti-IJ
uOpG%G
&Es
;
he,,
At
p& M > (27)
+ T~+M
For the mass balance the finite difference equationfor G i s (IV- l)and 1 sj < ais: yrj=-ridAt+yid
>
(fi,j-Q.A.i-“i+l,.i+~+lJ -
Q-lki-1Jh+1
(22)
using:
+ Yl+l.r %+1
=
a
a-1J
-
--M
P -At 6
‘i-l
‘i
)
- (28)
(23) 1.x 2 b
(24)
Using these two finite difference equations we obtain a new temperature and concentration distribution for the time level ttAl. These finite difference equations have the disadvantage of becoming instable under certain conditions. The numerical stability depends on the ratio of At/Ax2 and the convergence depends on the absolute values of At and Ax. The stability conditions[6] for the equations described above are:
and 'i-2
At
Ax
‘i*l
‘i+2
X
Fig. 2. The rectangular finite-difference mesh.
1489
D. VORTMEYER
and each term should be approximately of the same size. The values chosen by the authors were: 3 Q-lJ _ _. aij+l
afd
12’
-4.
a&j+1
&+I,.4 _
--I
12’ Uij+l
12
(29)
and stability was obtained in each example. The increments derived from the stated stability conditions are: For energy balance: A-&_=
A&
and
AtT=p.
2u@Cc,G
(30) az
For mass balance:
A&J!& 2ull
and
Ax”, At, = 3%;
(31)
Having calculated the new T(x) and y(x) distributions for the time step tij+l the values for the reaction rate r(x) and the radiation fluxes Z(x) and K(x) have to be computed with the algorithm described previously. The old values at the time tidy are now replaced by the new values calculated at the time tij+l and the whole process is repeated. In this way the time dependent properties of the reactor are obtained. 4. RESULTS
AND
DISCUSSION
4.1 Transient behaviour Having ignited the reactor it took about 300-400 set real time to establish a uniformly moving reaction zone. The above mentioned time naturally depends on the chosen ignition profile which was a sine curve. The transient behaviour is shown in Fig. 3 for two cases, which differ in the end temperatures of the ignition profiles. While in case (a) the end temperature at time level zero is too low and therefore the fixed bed has to be heated up in order to obtain the steady moving reaction zone, case (b) presents the opposite situation. The computations were done with the following data: gas inlet temperature
. . . TG = 473 K
and W. JAHNEL
effective solid density . . . pS = 450 kg/m3 gas density at inlet temperature. . . . . . . . . pG = 1.29 kg/m3 specific heat of solid. . . . . c, = 1423.6 J/kg grd specific heat of gas . . . . . c, = 1071.9 J/kg grd particle diameter. . . . . . . . d = 0.004 m emessivity of solid . . . . . . .E = 0.8 gas heat conductivity . . haso = O-0514 W/m grd reaction enthalpy . . . . . AH = - 1.4 X lo9 J/kmol activation energy. . . . . . . . E = 96.3 X 106J/kmol frequency factor. . . . . . . . . k,, = 9100 kmol/m3 set reaction order . . . . . . . . . . n = 0.35 radiation transmission number . . . . . . . . . . . . . . B==0*2-. The kinetic data are similar to those which were measured for the ethane oxidation on a Pd/,4l,O, catalyst by Wicke et al. [71. 4.2 Axial
dijksion processes and moving reaction zone In former papers[ 1,2] the movement of a reaction zone against the inflowing gas was explained by the axial heat conduction. It also was pointed out, that there are certain analogies to the theory of laminar flames. The theoretical treatment however was restricted to the energy equation only. Radiative transport processes and in particular axial mass diffusion processes were neglected. Because of the considerable mathematical complexity of the fully stated problem (Eqs. 3-7) an analytical solution seems to be impossible. Therefore the influence of the various parameters was investigated numerically. 4.3 The influence of axial heat conduction and radiation In a first paper[3] on a numerical treatment the energy equation was solved separately by using a model reaction rate function dependent only on temperature. When the resulting calculated velocity w of the reaction zone was plotted against the gas velocity u0 a complete qualitative agreement with experimental results (Eq. 1) was obtained. The same agreement was
1490
Moving reaction zones
d = 0.004
[m]
No radiation
65C
0.02 0.03 0.04 005 0.02 0.03 0.04
0.05 0.06 0.07 0.09 0.09 C
x Fig. 3. Time dependent
[ml
A& = haI,, + 0.8 u,,pcoGd.
(32)
As is shown in the following figures, these results remain unchanged if the complete set of simultaneous Eqs. (3)-(7) is solved. Figure 4 shows w as a function of the inlet velocity u. with the concentration parameter YOO. Because of reduced computer time the diffusion coefficient in the mass balance Eq. (4) was set to zero in these calculations since one of the following sections contains a detailed treatment of the influence of this transport property. If the results of Fig. 4 are plotted in double logarithmic form of Fig. 5 In (uo+i
*
temperature profile after ignition. (a) End temperature (b) End temperature of ignition profile too high.
observed when the inlet concentration cm was varied. In all cases the experimentally determined exponents 0.77 and 0.5 for u. and coo were confirmed when the following expression was used for the axial effective heat conductivity.
w) =f(ln
uo)
(33)
0.06 007
008
0.
[ml
of ignition profile too low
the gradients and the intercepts of the straight lines again lead to the exponents 0.77 for u. and 05 for cG. This result is valid for calculations with and without heat radiation. That radiative effects do not change the shape of the curves may be explained in the following way: The energy transport responsible for the moving reaction is the heat flux from the reacting to the non reacting part of the catalyst bed. The temperature of this cross section is about 500 K. Since this temperature is still rather low for larger radiation effects heat conduction is dominant. Our data for radiation effects led to the ratio QL/Qs = 3 if Q, and Q, are conductive and radiative heat fluxes. Therefore the radiative contribution in this cross section means an additive term in the Eq. (32) for A& and does not change its dependence on u. as a first approximation. That the functional relation w =f(uo) of Eq. (1) is important for the exponents can be seen from Fig. 6. This contains calculated results with a constant value of 0.538 W/m grd for z Al-
1491
D. VORTMEYER
and W JAHNEL 4 D,.= 0 [m */set] k, = 91000
4 t
-
No rodiotion
---
Rodiotion
[k mol /m3 set]
J--
;
p2ti $ -%
’
1
s
4%
O-
h,\
‘,
‘\
%R \
0006
\
\
-
‘\
’\
\
\
\
\
-I
0.05
\
\
0.007
0.006
\
\
\
\
0
0.009
\
\
\
\
0.15
0.10
0.20
0
Fig. 4. Calculated values for w as a function of u0 and cc0 with and without radiative contributions.
I
I
45
I
,,I,,
IO
1
I
15
20
"0
Fig. 5. Double logarithmic plot.
though no dramatic changes appear, there is an obvious difference in the exponent for u0 which now is 0.73 instead of the observed value 0.77.
0.04
0.08
0.12
u,
Fig. 6. w =f(uo)
[m/s4
0.16
for a constant and u,-dependent
-I
0 .2 0
value of
The injluence of axial mass d$usion Axial mass diffusion processes in isothermal fixed beds have attracted much attention during recent years and therefore many published experimental and theoretical results are available. These results will be applied to reaction zones with steep temperature gradients. In Fig. 7 the dimensionless P&let group (Pe’ = u,,d/ Da& is plotted as a function of the Reynolds number from a paper of Gunn[8]. The plot exhibits a slight maximum in the range 2 < Re < 20, where the P&let number is three. Since this range of Reynolds numbers agrees with the number occurring in our fixed bed simulation, it seems appropriate to evaluate D,,, from Pe’ = 3. However in order to show the generalised effect of axial mass diffusion, calculations were performed for Pe’ = 1,2,3 and Pe’ = 03.Numerical results for the last case with vanishing difiksion have already been presented in Fig. 4 and it has been pointed out that for this situation agreement exists between computation and 4.4
3
&x=7x Xo+O~i3u,pGcpGd
IL____
0.25
ug [m/secl
2
----
0
1492
Moving reaction zones
-I
&x43
0
msec-’ of Fig. 8 it is noted that increasing diffusive fluxes cause a decrease in the velocity W.
This is what one would expect because diffusion means a spreading of the concentration 2 profile and particularly a lowering of the concl.1 centration in the ignition cross section. Then the temperature gradient also becomes smaller which means smaller heat fluxes and smaller moving velocities w. 0.1 IO 100 1000 Surprisingly, however, at u. = 0.08 msec-l Fig. 7. Pe’ as a functionzf Reynolds number Re [8]. there is an intersection of the various curves, and beyond this intersection just the opposite experiment concerning the exponents of u,, and effects take place. Various tests have verified cc* this intersection so that numerical uncertainties Figure 8 presents the results for all four can be excluded. P&let numbers. Distinct differences are noted The explanation is difficult, but may be as particularly in the shape of the w =f(u,J curves, follows: The above argument for the decrease while the absolute values for w differ only by a of w with increasing diffusion are based on the few percent. Since the shape of the Pe’ -+ cc assumption that the ignition conditions remain curve is in agreement with experimental results, unchanged if D,,, changes. Since however it can be said that increasing diffusional prolarger D,,,, values lower the concentration in cesses distort the curve. the ignition cross section, there would well be an If we first consider the range O-02 < u0 6 0.08 effect on the ignition temperature which may be shifted to higher values. While this effect Y,_~O’Ol r-3 is small at low velocities with low diffusion k,=91000[kmolh processes, it may be more pronounced in the d -0.004 [m] range of higher gas velocities with higher diffuTd-0 [Kl sion processes. Then a range may occur where the temperature shift may be responsible for higher reaction rates near the ignition zone although the concentration is lowered in the ignition section. This then leads to an increase in w. However the previous remarks may be discarded as speculative if it is agreed that Pe’ = 3 is the right value to be used for the calculations in the above mentioned range of Reynolds numbers, since a comparison of the Pe’ = 3 curve with that of Pe’ + 03 shows only very slight differences which are within the limits of experimental accuracy. It may then be concluded that for the evaluation of the transit I I I I velocity of reaction zones axial mass diffusion 0.16 0.12 0.06 0.04 processes may be neglected, and that the w = u. Gn/secl f(u,,) is mainly determined by axial heat conducFig. 8. Reaction zone velocities w for diierent Pklet tion processes as assumed in [ 1,2]. numbers. I
o.o,l
1493
D. VORTMEYER 4.5
InJluence of particle diameter With increasing particle diameters the axial heat conductivity and the radiation contribution are increased. Consequently, the velocity of the reaction zone should increase if the particle diameters become larger, while all other parameters are kept constant. This is confirmed by the computed results of Fig. 9 which were calculated without heat radiation contributions for different diameters. An increase of d means larger transit velocities. For the inlet velocity u,, = 0.03 m/set the transit velocity w is presented as a function of d in Fig. 9, which also contains the calculations without heat-radiation contributions for the above value of uo. In order to obtain a functional relation between w and din the following form w
=
a’ruo@77cGO~5&
-
/juo
(34)
an exponent e = 0.069 was found for the curve neglecting radiation, while for the curve with
and W. JAHNEL
radiation the exponent (0.067 < e < O-075).
is no longer constant.
4.6 In..uence of reaction parameters Since the movement of the reaction is mainly governed by axial heat transport, there is no reason to believe that reaction order or frequency factor should have any effect on the general shape of the curves as long as the heat capacity of the bed remains constant. In fact the calculations performed in this work are done with the reaction kinetics of the ethane oxidation reaction and with fixed bed data from [l] for activated coal. No difference was noticed between the computed exponents and the measured ones. Naturally this is also valid for the frequency factor, which is changed in the computations of Fig. 10. The shape of the curves remains similar, only the absolute values of w being affected. Larger k. values imply an increase of w, since the temperature gradients become larger with increasing reaction velocity. 4.7 Heat capacity of the solid Although the exponents of u. and cG remain the same as before, the shape of the w =f(u,)
Do;0 [m2/sec]
y =0.009 d =0.004
[-I [m]
No rodiotion
-l -
-2
0
Fig. 9. Reaction
I 0.04
I 0.06 u0
I 0.12
Wsecl
/ 0.16
0,
-21
I
0
zone velocities for different particle diameters.
1494
0.05
I
1
0.10 u,
Fig.
10. w
0.15
,
0.20
0. 5
h/secl
as a function of u0 for three frequency factors.
Moving reaction zones
curves is influenced by the heat capacity of the solid as indicated in Fig. 11. This figure contains two calculated curves which differ in the heat capacities of the fixed bed only. As was pointed out in the introduction, changing heat capacities only effect the factor in formula (1). This again is in agreement with experiments.
--71_1
30
rmecl Fig. 11. w as a function of u0 for two heat capacities. 102xu,
5. CONCLUSIONS
After having developed a method for the numerical solution of the time dependent system of Eqs. (3)-(7) the dynamic behaviour of a fixed bed reactor with an exothermic reaction was studied with regard to moving reaction zones. Such an investigation has many advantages since it permits variation of parameters which are nearly impossible to change experimentally. The results concerning the influence of these parameters on the transit velocity of the reaction zone are as follows: (1) If the computed form
and if the axial effective taken to be equal to
heat conduction
was
the exponents p = 0.77 and v = O-5 are in complete agreement with experiment. (2) Heat radiation effects do not change the exponents, although the transit velocity w is increased by a few percent due to the larger heat flux. (3) Calculations with L& = const. resulted in a distorted curve with an exponent of p = 0.73 which is not in agreement with experimental evidence. (4) Results were obtained for the P&let number Pe’ = 1, 2, 3 and Pe’ + ~0.While the curves for Pe’ = 1 and 2 differ considerably from the shape of the experimental curves, Pe’ = 3 gives good agreement in the exponents. Pe’ = 3 is the P&let group recommended by Gunn for that particular range of Reynolds numbers. The P&let groups measured and derived under isothermal conditions seem to be applicable to non-isothermal cases. (5) Since there are only very small differences between the results for Pe’ = 3 and Pe’ + w (D,, = 0), it may be concluded in agreement with earlier work[ 1,2] that axial mass diffusion processes are negligible in comparison with axial heat conduction. (6) Because conduction and radiation depend on d, a slight influence of the particle diameter on the moving velocity is noted. (7) Different reaction kinetics do not change the exponents for u. and cc. (8) The above results are a further indication that the homogeneous mathematical model for a fixed bed catalytic reactor gives satisfactory results.
results are plotted in the
NOTATION
; B CGO
1495
coefficient of absorption, m-l coefficient of emission, m-l radiation transmission number, inlet concentration at NTP, kmol mm3
D. VORTMEYER
and W. JAHNEL
specific heat of solid, J kg-’ grd-l specific heat of gas, J kg-’ grd-’ cog radiation constant of the black cs body, Wm-l grd-’ d particle dia., m RU-0 effective axial diffusion coefficient, at NTP, m2 set-l coefficient of reflection, m-l reaction enthalpy, J. kmol-’ Z partial energy flux of radiation, Wmm2 K partial energy flux of radiation,
T TG z.4 u,, w
temperature, K gas inlet temperature, K gas velocity, msec-’ gas velocity at NTP, msec-’ velocity of the reaction zone, msec-’ x coordinate, m y mole fraction yGo mole fraction of reactant gas at NTP
C8
*Hg
Wm-2
2 Qs=Z-K
;t
Greeksymbols
--
constants in formula (1)
frequency factor, kmol mm3see-l specific heat flux by conduction, Wme2 specific heat flux by radiation, Wme2 reaction rate, kmol m3 set-’ time, set
PC
PS
emissivity of particle surfaces void fraction effective axial conductivity without radiative effects, Wm-’ grd-’ density of gas, kg rnw3 effective density of solid, kg me3
REFERENCES [l] WICKE E. and VORTMEYER D., Ber. Bunsenges. physik. Chemie 1959 63 145. 121 VORTMEYER D., Ber. Bunsenges. physik. Chemie 196165 282. [3] VORTMEYER D. andJAHNEL W., Chem. Znger. Tech. (CZT) 197143461. [4] VORTMEYER D., Ber. Bunsenges. physik. Chemie 1970 74 127; Fortschritt-Ber. VDI-Z. Reihe 3 Nr. 9, Dusseldorf 1966; Chem. Znger. Tech. 1966 38 404. [5] VORTMEYER D. and UORNER C. G., Chem. Znger. Tech. 1966 38 1077. [6] COLLATZ L., The Numerical Treatment ofDigerentiuZEquations. Springer-Verlag, Berlin 1966. 171 WICKE E., PADBERG G. and ARENS H., 4th European Symp. on React. Engng. 6th Session Briissel 1968. [81 GUNN D. J., Trans. Znstr. Chem. Engrs 1969 47 T3.56. [9] HARTIG H. and WICKE E., Z. phys. Chem. Neue Folge 1963 38 265. [lo] PADBERG G. and WICKE E., Chem. EngngSci. 1967 22 1035. R&sum&Les auteurs obtiennent des zones mobiles de reaction d’aprbs Ia solution dun groupe d’equations di@rentielles interdtuendantes et B l’etat instable de I’bneraie de I%auilibre de mat&e dans le cas dun reacteur a lit cataIytique fixe ayant une reaction exothermyque rapid;. Les auteurs font l’etude parametrique de l’inlhtence de la vitesse du gaz et de Ia concentration du corps reagissant, des diffusivites de chaleur et de masse, de la radiation de chaleur, des parametres cinetiques et du diambtre des particules. Les resultats experimentaux concordent tout a fait avec ces resultats. Zusammenfassung- Es werden bewegliche Reaktionszonen erhalten durch Losung des Satzes voneinander abhangiger und auf den nicht stationaren Zustand bezilglicher Energie- und Stoffbilanz Differential-gleichungen fur einen katalytischen Festbettreaktor mit einer schnellen, exothermen Reaktion. Es wurden parametrische Studien in bezug auf den Einfluss der Gasgeschwindigkeit und Konzentration der Reaktionsteilnehmer, der Warme und Stoffdilfusionsvermiigen, der W%rmestrahlung, der kinetisthen Parameter und der Teilchendurchmesser durchgefiihrt. Ein Vergleich mit experimentellen Ergebnissen envies sich als sehr befriedigend.
1496
Pergamon Press.
Printed in Great Britain
Moving reaction zones in fixed bed reactors under the influence of various parameters D. VORTMEYER and W. JAHNELt Technische Universit;it Miinchen, Germany (Received 28 September 197 1) Abstract-Moving reaction zones are obtained by the solution of the interdependent and unsteady state set of the energy and material balance differential equations for a fixed bed catalytic reactor with a fast exothermic reaction. Parametric studies concerning the influence of gas velocity and reactant concentration, of heat and mass diffusivities, of heat radiation, of the kinetic parameters and of the particle diameter were performed. A comparison with experimental work turned out to be very satisfactory.
IT IS established that in the case of fast exothermic reactions in fixed bed reactors the total conversion takes place over a few particle diameters range known as the reaction zone. Experimental investigations by Wicke and Vortmeyer [l] have shown that the reaction zone may move through the turbular reactor in either axial direction (Fig. 1). It was further demonstrated that direction and velocity w of the zone depend heavily on the gas velocity u. and on the reactant concentration cGo. The experimental relationship is:
where CYand p depend on other parameters of the system. Positive values of w denote movement against the gas flow while negative values indicate that the reaction zone is being blown out of the reactor. The exponents 0.77 of u,, and 0.5 of cc,, were later confirmed by further measurements by Wicke et al. [9, lo]. The movement against the flow is the most interesting feature and can only be explained by axial diffusion processes of heat and mass. Since it may safely be assumed that the transport of energy is the dominating process, a qualitative treatment of the energy equation in [l]
and an extended one in [2] by analogy to laminar flames gave very satisfactory results for the exponents of u. and cco. In addition, it led to interpretation of the coefficients OLand p in formula (1). ff is a coefficient of great complexity while /? turns out as the heat capacity ratio of fluid to solid. However many questions still remained open concerning the influence of various other parameters such as radiation, mass diffusion, kinetics, particle diameter etc. Since it is very difficult to change some of these parameters experimentally it seemed worthwhile to solve the complete set of the unsteady state differential equations numerically and thus simulate the dynamics of the system on the computer. Results by Vortmeyer and Jahnel[3] for the energy equation alone turned out to be very encouraging, since the exponents 0.77 of u. and 0.5 of cGo from the Computer solutions were in perfect agreement with formula (1). 1. THE EQUATIONS OF THE HOMOGENEOUS MODEL FOR A FIXED BED ADIABATIC REACTOR
Reactors with fast exothermic reactions exhibit substantial differences in temperature between solid and gas, particularly in the range of the reaction zone where the reaction rate is
tPresent address: Motoren-Turbinen-Union, Miinchen, Germany.
1485 CESVol27Na8-A
D. VORTMEYER
and W. JAHNEL
= a(g,aT/ax)/ax-
a(upGcpGT)/ax
.
- a(Z - K)/ax - AH . i.
T,=473
(2)
Although not necessary from the computational point of view, the average values A$, and c were introduced in the energy equation and the heat capacity of the gas was neglected in Eq. (2).
[K]
Energy balance: P,C~aTiat = ~&a2Tlax2 --uopco~pc aTlax - a( 1 - K)/ax - AH?. -5
-
-6
-
-7
1 5
0
1 IO
15
10%
0.
I 20
I 25
30
[m/secl
Fig. 1. Reaction zone velocity as a function of u0 and cGO,
very high. Nevertheless, the one dimensional homogeneous model was applied. Since the movement of the reaction zone depends on heat flux from the reactant to the non-reactant part of the bed, the application of the method appears to be appropriate in this range of relatively low temperatures with only minor temperature differences between solid and gas. The following additional assumptions were made: (i) A quasi homogeneous reaction is assumed to take place with no change in the number of moles. (ii) The axial conduction and diffusion processes under unsteady conditions may be evaluated from the Fourier and Fick equations. (iii) The radiative flux is approximated by the two-flux model as developed in [4] for fixed beds. The equations obtained were: Energy balance: (P,c,+
Although confirmed only for isothermal conditions, the effective axial diffusion coefficient was evaluated from Pe’ = udlD& = const. and was applied to the non isothermal reaction zone. Then the following material balance is obtained. E,aylat = D& a2ylax2 - u,,ay/ax - i
(4)
with (3 The transport equations for radiation [4] are: dZ/d_r+ al = bC,( TI 100)4 + gK
(6)
dZC/d_r- UK = bC,(T/ 1OO)4- gl.
(7)
The coefficients a, b, g depend on the emissivity and on the radiative transmission number of the bed. a=(1+B)2+(1+B)(1-B)(1-+J1-B)2 (1 +B)2-(1 -B)2(1-E)2
(l+B)d (8)
b= (1+B)2-(1+B)(1-B)(1-e)~E,(1-B)2 (1+B)2-(1-B)2(1-e)2
(l+B)d /a\ \7/
g=
lphd
(3)
am
(1+BY+(1+B)(1-B) (1+B)2(l-B)*(l-_E)*
(I_E)
(0;;;;. (10)
1486
Moving reactionzones
E is the emissivity of the particle surface and B is a dimensionless radiation transmission number for a fixed bed [5]. 2. THE
BOUNDARY
CONDITIONS
In order to calculate the dynamic behaviour of the tubular reactor and particularly the transit velocity of the reaction zone, boundary conditions have to be fixed and must be chosen in such a way that they are in agreement with the experimental system. The experiments [ l] were performed in such a way that the transit velocity of the reaction zone was measured at a great distance from the inlet and outlet of the bed. Mathematically this means an infinitely long reactor with the reactions taking place somewhere in the middle. However, only the finite reactor was treated. Consequently boundary conditions had to be introduced at the outlet and the second derivatives were chosen for this purpose. The conditions for the radiative transport equation at inlet and outlet are that there are no net fluxes of radiation. The following system of boundary conditions arises: t=O
--m
T,s
T(x)<
TE
YGO5 y(x) 2 0 (11) tz
0
x=-m
y=yGO T= TG I-K=0
t 2 0
x=+m
(12)
a2yla9=o a2TlaX2 = 0 I-K=O.
(13)
For ignition to take place, temperature and concentration profiles T(x) and y(x) have to be present somewhere in the reactor between plus and minus infinity at time level zero. At all other times the only gas stream conditions we know of are the inlet concentration and inlet temperature when the reactor zone is far away. Since the
reactor system in the infinitely long reactor is an initial value problem, the conditions at x = + m would come out by computation as aTlax = aylax = 0 or aTlaX
= a2yla2 = 0
for the adiabatic system. However for computational est condition
reasons the weak-
(14)
a2Tla.e = a2ylax2 = 0 was used at the now finite end of the reactor. Preliminary calculations for various reactor lengths have shown that in the case of particle diameters d = 0404 m, a length of O-1 m was sufficient to give results which were almost in exact agreement with those obtained when using very long reactors. Therefore all calculations were performed with this shorter reactor length. The influence of the boundary conditions at x = +m was also tested, with no change in the results. 3. THE
NUMERICAL
METHOD
The physical problem is described by a nonlinear system of two partial differential equations of the parabolic type (3), (4) and two ordinary first order differential Eqs. (6), (7). These four interdependent differential equations require two different types of solutions. Several methods for solving non-linear equations that result in linear finite difference equations, which can be solved by existing algorithms and which do not involve excessive iteration, have been developed. In this paper only forward Taylor Series projections are used. These methods do not need any iterations to treat the nonlinearity. Only in the case of ordinary differential equations does an iteration have to be made in order to satisfy suitable boundary conditions. The previously mentioned system of differential equations causes a combined initial-boundary-value problem. The well-known Runge-Kutta method was used to solve the ordinary differential equations
1487
D. VORTMEYER
which describe the transport mechanism of radiation. For this purpose it is convenient to write the equations in the following form: d+/dx = q($, x) + PVJ.
(1%
In this case I,IJrepresents the dependent variable while x is the independent variable. The nonlinear term q($,x) contains the coupling expression and p is a constant (see Eqs. 16, 17). If the radiation flux from x = - 00 to x = + CQis considered, I/J= Z(x) and p = --a
q=
bC,(T(~)/10O)~+gK(x).
Otherwise $ = K(x), the radiation opposite direction, and
and W. JAHNEL
After solving the differential equations (6) and (7) by the Runge-Kutta method we obtain a new distribution for Z(x) and K(x). The boundary conditions Z--K = 0 are now checked. If the condition is fulfilled and two successive distributions differ by less than a predicted error, the iteration is stopped and the correct Z(x) and K(x) distribution values are found. If the condition is not fulfilled the obtained distribution values are used to continue the iteration process. The radiation flux appears in the energy balance equation in the form
(16)
e(Z(x) - K(x))lax = 2X,( T(x)/ 100)4 - (g - a)(Z(x) - K(x)) (18)
flux in the
which can easily be determined by subtraction of Eq. (7) from Eq. (6). After having found the p = a q = - bC,(T(x)l loo)4 - gZ(x). ( 17) radiational heat flux distribution and the reaction rate, the parabolic equations for mass and The initial conditions for x = -CO and x = + 03 energy balance can be solved by an explicit are Z-K = 0 in both cases. forward method for the next time step using a Let us consider an arbitrary time step t at Taylor Series projection. This method has the which the temperature and concentration distriadvantage that no iteration is necessary because bution over the reactor length are known. Now the non-linearity only corresponds to the points the energy flux and the reaction rate which in the previous time step. The possibility exists depend on the distributions can be calculated that numerical instability arises when using this for this time step. method and therefore a stability analysis is The calculation of reaction rate Eq. (5) is easy necessary. to perform whereas the computation of radiaFor solving the partial differential equations it tional heat fluxes is much more difficult. In this is more convenient to write the equation in a case two interdependent ordinary differential general form such as the following: equations are to be solved. Because of the expression bC, ( T(x) I 100) 4, which contains the W) = a$ + P4t + Y& + 84, = A$& -&0. ( 19) interdependent variables, an iteration becomes necessary. The coefficients (Y,/3, y, 6 and f are determined In order to begin the iteration, Z(x) and K(x) by the special differential equations. In case of have to be chosen. For this purpose we use a energy balance, the dependent variable 4 is linear distribution of the form equal to the temperature T and the coefficients Z(x) = K(x) are in which the boundary values are set as Z(x = 0) = C,( T(x = O)/ 100)4 and
a=0 f=
ZW = x~) = C,(T(x = xiv)/ 1OO)4.
p=15,&.
y=l.lOpG&~
- S(Z(x) - K(x))bx
The mass balance 1488
- AHi.
equation
6=--h&
(20) is treated in the
Moving reaction zones
same way. In this case we set C#J = y and the coefficients are a=0
p=1
y=uo
a,,=-2&+6
(25)
a=-p&.o
(26)
(21)
f=-;.
For the energy equation is:
The boundary conditions can be written as
For
balance the finite difference
2
and
lsj
T IJtl=-H~*J(~)+Ti~(l-2~~)
It will be remembered that the reaction rate, the nonlinear term of the equations, contains the coupling elements as well. Because of the properties of the forward difference method, iteration is unnecessary. Both of the partial differential equations are to be solved independently from each other. After transforming the partial differential equation into a finite difference equation by using a forward Taylor Series projection over three discrete points shown in Fig. 2 we obtain the following equation:
4 lJ+l
=
+
--1 (2
Ti-IJ
uOpG%G
&Es
;
he,,
At
p& M > (27)
+ T~+M
For the mass balance the finite difference equationfor G i s (IV- l)and 1 sj < ais: yrj=-ridAt+yid
>
(fi,j-Q.A.i-“i+l,.i+~+lJ -
Q-lki-1Jh+1
(22)
using:
+ Yl+l.r %+1
=
a
a-1J
-
--M
P -At 6
‘i-l
‘i
)
- (28)
(23) 1.x 2 b
(24)
Using these two finite difference equations we obtain a new temperature and concentration distribution for the time level ttAl. These finite difference equations have the disadvantage of becoming instable under certain conditions. The numerical stability depends on the ratio of At/Ax2 and the convergence depends on the absolute values of At and Ax. The stability conditions[6] for the equations described above are:
and 'i-2
At
Ax
‘i*l
‘i+2
X
Fig. 2. The rectangular finite-difference mesh.
1489
D. VORTMEYER
and each term should be approximately of the same size. The values chosen by the authors were: 3 Q-lJ _ _. aij+l
afd
12’
-4.
a&j+1
&+I,.4 _
--I
12’ Uij+l
12
(29)
and stability was obtained in each example. The increments derived from the stated stability conditions are: For energy balance: A-&_=
A&
and
AtT=p.
2u@Cc,G
(30) az
For mass balance:
A&J!& 2ull
and
Ax”, At, = 3%;
(31)
Having calculated the new T(x) and y(x) distributions for the time step tij+l the values for the reaction rate r(x) and the radiation fluxes Z(x) and K(x) have to be computed with the algorithm described previously. The old values at the time tidy are now replaced by the new values calculated at the time tij+l and the whole process is repeated. In this way the time dependent properties of the reactor are obtained. 4. RESULTS
AND
DISCUSSION
4.1 Transient behaviour Having ignited the reactor it took about 300-400 set real time to establish a uniformly moving reaction zone. The above mentioned time naturally depends on the chosen ignition profile which was a sine curve. The transient behaviour is shown in Fig. 3 for two cases, which differ in the end temperatures of the ignition profiles. While in case (a) the end temperature at time level zero is too low and therefore the fixed bed has to be heated up in order to obtain the steady moving reaction zone, case (b) presents the opposite situation. The computations were done with the following data: gas inlet temperature
. . . TG = 473 K
and W. JAHNEL
effective solid density . . . pS = 450 kg/m3 gas density at inlet temperature. . . . . . . . . pG = 1.29 kg/m3 specific heat of solid. . . . . c, = 1423.6 J/kg grd specific heat of gas . . . . . c, = 1071.9 J/kg grd particle diameter. . . . . . . . d = 0.004 m emessivity of solid . . . . . . .E = 0.8 gas heat conductivity . . haso = O-0514 W/m grd reaction enthalpy . . . . . AH = - 1.4 X lo9 J/kmol activation energy. . . . . . . . E = 96.3 X 106J/kmol frequency factor. . . . . . . . . k,, = 9100 kmol/m3 set reaction order . . . . . . . . . . n = 0.35 radiation transmission number . . . . . . . . . . . . . . B==0*2-. The kinetic data are similar to those which were measured for the ethane oxidation on a Pd/,4l,O, catalyst by Wicke et al. [71. 4.2 Axial
dijksion processes and moving reaction zone In former papers[ 1,2] the movement of a reaction zone against the inflowing gas was explained by the axial heat conduction. It also was pointed out, that there are certain analogies to the theory of laminar flames. The theoretical treatment however was restricted to the energy equation only. Radiative transport processes and in particular axial mass diffusion processes were neglected. Because of the considerable mathematical complexity of the fully stated problem (Eqs. 3-7) an analytical solution seems to be impossible. Therefore the influence of the various parameters was investigated numerically. 4.3 The influence of axial heat conduction and radiation In a first paper[3] on a numerical treatment the energy equation was solved separately by using a model reaction rate function dependent only on temperature. When the resulting calculated velocity w of the reaction zone was plotted against the gas velocity u0 a complete qualitative agreement with experimental results (Eq. 1) was obtained. The same agreement was
1490
Moving reaction zones
d = 0.004
[m]
No radiation
65C
0.02 0.03 0.04 005 0.02 0.03 0.04
0.05 0.06 0.07 0.09 0.09 C
x Fig. 3. Time dependent
[ml
A& = haI,, + 0.8 u,,pcoGd.
(32)
As is shown in the following figures, these results remain unchanged if the complete set of simultaneous Eqs. (3)-(7) is solved. Figure 4 shows w as a function of the inlet velocity u. with the concentration parameter YOO. Because of reduced computer time the diffusion coefficient in the mass balance Eq. (4) was set to zero in these calculations since one of the following sections contains a detailed treatment of the influence of this transport property. If the results of Fig. 4 are plotted in double logarithmic form of Fig. 5 In (uo+i
*
temperature profile after ignition. (a) End temperature (b) End temperature of ignition profile too high.
observed when the inlet concentration cm was varied. In all cases the experimentally determined exponents 0.77 and 0.5 for u. and coo were confirmed when the following expression was used for the axial effective heat conductivity.
w) =f(ln
uo)
(33)
0.06 007
008
0.
[ml
of ignition profile too low
the gradients and the intercepts of the straight lines again lead to the exponents 0.77 for u. and 05 for cG. This result is valid for calculations with and without heat radiation. That radiative effects do not change the shape of the curves may be explained in the following way: The energy transport responsible for the moving reaction is the heat flux from the reacting to the non reacting part of the catalyst bed. The temperature of this cross section is about 500 K. Since this temperature is still rather low for larger radiation effects heat conduction is dominant. Our data for radiation effects led to the ratio QL/Qs = 3 if Q, and Q, are conductive and radiative heat fluxes. Therefore the radiative contribution in this cross section means an additive term in the Eq. (32) for A& and does not change its dependence on u. as a first approximation. That the functional relation w =f(uo) of Eq. (1) is important for the exponents can be seen from Fig. 6. This contains calculated results with a constant value of 0.538 W/m grd for z Al-
1491
D. VORTMEYER
and W JAHNEL 4 D,.= 0 [m */set] k, = 91000
4 t
-
No rodiotion
---
Rodiotion
[k mol /m3 set]
J--
;
p2ti $ -%
’
1
s
4%
O-
h,\
‘,
‘\
%R \
0006
\
\
-
‘\
’\
\
\
\
\
-I
0.05
\
\
0.007
0.006
\
\
\
\
0
0.009
\
\
\
\
0.15
0.10
0.20
0
Fig. 4. Calculated values for w as a function of u0 and cc0 with and without radiative contributions.
I
I
45
I
,,I,,
IO
1
I
15
20
"0
Fig. 5. Double logarithmic plot.
though no dramatic changes appear, there is an obvious difference in the exponent for u0 which now is 0.73 instead of the observed value 0.77.
0.04
0.08
0.12
u,
Fig. 6. w =f(uo)
[m/s4
0.16
for a constant and u,-dependent
-I
0 .2 0
value of
The injluence of axial mass d$usion Axial mass diffusion processes in isothermal fixed beds have attracted much attention during recent years and therefore many published experimental and theoretical results are available. These results will be applied to reaction zones with steep temperature gradients. In Fig. 7 the dimensionless P&let group (Pe’ = u,,d/ Da& is plotted as a function of the Reynolds number from a paper of Gunn[8]. The plot exhibits a slight maximum in the range 2 < Re < 20, where the P&let number is three. Since this range of Reynolds numbers agrees with the number occurring in our fixed bed simulation, it seems appropriate to evaluate D,,, from Pe’ = 3. However in order to show the generalised effect of axial mass diffusion, calculations were performed for Pe’ = 1,2,3 and Pe’ = 03.Numerical results for the last case with vanishing difiksion have already been presented in Fig. 4 and it has been pointed out that for this situation agreement exists between computation and 4.4
3
&x=7x Xo+O~i3u,pGcpGd
IL____
0.25
ug [m/secl
2
----
0
1492
Moving reaction zones
-I
&x43
0
msec-’ of Fig. 8 it is noted that increasing diffusive fluxes cause a decrease in the velocity W.
This is what one would expect because diffusion means a spreading of the concentration 2 profile and particularly a lowering of the concl.1 centration in the ignition cross section. Then the temperature gradient also becomes smaller which means smaller heat fluxes and smaller moving velocities w. 0.1 IO 100 1000 Surprisingly, however, at u. = 0.08 msec-l Fig. 7. Pe’ as a functionzf Reynolds number Re [8]. there is an intersection of the various curves, and beyond this intersection just the opposite experiment concerning the exponents of u,, and effects take place. Various tests have verified cc* this intersection so that numerical uncertainties Figure 8 presents the results for all four can be excluded. P&let numbers. Distinct differences are noted The explanation is difficult, but may be as particularly in the shape of the w =f(u,J curves, follows: The above argument for the decrease while the absolute values for w differ only by a of w with increasing diffusion are based on the few percent. Since the shape of the Pe’ -+ cc assumption that the ignition conditions remain curve is in agreement with experimental results, unchanged if D,,, changes. Since however it can be said that increasing diffusional prolarger D,,,, values lower the concentration in cesses distort the curve. the ignition cross section, there would well be an If we first consider the range O-02 < u0 6 0.08 effect on the ignition temperature which may be shifted to higher values. While this effect Y,_~O’Ol r-3 is small at low velocities with low diffusion k,=91000[kmolh processes, it may be more pronounced in the d -0.004 [m] range of higher gas velocities with higher diffuTd-0 [Kl sion processes. Then a range may occur where the temperature shift may be responsible for higher reaction rates near the ignition zone although the concentration is lowered in the ignition section. This then leads to an increase in w. However the previous remarks may be discarded as speculative if it is agreed that Pe’ = 3 is the right value to be used for the calculations in the above mentioned range of Reynolds numbers, since a comparison of the Pe’ = 3 curve with that of Pe’ + 03 shows only very slight differences which are within the limits of experimental accuracy. It may then be concluded that for the evaluation of the transit I I I I velocity of reaction zones axial mass diffusion 0.16 0.12 0.06 0.04 processes may be neglected, and that the w = u. Gn/secl f(u,,) is mainly determined by axial heat conducFig. 8. Reaction zone velocities w for diierent Pklet tion processes as assumed in [ 1,2]. numbers. I
o.o,l
1493
D. VORTMEYER 4.5
InJluence of particle diameter With increasing particle diameters the axial heat conductivity and the radiation contribution are increased. Consequently, the velocity of the reaction zone should increase if the particle diameters become larger, while all other parameters are kept constant. This is confirmed by the computed results of Fig. 9 which were calculated without heat radiation contributions for different diameters. An increase of d means larger transit velocities. For the inlet velocity u,, = 0.03 m/set the transit velocity w is presented as a function of d in Fig. 9, which also contains the calculations without heat-radiation contributions for the above value of uo. In order to obtain a functional relation between w and din the following form w
=
a’ruo@77cGO~5&
-
/juo
(34)
an exponent e = 0.069 was found for the curve neglecting radiation, while for the curve with
and W. JAHNEL
radiation the exponent (0.067 < e < O-075).
is no longer constant.
4.6 In..uence of reaction parameters Since the movement of the reaction is mainly governed by axial heat transport, there is no reason to believe that reaction order or frequency factor should have any effect on the general shape of the curves as long as the heat capacity of the bed remains constant. In fact the calculations performed in this work are done with the reaction kinetics of the ethane oxidation reaction and with fixed bed data from [l] for activated coal. No difference was noticed between the computed exponents and the measured ones. Naturally this is also valid for the frequency factor, which is changed in the computations of Fig. 10. The shape of the curves remains similar, only the absolute values of w being affected. Larger k. values imply an increase of w, since the temperature gradients become larger with increasing reaction velocity. 4.7 Heat capacity of the solid Although the exponents of u. and cG remain the same as before, the shape of the w =f(u,)
Do;0 [m2/sec]
y =0.009 d =0.004
[-I [m]
No rodiotion
-l -
-2
0
Fig. 9. Reaction
I 0.04
I 0.06 u0
I 0.12
Wsecl
/ 0.16
0,
-21
I
0
zone velocities for different particle diameters.
1494
0.05
I
1
0.10 u,
Fig.
10. w
0.15
,
0.20
0. 5
h/secl
as a function of u0 for three frequency factors.
Moving reaction zones
curves is influenced by the heat capacity of the solid as indicated in Fig. 11. This figure contains two calculated curves which differ in the heat capacities of the fixed bed only. As was pointed out in the introduction, changing heat capacities only effect the factor in formula (1). This again is in agreement with experiments.
--71_1
30
rmecl Fig. 11. w as a function of u0 for two heat capacities. 102xu,
5. CONCLUSIONS
After having developed a method for the numerical solution of the time dependent system of Eqs. (3)-(7) the dynamic behaviour of a fixed bed reactor with an exothermic reaction was studied with regard to moving reaction zones. Such an investigation has many advantages since it permits variation of parameters which are nearly impossible to change experimentally. The results concerning the influence of these parameters on the transit velocity of the reaction zone are as follows: (1) If the computed form
and if the axial effective taken to be equal to
heat conduction
was
the exponents p = 0.77 and v = O-5 are in complete agreement with experiment. (2) Heat radiation effects do not change the exponents, although the transit velocity w is increased by a few percent due to the larger heat flux. (3) Calculations with L& = const. resulted in a distorted curve with an exponent of p = 0.73 which is not in agreement with experimental evidence. (4) Results were obtained for the P&let number Pe’ = 1, 2, 3 and Pe’ + ~0.While the curves for Pe’ = 1 and 2 differ considerably from the shape of the experimental curves, Pe’ = 3 gives good agreement in the exponents. Pe’ = 3 is the P&let group recommended by Gunn for that particular range of Reynolds numbers. The P&let groups measured and derived under isothermal conditions seem to be applicable to non-isothermal cases. (5) Since there are only very small differences between the results for Pe’ = 3 and Pe’ + w (D,, = 0), it may be concluded in agreement with earlier work[ 1,2] that axial mass diffusion processes are negligible in comparison with axial heat conduction. (6) Because conduction and radiation depend on d, a slight influence of the particle diameter on the moving velocity is noted. (7) Different reaction kinetics do not change the exponents for u. and cc. (8) The above results are a further indication that the homogeneous mathematical model for a fixed bed catalytic reactor gives satisfactory results.
results are plotted in the
NOTATION
; B CGO
1495
coefficient of absorption, m-l coefficient of emission, m-l radiation transmission number, inlet concentration at NTP, kmol mm3
D. VORTMEYER
and W. JAHNEL
specific heat of solid, J kg-’ grd-l specific heat of gas, J kg-’ grd-’ cog radiation constant of the black cs body, Wm-l grd-’ d particle dia., m RU-0 effective axial diffusion coefficient, at NTP, m2 set-l coefficient of reflection, m-l reaction enthalpy, J. kmol-’ Z partial energy flux of radiation, Wmm2 K partial energy flux of radiation,
T TG z.4 u,, w
temperature, K gas inlet temperature, K gas velocity, msec-’ gas velocity at NTP, msec-’ velocity of the reaction zone, msec-’ x coordinate, m y mole fraction yGo mole fraction of reactant gas at NTP
C8
*Hg
Wm-2
2 Qs=Z-K
;t
Greeksymbols
--
constants in formula (1)
frequency factor, kmol mm3see-l specific heat flux by conduction, Wme2 specific heat flux by radiation, Wme2 reaction rate, kmol m3 set-’ time, set
PC
PS
emissivity of particle surfaces void fraction effective axial conductivity without radiative effects, Wm-’ grd-’ density of gas, kg rnw3 effective density of solid, kg me3
REFERENCES [l] WICKE E. and VORTMEYER D., Ber. Bunsenges. physik. Chemie 1959 63 145. 121 VORTMEYER D., Ber. Bunsenges. physik. Chemie 196165 282. [3] VORTMEYER D. andJAHNEL W., Chem. Znger. Tech. (CZT) 197143461. [4] VORTMEYER D., Ber. Bunsenges. physik. Chemie 1970 74 127; Fortschritt-Ber. VDI-Z. Reihe 3 Nr. 9, Dusseldorf 1966; Chem. Znger. Tech. 1966 38 404. [5] VORTMEYER D. and UORNER C. G., Chem. Znger. Tech. 1966 38 1077. [6] COLLATZ L., The Numerical Treatment ofDigerentiuZEquations. Springer-Verlag, Berlin 1966. 171 WICKE E., PADBERG G. and ARENS H., 4th European Symp. on React. Engng. 6th Session Briissel 1968. [81 GUNN D. J., Trans. Znstr. Chem. Engrs 1969 47 T3.56. [9] HARTIG H. and WICKE E., Z. phys. Chem. Neue Folge 1963 38 265. [lo] PADBERG G. and WICKE E., Chem. EngngSci. 1967 22 1035. R&sum&Les auteurs obtiennent des zones mobiles de reaction d’aprbs Ia solution dun groupe d’equations di@rentielles interdtuendantes et B l’etat instable de I’bneraie de I%auilibre de mat&e dans le cas dun reacteur a lit cataIytique fixe ayant une reaction exothermyque rapid;. Les auteurs font l’etude parametrique de l’inlhtence de la vitesse du gaz et de Ia concentration du corps reagissant, des diffusivites de chaleur et de masse, de la radiation de chaleur, des parametres cinetiques et du diambtre des particules. Les resultats experimentaux concordent tout a fait avec ces resultats. Zusammenfassung- Es werden bewegliche Reaktionszonen erhalten durch Losung des Satzes voneinander abhangiger und auf den nicht stationaren Zustand bezilglicher Energie- und Stoffbilanz Differential-gleichungen fur einen katalytischen Festbettreaktor mit einer schnellen, exothermen Reaktion. Es wurden parametrische Studien in bezug auf den Einfluss der Gasgeschwindigkeit und Konzentration der Reaktionsteilnehmer, der Warme und Stoffdilfusionsvermiigen, der W%rmestrahlung, der kinetisthen Parameter und der Teilchendurchmesser durchgefiihrt. Ein Vergleich mit experimentellen Ergebnissen envies sich als sehr befriedigend.
1496