Generalized Criteria for the Stability of Catalytic Reactors C. McGREAVY and J . M . THORNTON Department of Chemical Engineering, Leeds University, England

It is well known that particles catalyzing highly exothermic reactions can have multiple steady states, which are characteristic of unstable behavior. Previous work has shown that under practical conditions the pellets are essentially isothermal, all the temperature rise being across the surrounding fluid film. A technique has been developed for such cases which enables the regions of potential local instability, if present, to be defined for arbitrary operating conditions. It is thus possible to monitor the calculation of concentration and temperature profiles in a fixed bed reactor to ascertain whether the conditions would represent multiple states and therefore non-feasible operating conditions.

I

C’est un fait bien connu que les particules qui catalysent des rkactions fortement exothermiques peuvent avoir de multiples regimes permanents, lesquels caractkrisent un comportement instable. On a demontre dans les travaux prkckdents qu’en pratique les pastilles sont essentiellement isothermiques et que l’klkvation de tempkrature se fait sur la pellicule du fluide environnant. On a mis au point une mkthode relative des cas de ce genre, laquelle permet de dkfinir, pour des conditions arbitraires d’opkration, les regions 0; une instabilitk locale peut se produire, le cas kchkant. On peut ainsi contraler le calcul des profils de concentration et de temperature dans un lit fix6 pour dkterminer si lee conditions qui y existent representent des regimes multiples oii la marche des operations est impossible.

t is well known that the solutions of the equations for the

catalyst pellets supporting highly exothermic reactions can exhibit multiple solutions for certain conditions. Since such regions indicate instability, considerable attention has been given to the analysis of such problems as a first step in a more general and rigorous study‘l). . Attempts to develop criteria for establishing the regions of non-unique solutions have met with only limited success and no satisfactory treatment appears to exist relating local and global instability in a fixed bed reactor. Primarily, the greatest source of difficulty has been the fact that the usual form of effectiveness factor chart expressed in terms of the Thiele modulus 4 as a function of the thermicity factor /3 and activation factor y only applies to a single pellet under one set of conditions. McGreavy and Cresswcll(2),however, have shown that by taking account of the finite rates of interphase heat and mass transfer, and hence rewriting the boundary conditions as a Neumann rather than a Dirichlet problem, considerable simplifications may be made. Primarily, this is because of the restricted range of values that the characteristic dimensionless groups can take in a physically realizable system. I n particular, the relatively large thermal film resistance at the surface of the catalyst and low solid diffusion coefficient generate large temperature differences across the film but result in the solid being essentially isothermal. Typical effectiveness factor curves become a function of modified Sherwood and Nusselt groups, as shown in Figure 1, where the main regions have been identified. Multiple solutions exist between the values of Thiele modulus given by 41and &, but not every set of conditions produces a curve with multiple solutions. Further information is required on the permissible magnitude of perturbations which would not cause the pellet to change its state. T o analyse a whole reactor in detail an infinite set of such charts would be required i.e. one for each point in the bed. This is because /3 and y vary throughout the reactor. T h e problem must be reformulated to enable a single graph to be used for the whole reactor. Preferably, the effectiveness factor

curve should be reduced to a single point in the new representation, and the regions of multiple solutions i.e. potential instability, identified. I t should then bc possible to represent a trajectory through the reactor on the diagram to see if multiple solution conditions will be generated.

Formulation of the equations For spherical geometry, the heat and mass balances in the pellet are described by the equations: d2C, 2 dC, + -ds2 s ds ~

k -c,

D

= 0

-d +2 T- ,- + (2- dATH , )-Ct=O k ds2 s ds K

where

k = A exp (- E/RT,)

In dimensionless form these may be written d2c -

dxZ d2t

dx2

dc + x2 dx -

dt + x2 d+ x -

e2

