convection model

Compurers them. Engng, Printed in Great Britain.

Vol. 16, No. 8, pp. All rights reserved

721-751.

1992 Copyright

0098-I 354/92 $5.00 + 0.00 6 1992 Pergamon Press Ltd

DYNAMIC BEHAVIOR OF FIXED-BED REACTORS WITH “LARGE-PORE” CATALYSTS: A BIDIMENSIONAL HETEROGENEOUS DIFFUSION/CONVECTION MODEL R. M. QUINTA

FERREIRA,'

A. C. COSTA’ and A. E. RODRIGUES~~

‘Department of Chemical Engineering, University of Coimbra, 3000 Coimbra, Portugal 2Department of Chemical Engineering, University of Porte, 4099 Porto Codex, Portugal (Received 5,February

1992; received for publication 18 February 1992)

Abstract-A complete bidimensional, heterogeneous model taking into account intraparticle diffusion and convection (HT,) was developed and used to analyze the dynamic behavior of a fixed-bedcatalyticreactor

containinglarge-porecatalysts. Results were compared with those obtained by simpler bidimensional models: HT, , a heterogeneous model which considers only intraparticle diffusion and the pseudo-homogeneous (PH) model. Two situationswere tested:the pseudo-steady-state approximationfor the catalyst and the transientregime inside particle.The integra$ionof model equationswas carriedout using the method of lines. Orthogonal collocation in finite elements was employed to discretize space coordinates: reactor axial coordinate (accomplished by available software), reactor radial coordinate and particle space position (handled by the authors). Hot spots obtained with bidimensional HT models are lower than those predicted with unidimensional models. However, the maximum temperature rise predicted by the PH model is higher than in the corresponding unidimensional PH model.

INTRODLJCTION

Since Nir and Pismen (1977) pointed out the improvement in efficiency of large-pore catalysts due to the additional mass transport by intraparticle convection, several studies have been carried out in this area. This pioneering work dealt with the case of isothermal particles with irreversible first-order reactions. It has been followed by other papers dealing with irreversible sequential reactions (Nir, 1977), zeroorder reactions (Rodrigues et al., 1984; Stephanopou10s and Tsiveriotis, 1989), consecutive/parallel reactions (Cresswell, 1985) and non-isothermal reactions (Rodrigues and Quinta Ferreira, 1988). Some experimental results on the importance of intraparticle convection are now available (Rodrigues et al., 1982; Cogan et al., 1982; Cresswell, 1985; Cheng and Zoulalian, 1989; Dogu et al., 1989). Therefore it is not surprising to find in progress the manufacture of large-pore materials: catalysts (pellets, monoliths), ceramic,membranes, HPLC supports (Carta, 1988; Horvath, 1990; Rodrigues et al., 1991a) in order to take advantage of the intraparticle convective phenomenon. It is then important to clarify the behavior of catalytic reactors which make use of such non-conventional catalyst particles. In a previous work, we presented the importance of the intraparticle convective flow on the behavior of ~To whom all correspondence should be addressed.

fixed-bed catalytic reactors operating in the steadystate regime (Rodrigues and Quinta Ferreira, 1990). The critical conditions that separate the operation of the reactor in a stable region from the parametric sensitivity zone, where the temperature runaway may occur, have been established. The dynamic behavior of such reactors was studied further (Quinta Ferreira and Rodrigues, 1991). Here, the simulation results have been achieved by using three unidimensional (1 D) models: the pseudo-homogeneous model (PH); the heterogeneous model with diffusion as the only mechanism of transport inside the solid (HT,); and the heterogeneous model which considers the additional intraparticle transport by convection (HT,,). The main objective of this paper is to analyze the effect of intraparticle convection in large-pore catalysts on the dynamic behavior of packed-bed reactors. In order to reach this goal a bidimensional(2D), heterogeneous model which takes into account irztraparticle diffusion and convection (HT,,) was developed. Results from this complete model will be compared with simpler 20 models: heterogeneous with intraparticle diffusion only (HT,) and pseudohomogeneous (PH). The synthesis of the phthalic anhydride by o -xylene oxidation was used as an example of an exothermic reaction. In such highly exothennic reactions, radial thermal gradients lead to non-uniform temperature and concentration profiles. Therefore, we have used now the 2D models: PH, HT, and HT,,. in order to 721

R. M.

722

QUINTA

establish the influence of the intraparticle convective flow on the transient behavior of fixed-bed catalytic reactors (Quinta Ferreira, 1988). For the HT models the radial heat transfer has been considered separately in the solid and in the fluid phases. For the PH model we have used lumped parameters which represent the contribution of both phases. The dynamic responses have been established for the start-up of the process. We have also tested the validity of the model which considers the steady-state approximation for the non-conventional large-pore catalysts, where the convective transport cannot be neglected. PREVIOUS

WORK

By using the simplest 1D PH model, Hansen and Jorgensen (1976a, b) (for the reaction between H, and O2 over Pt/Al,O,) and Von Doesburg and De Jong (1976) (based on the hydrogenation of CO to CO2 and methane over Ni/A120s) found good agreement with the experimental transient results. Considering the radial gradients in the PH model, Jutan et al. (1977) analyzed the dynamic response of the hydrogenation of butane catalyzed by NijSiO,. Lee and Agnew (1977b) studied the experimental transient behavior of the reactor for the synthesis of the vinyl chloride over Cl,Hg/activated carbon. The results were reasonably adjusted by the 1D HT model with the pseudo-steady-state approximation for the particle. Nevertheless, this approximation leads to transient profiles predicted by the model which developed at higher velocities than the ones obtained experimentally. Sharma and Hughes (1979) concluded that the simple HT model was not adequate for the modeling of the experimental transient response of an adiabatic reactor where the catalytic oxidation of CO to CO, occurred. The complete 2D HT model was applied by LopezIsunza (1983) to a fixed-bed catalytic reactor for the oxidation of o-xylene to phthalic anhydride over V,05/TiOz. While for the steady-state regime model predictions were in good agreement with the experimental response, it failed in the case of the transient behavior. According to the author this was due to the fact the kinetic model did not take into account the catalyst deactivation observed during the experiments. Castro (1983) compared the theoretical results for the dynamic behavior of a fixed-bed catalytic reactor (synthesis of maleic anhydride from partial oxidation of benzene). By using the 2D pseudosteady-state model, which considered the solid temperature as the only transient process (neglecting the accumulation of mass and heat in the fluid phase), he concluded that one does not save computing time when compared with the model which included the

FERREIRA et al.

dynamic features of the fluid and the solid; in both cases the pseudo-steady-state approximation inside particles has been considered. Several studies in the literature compare responses predicted by ID and 2D models, trying to reduce the complexities of the mathematical models. For liquid-phase reactions, Crider and Foss (1966, 1968) and Sinai and Foss (1970) concluded that the radial gradients could be neglected; they have studied, experimentally, the homogeneous reaction of sodium thiosulfate with H,O, over an inert packing. Based on the experimental results of a system with the reaction between H, and 0, catalyzed by Pt/SiO,, Hoiberg et al. (1971) stated that the 2D models (where local linearization has been applied) should be used when the radial temperature gradients are higher than 20%. Lee and Agnew (1977a) concluded that, for the steady-state regime, the maximum temperature predicted by the 2D HT model (isothermal particles) was lower than that obtained with the ID model. However, 2D HT model did not take into account the radial heat transfer through the solid phase. According to Froment (1967), the 2D PH model led to higher hot spots when compared to those predicted by the 1D model. For the HT models, without intraparticle resistances and taking into account the heat transfer separately in the solid and fluid phases, De Wasch and Froment (1971) also reached the same conclusions. These two simulation studies have been limited to the steady-state behavior of the oxidation of the o-xylene in the phthalic anhydride over V,O, . Holton and Trimm (1976) studied the case of mixed heterogeneous and homogeneous reactions in the gas phase, in the steady-state regime. They observed that for the 1D PH model the temperatures obtained were higher than those predicted by the 2D model due to the increase in the heat released in the homogeneous reaction. Pereira Duarte et al. (1984a,b) compared the steady-state response for exothermic reactions using different mathematical models. They used 1D and 2D PH models. In the case of the HT models they analyzed two situations: the use of lumped heat transfer parameters and distributed ones, by considering the heat transfer separately in the fluid and solid phases. Froment and Bischoff (1979) pointed out the importance of the thermal transport through the stationary phase; for typical industrial conditions they stated that the solid and interparticle film account for 25% of the radial heat flow. Several aspects involved in the steady-state modeling of fixed-bed reactors were analyzed by Martinez et al. (1985), who established general criteria to

A

20 heterogeneous diffusion/convection model

choose the model that fulfills the objectives of the simulation. They also discussed different methods for the evaluation of the model heat transfer parameters and analyzed the effect of data dispersion on the model predictions. In order to decrease the computing times required to solve the 2D HT models, MC Greavy and Naim (1977) implemented a strategy to reduce these modeis to the corresponding unidimensional ones, representing the mass and thermal radial gradients by simple algebraic expressions. Mascazzini and Barreto (1989) developed analytical expressions for the evaluation of heat transfer coefficients for the 1D HT fixed-bed reactor model. Those expressions were based on the 2D HT model at the hot spot conditions which occur for highly exothermic reactions. By comparing the results of both 1D and 2D models for an important range of typical industrial conditions, they concluded that very precise estimates can be achieved with the former, if conditions for parametric sensitivity do not occur. Using a two-phase HT model of a fixed-bed catalytic reactor, Odendaal et al. (1987) studied the sensitivity of thermal transport in the axial and radial directions, in the highly exothermic case of naphthalene oxidation. By analyzing the reactor performance with respect to the axial and radial P&let numbers, Biot numbers and local intra- and interphase transport coefficients, the authors came to important conclusions concerning the relative importance of each of these parameters. They also suggested correlations and presented limiting values of thermal transport coefficients in order to establish a “window of reliability” of these key parameters. The parametric sensitivity of the model was then examined regarding the feed concentration, the feed and wall temperatures and the tube to particle ratio, d,/d,, for reasonable values of the main thermal transport parameters. A packed-bed reactor model with axial dispersion in the fluid phase and internal diffusion and reaction within the solid, was used by Ramkrishna and Arce (1989) for a first-order reaction under isothermal conditions. By analyzing the properties of the linear operator implied by the equations of the 1D model, they established an analytical criteria which allows the determination of the range of the parameters values for which the use of the simple PH model is valid. Windes et al. (1989) studied the steady-state and dynamic behavior of a packed-bed reactor for the partial oxidation of methanol to formaldehyde using both 1D and 2D PH and HT models. The PH models used in this work included a non-

723

