Use of the Heaviside function as a substitute for Arrhenius' law in mathematical models of distributed parameter reactors

The Chemical Engineering Journal, 9 (1975) 207-213 0 Elsevier Sequoia +%A.,Lausanne. Printed in the Netherlands

Use of the Heaviside Function as a Substitute for Arrhenius’,Law in Mathematical Models of Distributed Parameter Reactors Part I: Zeroth Order Reaction ITALO H. FARINA Departamento

de Ingenieria Quimica, Facultad de Ingenieria,

Universidad National de La Plata, 1 esq. 47-Lu Plata (Argentina)

(Received 3 June 1974; in revised form 10 March 1975);

Abstract

the exponential term of Arrhenius’ law by the Heaviside step function. In other words, the function

The possibility of using the ffeaviside step function as a substitute for Arrhenius’ law is investigated in order to avoid the large amount of computation involved in the exponential function. The results for a very simple mathematical model of a distributed parameter reactor are compared when a zeroth order reaction is taking place. The results are found to be very satisfactory although some discrepancies are observed.

rA = ko exp(- E/R W’) f(C) where W’ is the catalyst temperature of concentrations, is replaced by rH =

and C the vector

koH(W’ - TL) f(C)

where 0 H(W’ - T:) = 1

if if

W’< T; W’>T;

and Tf is the critical or ignition temperature. Since a costly computational term is replaced by one which consists of a logical decision and an assignation, it is reasonable to ask how good the Heaviside step function is as an approximation to Arrhenius’ law in realistic models for fixed-bed catalytic reactors. Burgess and Lapidus3 have made an important contribution in the investigation of extrapolation methods that reduce computing times by several orders of magnitude. The combination of both proceduresextrapolation and use of simple functions like Heaviside’s-yields the possibility of implementing distributed parameter models in real time. The objective of the present work is to evaluate rH as an approximation to rA in a very simple model for a catalytic reactor when f(C) E 1, i.e. for a zeroth order reaction. In a second part the analysis will be extended to first order reactions.

1. INTRODUCTION

The present state of development of mathematical models and numerical methods for lumped parameter systems is such that it allows a direct and relatively easy implementation in real time computers. However, the same conclusion is not true for distributed parameter systems described by partial differential equations. A dramatic example is afforded by the mathematical models for heterogeneous catalytic reactors. In these cases the long computer time spent in solving the equations of the models forbids the calculations even by off-line computers. It is therefore of great importance to search for numerical approximations which are exact enough and of short computing time. Aris et al. L 2 have used an idea that is equivalent to replacing the exponential term derived from Arrhenius’ law in order to study the transient behaviour of distributed parameter systems. They state that I: “In looking for principles by which the final steady state may be determined it seems desirable to look first at a much simpler model which retains the characteristic features of the packed bed”. The simplified model, among other things, replaces

2. STATEMENT OF THE MATHEMATICAL MODEL

Our mathematical model consists of a tube through which the reactive fluid flows. The zeroth order reaction -Al + Aa = 0 takes place at the inner wall of 207

I. H. FARINA

208 the tube in such a way that no mass accumulation in the solid takes place. The system is adiabatic and the tluid is in plug flow. Furthermore, we assume that there is no axial conduction or dispersion. With these assumptions the conservation balances are as follows. (a) Mass balance

at the wall:

PF=-r+k,(Cr-C,)=O

(1)

where Q = Pke (-AH); TA is the critical temperature whose meaning is considered below. Introducing these variables into the balance equations we obtain --+-_=

ar aT -+-=W at az

(b) Mass balance A yt;

- k,(Cf

(2)

For a zeroth order reaction we have r=

ifCfZ0

kodw’) (0

ifCf=

We have introduced

I

~ C,y)

(3) 0

where g is a given function of the wall temperature. Combining eqns. (1) and (2) we get

ifX=

1

(7)

(8)

W+g(W)

=T-IV l

fl = k+;:

for the fluid:

= ~ VA g+

Tp

“at

1

T