(+ (:) (- :-> 82

exp

exp

0. . . . . . . . . .

c =

0. . . . . . . . (2)

where

The Canadian Journal of Chemical Engineering, Vol. 48, April, 1970

187

Rather than using the Thiele modulus I$ which would be a function of the position in the bed, a ncw parameter 0 is used, based on the pre-exponential factor and not on the rate constant. It is therefore independent of position in the bed and characteriscs the reactor system. Equations (1) and (2) are subject to the boundary conditions:

t

7)

-dc_-- dt dx dx Sh' --

2

Figure 1-A typical effectiveness factor chart showing the main controlling regions: (1) Kinetic control; ( 2 ) Diffusion control; ( 3 ) Unstable steady states; ( 4 ) Interphase mass transfer control.

=o

"1

atx=O

(Iwc) = dx } n t x =

1

Nu' -(T-t) = dx

2

where

Equations ( 1 ) and (2) together with their boundary conditions would normally be solved numerically. It has bccn found that, for the whole practical range of heat transfer cocfficients and thermal conductivities, the pellet is essentially isothermal - i.e. the temperature rise is almost exclusively across the fluid film surrounding the catalyst(*). Equation ( 2 ) may therefore be replaced by an overall heat balance on thc pcllct

-02'i--l

which in dimensionless form becomes t - T = B Sh' (1

- 04

0

02

M

06

OB

- ~ ~ - 1 ) .. . . . . . . . . . . . . ( 3 )

where .I0

.I2

14-

16

I?= (-AH)C,DK = - P

I8

Nu'y

2bhE

-t

Figure 2-The relationship between T and t obtained solving equation (4) for typical values of 8 and B.

R is a convenient combination of several dimensionless groups, and enables a reactor system to be represented on a single chart, as will be shown later. For a given system, B depends only on the concentration in the fluid and can be regarded as equivalent to reactant concentration. It is, therefore, an implicit function of axial position. Since t is now constant throughout the pellet, Fquation (1) is linear and can be solved analytically, after making the substitution 2; = xc, giving: Sh'

02

10-5

Figure 3-The

unique and non-unique regions predicted by equation (5).

whcrc r = eexp(-

and 1

Sh'- 500

k)

g = tanh r

Substituting this in ( 3 ) givcs

-8

Figure &The

I88

effect of 8 on the non-unique region.

This is the condition for a steady state to exist. t may t)e obtained by any of the normal root finding techniques. Figure 2 shows graphically a plot of the right hand side of Equation (4) for three typical sets of reaction parameters. As would be expected, a region of three steady states may exist in some cases, the middle one, of course, being unstable.

The Canadian Jnurnal of Chemical Engineering, Vol. 48, April, 1970

This equation is easily solved using a Newton-Raphson procedure.

B -104 Y

Local stability

-discusdon

of

the rerults

T h e general form of the solution to Equation (5) is shown in Figure 3 . T h e lines ( I ) , (2) and ( 3 ) correspond to the curves in Figure 2. In case ( I ) , starting at A , the pellet temperature will slowly increase as T increases until the line X-Y is crossed, when the pellet temperature will jump to its maximum possible value. This is the region of mass transfer control between the bulk fluid and the catalyst surface, where t

fi

-B

effect of S h on the non-unique region.

Figure 5-The

14

-8

limits of non-unique solutions for various values of S h and 8.

Figure &The

h

I k

10-5

I

I

t

1

1

10-4

1

1

1

1

10-3

6

upper curves for Sh’ of 8.

Figure 7-The

=

M

at various values

Calculation of the bou.nds on the non-unique region For a given system B, .%’and 0 are constant. If non-unique solutions exist, then the bounds on T may be found by solving the equation dT _

dt

--

0 to give local extrema (see Figure 2)

Differentiating (4) 1

Sh’* - B--

t*

+g*

- I’

+ rg

= 0.. . . . . . . . ( 5 )

T

+ B Sh‘.

,

....

,

. . , . . , . , . . . .(6)