isothermal effectiveness factor in order to compensate for the differences between the reaction rates in both models, making them equal. These authors still refer to some methods for parameter adjustment in order to approximate the results of both models. The parametric sensitivity analysis they performed pointed out that the 2D models were more sensitive at predicting higher transient hot spots than the 1D models. Several reviews on the dynamic behavior of fixedbed catalytic reactors have been presented (Hlavacek and Von Rompay, 1981; Froment, 1967). MODEL

EQUATIONS FOR TRANSIENT OPERATION OF FIXED-BED REACI.ORS WITH LARGE-PORE CATALYSTS

Since we are dealing with a highly exothermic reaction in a non-adiabatic reactor the concentration and temperature gradients in the radial direction must be taken into account. The mass and heat fluxes in the radial coordinate have been established based on the concept of the effective transport. As a consequence, the radial dispersion coefficient of mass, O,, , has been used. The contribution of the solid to the thermal transport is an important factor, even for typical industrial conditions (Froment and Bischoff, 1979). Therefore, for the HT models the radial heat transfer through each phase has been considered separately. We have then calculated the effective thermal radial conductivities for the fluid, A%, and for the solid, A:, . The heat transfer to the wall has been taken for the two phases by using coefficients for the fluid, hk, and for the catalyst, k;. For the PH model, where concentration and temperature gradients in the solid are neglected, we have used lumped parameters (which include the contribution of both phases): effective radial thermal conductivity, A,, , and wall heat transfer coefficient, h,. Table 1 shows the correlations used for the estimation of these parameters and their values at the reference temperature. It can be observed that in our case L:, is about 27% of the A,, value. These HT non-isothermal systems have, in general, associated to them a quite large gap between the velocities of propagation of the mass and heat waves. This is mainly due to the high thermal capacities of the beds, in which the solid phase has the most important contribution. In fact, in the present study, the time constant the wave for thermal (h:zh=,fL(I - ~hc,,Il(w~~c,~) = 2540 s) is much than the space time for the fluid (r = L/z+ = 0.55 s). Therefore, along the catalytic bed the concentration wave will propagate with a velocity that is 4800 times higher than that of the heat wave.

y[&l]]&(r*,z*,t*)

AWe: whensystem properties donotchange withtemperature: a; =. . . - t$,= I anda{9= 0.

a;,Bi,[4(l,z’,“)-BY]=-

d&k*,z’, t*) =o 8r* fl.1

a&r*,z*,1’) = o de&‘,z*,t’) =o lJr* ~~0 Jr* p-0

z*>o, r’=O:

r*= I:

fb(r*,O,t*)=i

z*=O, r*)O:

fl&*,O,I*)=I

~'20, r*BO: f&r*,z*,O)=O Bh(r*,z*,O)=fJw

~2 I a~h(r*,2*,t*)+la~(rl,2',t*) tailf-69 drs2 r* Jr’ I R0Pe, [

alb(r*,z*,r*) al, =-a; ‘~(r~~~“*)-R,(r*.;‘.l’)Daup[-

Boundary conditions, I*> 0

Initialconditions, t’ = 0

Thermal balance

Massbalance

Table2. Dimensionless modelequations forthetransient, 2DPHmodel

(3)

(2)

(1)

R. M. QUINTA FERREIRAet al.

726

Table 3. Model parameters for the 2D PH model E

Arrhenius number DamkhSlcr

Da =-zk(G)

number

Dimensionless

21.8

y=m

B

adiabatic temperature riac

+&G Toprep l-%,=3

Radial mass PCclet number

0.7 0.7 7677

m

Radial heat P&let number Bi =h,% w ,% El

Wall heal Biot number

1.5 2.08 x 1O-4

Space time for fluid/space time for the thermal wave ‘Radial

P&let

into account the contribution phases. Bidimensionai (HT,,)

(2D)

convection

numbers

based

on

the

of the fluid and solid

particle

diameter:

Bidimensional

(20)

Pe,(d,)=~=20.5;

heterogeneous

dz&kion/convection

model

heterogeneous

For the HT models we have also carried out simulations by using a simpler model where the accumulation of mass inside the catalyst has been neglected. In this case, the intraparticle concentration profiles can be calculated analytically:

X

(IT,)

, a’f,(r,c,

r*, z*, t*) i7rz2

-4+z;f,(r,*,

r*, z*,

t*)

1 .

(18)

The effectiveness factor is still calculated by equation (15). For the pseudo-steady solid, the concentration through

state approximation in the profiles can be calculated

z*, t*)

The effectiveness factor is then:

z+, t*)

sinh Q, exp[a,(2r,* - 1)] - sinh a2 exp[a, (2r,* - 1)] > (16) sinh(a, - a2)

where

and the effectiveness factor is given by 1 ?a= =

d@iiion

For the HT, model only the mass balance for the particle differs from the one established for the complete model HTdc; now & = 0 and so equation (12)--Table 4--is replaced by

f,(r,*, r*, z*, t*) =f.(r*,

r*, z*, t*) =A@*,

CT

model

Table 4 shows the dimensionless equations when the additional transport by convection inside the solid is taken into account. The model parameters given in Table 5 are calculated at T, = 625 K and PO = 0.01 atm. We have considered isothermal particles since the thermicity factor is very low: Prater /I = (- An)o, C.J& T, = 0.005. The effectiveness factor of the slab catalyst referred to the surface conditions is calculated at each time f * along the radial and axial coordinates of the reactor:

f,(r,*.

du,

al

1 a2

coth a, - coth a2’

(17)

For all models the influence of the temperature on the system properties can be observed in Table 6. NUMERICAL.

METHODS

The integration of the system of partial differential equations (PDEs) for the three models has been carried out using the method of lines by makibg use of the package PDECOL (Madsen and Sincovec, 1975). This code discretizes automatically only one space variable that we have chosen to be the axial coordinate along the reactor z*. Therefore, we have

Jr*

ar*’

r*

,

Z*

1

9

r*aO; r;=O

and $=I:

-

ar*

P

r’=l

F-l

fp(r:,r*,z*,O)=O

4cl;,Bi,[r,(r*,z*,r*)-f,(r*,z+,t*)]=-

r*>O, $30:

~2j~(rp*,r*,z*,~*) df,(r;,r*,z+,r*) -‘An ar* ar;z

catalyst Partlcle

K;sBi:[e,(l,r*,r*)-e,]-

aO,(r*,z*, I*)

ae,(r*, Z* I *) ar* ’

= dO,(r*,z*,I*) afJ(r*, 2*,t*) ar* ar* ._s’O 1 =s

=O

=

Me: when systempropertiesdo not changewith temperature:ai =. . = ais = 1 and aiq= 0.

2.30,

Boundaryconditions,t* > 0

,_,

e-0

a;,Bii[eb(l,z*, r*) -e,] = -

ar*

Jfb(r’, z*,t’)

i%*

am*.z*,I*)

~‘30,

r*= I:

z*>O; r*=O:

r*=O: r*)O: ,fb(r*,O,t*)=l B,(r*,O,I*)=l

Boundaryconditions,I* > 0

Initialconditions,I* = 0

1

r*aO: fb(r*,z*,o)=oe,(r*,z*,o)=e, e,cr*,z*,o)=e,

,

~‘20;

Mass balance

I

ar*

t*)-ebv Z* I*)] L* I a2e,(r*,2*,r*)+lae~(r*,2*, I*) +Kl~~ope;n Jr*’ P dr* II [ 7

IWl/Partlcle Interface

[

ta;N~[e,(r*,r*,l*)-Bb(r*,z*,f*)]

-a;N,[f,(r*,z*,I*)-f,(r*,t',f*)]

a2e,(r*,z*,r*)+ideb(r*,Z*,(*)

,a&&*,z*,r*)

&’

,dfbv>z*,f*)

, L* I +a”j$s

=-al

=-a1

-a;N,[B,(r*

z*,I*)

a,*

J&Jr*,

a,*

%(r*,z*> 1’)

Initialconditions,f l = 0

Thermalbalance

Thermalbalance

Mass balance

FIuldPllase

Table4. Dimensionlessequationsfor the transient,2D HT, model

(I2)

(11)

(9)

03)

(7)

(6)

R. M.

728

QUINTA FERREIRA et al.

Table 5. Model parameters for the HT,

model

Arrhenius number Damkhbler number Dimensionless adiabatic temperature rise 2.08 x lo-*

Space time for fluid/space time for the thermal wave Number of film mass transfer units Number of film heat transfer units

143.6

Radial mass P&&t number

7677=

Radial heat P&let number for the fluid Radial heat P&let number for the solid

14,717.

Wall heat Biot number for the fluid

2.4

Wall heat Biot number for the solid

3.8

Mass Biot number

66.9

Number of intrapartick diffusion units

1.2

Thicle modulus

0.8

Intraparticlc mass P&let number

10

‘Radial

P&let

(dpuoPrcp,)/(&r)

numbers =

15.4;

based

PG(d,)

on

the particle diameter: = 39.2.

R,(d,)

= (dpt+,)/(&)

- 20.5,

Pek(d,)

-

= (d,uoP,c,,)l(G)

applied the orthogonal collocation in finite elements (Carey and Finlayson, 1975) to discretize those equations with respect to the reactor radial variable r l and the particle space coordinate rz . The discretization in the solid has been detailed in a previous paper (Quinta Ferreira and Rodrigues, 1991).

functions. The radial coordinate of the reactor was divided into NER subintervals, inside which r* has been normalized between 0 and 1 through a change of variable:

Discretization

where, rbint(kR) is the value of r * at the first node of the kR finite element and hR,, is its length

of the reactor

radial coordinate

The method of orthogonal collocation in finite elements divides the total domain to integrate into subintervals where the exact variables are approximated by trial functions. We have chosen Hermite polynomials as trial functions, since they ensure the continuity of the functions and their first derivatives at the nodes of the finite elements (Finlayson, 1980). In the model equations only the terms corresponding to radial dispersion present derivatives with r * as an independent variable. Therefore, following the methodology used for the particle discretization, only these derivatives are going to be replaced by the trial

uR =

r* -

rbint(k, hR,,

R)

kR = 1,2,.

’

__ , NER,

(21)

(kR=l,...,NER).

Finlayson (1980) states that for high-dimension problems it is not necessary to use high-order polynomials. Using cubic polynomials there are then four unknown parameters in each finite element. We considered the nodes of each finite element and two more interior collocation points, that have been taken as the zeros of the second degree Jacobi polynomial: uR, = 0.2113248654 and uR2 = 0.7886751346. This discretization process will be illustrated with the HT models. Here, the exact solutions fb, 6, and

729

heterogeneousdiffusion/convectionmodel

A 2D

The first derivatives are then discretized as follows:

j&-j, -$+,+2kR--2,

= ae,(r*,

z*,

__ae; 1

r*)

&*

h&c, auR

rf

kR=l,...,

us,

NER;

j=l,2.

(23)

The second derivatives are given by:

1

a2r3,(r*, z*, t*) are2

al&

9 = hR:R

duR2 u.,

=&

a2e,(r*, z*, t*) ara2

1 = hR:, r:

$,

a2e, aUR2”,, I

kR=l,...,

0, are approximated at those collocation points of each kR subinterval by:

The matrixes through:

Hi,,

NER;

z*, t *)

1

LQUR,, 2

*,

t*) =

5 ffifat,+,,_,,

es+/? >z*, t*)