aw This balance simply states that the number of moles of Ai that disappear by chemical reaction is equal to the number transferred from the fluid to the wall.

ifX<

ifX<

1

ifX=

1

the dimensionless

(9)

parameters

W=p cwAw cfAf

Parameter fl is the inverse of the dimensionless adiabatic temperature rise and w is the relation between the heat capacities of the fluid and the wall. Initial and boundary conditions for eqns. (7)-(9) are T(0, t) = G(r)

W(z, 0) = We(z)

X(0, t) = 0

T(z) 0) = To(z) X(z, 0) = X0(z)

A$=__

VA$_&

(c) Heat balance cwAw $

cfA $

for the wall:

= Ph(T’

(d) Heat balance

(4)

(5)

for the fluid: + Ph( W’ -~ T’)

These equations can be put in dimensionless form by defining the following dimensionless variables’ : T=;(T’-

T;)

hPz’ z =-VAcf

x=

ci ~ Cf Ci

W=h;(W’p - hPr’ * =-

Acf

is that 0
3. THE CRITICAL TEMPERATURE

~ W’) + (-AH)&

= - VcfA ‘2

An auxiliary condition

T;)

A parameter artificially intr.educed into the model is the critical temperature TA. It is defined as that temperature which results in the ignition of the reaction. It is precisely defined when the step function is used as the function between the rate of reaction and the temperature, but its meaning is not clear when the rate of reaction follows Arrhenius’ law. In this case TL is another parameter that must be adjusted in the model. Let us choose Ti such that e = exp(-q/T;)

(10)

where E is analogous to the “boundary layer thickness’ used in fluid dynamics. The choice of E is equivalent to the choice of Ti except that it is easily visualized as a fraction of the heat generation curve. E (or Ti) is a parameter which has to be identified in a real case according to the criteria determined by the identification method to be used.

THE HEAVISIDE FUNCTION AS A SUBSTITUTE FOR ARRHENIUS’ LAW. I

209

In Arrhenius’ law we can write

L*(gw+Tq-l=(y-)-l W'

(11)

where (Y= qhP/Q is a new parameter representing the activation energy. Then Arrhenius’ law is expressed by

a=01 E ze-2

(12) It should be observed that for W’ = O”K, W = cy/ln E. Since it is meaningless to work with temperatures below O”K, we have the following restriction: T,

W$

(13)

This restriction was not considered in the model of Aris. Therefore his model is invalid as an approximation to the case of zero activation energy since in that case Th would be the absolute zero and the only possible solution would be that of an ignited reactor. Furthermore, negative dimensionless temperatures would be negative absolute temperatures.

4. STEADY STATE

If ri is constant Ws - Ts =

the steady state is defined by

d%)

ifX,<

0

ifX,

dTs = W, - Ts

L

dX,_ Pd Ws) ifX,< 0

= 1

Ts(0) = Ti

a.?.

dz-

1

ifX,

(14)

(15)

1 = 1

(16)

Fig. 1. Wall temperature W, corresponding toXs=lasa function of Ti; 01= 0.1, E = e-2.