If the temperature T is now reduced, Equation ( 6 ) will apply until the line X Z is crossed. The pellet will again be in a region of unique solutions. In case (2), it is not possible to return to the initial steady state by reducing only the temperature, while case (3) is unique at all temperatures. Figures 4 and 5 show the effect of varying 0 and Sh’ respectively. For a given reactor, 0 is constant and Sh’ can only vary between fairly narrow limits, since the mass transfer coefficient k, (and hence Sh’) is proportional to the square root of the flow rate ( 3 ) . A set of curves similar to those in Figure 5 will cover all possible operating conditions, but, because of this small variation in Sh’, one curve may often be sufficient, especially as the upper arm of the curve is the most important and is relatively insensitive to changes in Sh’. T h e insensitivity of the upper curves of Figure 5 to changes in Sh’ was expected since any point on an upper curve corresponds to a position such as that of +2 in Figure 1. This point is in a region of intraparticle diffusion control and provided the interphase mass transfer coefficient is not low it will not be expected to have much influence. Figure 6 is formed by plotring lines similar to the dotted lines in Figures 4 and 5 for various values of 0 and Sh’. This graph gives an immediate indication of a h i i t on the concentration (i.e. on B ) bclow which non-unique solutions cannot exist at any temperature. It is interesting to note that for almost the whole practical range of B, no multiple solutions can exist for T > 0.125 i.e. y < 8. In practice it is nor only desirable to avoid regions of multiple solutions, but also regions of high sensitivity. In Figure 1, for instance, curve (3) corresponds to a point just beyond the cusp (see also Figure 3), where, between V and W ,the pellet temperature t increases very rapidly for small changes in the fluid temperature T . In Figure (S), it was found that beyond the cusp for any particular value of Sh’, the curve for Sh’ = 03 continued to predict regions of high sensitivity. It may well be desirable in practice, therefore, to use the curves for Sh’ = 00 in all cases, as shown in Figure 7 , since these enable all regions of potential difficulty to be avoided. In effectiveness factor charts, the bounds on 4 may be determined by taking points at constant 7 from graphs similar to Figure 4. A typical set of curves is plorted in Figure 8 for Sh’ = 500.

Accuracy of the method The results obtained from Figure 8 may be used tn check the accuracy of the proposed method, since 41 and 42 (see Figure 1) may also be found by solving Equations (1) and (2) by a suitable finite difference method and plotting 9 as a function of 4. Great care is needed in the numerical solution to ensure that the values of effectiveness factors are accurate. It is parricularly important that a sufficiently fine finite difference network be used since very steep gradients in the pellet are common. For example, in region (4) the concentration often falls to zero in the outer 1% of the radius. A convenient way of dealing with this problem

The Canadian Journal of Chemical Engineering, Vol. 48, April, 1970

189