=

fT&4Rj,z*, r*) =

5 H,,aS,+&-2, I=

kR = I,.

..,NER;

(22)

where, ab, at and as are the parameter vectors for the approximation of the corresponding functions: bulk fluid concentration, fb; bulk temperature, t?,; and solid temperature, 8..

are calculated

“+’

auR2

uR/y

1=1,...,4;

and the trial functions H,(uR) j=1,2,

(24)

a2H,(uR)

j=1,2;

1

Bj,

j=l,2.

Hji = H,(uR/),

BjI =

I-I

and

Ail

A, _ an&R) JJauR

%(r;L,

Bj@lt2kR-2,

If, = (1 - uR)*(l H, = uR(1 H, = uR’(3

(251

are (Finlayson, 1980): + 2uR),

uR)‘hR,,, -

H., = uR2(uR

2uR), -

l)hR,,.

(26)

R. M.

730 Table

7. Discretiaed

equations

on the reactor

QUINTA

FERREIRA et 01.

radial coordinate

and on the particle

spatial

coordinate

for the HT’,

model

FluldPbasc (KR-I,....NER; Mass

j-1,2)

balance

(32)

(33) Fluld/Partlcle (KR=I,....NER; Thermal

Interface j=l,2)

balance -1

z*)q DaB exp ~

a;N,&,(uR+

+Gz~e

Initial conditions,

Boundary

r* = 0

conditions,

t* > 0

Z*

L*

5 0:

=* = 0:

’ 1 IhI-

I*)

(34)

‘=‘ATB,,a%+,,-,

j&R,.

z-,0)

g&R,,

z*,o) = 8,

= 0

Jb(uR,. 0, t’) = 1

Catalyst

e;(uR,,r*,O)=6’,

(35)

&(uRj, 0, r*) = l

(36)

Particle

(k=I,...,NE;j=l,2) [for each uR,(j = 1.2) of each subinterval

Mass

z*,

z*, r*) - &(uR,, z*, r’)]

1

0

x’“RJ*

/CR = 1,

, NER]

balance

Initial conditions,

(37)

I* = 0

~&,uR,z+.O)=O

2.20:

(38)

where 1 I CMA TB,, = B,, + LA,, [rbint(kR) + kR,,uR,] hR,, kR,,

(kR=l,...,NER;

j=1.2)

ab,=ag=a&=O a? + mm = 0 at1+2NER = ~;7Bi~(%

- atI +ZNER) (39)

IL12t2NE.=a;~Bi:,(e~--L11+2NEll) I 1 CUAT~,=a;,~B,,-22rZ,-A,,--OH,, kr kr P, = I, +2NE=f,(uR,

Inside each finite element we have defined a local numbering, Z = 1, ___,4. For the entire radial domain a globa numbering system has been used, Z+2kR-2, (kR=1,2 ,_.., NER). In each kR element and for each function there are four unknown _ r(Z=1,...,4), parameters: ab,, 2kR for the bulk fluid concentration; et,+x,R_2(Z = 1, . . . ,4), for the bulk temperature; and as,, 2*R_ 2(Z = 1, . __,4), for

I*, I*) -

(k=l,....NE;j=1.2) -‘A-,.,,=,

m

the solid temperature. However, it can be easily shown that these parameters are given in each case by the corresponding functions and their first derivatives at the knots. Thus, along the reactor radial coordinate there are now 2NER + 2 unknown parameters for each one of the three functions: fb, 0,, and 0,. By using the first boundary conditions for r* = 0 and r* = 1 [equations (1 l), Table 41 and handling

A 2D heterogeneousdiffusion/convectionmodel equations (21)-(26), those parameters:

we get for each function two of

731

and temperature of the solid phase fj = I, 2),

kR=I(r*=O;uR=O); ab2 = 0

at, =0

as, = 0;

(27)

5

I-1

as,Hj,=

&(uR,, z*, t*),

I#2

and kR = NER(r*

= 1; uR = l),

Ii asl+xR-2Hif=

e”,
I-1

at2+2N&R

=

0:;~BiXL

-

atI

+z~~R),

~2+2mR

=

a;8Bi:,(@w

-

=I

+2rmR)-

(28)

The remaining 2NER parameters of each function are related with the corresponding functions: bulk concentration, bulk temperature and solid temperature through equations (22). One then gets the following systems of linear algebraic equations: concentration in the bulk fluid phase (i = 1,2),

5 ab,HjI=x(uRj,

I- I

z*, t*),

I#2

,$,ab,+,R-2

Hj, =j=..(uR], z*, t *),

kR=2,...,NER-11;

... ... . .... .. . ..... .

ab,+mm-2~~=_%(uRj9 z*, t*h

,i,

(29)

temperature of the f7uid phase (j = I, 21, 4

1

I= I

at,Hj,=

&(uRj,z*,

t*),

If2

4

c I- I

atI+2kR--2

H,

= &,(uR,, z*, t*),

kR=2,...,NER-11;

kR=2,...,NER-1.

%+2,-2H,,

=O,(uRj,z*,

+


-

~;sBi~ff_,.a)=,

+INER

t*)-a&,BitB,Hj4,

(31)

For the concentration, the matrix of the coefficients is constant for all the reactor axial positions. Thus, we have calculated once its inverse and the parameters ab have been obtained along the reactor multiplying it by the concentration vector x. However, in the case of the bulk and solid temperatures, the matrixes of the coefficients are dependent of these temperatures (through the parameters CC;,and a;*). Therefore, these systems of equations have been solved in each axial point by using the Gauss elimination method. Substituting now equations (21)-(28) in the terms of the model equations accounting for radial dispersion (Table 4), we will get the discretization in r *. Table 7 shows the final equations for the complete model HT,, . Also the mass balance for the particle is discretized in the spatial position r: ; this discretization process has already been presented previously (Quinta Ferreira and Rodrigues, 1991). Here, NE represents the number of finite elements in the solid, hk is the length of each k subinterval (k = 1, . . _, NE), a,, 2k_ 2 are the parameters which define the approximation of the particle concentration (using also cubic Hermite polynomials) and u is the normalized spatial coordinate inside the elements. For each time t* and each axial position z+, the catalyst effectiveness factor is calculated in the radial collocation points of each NER subinterval through 1 ’ =&uR,

ff

z*, t*) t_,

+,

ar+,-2C&),

with = &(uRj, z*, t*) - a;,Bi~B,H,,,

(30)

CH, =

H,(u) du.

(40)

732

R. M.

Number

of ordinary

differential

equations

QUINTA

(ODES)

As presented above we have handled the first steps of the discretization process using cubic Hermite polynomials for the trial functions of the orthogonal collocation method in finite elements. The discretization in the axial coordinate of the reactor z* is automatically done by the PDECOL package which also uses the same method with B-splines as trial functions (De Boor, 1977). Defining NE as the number of finite elements in the space coordinate of the catalyst particle, NER as the number of subintervals in the radial coordinate of the reactor and NINT the number of elements in the axial coordinate of the reactor, we will obtain first for each of the interior collocation points of the axial position of the reactor the following number of PDEs discretized in r* and rz: Fluid phase

Number of equations

Concentration Temperature

2NER 2NER

Particle (isothermal) Temperature Concentration

2NER 2NE x 2NER

Taking into account that in each of the axial subintervats there are also two interior collocation points (cubic B-splines), the final numbers of ODES are: PH model

[2 x 2NER][2

HT model

[(2NE

NINT],

+ 3) x 2NER][2

NINT].

These ODES form an initial value problem and are then solved by an implicit time integrator GEARIB (Hindmarsh, 1976), included in the PDECOL package, specific for problems with banded Jacobians. The use of the 2D models strongly increases the numerical dimension of the problem, increasing its stiffness. Thus, the strategy we have followed of dividing the reactor into sections for the case of the 1D models (Quinta Ferreira and Rodrigues, 1991) is even more adequate in this case. Therefore, we have also used a large number of sections and a low number of axial finite elements inside each section in order to minimize the number of ODES to solve simultaneously. The calculations in one section were carried out just after all the dynamic response of the previous section had been achieved (i.e. after the determination of the concentration and temperature profiles for the final steady state). For each integration time the boundary conditions for a given interior section were determined by an interpolation process using cubic splines (De Boor, 1978) of some of the variables calculated at the node between the last two finite elements of the previous section. Therefore, the con-

FERREIRA et al.

tinuity of the functions and their first derivatives between sections is guaranteed. COMPUTER

RESULTS

The dynamic behavior of the reactor was studied in all cases for the startup of the process. Initially the reactor was heated at the wall temperature and was free of reactant mixture. At the time t* = 0 the feed stream entered the reactor; thus, a step change in the inlet concentration was introduced in the system. This situation corresponds to the most difficult numerical problem to deal with, since the discontinuity introduced at the boundary of the system propagates along the reactor in sharp fronts. It is a hard task to obtain precise computer solutions for this kind of problem, which has the associated phenomena of numerical dissipation and numerical dispersion (Lapidus and Pinder, 1982). When dealing with 1D models, the influence of some parameters on the numerical transient response has been studied in order to minimize those phenomena (Quinta Ferreira and Rodrigues, 1991). According to those results, we have the same final parameters here: 50 sections along the reactor (N, = 50), two finite elements in each section (NINT = 2) and 398 interpolation points between sections (to calculate the boundary conditions of each section inside the reactor). These parameters also led to dynamic responses of the present system, almost free of numerical fluctuations, and just a little numerical dissipation associated with the sharp front. Pseudo -steady -state

regime

inside the particles