From eqn. (14) we see that to avoid multiplicity of the steady state it is sufficient that the slope of ge(Ws) at the inflection point be less than one. This occurs if CY > 4 eM2 = 0.54144; otherwise we will have the situation shown in Fig. 2. Thus if CY< 4 eP2 there will be multiple solutions if Ts is in the interval (71, yr). The values of 71 and yr are obtained by solving dge(W)/dW = 1. They are represented as a function of (Yin Fig. 3. In the same figure we show the location of the absolute zero of temperature, i.e. ol/ln E. It can be seen that, for values of a less than 0.36, y1 < cy/ln E. Therefore in these cases we will have multiplicity from the inlet with a greater probability of an ignited or upper solution. If T, > yr the only possible solution is the ignited one. Let W2 be the positive solution of the equation

Hence Xs = P(Ts - Ti)

(17)

Because a zeroth order reaction is considered, total conversion will be obtained at a finite distance z 1. At z1 chemical reaction ceases and the temperatures reach the values Tl and WI . From eqn. (17) 1 =fl(Tl - Ti), and from eqn. (14) WI - T1 =g(Wl). Therefore W1 =$ + Ti +g(Wl)

(18)

The root of this equation gives us T1 and W 1. Some of the values of WI as a function of ri when g sg, are shown in Fig. 1 for illustration.

Fig. 2. Possible steady state solutions in the heat generation curve, eqn. (14).

210

I. H. FARINA

x

s

=

(3ge(Wdz

ifX,<

I

(7-l)

i1

This solution is obtained whenever Ti 2 7, and Ti> Wa - g,(Wl). It corresponds to a type A solution and is represented in Fig. 4.

5. USE OF THE HEAVISIDE

FUNCTION

In order to approximate Arrhenius’ law by the Heaviside step function we must take

6-

(22) \

c

L-1

Thus, in the steady state,

Ws- Ts =ge(W)H(Ws) =

Fig. 3. Variation of the range of non-uniqueness of steady state solutions and the absolute zero as a function of LY.

if W>OandX<

&eWYdW= I.

if WGO

Ifge(W1) or -- ge(Wr), the ignited solution will be ifX,<

T

ws=

ri

(

+ge(wl)(z

+

1)

ifX,

TS

:I/“-

(19)

= 1 ifX,<

ifW
1 = 1

CO)

Fig. 5. Type B and C steady

(24)

X=1 1

ifW
1

_ state solutions.

(23)

1

if W>OandX<

2

Fig. 4. Type A steady state solution.

or X = 1

if W>OandX<

1

S

ifX,

1

or X=

(25)

THE HEAVISIDE FUNCTION AS A SUBSTITUTE FOR ARRHENIUS’ LAW. I

In this case we have yt = -g,( Wr) and yr = 0. Therefore, if Tr > 0 the only possible solution is the ignited one and it is the same as the one given by eqns. (19)(21). If Ts< -ge(W1) the solution will be the quenched one, i.e.

x, =o

(26)

Ws = Ts=Ti

This solution is always obtained for fi < -ge(W1). Within the interval -g,( WI)< Ts< 0 nothing can be said without considering the transient state. In this situation a lot of intermediate solutions are possible, eg. the type B and C solutions represented in Fig. 5. This point has been discussed in refs. 1 and 2.

TYPO

6. TRANSIENT STATE

T(0, r) = Ti (27)

X(0, t) = X(z, 0) = 0 With these conditions, and taking into account the step function, the wall temperature at the inlet will evolve according to the equation

aw,

r)

=

2-i- MO, t) •t g,( W,)H{

boundary

and

W(z,O)= wo

cdJ___at

0°K

Fig. 6. Division of the plane W(0, 0), Ti according to the type of solution that is reached at the steady state.

We will restrict ourselves to constant boundary initial conditions of the type T(z, 0) =

0/

in practice is W(0,0)> 0 for all Tiand r, > 0 for all YO, 0). When the reaction rate is given by Arrhenius’ law it is no longer possible to solve eqn. (28). For this case conditions (27) were used for the numerical solution of the transient state (see Appendix). It was found that the zone of the IV(O,t), Ti plane in which quenched solutions are obtained is very small. Thus the ignited solution is obtained for WO = -0.025 and Ti= -0.044. On making Tt smaller, type B solutions begin to

W(0, t)} (28)

This equation is analogous to eqn. (19) of ref. 1 and using a similar procedure to Aris and Schruben we can see that K(O) = 2-i + ge(W1) will be reached only if T, > 0 or W(0, 0) > 1-i > -g,( WI),while the quenched solution W,(O) = Tr

(29)

0 when (30)

will be obtained if W(0,0) < 0 and Ti< 0.This will always be so for Ti
Fig. 7. Evolution of the transient state for WO = -0.025 and Ti = -0.0499. Exponential function; OL= 0.1; P = 1; e=e-2.

1. H. FARINA