TABLE1 COMPARISON OF APPROXIMATE A N D RIGOROUS VALUES OF ROUNDS OF T H E NON-UNIQUE REGION (Sh' = 500)

P

y

Nu'

B

10 10 20 20

1 0.2 1 10

0.005 0.005 0.0005 0.0001

___ 0.05 0.01 0.01 0.02

_ - - ~ _ _ _ - ~ -

__

I

0.81 0.81

0.82 0.81 0.19 9.0

0.20 8.9

1.32 1.32 2.45 12.5

SUMMARY OF PHYSICAL

1.31

2.45 11.8 ____

--

~

1.26

TABLE2 DATAFOR

T H E REACTOR SYSTEM

&#-

Figure 8-The bounds on q5 for non-uniqueness of effectiveness factor charts ( S h = 500).

Tube radius Pellet radius Bed porosity Pellet porosity Tube length Peclet number for radial heat transfer Peclet number for radial mass transfer

Fluid density Superficial velocity Specific heat of fluid Overall wall heat transfer coefficient Pellet to fluid heat transfer

10 10

6.83 X g/cm3 65.7 cm/sec 0.585 cal/g OK 1.05 x 10-2cal/cm2 see O

K

coefficient

1.88 x 10-2cal/cm2sec "K

coefficient

4.36 cm/sec

pellet

3.66 X cm*/sec 500 K cal/g mole 26.6 K cal/g mole 8.29 X 10losec-l 4.45 x 10-8 g.molt:s/cms 500 "K

Fluid to pellet mass transfer Effective radial diffusivity in -

2.1 cm 0.21 cm 0.4 0.4 125 cm

Heat of reaction Activation energy Arrhenius pre-exponential factor Inlet concentration Inlet temperature

and 42can be predicted to better than 6% for all values of Sh'. However, the accuracy of predicted fluid temperatures leading to multiple solutions is much better than this since $ is exponentially dependent on ?'.

Local and global stability

Nofwnipor

I

Analyses of multiple solutions in connection with stability in tubular reactor systems have tended to dcal with either the quasi-homogeneous (adiabatic) reactor or with the behavior of single particles. In the heterogeneous reactor, there is an interaction between the two which inevitably restricts the degrees of freedom in specifying the state variables, and which may tend to limit the development of instabilities. This interaction is particularly difficult to investigate if the usual dimensionless groups for the pellet - 4, @ and 7, are used, since each point must be represented on a different chart. Consider the use of the modified groups proposed here, which arise out of the simplifications possible in real systems (*). A plot of the fluid temperature, T , against the group B is characteristic of a ser of operating conditions for a given system as shown in Figure 3, where the region of non-unique solutions is indicated. If, for example, equations describing the heterogeneous, two dimensional catalytic reactor are solved, it is possible to plot longitudinal trajectories for particular radial positions on the same chart. Only if any of the curves pass through the multiple solution region will the reactor tend to have multiple solutions in some region, and therefore be potcntially unstable.

m . u

.02

I

2

3

4

5

6

7

8

9

I0

ax&

Figure 10-Profiles of T and B at various radial and longitudinal positions for a coolant temperature of 488°K.

is to use uncvenly spaced grid points, the majority being in the region x = 0.99 to 1.O. Failure to recognize rhis seems to have caused snme errors in published charts. I n region (2) the gradients are not so steep and a uniform grid can be used. A comparison between the valucs of $1 and 42 calculated by both techniques is given in Table 1. It was found that if N7t < 20, which covers all cases likely to arise in practice, then 41

190

The Canadian Journal

of

Chemical Engineering, Vol. 48, April, 1970

Typical axial trajectories for the reactor system defined by Table 2 are shown in Figure 9. T h e influence of coolant temperature is demonstrated, values greater than about 480% indicating possible instability. Figure 10 shows longitudinal trajectories for various radial positions for a coolant temperature of 488'K. No complete radial profile lies in the multiple solution region, so it is possible that instabilities will be damped down. Besides indicating where a reactor will tend to be unstable by virtue of trajectories passing through the non-unique region, it is possible to obtain some idea of how the reactor will behave outside, but close to this region. In the coiirsc of numerical solutions, the criteria provided by these charts makes it fairly simple to assess whether undesirable operating conditions have been presented. Where the reactor model is part of an optimization procedure, this is very useful.

Conclusions

1

= oexp(-k)

I'

=

R R'

= =

S

=

Sh'

=

1

= dimensionless temperature

T

=

TI T,

=

X

=

.

=

I

distance from tube axis gas constant tube radius distance from pellet centre 2bk, Sherwood number D R Ts E RT, dimensionless temperature E fluid temperature solid temperature dimensionless radial coordinate s / b

Greek Letters

A method of analysing multiple solutions in reactor stability problcnis has been proposed which enables local and global stability to be related on a simple diagram. It also indicates when thc reactor is in the region of potential instability. This information is useful in design, optimization and control studies. Nomenclature pre-exponential factor in Arrhenius rate expression pellet radius ( - A H ) CfDR 2bhE dimensionless concentration C,/C, reactant concentration in fluid reactant concentration in solid effective radial diffusivity inside the pellet activation energy tanh(r) heat transfer coefficient a t pellet surface heat of reaction rate constant per unit volume of catalyst mass transfer coefficient a t pellet surface thermal conductivity of pellet 2bh Nusselt number K

4

=

b.&

e

=

b dD d

References ( 1 ) Aris, R., Chem. Eng. Sci. 24, 149 (1969). (2) McGreavy, C. and Cresswell, D. L., Can. J. Chem. Eng. 47, 583

(1969). (3) Carberry, J. J., AIChE Journal 6, 460 (1960).

Manuscript received September 29, 1969; accepted January 5, 1970.

The Canadian Journal of Chemical Engineering, Vol. 48, April, 1970

* * *