Figure l(a) shows the transient axial profiles of the mean radial values of the bulk-phase temperature and concentration for the HT,, model (PO = 0.017 atm, T, = T, = 625 K and uo(SPT) = 1 m/s) when the accumulation term in the solid is neglected. The properties of the system were allowed to change with temperature along the reactor. Due to the high ratio between the propagation velocity of mass and thermal waves (zh = 4800~), the first stage corresponds to a pseudo-isothermal behavior of the reactor which occurs in a very short time (t * -c 10). Here the process is controlled by the propagation of the concentration wave and the discontinuity introduced at the inlet propagates along the reactor Lwtihdecreasing amplitude due to the chemical- reaction. For t * > 10, the heat wave dominates the overall response during approx. 99.5% of the total time needed to obtain the final steady-state profiles. Figure l(b) shows the steady-state radial temperature profiles for the fluid phase for different axial positions in the reactor, (T,, is the mean radial

733

A 2D heterogeneous diffusion/convection model

t* A

. .

.

..

z----D -.-__-._E-mF __-_G--H

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

;;: 10.0 100.0

---I -------

J --.-L -

0.1 00::

-

400.0 1000.0 2500.0

1.0

(a)

1.04

1.03

t

F

2500

I=4 1000, ax*

l1.01 t: 1.00 I 0

0.5

I.

3

1.02 1.00

L

100

I

1.0

r*

0

0.5

.,

,

1.0

r*

(b) (c) Fig. 1. Model HT, with the pseudo-steady-state approximation for intraparticleconcentrationprofiles (A, = 10; PO = 0.017 atm; To = T, = 625 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different 1* values: A, 0.1; B. 0.2; C, 0.4, D, 0.6; E, 0.8; F, 1.0; G. 10.0; H, 190.0; I, 490.0; J, 1000.0; I,, 2500.0. (b) Radial temperature profiles: steady-state solution (t* = 2500) for different axial positions z*. (c) Radial temperature profiles: transient profiles for I* = 0.1.

temperature)_ For Z* = 0.1 (near the hot spot), Fig. l(c) represents the transient radial profiles of temperature; we have used two finite elements in the radial coordinate, NER .- 2. Our results show that the radial gradients are severe, mainly at the hot spot. Here, the difference between the centerline and wall temperatures is 13 K, which corresponds to 60% of the maximum temperature rise in the center of the _ reactor (-22 K). For a higher inlet and wall temperatures, the magnitudes of the hot spots and the final reactant

conversion are also higher, as seen in Fig. 2(a): which shows the bulk temperature and concentration profiles along the catalytic bed (mean radial values) when T, = T, = 660 K and P,, = 0.017 atm. We have also plotted the steady-state profiles obtained with the steady-state models (curves M) that are in quite good agreement with the steady-state calculated through the startup transient responses. The radial thermal gradients shown in Figs 2(b,c) indicate a maximum radial difference of 56 K, much higher than in the previous case.

R. M. QUINTA FERREIRAet al.

734 1.12 1.1 l.OB

I-” .

1.06

E

1.04

i= 1.02 1

0-0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.0

0.9

1.0

1 0.8 ’ ii 8

06. 0.4 0.2

/J.L

0 0.0

I

I

I

0.1

0.2

0.3

I

OA_

I

zti,

1

f.6

L

1

I

0.7

0.6

0.9

1.0

a

A-0.1

a

B-O.2

0

c-o.4

I3

D-0.6

m

E- 0.8

e

F-l.0

a

G- 10.0

m

H- 100.0

a

L-400.0

0

J-1ooQ.O

8

L-25w.o

(a) 1.16

c

1.16

1.14

L

1.14 1.12 1.10 1.06 1.06

0.6

0.5 r*

(b) Fig. 2. Model HT,, with the pseudo-steady-state approximation for intraparticle concentration profiles (A,,,= 10; P,, = 0.017 atm; r, = T, = 660 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different t* values: A, 0.1; B, 0.2; C, 0.4, D, 0.6, E, 0.8; F. 1.0; G, 10.0; H, 100.0; I, 400.0; J, 1000.0, L, 2500.0; M, steady state. (b) Radial temperature profiles: steady-state solution (I * = 2500) for different axial positions z l. (c) Radial temperature protiles: transient profiles for z* = 0.07. The computing time on a Cyber 830 was 10 h 38 mm, about 17 times higher than the time needed to obtain the transient numerical solution of the 1D model. This drastic increase is due to the higher number of equations to be solved simultaneously when radial gradients are taken into account, increasing the stiffness of the problem. Since the radial

profiles are not very different from the ideal parabolic shape we have tested the numerical solution of the system by using just one finite element in the reactor radial coordinate (NEZ? = l), which led to the same results achieved with iVER = 2 for both situations analyzed previously (T, = 625 and 660 K). However, the CPU time is significantly reduced by

A 2D heterogeneous diffusion/convection model

735

,,oo ~_.~_---.--_.~._._._._._.~-. A-F 0.99

t* A

0.1

. . . ..m.........

;-----

::a

E-wF _-__ GH-s._I --_--

J L -

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

-

::t &k: 1oo:o 400.0

-

1000.0

2500.0

1.0

2' = Z/L (aI

Fig. 3. Model HT.,, with the pseudo-steady-state approximation for intraparticle concentration profiles (A, = 10; PO= 0.017 atm; T, = T, = 625 K). System properties constant with temperature. (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different t* values: A, 0.1; B, 0.2; C, 0.4; D, 0.6; E, 0.8; F, 1.0; G, 10.0; H, 100.0; I, 400.0; J, 1000.0; L, 2500.0. (b) Radial temperature profiles: steady-state solution (t * = 2500) for different axial positions z*. (c) Radial temperature profiles: transient profilesfor I l = 0.1.

70%. Moreover, the results obtained for the 1D models have shown that most of the CPU time is employed in obtaining the numerical solution of the first pseudo-isothermal state (N 85% of the total CPU time). This shows the great numerical difficulties associated with such problems, where the discontinuities introduced at the boundary of the system propagate along the catalytic bed. Therefore, further significant reductions in the computing times can be achieved if one takes the pseudo-isothermal solution as the initial condition of the problem.

Figure 3(a) shows the longitudinal profiles and Figs 3(b, c) represent the radial temperature gradients when constant system properties (calculated at the feed temperature, T, = 625 K) are considered. One can observe a slight increase in the hot spot and the radial gradients, when compared with the previous case [Figs l(a-c)]. Even for the highest inlet temperature Lr’,, = T, = 660 K, the results obtained with constant system properties [Figs 4(a-c)] are not very different from those obtained when properties change with temperature [Figs 2(a-c)].

R. M. QUINTAFERREIRA et al.

736 1.12 1.1 1.06 1.06 1.04 1.02 1

I

I

I

0.0

0.1

0.2

I

0.3

I

0.4

1

I

0.5

0.6

I

0.7

I

0.8

I

0.9

I 1.0

ia

I

1 0.8 s

0.6

-is 8

0.4 0.2 0

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.6

0.9

1.0

t* A-0.1

E

B-0.2

E3 I3 I3

c-o.4 D-0.6 E-O.8

8

F-l.0

0

G - 10.0

El

n-100.0

r3

I-400.0

0

J-IWO.0

8

L-2500.0

Z’=Z/L

(a)

1.16

1.16

1.14

1.14

k_

1.08

t1.06

1.06

1.04 1.02

0.6

(b) Fig.

4. Model HT,, with (A, = 10; PO = 0.017 atm; axial protiles of the radial A. 0.1; B, 0.2; C, 0.4; D,

the pseudo-steady-state approximation for intraparticle concentration profiles To = T, = 660 K). System properties constant with temperature. (a) Transient mean temperature and concentration in the bulk phase for differentI * values:

0.6; E, 0.8; F, 1.0; G, 10.0; H, 100.0; I, 400.0; J, I000.0; L, 2500.0; M, steady state. (h) Radial temperatureprofiles:steady-statesolution(t * = 2500) for differentaxialpositionsz*. (c) Radial temperatureprofiles: transientprofilesfor I* = 0.07.

The software used to obtain the transient nnmerical responses of such reactor models is quite complex. Taking into account the influence of the temperature changes on the system properties means that they must be recalculated at each step, which requires an

enormous effort in writing and testing the overall program. Moreover, that influence must be important for more drastic operating conditions, with higher temperatures_ But probably, these cases will correspond to undesired conditions, i.e. close to the

737

A 2D heterogeneous diffusion/convection model

0.60 s

2 0

0.60

.-

10.0 100.0 400.0 1000.0 2500.0

0.40 0.20

1.04 t 1.03

t= > l-

1.02 1.01 1 .oo

I

I

,

,

,

0.5

I

1.0

b

1’ 7

1000~2500

400

100 10

0.5

0

r*

1.0

r”

W

(c)

Fig. 5. Model HT, with the pseudo-steady-state approximation for intraparticle concentration profiles (A,,, = 10; PO = 0.017 atm; To = T,, = 625 K). (a) Transient axial profilesof the radialmean temperature and concentrationin the bulk phase for different t* values: A, 0.1; B, 0.2; C, 0.4, D, 0.6, E, 0.8; F, 1.0; G. 10.0; H, 100.0; I, 400.0; J, 1800.0; L, 2500.0. (b) Radial temperature profiles: steady-state solution (t l = 2500) for different axial positions .z*. (c) Radial temperature profiles: transient profiles for z * = 0.1.

parametric sensitivity of the reactor that will lead to the dangerous situation of the temperature runaway. Therefore, it is our belief that it is not worthwhile to consider the properties’ changes with temperature, taking into account the highly significant simplifications introduced by this approximation. This is still supported by the fact that the critical conditions achieved for this same system (using the steady-state regime) were more conservative, and therefore safer, when that inl?uence was neglected (Rodrigues and Quinta Ferreira, 1990). For as

the

the

HT,

only

model,

which

intraparticle

considers

mechanism

diffusion of

mass

transport, Fig. 5(a) represents the axial profiles of the bulk-phase concentration and temperature (radial mean values&-P, = 0.017 atm, T, = T, = 625 K and q,(SPT) = 1 m/s. The radial gradients of the temperature are still significant-Figs 5(b,c). During the pseudo-isothermal behavior of the reactor t * < 10, the concentration decays sharply to zero in each axial position of the catalytic bed. This is independent of what is going on inside the particles. Therefore, the sharp front arrives at the same time at a given point of the reactor for both models [Figs l(a) and 5(a)]_ However, since the reaction rate is higher for the HT, model, owing to the

R. M. QUANTAFERREIRAet al.

738

0

I

0.0 I

0.1

0.2

0.3

0.4

I

I

I

I

I

I

1

I

I

0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.5

0.6

0.7

0.0

0.9

1.0

0.8 E

0.6 0.4

B

0.2

I

1.0

a

J-1000.0

8

L-25co.o

Z’=Z/L (a)

1.06 1.04 c 2

1.02 1.00

1 .oo

0

0.5

1.0

r*

I 0

I

I

10

I

I

I

0.5

I

I

I 1 .o

r*

Fig. 6. Model HT, with the pseudo-steady-state approximation for intraparticle concentration profiles (A, = 0, PO= 0.017 atm; To = T, = 660 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different I* values: A, 0.1; B, 0.2; C, 0.4; D, 0.6; E, 0.8; F, 1.0; G, 10.0; H, lOO.0; I, 400.0; J. 1000.0; L, 2500.0, M. steady state. (b) Radial temperature protiles: steady-state solution (t * = 2500) for different axial positions I*. (c) Radial temperature profiles: transient profiles for z* = 0.09.

higher catalyst efficiency, the bulk concentration before the sharp front is now lower. When the propagation of the thermal wave is dominant, all the axial concentration profiles have lower values. Then, for the HT, model the maximum temperatures are not so high, since the lower reaction rates also lead to lower reaction heat generation. Neglecting the mass convection inside the solid one gets a lower reactant conversion, as expected. These differences are even more evident when we compare the results obtained for both models with more severe operating temperatures (TO = T, = 660 K): Figs 6(a-c) for the HT,

model and Figs 2(a-c) for the HT,, model. The conversion predicted from this HT, model is higher (9%) than that predicted from HT,. This fact is identified as one of the main objectives in promoting the research for large-pore catalysts where convection can be implemented. Intraparticle convection drives catalyst effectiveness between diffusion-controlled and kinetic-controlled limits. We will then have higher final conversions without increasing the severity of the operating conditions; in particular, it is important to reduce the wall/feed temperatures in large-production systems to save energy costs.

739

A 2D heterogeneousdiffusion/convectionmodel 1.04 1.03

I -

5m

-L +---J<_L.

1.02 -

TO lw

-‘---~4_ / .I-01 - r _ H _ ----__

1 .O (m/s)

UO (SPT)

----_-____

G

--_

__~_._~___.___._._____.___

1.00

10

(atm) 625 (K) 625 (K)

0.0170

PO

.-.

A-F

t= * .............._

0.1

B0

D-E--fz -_-__ G----H ---, -__-___ J L 0.1

0

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

ES 0:6

-w--m

-

?i 10:0 100.0 400.0 1000.0 2500.0

1.0

if* = Z/L (a)

”

1 .o

“2

0

0.5

r’

1

.o

r*

(b)

(c1

concentration profiles (A,,, = LO; HTao with the transient regime for intraparticle T, = T, = 625 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different t* values: A, 0.1; B, 0.2; C, 0.4; D, 0.6; E, 0.8; F, 1.0; G, 10.0; H, 100.0; I, 400.0; J, 1000.0; L, 2500.0. (b) Radial temperature profiles: steady-state solution (r * = 2500) for different axial positions z *. (c) Radial temperature profiles: transient profiles for z * - 0.1. Fig. 7. Model PO = 0.017 atm;

The

hot

spots

predicted

lower than those obtained models

(Quinta

means

that

improved persion

Ferreira

the heat by

system

are 1D

HT

and Rodrigues, the

the catalytic

fluid and solid phases. same

by 2D

transfer

considering

through

models

with the corresponding 1991).

to the wall radial bed

in the steady-state

thermal

separately

In a previous

work,

regime

This

has been disin the

using the

and the 1D

which gives the prediction of the operating conditions of the stable and unstable operation of the reactor (Rodrigues and Quinta Ferreira, 1990). The present

models,

we established

the runaway

study shows that if the 2D modeIs

diagram

had been used, we

would get less severe conditions. Therefore, the conditions obtained before are conservative, which enables us to have confidence in those results. Transient regime inside the catalyst particle Figure 7(a) shows the dynamic response of the complete transient model HT,. Figures 7(b,c) represent the radial temperature profiles. At the initial times, one can observe the competition between the two mechanisms of transport inside the solid. When the intraparticle convection flow is important the concentration wave in the bulk phase is delayed, due to a higher penetration of the reactant inside the

R. M.

QUINTA

FERREIRA

et al.

1.12

1.1 1.06 g

1.06

-z e

1.04 1.02 1

0.0

0.1

0.2

0.3

0.4

0.6

0.5

0.7

0.9

0.8

1 .O

I3

I

1 0.6

1’ A-0.1

e

6 - 0.2

0

c - 0.4

s

I3

0 - 0.6

0.6

I3

E-O.8

2

0.4

6

F-1.0

0

G-10.0

B

H-100.0

8 0.2 0 I

0.0

0.1

I

I

0.2

I

I

0.3

I

0.4 0.5 Z’=Z/L

I

0.6

0.7

I

I

0.9

0.0

1.0

I3

I-400.0

0

J-1000.0

6

L-2wo.o

(a) 1.16

1.16

1.14

1.14

t

1.12 1.10 t” I=

1.08 1.06

0.6

0.5

0.5

r’

r* 03

@‘I

with the transient regime for intraparticle concentration profiles (A, = 10; Fig. 8. Model HT, PO= 0.017 atm; To = T, = 660 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different I* values: A, 0.1; B, 0.2; C, 0.4, D, 0.6; E, 0.8; F, 1.0; G, 10.0; H, 100.0; I, 400.0; J, 1000.0; L, 2500.0; M, steady state. (b) Radial temperature profiles: steady-state solution (I * = 2500) for different axial positions z*. (c) Radial temperature profiles: transient profiles for z* = 0.07.

particles When

[upper

convection, phase

part of curves

intraparticle

the reactant

is higher.

separated

diffusion

These

Et and C,

dominates

concentration

two concentration

by the plateau

shown

later times, the intraparticle

Fig.

be important

throughout

tant conversion

in the bulk

conventional

fronts

in those curves.

convective

7(a)]_

intraparticle are For

flow seems to

HT,

This behavior

The final reacfor the

model.

can be better understood

ing the time constant of transport

the reactor.

is higher than that achieved

by compar-

values for the two mechanisms

in the catalyst.

The time constant

for the

A ZD heterogeneous diffusion/convection

741

model

1.06 I-O -i% kf

1.04 1.02 A-F

0.98 A

0”

0.80

-= 0”

0.60

. . . . . . ..*......

::: 0.4

Bpa---_E---

0.8 In 10.0 100.0 400.0

.._

p -_-”

0.40 E

-_.....L_ 0

0.1

0.2

‘_ 0.3

2’

0.5

0.6

0.7

J -

‘\

i ..__

0.4

x

0.6

m-s, -e-e_-

1000.0

-.

0.9

1.0

= Z/L

(a)

1.08

-

.O

(cl Fig. 9. Model PH (PO= 0.017 atm; To = T, = 625 K). (a) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for different I* values: A, 0.1; B, 0.2, C, 0.4; D, 0.6; E, 0.8; F, 1.0; G, 10.0; H, 100.0; I, 400.0; J, 1000.0; L, 2500.0. (b) Radial temperature profiles: steady-state solution (t+ = 2500) for different axial positions z*. (c) Radial temperature profiles: transient profiles for z* = 0.2. to conclude

that the results just differ for the initial

pore diffusion is T* = tR:JD, = 0.5 s and that for intraparticle convection is tC= ER,/v,= 0.05s (the intraparticle P&let number represents the ratio

isothermal

between

ent behavior

these

Therefore, internal goes

two

when

this phenomenon

resistances

inside

the

concentration

steady-state

are lowered

catalyst--causing

R, = q,/r= = 10). is considered, and

more a

delay

the

reactant in

the

wave for the fluid phase.

A comparative the system

parameters,

study

of the dynamic

response

of

for this case and that using the pseudoapproximation

in the catalyst

enables

us

times (z* < lo),

when the overall

[Figs l(a)

process

is pseudo-

and 7(a)]_ Therefore,

the transi-

of both models

is similar almost

all the

Nevertheless, for the solid transient regime the computing times are more than twice those required by the simplified model. Hence, we suggest performing system modeling with HT models, which neglect the accumulation term inside the particles. These conclusions are also supported by Ferguson and Finlayson (1974), who con&m the time

(t*w

10 to 2500).

742

R. M. QUINTA FERREIRA et al. PH model TO=Tw=625 K

700

-

2 500 ---ii

-

+6

300

-

100

-

700 g s

500 ii

300

__-__-*

=_--

t ,,; Bi,

HTdo ,&

100

5

= 1.5

I 10

15

20

25

-

80

-

60

-

40

‘;j6 $

,____-._______________ Xf AT,,,

a

______-a----

30

, 5

!! ‘d 300

-

100

t :

L

*-._ Xf ---*-_-___.___________________________ AT,,, 5

10

15

-

60

-

40

3

.A--‘-‘Tdc4 ,_/’ 10

= 3.0

Bi,

15

I 20

25

,-

30

(b)

700

_:

60

P%h (W

(a)

500

-

_..*

Perh (WI

8

- 100

20

25

-

100

-

60

-

60

-

40

‘;; 6

z

700

>i

500

-_ ---._____-__ _.______ ---.---.-.-__-._______ *t II

100 -

I-

30

5

Si,

10

15

20

25

30

Per,,, Wp)

(cl

(W

Fig. 10. Parametric sensitivity of the PH model (P,, = 0.017 atm; T, = T, = 625 K). Hot spot magnitudes (AT,,,,,) and final conversions (q) as a function of: (a) Pe,(d,)(Bi, = 1.5); (b) Pe, (d,)(Bi, = 3.0); (c) Bi,;

(d) I%,($,). validity of such an approximation when the ratio between the time constant for heat transfer to the catalyst (T,_) and the time constant for intraparticle mass transport (thd) is higher than 1. In fact, rhs= p,c,/h = 18.5 s in our system and rd = e,R;/D, = 0.5 s, so t,,Jzd = 37 2~ 1. The same conclusions can be achieved for the higher conversion situations if we compare the dynamic response of the system when the transient regime in the particles is also taken into account: Figs S(a-c) (r, = T, = 660 K) and the corresponding profiles obtained with the pseudo-steady-state approach for the solid, Figs 2(a-c). Pseudo -homogeneous (PH) model For the simplest 2D model, Fig. 9(a) represents the longitudinal temperature and concentration profiles when T, = T, = 625 K. Figures 9(b,c) show the radial temperature gradients. Using the same operating conditions as before, one can observe much higher hot spots than for the HT, model, Fig. l(a). These significant differences between the PH model and the HT, model with intraparticle convection observed in the transient and steady-state responses are even more pronounced for higher inlet temperatures.

When T, = T, = 660 K, one cannot even obtain the numerical solution of the PDEs of the PH model. The high hot spots and great parametric sensitivity in this case leads to strong fluctuations in the temperature profiles with further divergence of the numerical method. Therefore, these results suggest that the PH model is not adequate for modeling our system. As mentioned earlier, we have used here lumped parameters for the radial thermal conductivity (a,,) and for the wall heat transfer coefficient (h,), which include the influence of the thermal transport in the fluid and solid phases. In the next section we will address in more detail the problem of equivalence between this simple PH model and the HT one which accounts for the inter- and intraparticle resistance by analyzing the values of the lumped radial parameters which would allow a process representation similar to that obtained with the full model. The maximum temperature rise in the catalytic bed is now higher than predicted by the 1D PH model (Quinta Ferreira and Rodrigues, 1991). This suggests that the parameters we used in the first case (Pellet numbers and Biot numbers) led to higher radial thermal resistances, located within the bed and near the wall, when compared with the 1D model where

743

A 2D heterogeneousdiffusion/convectionmodel PH model TO=Tw-660 K

700

8

_

500

-

300

-

__--

-

60

-

60

2 lq

-r _______._________.__--.~~-~~-~-----~

- 100 700 *________‘_~~‘______~~~_._____________--_ mar Bi /

loo-

5

1.5

w=

3 Y x-

HTdc 10

g ‘d

ATmax

500 2

300

Bi,-3.0

-

100

-

60

-

60

3z

100 15 %h

20

25

5

10

15

20

25

30

,

,

,

,

,

,

5

10

20

25

30

30

(‘NJ

(a)

Xf *------------_______-_-----___________.____.

700 g

500

f Q

300

100

100

1001 5.

10

15

20

25

30

Si,

15 pe,,,,

(cl

1

Wp)

(d)

Fig. 11. Parametricsensitivityof the PH model (P,, = 0.017 atm; T,,= T, = 660 K. Hot spot magnitudes (AT_,) and finalconversions(x,) as a functionof: (a) Pe*(d,,)(Bi, = 1.5); (b) Pe, (d,)(Bi, = 3.0); (c) Bi,; (4 Pe, (4 ). the radial heat transfer is considered to he located only near the wall. So, it seems that using the 2D model, one will obtain operating critical conditions more severe than those suggested by 1D model (Rodrigues and Quinta Ferreira, 1990). However, we concluded in that study that the PH model is not appropriate for the prediction of the runaway conditions of such HT systems. Equivalence sion fconvection (PH)

between (NT,,)

the

heterogeneous

dzfi -

and the pseudo -homogeneous

models

In order to analyze the validity of using a simple PH model to describe the behavior of highly exothermic reactors with large-pore catalysts, we tested the sensitivity of the steady-state response for P&let and wall Biot numbers with respect to the hot spot magnitude (AT,,,,,) and the final reactant conversion (q), using the HT,, and PH models for a large range of these radial parameters values. For the PH model, Figs 1O(a4) and 11(a-d) show the influence of the effective radial parameters we,&), Bi, and Pe,(&); the P&clet numbers used in this study are based on the particle diameter] on AT,, and xr for the two inlet temperature situation discussed previously in the transient analysis: T, = T, = 625 K, Figs lO(ad); and

T, = T, = 660 K, Figs 11 (a-d). An increase in Pe,, i.e. a decrease in the radial thermal dispersion, leads to an increase in the maximum temperature and final conversion due to the higher reaction rates occurring in the process. In the cases analyzed, there is a Pe,, value which corresponds to a strong sensitivity of the system leading to the temperature runaway behavior, with extremely high hot spots (AT,, s 400 K), as shown in Figs 1O(ad) and 11(ad) [in parts (a) we have tixed Bi, = 1.5 and in parts (b) Bi, = 31. For higher Bi, values, this situation will also take place for higher P&let numbers, since the thermal stability of the process can be maintained with lower heat radial dispersion providing that the heat ransfer to the wall is increased. For fixed Pe, values (obtained through the correlations in Table l), Figs 10(c) and 1 l(c) show a pronounced decrease in AT,, and xr for low Bi, values; furthermore, the system response is almost insensitive to changes in this parameter. The effect of the radial mass dispersion can be neglected as seen in Figs 10(d) and 1 l(d), where AT_ and x, are not influenced by Pe,. We can now observe that the thermal radial parameter values we used for the dynamic analysis,

R. M. QUINTA FERREIRAet al.

744

HT.,, model TO-Tw-625

K

1

110

110

100

llO0

90

90 80

8

8

70

1

60

+Q

2

50

XI ________________ ___... ___.- _--.-

70

s

50

---.---..A--

___._---

-

60

3

-

60

x-

-

40

40 30

30

*T,,,

10

10 10

20

40

30

50

60

10

70

20

30

pefh(dp)

40 P&,

(a)

50

60

70

WI (b)

110 90

90 8

100

8

60

3 a

70 ii

$

50

100 60

7o 50

60

60 30

30 40

*Trnax

10 5

IO

15

20

10 25

I

I 5

30

10

15

20

25

30

Bi;

Bit,

(d)

(Cl

110 100 90 8

60 I

70

l-E a

50

60

g x-

40 30 10 10

20

30

40

per,,,

50

60

70

(dp)

08

Fig. 12. Parametric sensitivity of the HT,, model (PO= 0.017 atm; To = T, = 625 K). Hot spot magnitudes (AT_) and final conversions {x0 as a function of: (a) Pe:,,@); (b) (Pek)(d,); (c) Bi:; (d) Bi&; (e) Pe,(%).

Bi, = 1.5 and Pe,,(d,) = 10.5 (Table 2), correspond to a situation either close to the runaway phenomenon when T, = T, = 625 K, Figs. 10(a), or even to a thermal instability behavior, for T, = T, = 660 K, Fig. 11(a) (these particular cases are indicated by solid circles). This explains the numerical difficulties

associated with the determination of the transient solutions of the PH model, in particular, when r,=2-,=66OK. The previous results allow the prediction of the lumped radial heat parameter values which should be used for the PH model in order to approximate its

A 2D heterogeneous diffusion/convection HT,,

745

model

model

TO=Tw=660

K

150 110 100 90

110

8

90

2 IO

70

//___..-.-~----

..____.

80

z

60

g

:2--I

B 70 !z 5

Xf _________._._....--_..

50

ATlllax

50

30

30 10

10 10

20

30

40

+,

50

60

10

70

20

30

40

P&t

(dp)

50

60

70

(@I

W

(a)

I 110

110

90

90

-

8 2 $

7o

-

50

80

-

30 10

100

60

-

8 2

z C-

7o

$

80

g

60

=

50 30

40

IO

I

I

5

10

15

20

25

5

30

10

15

20

25

30

si”,

Sit,

00

(C)

110

-

100

-

80

-

60

90

zp kE u

70

2 az

50

10 10

20

30

40

P%

50

60

70

RW @I

Fig. 13. Parametric sensitivity of the HT, model (PO = 0.017 atm; To = T, = 660 K). Hot spot magnitudes (AT,,,) and final conversions (q) as a function of (a) Pe$,(d,); (b) Pet@,); (c) Bi:; (d) Bi;; (e) Pe,(4).

to the value obtained with the HT, model. For this complete model, the steady-state results we achieved were: AT_, = 16 K, xr= 61% for To = T, = 625 K; and AT,,,,, = 64 K, xr = 94% for To = T, = 600 K. The different parameter values esti-

response

mated through different correlations allow the possible use of various sets of model parameters. Therefore, according to the values for the first case, TO= T, = 425 K, one could use for the PH model (among other possible combinations) Bi,,. = 1.5 and

R. M. QIJINTA FERREIRAel al.

746

1

0.2

0.1

0.0

0.3

0.4

OS

0.6

0.7

0.6

O-9

1.0

0.6

8

l+Tdc mocief

A -

0.6

PH moded B _____- -

1.5

C m-m.._ 0.2

,

0

I

0.1

0.0

1

0.2

I

0.3

I

0.4

0.5

I

0.6

I

0.7

0.6

0.9

Pe&

WW

3

D.--m

1.5

E-

1.5

6.5 $0 a 10.5

1.0

Z’=ZlL

(a)

1.03

,

1

E3

0.99

1 0.0

I

0.1

I 0.0

I

0.2

I

0.1

0.2

I

1

0.3

I

0.3

0.4

I

0.5

I

0.4

I

0.5

I

I

0.6

I

0.6

I

0.7

0.8

I

I

0.7

0.6

I

I

0.9

,

0.9

1.0

t+ A-O.1

e

B - 0.2

E3

c-o.4

m

D-O.6

I3

E-0.8

Ei

F-l.0

E3

a-

I3

n-100.0

10.0

I3

l-400.0

0

J-IWO.0

E

L-2500.0

I 1.0

Z’=Z/L (b)

Fig. 14. (PO = O.Ot7 atm; T, = T, = 625 K). (a) Steady-state axial profiles of the radial mean temperature and concentration in the bulk phase for the HT,, and PH models with different values of Pe, and B&.. (b) Transient axial profiles of the radial mean temperature and concentration in the bulk phase for the PH model with Pe,(Cr,) = 6.5 and Bi, = l.S for different I * values: A, 0.1; B, 0.2; C, 0.4; D, 0.6; E, 0.8; F, 1.0; G. 10.0; H, 100.0; I, 400.0; J. 1000.0; L. 2500.0.

= 10, as shown Pe,(ci,) - 6.5 or Bi, = 3 and Pe,@,) in Figs IO(a,b) (these two situations are indicated by a white rectangle). When To = T, = 660 K, there are also two possible sets of lumped parameters which = 3.9 or could he chosen: Bi, = 1.5 and Pe,(4) Bi, = 3 and Pe,(d,) = 5.8 [Figs 1 l(a,b)].

To take advantage of the utilization of the simplified PH model, the lumped parameters calculated through the correlations [Bi, = 1.5, Pe,&) = 10.51 have to be adjusted in order to increase the heat transfer to the surroundings, either by increasing the wall heat transfer coe.fEcient and/or increasing the

747

A 2D heterogeneousdiffusion/convectionmodel 1.16 1.13 1.11 g

1.08

E

fk+ i

A

1.07 1.03 1.03 1.Ol

c&v.“____

O.SQ ,o.o

0.1

0.2

0.3

0.4

0.5

4

0.3

0.7

0.8

0.3

1.0 I

A-

0.8

2 8

HTdc model

0.6 0.4 0.2 0 0.0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Z’=Zfl (a)

E3 B 0

m m 1.Ol

tA-O.1 B-O.2 c-o.4 D-O.6 E-O.8

6

F-1.0

E3 I3

a-10.0

H-100.0

m

I-400.0

0

J-1000.0

E

L-z500.0

(b)

Fig. 15. (PO = 0.017 atm; To = T, = 660 K). (a) Steady-state axial profiles of the radial mean temperature and concentration in the bulk phase for the HT, and PH models with different values of Pe,,, and Bii. (b) Transient axial profiles of the radial mean temperature and concentration in the bulk phase PH model with Pe,(d,) = 3.85 and Bi, = 1.5 for different I * values: A, 0.1; B, 0.2; C, 0.4: D. 0.6; E. 0.8, F. 1.0: G, 10.0; H, 100.0: I, 400.0; J, 1000.O; L, 2500.0.

radial thermal dispersion (which will give higher Bi, values and/or lower Pe* values). In the cases analyzed, the convenient Bi, and Pe, values can still be included in the normal range of those effective par-

ameters. Thus, apparently, the PH model could t-eplace the full HTd, model after an adequate choice of those radial parameters. However, this is not straightforward and some problems can arise from this

748

R. M. QUINTA FERREIRAet al.

approach. For instance, if we chose the above-mentioned first set of parameter values for the studied case To = T, = 625 K [Biw= 1.5, Pe,(d,) = 6.51 they would lead to a thermal instability situation when used for To = T, = 660 K, as shown in Fig. 11(a). Therefore, in spite of being able to find, for each particular situation, a convenient set of lumped parameters for the PH model which could fit the global results of the more complex HT.,, model, one must be very cautious when doing so. This is due to the fact that in our system, the PH model shows an extremely high sensitivity to such radial parameters, leading to an abrupt increase in the hot spot magnitude for relatively low P&clet numbers. For the HT models our results show a more gradual effect of the thermal radial dispersion on the maximum temperature of the process, as seen in Figs 12(a-e) and 13(a-e), which show the effect of Pe,, , Per,,,and Bi, on the process parameters analyzed above: the hot spot magnitude, AT,,,, and the final reactant conversion, xfr when intraparticle convection is taken into account. This study has also been situations: for two steady-state carried out To = T, = 625 K and T, = T, = 660 K. The system is more sensitive to changes in the thermal parameters for the fluid phase (Pef,,, Sic) than to changes in those associated with the catalyst particles (Pe:,,, Bi&) as seen in Figs 12(a-d) and 13(a-d), by comparing parts (a) with (b) and parts (c) with (d) (the situations corresponding to the parameter values used in the dynamic simulations presented above are indicated by solid circles). The effect of the radial mass dispersion can also be neglected in the HT models, Figs 12(e) and 13(e). These results are in agreement with the work of Odendaal et al. (1987), who examined the relative importance of the parameters associated with the 2D HT model of packedbed reactors where the highly exothermic naphthalene oxidation occurred. They concluded that the system was insensitive with respect to the radial mass dispersion and that there was a significant dependence on Pe,, and wall heat transfer coefficients. Windes et al. (1989) also pointed out that the heat transfer within the bed is the most important thermal transport mechanism, in particular the radial heat dispersion. They also concluded that the radial mass dispersion could be neglected. For the steady-state regime, when To = T, = 625 K, Fig. 14(a) shows the comparison between the radial mean temperature and concentration profiles obtained with the complete HT,, model (curves A) with those predicted by the PH model when the previously chosen set of parameter values are used [curves B with Bi, = 1.5 and Pe,,(d,) = 6.5; and curves C with Bi, = 3 and

Pe,,,(d,) = lo]. One can observe quite good agreement between these steady-state results. In order to see the sensitivity of the PH responses, we have also included the profiles obtained with Bi, = 1.5 and Pe,, @,) = 8 (curves D) and those corresponding to the initial set of PH parameters, Bi, = 1.5 and Pe,(dp) = 10.5 (curves E), which indeed deviate from the correct ones. For the same operating conditions and using the parameter values corresponding to curve B [Biw= 1.5 and Pe,,(d,,) = 6.51, the startup transient profiles for the simple PH model are shown in Fig. 14(b). These results must now be compared with those shown in Fig. l(a), for the HT,, model, which predicts slower transient responses, i.e. for the same t* values, the PH model indicates higher temperatures. For To = T, = 660 K and when the process is running in steady-state conditions, Fig. 15(a) shows the HT.,, results (curves A) and the PH profiles, with Bi, = 1.5 and Pe,,(d,) = 3.9 (curves B) and those obtained with Biw= 3 and Pe,(d,) = 5.8 (curves C). In spite of the relatively close values for ATmm and x, in the three cases (this has been the basis for the selection of these sets of PH parameter values), the internal temperature and concentration gradients show some differences, in particular, the catalytic-bed zone with high temperatures is larger when predicted by the HT,, model. The system is now operating with more severe thermal conditions, leading to greater differences between the processes which occur in the fluid and solid phases. The high thermal parameter sensitivity of the PH model can be observed if we maintain the Bi, value corresponding to curve B and slightly increase the Pe,, value: the temperature runaway phenomenon will then take place as shown in curve D, obtained with Bi, = 1.5 and Pe,(d,) = 4. The parameter values calculated through the correlations used in the transient analysis are still more inadequate than the representation of the system behavior by the simple PH model (curves E). The startup PH dynamic profiles obtained with Bi, = 1.5 and Pe,(d,) = 3.85 are shown in Fig. 15(b). The hot spot magnitude of the final steady-state (for t* = 2500) is now lower than that calculated with the steady-state equations [Fig. 15(a)]. Also, the dynamic results indicate faster transient responses when compared with those obtained for the same operating conditions using the complete HT,, model-Fig. 2(a). These results are consistent with the analysis presented by Windes et al. (1989), which found a good qualitative agreement between the PH and HT models for steady-state, in particular for moderate conditions; however, the dynamic responses could be very different. Ramkrishna and Arce (1989) considered the rational use of a common PH model for

749

A 2D heterogeneous diffusion/convection model

first-order isothermal reaction systems and their results also point out that this simpler model is not adequate to describe the transient behavior. It is our belief that the HT models should still be used for a more thorough analysis of the dynamic behavior of highly exothermic catalytic reactors operating in severe conditions, since the fluid and solid phases have mass and heat gradients that are very different. However, we think that the equivalence of the mathematical models can be investigated further in order to explore all the possibilities of using the simpler models which, neglecting the resistances in the solid phase, will still allow a convenient representation of the behavior of such catalytic systems; the CPU times involved in the numerical calculations will then be reduced substantially; therefore these models, will be more adequate for process control. CONCLUSIONS

The analysis of the influence of the intraparticle convection on the transient behavior of fixed-bed catalytic reactors has been carried out by considering radial gradients in three reactor models: a heterogeneous model which takes into account intraparticle convection and diffusion inside the catalyst (HT,,); a heterogeneous model where intraparticle diffusion is the only mechanism of transport (HT,) and the pseudo-homogeneous model (PH). The simulation work was based on the oxidation of o-xylene to phthalic anhydride, which is a highly exothermic reaction-wherein radial gradients are important. By adjusting the lumped parameters of the PH model, it is possible to fit the steady-state profiles to those predicted by the complete HT,, model. However, the transient startup responses show some significant differences. For large-pore catalysts the intraparticle convective flow contributes to the increase in final reactant conversion; it is then possible to reduce the wall/feed temperatures, which is of interest for largeproduction processes in order to save energy costs. In any case, results from the HT,, model lie between those from the HT, and PH models; in fact, the effect of intraparticular convection is to “augment” the effective diffusivity, and so steady-state conversions from HT,, are between “diffusion’‘-controlled and kinetic-controlled cases. Even for the complete HT,, model, the results favored the use of models which neglect the accumulation term inside the particles, which have lower computing times. The maximum temperature rises inside the catalytic bed predicted by these 2D HT models were lower than

those

achieved

previously

with

1D

models

(Quinta Ferreira and Rodrigues, reverse is true for PH models.

1991). However,

Acknowledgement-Financial support from JNICT is gratefully acknowledged.

INIC

the

and

NOMENCLATURE

Bi, = Bi, = C = C,,,,, = C,, = C, = clr = cil= Da = D, = D, = dp d, E f

= = = =

tB = (-AH)G/~ow,r1 Mass Biot number (k,R,lD,)

Wall heat Biot number (Bi, = h,&jl,,) Reactant concentration, mol/m-’ Mean radial concentration in the fluid phase, mol/m” Reactant concentration in the feed, mol/m” Reactant concentration at the catalyst surface, mol/m3 Heat canacitv of the fluid, Cal/kg K Heat capacit; of the catalyst, ca%kg K Damkiihler number [Da = k(TO)r/c] Effective diffusivity, m2/s Effective diffusivity in the reactor radial coordinate, m2/s Equivalent diameter of the particle, m Reactor diameter, m Activation energy, cal/mol Dimensionless reactant concentration

(f = C/Co) f= Approximated -AH

=

h= h, = h, = hR,,

n.,

=

k= k, = L = Nd = N, = N,,, =

dimensionless reactant concentration Heat of reaction, cal/mol Film heat transfer coefficient, Cal/m2 s K Length of the finite elements for the space coordinate inside the oarticle Wall heat transfer co&icient, Cal/m’ s K Length of the finite elements for the reactor radial coordinate Reaction kinetic constant, s-’ Film mass transfer coefficient, m/s Reactor leneth M Number of pore diffusion units (Nd = D,t /Q.R.E) urnts film mass transfer Number of [Nf = (1 - L)kyCIpr.16] heat transfer units Number of film

[Nn, = (1 - ~W,=l~~rc~rl

N. = Number of sections in the reactor NE = Number of finite elements inside the particle NER = Number of finite elements in the reactor radial coordinate NINT= Number of subintervals in a section PO= c-Xylene partial pressure at the reactor entrance. atm Per,,= Radial heat P&let number(Pe, = Lu,p,~,~/l,,) Pe,(d,) = Radial heat P&let number based on particle r diameter [Pe,,(d,) = dpuOpfcpf/A_] Pe, = Radial mass P&let number (Pe, = LqJD,) Pe,(d,) = Radial mass P&let number based on particle diameter [Pe, (4) = dp t+,/D.,] Pr = Prandtl number R = Ideal gas constant & = Reactor radius, m R, = Half thickness of the slab catalyst Re = Revnolds number _ r = Reactor radial coordinate, m r * = Dimensionless reactor radial coordinate

CT* = rl%l

r = Spatial coordinate in the particle, m P = Dimensionless particle coordinate (r,’ = rr,/2&.) rp rbinr = Nodes of the finite elements in the reactor radial coordinate

R. M. QUINTA FERREIRA et al.

750

rint = nodes of the finite elements in the particle spatial coordinate T - Absolute temperature, K Tbm = Mean radial temperature in the fluid phase, K T0 = Feed temperature, K T, = Surface temperature, K T, = Wall temperature, K AT,., = Hot spot magnitude, K t=Time, s I+ = Dimensionless time (t* = t/r) (I = Wall heat ransfer global coefficient, Cal/m2s K u = Normalized coordinate inside a finite element for the particle (u = r: - rint(k)/h,) ui = Interstitial velocity, m/s. u, = Superficial velocity, m/s uR = Normalized reactor radial coordinate inside a finite element (uR = r* - rbint(kR)/hR,,) Do = Intraparticle convective velocity, m/s x, = Final reactant conversion, % L = Reactor axial coordinate, m reactor axial coordinate .z* = Dimensionless (I* = z/t) Greek

symbols a’ = /3= y= t = cp = q= 8= ff = & = A,, = 1, = ob = pr = ps = pSp =

Coefficients defined in Table 6 Prater thermicity factor p = (-AU)C,D,/I,T,] Arrhenius number (y = E/RT,) Bulk porosity Catalyst porosity Effectiveness factor referred to the catalyst surface Dimensionless temperature (6 = T/T,) Approximated dimensionless temperature Effective thermal conductivity, Cal/m s K Radial effective thermal conductivity, cal/m s K Mass intraparticle P&et number (1, = u,,R,/D,) Bulk density, kg/m3 Fluid density, kg/m3 Catalyst density, kg/m3 Average heat capacity of the bulk

5 = 5% :%zlt;e%i?& = L/a.) Space time for the thermal waie for the HT models [r,, = (1 - ti)p,c,r/~p~c,,] rhl = Space time for the thermal wave for the PH model (r,,, = &,s /c~~Q,,) T,, =

@,=Thiele

modulus [&=RpJz]

Superscripts f = Fluid s = Solid Subscripts b d dc p s w

= = = = = =

Bulk conditions in the fluid phase Diffusion Diffusion/convection Particle Particle surface Wall

REFERENCES Carey G. and B. Finlayson, Orthogonal collocation on finite elements. C/rem. Engng Sci. 30, 587-596 (1975). Carta G., Personal communication, Univ. of Virginia, (1988). Castro J. A., Aspects of modelling chemical processes for adaptative control. Ph.D. Thesis, Univ. of Leeds, U.K. (1983).

Cheng S. and A. Zoulalian, Influence de la &action chimique et de la convection interne sur la determination de la diffusivite apparente d’un catalyseur. Chem. Engng J. 41, 91-104 (1989). Cogan R., G. Pipko and A. Nir, Simultaneous intraparticle forced convection, diffusion and reaction in porous catalyst-111. Chem. Engng Sci. 37, 147-151 (1982). Cresswell D., Intra-particle convection: its measurement and effect on catalyst activity and selectivity. Appl. Cutalysis 15, 103-116 (1985). Crider J. and A. Foss, Computational studies of transients in packed tubular chemical reactors. AIChE J1 lZ(3). 514-522 (1966). Crider J. and A. Foss, An analytical solution for the dynamics of a packed adiabatic chemical reactor. AZChE J/ 14(l), 7784 (1968). College R. and W. Patterson, Heat transfer at the wall of a packed bed: a j-factor analogy established. In Proc. flth A. Mtg on Heat Transfer and Catalysis and Catalytic Reactors, Inst. Chem. Engrs, p. 103 (1984). De Boor C., Package for computing with El-splines, SIAM JI namer. Analysis 14, 44-472 (1977). De Boor C., A Practical Guide to Splines. Springer-Verlag, New York (1978). De Wasch A. and G. Froment, A two dimensional heterogeneous model for fixed-bed catalytic reactors. Chem. Engng Sci. 26, 629-634 (1971). Dixon A., Thermal resistance models of packed-bed effective heat transfer parameters. AIChE Jl 31(5), 826-834 (1985). Dixon A. and D. Cresswell, Theoretical prediction of effective heat transfer parameters in packed beds. AIChE JI 25(4), 663-676 (1979). Dogu G., A. Pekediz and T. Dogu, Dynamic analysis of viscous flow and diffusion in porous solids. AIChE JI 35(8), 1370-1375 (1989). Ferguson N. and B. Finlayson, Transient modeling of a catalytic converter to reduce nitric oxide in automobile exhaust. AIChE J/ 20(3), 539-550 (1974). Finlayson B. A., Nonlinear Analysis in Chemical Engin wring. McGraw-Hill, New York (1980). Froment G., Fixed bed catalytic reactors. Ind. Engng Chem. 59(2), 18-27 (1967). Froment G. and K. Bischoff, Chemical Reactor Analysis and Design. Wiley, New York (1979). Gunn D., Theory of axial and radial dispersion in packed beds. Trans. Znstn them. Engrs 47, 351-359 (1969). Handley D. and P. Heggs, Momentum and heat transfer mechanisms in regular shaped packings. Trans. fnstn them. Engrs 46, 253-264 (1968). Hansen K. and S. Jorgensen, Dynamic modelling of a gas phase catalytic fixed bed reactor-I. Chem. Engng Sci. 31, 579-586 (1976a). Hansen K. and S. Joraensen. Dvnamic modellinn of a pas phase catalytic fix&bed reacior-II. Chem. &tgng &i. 31, 587-598 (1976b). Hindmarsh A., Preliminary documentation of GEARIB: Report UCID-30 130, Lawrence Livermore Lab., Livermore, CA (1976). Hlavacek V. and P. Von Rompay, Current problems of multiplicity, stability and sensitivity of states in chemically reacting systems. Chem. Engng Sci. 36, 1587-1597 (1981). Ho&erg J., B. Lyche and A. Foss, Experimental evaluation of dynamic models for a fixed-bed catalytic reactor. AIChE Ji 17(6), 1434-1447 (1971). Holton R. and D. Trimm, Mathematical modelling of heterogeneous-homogeneous reactions. Chem. Engng Sci. 31, 549-561 (1976). Horvarth C., Preparative HPLC. Lecture given at the iVA TO ASI on Chromntographic and Membrane Processes in Biotechnology, Azores (1990).

A 2D heterogeneous

diffusion/convection

Jutan A., J. Tremblay, J. Mac Gregor and J. Wright, Multivariable computer control of a butane hydrogenolysis reactor I-State space reactor modelling. AIChE JI B(5), 732-742 (1977). Kunii D. and J. Smith, Heat transfer characteristics of porous rocks. AIChE Jl a(l), 71-78 (1960). Lapidus L. and G. Finder, Numerical Solution of Partial Dtfirentiai Equations in Science and Engineering. Wiley, New York (1982). Lee R. and J. Agnew, Model studies for a vinyl chloride tubular reactor I. Ind. Engng Chem. Process Des. Dev. M(4), 490-495 (1977a). Lee R. and J. Agnew, Model studies for a vinyl chloride tubular reactor II. Znd. Engng Chem. Process Des. Dev. 16(4), 495-501 (1977b). Lopez-Isunza H., Steady-state and dynamic behaviour of an industrial fixed-bed catalytic reactor. Ph.D. Thesis, Univ. of London (1983). Madsen N. and R. Sincovec, PDECOL: general collocation on finite elements. Chem. Engng Sci. 30, 587-596 (1975). Martinez O., S. Pereira Duarte and N. Lemcoff, Modeling of fixed bed catalytic reactors. Computers them. Engng 9, 535-545 (1985). Mascazzini N. and G. Barreto, Development of overall heat transfer coefficients from the heterogeneous two dimensional (2-D) model for catalytic fixed bed reactors. A comparison between 1-D and 2-D models. Chem. Engng Commun. 80, 11-32 (1989). Melanson M. and A. Dixon, Solid conduction in low d,/d, beds of spheres pellets and rings. Heal Mass Transfer 28(2),

383-394

(1985).

MC Greavy C. and H. Naim, Reduced dynamic model of a fixed-bed reactor. Can. J. them. Engng 55, 326-332 (1977).

Nir A. and L. Pismen, Simultaneous intraparticle forced convection, diffusion and reaction in a porous catalyst I. Chem. Engng Sci. 32, 35-41 (1977). Nir A., Simultaneous intraparticle forced convection, diffusion and reaction in a porous catalyst-II. Chem. Engng Sci. 32, 925-930 (1977). Odendaal W., W. Gobie and J. J. Carberry, Thermal parameter sensitivity in the simulation of the non-isothermal, non-adiabatic fixed-bed catalytic reactor-the twodimensional heterogeneous model. Chem. Engng. Commun.

58, 37-62

(1987).

Olbricht W. and 0. Potter, Mass transfer from the wall in small diameter packed beds. Chem. Engng Sci. 27, 1733-1743 (1972). Pereira Duarte S., G. Barreto and N. Lemcoff, Comparison of two-dimensional models for a fixed-bed catalytic reactor. Chem. Engng Sci. 39, 1017-1024 (1984a).

model

751

Pereira Duarte S., 0. Ferreti and N. Lemcoff. A heterogeneous one-dimensional model for non-adiabatic fixedbed catalytic reactors. Chem. Engng. Sci. 39, 1025-1031 (1984b). Quinta Ferreira R., ContribuiPgo para o estudo de reactores cataliticos de leito fixo: efeito da convec@o em catalisadores de poros largos e cases de catalisadores bidisperSOS. Ph.D. Thesis, Univ. of Porto, Portugal (1988). Quinta Ferreira R. and A. Rodrigues, E&ct of intraparticle convection on the transient behavior of fixed-bed mactors: unidimensional models. Submitted (1991). Ramkrishna D. and P. Arce, Can pseudo-homogeneous reactor models be valid? Chem. Engng Sri 44, 1949-1966 (1989). Rodrigues A., B. Ahn and A. Zoulalian, Intraparticle-forced convection effect in catalyst diffusivity measurements and reactor design. AZChE JI 28(4), 541-546 (1982). Rodrigues A., J. Orfao and A. Zoulalian, Intraparticle convection, diffusion and zero order reaction in porous catalysts. Chem. Engng Commun. 27, 327-337 (1984). Rodrigues A., J. Loureiro and R. Ferreira, Intraparticle convection revisited. Chem. Engng Commun. 107, 21-33 (1991). Rodrigues A., L. Zuping and J. M. Loureiro, Residence time distribution of inert and linearly adsorbed species in a fixed bed containing “large-pore” supports: applications in seoaration eneineerine. Chem. Emma Sci. 46, -2765-2773 (1991). Rodrigues A. and R. Quinta Ferreira, Convection, diffusion and-reaction in a large-pore catalyst particle. AZChE Symp.

Ser. 84, 80-87

(1988).

Rodrigues A. and R. Quinta Ferreira, Effect of intraparticle convection on the steady-state behavior of fixed-bed catalytic reactors. Chem. Engng Sci. 45,2653-2660 (1990). Sharma C. and R. Hughes, The behaviour of an adiabatic fixed bed reactor for the oxidation of carbon monoxide II. Chem. Engng Sci. 34, 625-634 (1979). Sinai J. and A. Foss, Experimental and computational studies of the dynamics of a fixed bed chemical reactor. AIC’hE J/ 16(4), 658-669 (1970). Stephanopoulos G. and K. Tsiveriotis, The effect of intraparticle convection on nutrient transport in porous biological pellets. Chem. Engng Sci. 44, 2031-2039 (1989). Von Doesburg H. and W. De Jong, Transient behaviour of an adiabatic fixed-bed methanator I. Chem. Engng. Sci. 31, 45551 (1976). Windes L., M. Schwedock and Harmon Ray, Steady-state and dynamic modelling of a packed bed reactor for the partial oxidation of methanol to formaldehyde. I. Model development. Chem. Engng Commun 78, l-43 (1989). Yagi S. and D. Kunii, Studies on effective thermal conductivities in packed beds. AIChE JI 3(3), 373-381 (1957).