evolution of the profiles is qualitatively the same. Furthermore, the maxima of Wand Tin the steady state are smaller for the case of Arrhenius’ law. This is because the heat dissipation line will intersect the heat generation curve at values less than unity (see Fig. 2).

7. CONCLUSIONS

Kg. 8. Evolution of the transient state for Wo = 0 and ri = 0.01. Step function: 01 = 0.1; (3= 1; E = e-2. emerge and they move towards z = 00 as Ti -+ a/in E = -0.05. Since in practice the reactor has a finite length, it is obvious that it will quench for a value -0.044 > ri > -0.05. The evolution of the transient state when Wo = -0.025 and F, = 0.0499 is shown in Fig. 7. In order to compare solutions, we show the transient state using the step function for Wo = 0 and Ti= 0.01in Fig. 8 and the same state using the exponential in Fig. 9. It is observed that the response of the system using the step function is a little faster than the one using the Arrhenius’ law. However, the

From the preceding study we can conclude that the replacement of Arrhenius’ law by the step function gives fairly acceptable results for the model under consideration. The profiles obtained for the steady state using the step function accurately predict the profiles when the solution is ignited. When the solution is quenched the step function gives Xs = 0, W, = Ts= Ti. In this case the heat generation curve intersects the heat dissipation curve at a lower point. However, the profiles will be practically the same as before since the heat generated is negligible. In the transient state the step function gives a time evolution comparable with the Arrhenius’ law solution if initial and boundary conditions are not in the third quadrant of the W(0, t), Tiplane. The zone of that plane in which a quenched reactor solution is obtained is very small for Arrhenius’ law but covers the whole quadrant for the step function. However, the physically valid region of the third quadrant is very small, and therefore the disagreement is negligible. Type B solutions are possible with Arrhenius’ law but have not been observed with the step function. Obviously they are theoretically obtainable with adequate initial conditions.

NOMENCLATURE

A c C

AH h k, ko P Q 4 Fig. 9. Evolution of the transient state for W,J = 0 and 7’i = 0.01. Exponential function; (Y= 0.1; p = 1; B = e-2.

T

cross section concentration volumetric heat capacity heat of reaction heat transfer coefficient mass transfer coefficient pre-exponential factor in Arrhenius’ law internal reactor perimeter = Pko(-AH) = E/R,the energy of activation divided by the gas constant fluid temperature

THE HEAVISIDE

FUNCTION

AS A SUBSTITUTE

maximum fluid temperature chemical reaction ceases on time coordinate reactor wall temperature maximum wall temperature chemical reaction ceases on conversion length coordinate

Tl t W Wl

X .a

FOR

obtained when total conversion

obtained when total conversion

Greek symbols (Y = WIQ parameter defined in Fig. 2 -f1 parameter defined in Fig. 2 ?i E parameter defined by eqn. (10) Indices f i S W

I

fluid entrance steady state wall dimensioned

1 R. Aris and 0. Schruben, Chem. Eng. J., 2 (1971) 179. 2 I. H. Farina and R. Ark, Chem. Eng. J., 4 (1972) 149. 3 W. P. Burgess and L. Lapidus, Chem. Eng. Commun., I

213

LAW. I

APPENDIX

Numerical solution of eqns. (7), (8) and (9) Let r be a direction along the characteristics of eqns. (7) and (8):

(A.1) dT

aT aT -x+z=W-T

z-

(A.3

We use the net shown in Fig. A. 1. Given the values of T, W and Xat time t for all z we can use eqn. (9) to compute the values of W for t + 6 using a modified Euler or any Runge-Kutta method. Now, using the value W = ( WJ + W1)/2 in the right-hand side of eqns. (A. 1) and (A.2) they can be solved analytically to compute Xs and Ts . Several runs halving the step 6 led to a value of 6 = 0.02, as in ref. 1.

physical value

REFERENCES

(1973).

ARRHENIUS’

‘f t+s

3 4

3

--_$

t

2

I z

Fig. A.l.

Ii8

> z

The net used for the numericalcomputation.