Chaotic behavior in the dynamical system of a continuous stirred tank reactor

Physica 25D (1987) 131-154 North-Holland, Amsterdam

CHAOTIC BEHAVIOR IN THE DYNAMICAL SYSTEM OF A CONTINUOUS STIRRED TANK REACTOR D.G. RETZLOFF, P.C-H. CHAN, C. CHICONE, D. OFFIN and R. MOHAMED Unioersity of Missouri-Columbia, Columbia, 340 65211, USA Received 29 November 1985 Revised manuscript received 24 June 1986

The dynamical system describing a continuous stirred tank reactor (CSTR) for the reactions A ~ B ~ C and A--, C, B --, D is considered. A circulating attractor with accompanying circulating orbits is shown to exist when the critical point of the system is unique and unstable. The orbit structure has been numerically found to consist of periodic orbits and chaotic behavior.

1. Introduction It was recently discovered that the dynamical equations describing the behavior of simple first order chemical reactions in a continuous stirred tank reactor (CSTR) exhibit both periodic and chaotic phenomena [1, 2]. Lynch, Rogers and Wanke [1] numerically observed chaos for the first order reactions A --* C, B ~ D by varying one of the nine lumped parameters that occur in the mathematical description of the CSTR. The chaos was studied using autocorrelation analysis and next-amplitude plots. Recently Jorgensen [2] in her thesis also observed chaotic behavior for the first order reaction sequence A ~ B ~ C, which included periodic doubling of at least the fundamental period. At abou~ the same time Retzloff and Chicone [3, 4] showed that for the adiabatic case of A ~ B ~ C there exists a positively invariant plane containing the critical points such that all orbits are ultimately attracted to this plane. Once on this plane the dynamical behavior becomes two dimensional and therefore chaotic behavior is forbidden by the Poincart-Bendixson theorem. In this paper we will analyze the dynamical behavior of a CSTR for both first order reaction processes A ~ C, B --* D and A ~ B ~ C. We will

show that both systems have similar dynamics, chaotic structures and routes to chaos. In section 2 we will show that all critical points lie on a single plane and that a positively invariant trapping block exists for the dynamics of these systems. It will also be shown that the trapping block collapses to two-dimensional planar regions and that chaos will be forbidden in two limiting situations (one being the adiabatic case). In section 3 we will consider the stability behavior of the critical points and establish the existence of a circulatory attractor which is not associated with any critical point for some situations. In section 4 we will discuss their periodic and chaotic behavior and identify the route to chaos for these systems. We present our conclusions in section 5.

2. Invariant regions and critical point planes There are a number of equivalent formulations of the equations describing the dynamics of a CSTR for the reactions under consideration. We will follow the notation of [5] as it lends itself more readily to the analysis we will present in this paper. Therefore for the sake of clarity we state below the relevant equations for the two reaction systems.

0167-2789/87/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

132

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

a) Dynamical equations for A ~ B ~ C: dX x

dl

dx 2 dt

= 1 - x 1 - 0/axl e -v~/x3, ----

m2 -- x2 -- 0/2X2 e-Y2/x3 + 0/lXl e - ~ l / x 3 '

(2.1) dx 3

+ 0/2 f12X2 e --)'2/X3

b) Dynamical equations for A ~ C, B ~ D." dx 1

dt

dx2

dt

= 1 - x 1 - 0/lXl e - ' n / x 3 ,

= 1 --X2--0/2X2 e-;'2/x3

(2.2)

dx3 = q( P dt q - x3) + aaflxXl e-V#x3 q- 0/2~2X2 e -3'2/x3.

It has been shown [5] that nonphysical critical points can occur for (2.1) or (2.2). We therefore note that all physically realizable dynamics and critical points must be restricted to the set S c R 3 defined as S = { ( x 1,x 2,x3)lx 1, x 2 > 0 , x 3>0}. Furthermore, the parameters in (2.1) and (2.2) are constrained by the following relationships: 0 < 0/1~ 0/2 < ~ ' 0 ~ Y1, Y2' m2 < o0,

l
oo,

(2.3)

-- OO < ]~1, ~ 2 < O0.

In what follows, we employ the standard notation that d x J d t = V, and let V denote the vector field V = (V1, 112,1/3). For simplicity, only (2.1) will be considered in detail. The analysis for (2.2) is identical, and the corresponding results are presented at the end of this section. For the adiabatic case of (2.1) (i.e. p = q = 1), it has been shown [3, 4] that S is positively invariant so that all orbits are trapped in the physically realizable domain. Furthermore, all the critical

points lie on a positively invariant plane and therefore chaotic behavior cannot occur. In this section, we extend the analysis to the nonadiabatic case. We will show that all the critical points still lie on a plane, henceforth called the critical point plane. We then show that there exists a subset of g which is a globally attracting sink. This subset is seen to collapse to the critical point plane in the adiabatic limit. Furthermore, for high temperatures (e.g. p / q - - , oo), it also collapses to a plane (albeit different from the critical point plane), thereby forbidding chaos as well. In particular, we will prove that chaos must ultimately disappear for a sufficiently large but finite value of the bifurcation parameter p / q . Proposition 2.1. All the critical points for (2.1) lie on a plane which we call the critical point plane. Proof. The critical points are given by V = 0. This can easily be reduced to the following equation for the critical points:

~'c(D¢):

(fla+fl2)xl+fl2x2+qx3=Dc,

(2.4)

where D c = B1 + f12(1 + m2) +P- Therefore all physically realizable critical points are contained in ~ ( D c ) A S. [] Next, we show that a set R can be defined which contains all the critical points. R will be shown to be a global sink for physically realizable dynamics. In preparation for identifying R we consider the foliation of S by the class of planes ~b(Db) :

(I~I+~2)XI"FI~2X2+x3=Db

(2.5)

with positive normal N b given by N b = ((ill + fl2),#2,1}. Our objective is to obtain a priori bounds for D b (D,,~ and /)max) such that the dynamics of (2.1) will always be trapped within R, i.e. we require the bounds to satisfy < 0, ND" Vl'~b(Ob) > 0,

for D b >_Dmax, for D b
(2.6)

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

:3-P/q

Table I Dmin and Dmax for different signs of fll and f12.A --0 B --0 C. ill, B2

Dm~

Dmi~

>0, >0

fll+f12(l+m2)+ p-q

fl~l p + ~ (q1 + m2) + q

<0, <0

f pl ~ L + ~q( l + m 2 ) +

#1+#2(1+m2)+ -pq

>0, <0

fll+~(l+m2)+

pq

fl~+fl2(l+m2)+ pq q

<0, >0

fll + f l : ( l + m 2 ) + p-~ ~-

fll+~(l+m2)+p q

By (2.1) and (2.5), it can be shown that

133

×t Fig. 1. Trapping planes ~rb(Dmin) and ~'b(Dmax) and critical point plane %(De).

Nb" Vl#b(Ob)= fll + f12(1 + m2) + P + ( q - 1 ) f l , x 1 + ( q - 1)fl2(x 1 + x2) - qDb,

(2.7) so it is evident that the calculation of Dmi~ and D ~ x depends strongly on the signs of fl~ and f12. For illustration, we consider only exothermic reactions (ill > 0, f12 > 0). The results for the other cases are tabulated in table I. 1 Drain = q [fll + f12(1 + m2) + P ] ,

(2.8)

positively invariant and all trajectories of (2.1) initially in S \ R eventually enter R.

Proof. We show that for points in S \ R

and on 3/~ the vector field II always points into R. By (2.1) Vllxl~0-- 1,

Dmax = fll + f12( 1 + m2) + P / q , Vl l x~ ~ l < - al e - ~ / x 3 < O,

whereupon R is defined as V2[x2, o = m : + a l x l e-r~/~3 > O, m

R:

( ( x l , x 2 , x 3 ) l ( x l , x2, xs) ~ ~rb(Db) ;

V3 [o < ~ <_p/q >- a l f l l x l e - ~/(P/q)

Omin -~ D b < Dm~x; 0 < x 1 < 1;0 < x 2 < 1 + m2; p / q < xs < fll + fl2(l + m 2 ) + p / q ).

(2.9)

The trapping planes ~rb(Dmi~) and ~rb(Dmax) together with the critical point plane ~rc(Dc) are shown in fig. 1. We now establish that R is a trapping region for the dynamics of (2.1).

"J'-0t2fl2X 2 e - v2/(p/q) > O.

To prove the upper bound for x2, it is sufficient to establish that x 1 + x 2 < 1 + m E. This follows from the observation that ( V 1 + V2)lxl+x2 >_l+m2 < - o t 2 x 2 e - r 2 / x 3 < O.

n

Proposition 2.2. The set R is a globally attracting

sink for the reaction system (2.1). That is, R is

The upper bounds for x 1 and x 2 are also evident from the principle of mass conservation, For D b

134

D.G. Retzloff et a l . / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

>_ Dmax, by (2.7) Nb"

c R is n o w obtained as follows:

Vl%(o~)--< ~1 "-1-f12( 1 +

m2) + p + (q - 1)/31

Proposition 2.3. F o r A --* B ~ C if both reactions

are exothermic, then the set R can be replaced by R * given b y

+ ( q - 1)f12(1 + m 2 ) - qD~a x = O, so that the m o t i o n is towards the b o u n d e d region in fig. 1. Similarly, for D b < Dmi.,

R*:

{ ( X l , X 2 , X 3 ) [ ( X x , X2, X 3 ) ~ r b ( D b ) ;

D ' i n _< D b < D*a~; r < x I < s; u < x I + x 2 < v;

Nb" V[%(Db) ~--"fll "~ f12( 1

"b m 2 ) + p - q D ~

= O,

so that the m o t i o n is towards the u n b o u n d e d region in fig. 1. Therefore/)max and Dmi, together t r a p the dynamics. Finally the upper b o u n d for x 3 is obvious f r o m (2.5) and (2.8). Thus R is a global sink for (2.1). [] W e conclude that all dynamics of (2.1) will ultimately b e c o m e trapped in the region R. Alt h o u g h R does not completely contain %(De), it is evident f r o m the p r o o f of proposition 2.2 that any critical points for (2.1) cannot lie outside of R, i.e. they all m u s t lie within ~rc( D e) n R. F o r p = q = 1, it is seen that Dmi~ = Dmax = D c a n d the sink region R collapses to the critical p o i n t plane. Furthermore, (2.7) goes to zero, so that the critical point plane itself is positively invariant in this limit. F o r e n d o t h e r m i c reactions, the same analysis yields results that are similar to (2.8). F o r all cases, 0 _ < x l _ < l , 0 _ < x 2_
Dai,

max ~ °ti e - ~'i/[B~ +f12(1+ m2)+ P /q],

Dai, min = ai e - Y J [ P / q ] , where i = 1, 2. A new globally attracting sink R *

P / q -< x3 < ~1 q- f12( 1 + m2) + P / q } ,

where 1 Dma x* = -~ [~81 + ~z(1 + m2) + p + (q - 1)ills + (q - 1)fl2v],

1 Din*in = q [ill q- f12( 1 + m2) -{-P q- (q - 1)fllr + (q - 1)fl2u], and r = 1 / ( 1 + Dal,max), s = 1 / ( 1 + Daa,~in), u = [1 + m 2 + Da2, max/(1 -4- Dal, max)] / ( 1 + Da2, max),

V = [1 + m 2 + D a 2 , = J ( 1

+ Dal,min)]

/ ( 1 + Da2, mi.). F u r t h e r m o r e in the limits q = 1 and p / q - - * oo, chaotic b e h a v i o r is not possible. P r o o f . We note that the bounds on x 3 are still

valid f r o m proposition 2.2. T o verify the bounds on Xl, assume that x I > s. Then V1 = 1

-

x 1 -

alX

1

e -~1/x3

< 1 - x 1 - oqx I e -r~/(p/q) < O.

T h u s x a will definitely decrease until x 1 < s. Similarly, if x I < r, it will definitely increase until x 1 > r. T o verify the bounds on x 1 + x 2, assume x 1 + x 2 > v. T h e n eventually x 2 > v - x~ > v - s.

D.G. Retzloff et a L / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

Therefore

the change of coordinates given by

V1 + V2 = 1 + m 2 -- ( x 1 + x 2 ) -- a 2 x 2 e -v2/x3

Yl = x l ,

< I + m 2-

v - (v --s)Da2,

min=0.

So x I + x 2 will decrease until x t + x 2 < v. Sinailarly, if x 1 + x 2 < u, it will definitely increase until

Y2 ~- X 2 '

Y3 = eX3

135

(2.10)

transforms (2.1) into the system dy 1 - -dt

= 1 - Yl - a l y l e - n ' / Y L

X1-.t.- X2 >__U.

We have established that all orbits starting in S eventually (say at time t*) satisfy the constraints r<_xt t * will continue to satisfy these constraints. Therefore for t > t* we have by (2.7) < 0, Nb° Vl%(Ob) >-- 0,

d Y2 ----"m 2 - - Y 2 -- °12Y2 e-v2e/y3 q'- °tlYl e-Vx~/Y3, dt

--

(2.11)

-d- Y3 = q ( 1 - Y 3 )

dt + e ( a t f l t y x e -vl~/y3 + a2fl2y2e -v2~/y3 ).

,

for D b >_ Dmax, for D b_< D ~ .

We conclude that R* is a global sink. The results presented in section 4 indicate that Din,x and Dm~ are excellent bounds for the dynamics of (2.1). [] For p = q = 1 the region R* collapses to a single plane which must be a global sink in agreement with previous results. Furthermore, for p / q o~, it is easily shown that r--*s, u--* v whereupon D%x ~ D ~ and R collapses to a single plane also. Therefore in these limits chaos cannot exist by virtue of the Poincare-Bendixson theorem. Since the parameter p / q characterizes the dynamics, it will be chosen as the bifurcation parameter. For convenience, we will define B = p / q to be used in our discussions. In the results to follow we establish an even stronger relationship between the bifurcation parameter B and the disappearance of chaotic behavior for the system (2.1). For simplicity, we will again restrict attention to exothermic reactions where fix > 0 and f12 > 0. The constants al, a2, ill, f12 and m 2 remain nonnegative as before and q >_ 1 which follows from its definition. These parameters are considered fixed for the following discussion. Since we will be interested in the behavior of (2.1) for large B it is convenient to define e = 1 / B as a new parameter. Then 0 < e < 1 and, of course, e = 0 as B ~ oo. With these choices,

After an extension to Y3 = 0 by continuity it is clear that the dosed positive octant ~2 = ((Yx, Y2, Y3) E R3[yl, Y2, Y3 >- 0) is a positively invariant set for the flow of (2.11). If e = 0 in (2.11), we obtain the linear system dYl d----T= 1 - (1 + a l ) y 1, dy2 dt = m 2 -

( l + a 2 ) Y 2 + a l y 1,

dy3 dt = q ( 1 - Y 3 ) , which has a hyperbolic stationary point )3 with coordinates

m2(1q- Oll) "]- a 1

1

i l l = l+a-----~' )32= ( 1 + a l ) ( 1 + a 2 ) ' ) 3 3 = 1 " This stationary point is a global sink for the linear flow. We now identify a compact positively invariant set for (2.11) which contains ~. To this end, we define L: R 3 --, R by L(y)

= (ill + fl2)Yl + f12Y2 +Y3

and the plane w(D) as

L(y)=D, where D = fll + f12(1 -t- m2) + q. The normal to

136

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

rr(D) is denoted as N. Then if we define 2 to be the closed subset of $2 bounded by the coordinate planes and ~r(D), we have the following lemma:

Lemma 2.1. The set X is compact and for 0 ~ e < 1 in (2.11) satisfies the following properties: (i) X is positively invariant; (ii) the point .f is in the interior of X; (iii) if Y0 ~ 82\2, then the trajectory of (2.11) starting at Yo enters the interior of 2 at some finite time. Proof. Clearly ,~ is compact and is positively invariant, as is easily checked by computing the vector field along the boundaries of X. Also, all the coordinates of )3 are nonzero, and L ()3) < D, so that 33 lies in the interior of X. Finally dL

N- V = ~

= H I +/32(1 + m 2 )

+ (1 - q)Y3 + (e

-

+q-L(y)

1)(o~1/~1y 1 e -'yle/y3

+ a2flzy2e-V2E/Y3 ).

For all points which lie in [ 2 \ 2 , N . V< 0 since L ( y ) > D, which means that after a finite time the trajectories starting at these points enter the interior X. []

Theorem 2.1. There exists a finite Bmax of the bifurcation parameter B for the dynamical system (2.1) such that for B > Bm~x (2.1) has a globally attracting hyperbolic sink. Proof. The dynamical system (2.11) is linear when e = 0 and thus admits a ball Br(.13) c int 2 whose boundary is transverse to the vector field of (2.11). The flow within Br()3) at e = 0 is structurally stable by virtue of the fact that )3 is the sole invariant set and is hyperbolic as a critical point. Thus if we perturb e slightly from e = 0, the vector field remains transverse to OB~()3) and trajectories originating on the boundary remain within Br(.13) for positive time. Furthermore they ultimately collapse to the critical point within Br()3). The argument will be complete if we show that all trajectories must enter Br(.13) if e is sufficiently close to

zero. To this end we note that by lemma 2.1 for e = 0 trajectories originating at Z ~ OX cross the boundary of J~r()3) at some finite time To(Z ) which is continuous in Z. Because of the continuous dependence of the flow on e, for e > 0 and small, the trajectories of (2.11) originating on 0 2 must also cross OBr(y) at some finite time T~(Z) which is jointly continuous in (e, Z). By lemma 2.1 all trajectories outside of X must ultimately cross 02 and hence 0Br()3 ) transversely. Upon entering B~(j3), a trajectory is ultimately attracted to the unique critical point within Br(.V)" This establishes the globally attracting character of the unique critical point for the system (2.11) when e is near zero. However, for e > 0, the dynamical system (2.1) is conjugate via (2.10) to the dynamical system (2.11). From this the conclusion of the theorem follows. [] Similar arguments can be used to show that the conclusion of theorem 2.1 is true for the system (2.2). It is an immediate consequence of theorem 2.1 that chaos cannot occur for B > Bm~. Thus the domain of the bifurcation parameter for which chaos occurs is a finite interval. This dynamical behavior is exactly what is observed in the numerical results discussed in section 4 and summarized in table II. It is apparent from this table that the trapping block collapses, as measured by the distance between the bounding planes ~rb(D*~, ) and rrb(D*n), when the observed chaotic behavior disappears for increasing values of the bifurcation parameter. For the numerical example shown in table II, Bm~, 1.2. We now state the significance of the critical point plane for the dynamical behavior under consideration. This is the content of the following proposition:

Proposition 2.4. In the adiabatic case (p = q = 1) the critical point plane is a global sink. For the nonadiabatic case orbits that intersect the critical point plane have a nonzero component of V along the normal direction N c and hence will pass

D.G. R e t z l o f f et a l . / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

137

Table II Collapse of the trapping block R* with increasing bifurcation parameter B for the dynamical system (2.1). The parameter values are eq = 7.20 × 1011, a 2 = 1.98 x 10 u , 131 = 10.0, 132 = 3.725, 35 = 72 = 25, q = 250.0, m 2 = 0.0.

Bifurcation parameter B

d(~rb(Dm*a~) - ~rb(D*i,)) = D~a * x - Dmi *n Comments

0.004

13.67

0.2

13.67

0.8

13.48

0.95 0.9999

single stable critical point single stable critical point single stable critical point periodic orbit * periodic orbit *

5.08 2.146

0.9999-1.0000

chaos

1.0003 1.1 1.2

2.018 0.274 0.042

3

1.6 x 10 7

periodic orbit * periodic orbit * single stable critical point single stable critical point

*Single unstable critical point.

Table III Invariant regions and critical point planes for (2.2). 131 > 0, 132 > 0 Critical point p l a n e

X5

~rc(Dc) = 131Xl q- 132x2.4_ q x 3 = D c

t

Dc ~ flt + fl2 + p Invariant region

!

¢rb(Db):

131X1 4. 132X2 4- X 3 = D b

1

Omin = -~ (131 + 132 + P ) Xz

£ Dmax = 131 + 132 4- q

am*in = q[131 + 132 + P + 1 3 1 ( q - 1)/(1 + DaLmax ) + 132(q - 1)/(1 + Da2,m~,) ]

1 Om*ax= q [ f l t 4- 132 4- P 4- 13x(q - 1)/(1 + Oa,,mi,) ×l

Fig. 2. Critical point plane for large activation energies. - - - , loci of points where V3 = 0.

+ 132(q - 1)/(1 + Da2.mi,) ] where Dai,max = a t e Y'/(fll+B2+P/q) Dai, min = oti e Yi/( P/q)

13 8

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

through the critical point plane, except at points where V3l,~c(oo) = 0. Furthermore, for A ~ B ~ C when either

(i)

131 + f12(1 + m2) < O,

1)

or (ii)

3'1 ' >4

l[qP

q

BI+B2(I+mE)+p

1'

the set No. V = 0 is a continuous nonclosed curve F parameterized by x 3. In particular, this curve intersects the coordinate planes given by x 1 --0 and x 2 = 0. It therefore separates the steady state plane into two separate regions where Arc. V > 0 and N c • V < 0 respectively.

It is therefore evident that for any x 3 > 0, we can calculate unique values of x 1 and x 2 by (2.12) and (2.13). However, only %(De) N S is of interest for the present problem, so that within some range of x 3, there may not be any Xl, x 2 ~ S which contribute to the curve F. For example, as x 3 ---, 0, V3 ~ p > 0, so that F cannot intersect this plane. It only remains to show that under the hypothesis F intersects the x I = 0 plane and the x 2 = 0 plane once each. We therefore look for zeros of (2.12) and (2.13). For (2.12), we consider g(to) + h(to) = 0,

(2.14)

where g(to) = {[ill + f12( 1 + rn2) + P ] 0.1- q} a2 e-Y2",

Proof. T h e normal N¢ to the critical point plane %(De) is

Nc = ( ( i l l nt- f12), f12, q}, so that N c • VI~,(D0 = ( q - 1) V3 I.(o0" Therefore, in the adiabatic case, the critical point plane is a global sink. For q > 1 (nonadiabatic case), the motion is towards the bounded region lying below %(De) in S if V3 < 0 and towards the unbounded region lying above %(Dc) if V3 > 0 . Thus, N c . V = 0 i f and only if V3 = 0 . The loci of all the points on %(Dc) where this is satisfied can therefore be obtained by solving V3 = 0 and (2.4) simultaneously to give

X 1 ~---

x2=

P-

qx3 + [ill + f12( 1 + m2) + P -

-- 0till e-V'/x'

+

h ( to ) = p to - q ,

and 0.1= 1 / x 3. Clearly, g(to) has a zero at 0.11= q / [ f l l + f12(1 + m2) + P ] and h(to) has a zero at 0.12 = q / p . The analysis of (2.14) follows identically a procedure outlined in [3], which will not be repeated here. The results are: i) For fll Jr- f12(1 + m2) + p > 0, a) 6.12_~ 0.11 (i.e. fiX + f12(1 + m2) < 0): (2.14) has one solution. b) toE > 0.11 (i.e. fix + fiE(1 + mE) > 0): (2.14) has one solution if 4/3'1 > q / p - q / [ f l l + fiE(1 + m 2) + P], and may have three solutions otherwise. ii) For fll ÷ f12(1 + m2) + P < 0, (2.14) has only one solution for all parameter values. For (2.13), an identical analysis gives the same results but with 3'1 replaced by 3'2[]

qx3]a2 e-v2/x'

°/2(]~1 ÷ ~2)

(2.12)

e-3'2/x3

--~t1~1 e-Y1/x3 [~1 ÷ J~2( 1 ÷ m2) ÷ P -- q x 3 ] -- (~1 ÷ ~ 2 ) ( P f12 [ - alfll e-:'I/x' + ~2(fll + f12) e-:'2/x']

-- q x 3 )

(2.13)

139

D.G. Retzloff et al. / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

The significance of the critical point plane in the nonadiabatic case lies in the fact that it is a natural place to find a Poincar6 section as a subset, since it divides the invariant region R* into two subsets neither of which are themselves invariant. Furthermore, we have established in proposition (2.4) that under the assumption of sufficiently small activation energies the critical point plane %(D~) itself is divided into two regions by the curve F. These two subsets of %(Dc) will each be a Poincar6 section. That this is the case together with an analysis of the orbit structure will be the subject of section 4. For large activation energies, however, the hypothesis of proposition (2.4) is violated whereupon more exotic dynamic behavior may arise. A possibility is depicted in fig. 2. It may be of interest to note that the assumption of infinitely large activation energies was used by Jorgensen [2] in her thesis. An analysis of the Poincar6 sections for large activation energies will be presented in a future communication. The invariant regions and critical point plane for the dynamical system (2.2) are listed in table III.

be odd. We shall denote the critical points as x~i~), where i = 1 , 2 . . . . . 2 n + l and xt~
X3 _ A1X2 +

A 2 ~ _ A3 =

0,

(3.2)

where A 1 = trace of J, A 2 = sum of principal minors of J, A 3 =

determinant of J.

Representative examples of the critical points that can occur and their eigenvalues are given in table IV. It is significant to note that the eigenvalues follow well-defined patterns for all the cases that are considered. This is the subject of the next two propositions. Again, for simplicity, we shall only consider the consecutive reactions A ~ B ~ C in our analysis, since the parallel reactions A ~ C, B ~ D can be handled in identical fashion. Proposition 3.1. If there are no degenerate roots of

3. Stability analysis and the Hopf bifurcation The number of critical points that exist for the dynamical systems A --* B --* C or A --* C, B -o D can be obtained from (2.1) or (2.2). For A -o B -o C, this is done by calculating the zeros of x 3 ) + fll

Proof. We observe that F(0) > 0, F(oo) < 0, and therefore a root must exist. When F ( x s ) crosses

the x 3 axis, the derivative of (3.1) must be such that

otI e - yl/x3 F(x3) = q( B-

(3.1), then the number of positive real eigenvalues at x~/) is even when i is odd, and odd when i is even.

1 + Ot1 e - T t / x 3

dF f <0, dx3 x~3~ > 0 , +f121+a2 e-~/~

m 2 + l + a l e -~/~3

for i = 1,3 . . . . . 2 n + l , fori=2,4,...,2n.

(3.3)

"

(3.1) It has been shown [3, 4, 5] that as many as seven critical points can exist for (3.1). Furthermore, it is well known that the number of critical points (assuming that there are no degenerate roots) must

Furthermore, it can be shown by a straightforward calculation that at each critical point x ~ ) A~ i ) = (1 + cqe-VJx~'))(a + a2e v2/x~) d___ff__F

-

'

dx 3 xt3~'

(3.4)

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

140

Table IVa Representative

list o f s t e a d y s t a t e s ,

their eigenvalues

a n d s t a b i l i t i e s f o r r e a c t i o n A -* B ---, C

Steady state 1

x1

x2

1.000

x3

0.10097

1.3

2~l - 1.00

X2

?t 3

- 1.00

-0.999

Stability 3-D stable

Parameters a 1 = 2.279; a 2 = 4.2 × 101° ,81 = 0.79; ,82 = 43.6

2 3

1.000 0.9643 "

0.08714 0.00

1.9007 7.2800

-1.00 -4.37

4

0.9309

0.00

5

0.6313

0.00

8.758

-1.393

× 10 s

-1.00

-4.342

× 109

-1.00

-0.575

3-D stable

1

0.9957

0.1134

3.04

2 3

0.6285 0.2531

0.01049 0 . 5 2 4 × 10

5.457 7.35

-1.0206

- 1.00

-0.996

3-D stable

a 1=300;

- 41.445 - 1627.36

- 1.00 - 1.64

0.75 1.00

1-D unstable 3-D stable

/31 = - 1.05; f12 = 5.95

22.00

3

- 1.00 - 1.00

× 107

7.156 -0.142 0.122

1-D unstable 3-D stable

3'1 = 30.0; 3'2 = 50.0 m 2=0.1;p=1.3 q = 1.0

1-D unstable

3'i = 34.0; 3'2 = 76.0 m z = 0.109; p = 3.05 q=

1

0.998

3.001

1.004

- 1.00

- 1.00

2

0.833

3.167

1.50

- 1.00

- 1.00

- 99.06

- 1.00

-0.937

a 2=5×107

1.0

3-D stable

a 1 = 4405.293;

1-D unstable

,81 = 3.0; f12 = 0.75

3-D stable 1-D unstable

3'1 = 15; 3'2 = 100 m 2=3.0;p=1.0 q=

a 2 = 36346337.88 2.133

3

0.981 X 1 0 - 2

3.988

3.97

4

0.278 × 10-2

1.35

5.97

- 357.4

- 1.00

-0.992

5

0.205 × 10-2

0.248

6.806

- 486.3

- 10.03

- 1.00

3-D stable

1

0.998

3.00

1.004

-1.00

- 1.00

-0.937

3-D stable

2.624

1.0

a I = 4405.293; a 2 = 3634337.88

2

0.833

3.167

1.5

1-D unstable

,81 = 3.0; 132 = 0.75

3

0.982 × 10-2

3.99

3.97

- 98.94

1.00

-1.00 - 1.00

- 1.00

2.133

3-D stable

3'1 = 15.0; 3'2 = 100.0 m 2=3.0; p=1.0 q=l

1

0.0179

0.0612

1.073

- 22.546

- 5.711+97.38i

3-D stable

a 1 = 7 . 2 0 0 4 9 x 1011; a 2 = 1.98 × 1 0 1 1

/31 = 10.0; ,82 = 0.75 3'1 = 3 ' 2 = 2 5 ; m 2 = O p = 254.98; q = 250 1

0.0283

0.093

1.052

-14.374

6.07 +

77,28i

2-D unstable

a I = 7 . 2 0 0 4 9 X 1011; a 2 = 1.98 × 1 0 1 1 ,81 = 1 0 , ,82 = 3.725,

3'1 = )'2 = 25 m 2 = 0; p = 249.751; q = 250.0

whereupon by (3.3) A(3i) = { < 0, > 0,

for i = 1 , 3 , . . . , 2 n + 1, for i = 2 , 4 , . . . , 2 n .

(3.5)

Since A~0 is also equal to the product of the three eigenvalues, the proposition follows i m m e d i a t e l y . []

Remarks. 1) Our results do not depend on whether there are complex eigenvalues, since any complex eigenvalues will always exist in conjugate pairs and therefore do not change the odd/even properties. 2) There are no restrictions on whether the reactions are endothermic or exothermic. However, the dynamical behavior can be restricted

141

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor Table IVb Representative list of steady states, their eigenvalues and stabilities for reaction A ---, C and B ~ D Steady state

x1

x2

x3

1 2 3 4 5

0.999 0.999 0.859 0.836 0.114

0.955 0.687 1.400 × 10- 2 1.17 × 10-2 3.0027 × 10 -4

1.0230 1.157 1.6341 1.658 2.386

1 2 3

0.986 0.7075 7.731 × 10 -2

0.9363 0.326 1.648 × 10 -2

1.00308 1.00386 1.0762

- 1.0246 - 2.0989 -57.505

1

0.0223

0.0766

1.063

- 16.39

0.50708 + 92.42i

1

0.01785

0.06201

1.074

- 20.378

- 6.357 + 103.253i

~kl

)k2

- 1.0036 - 1.0048 - 67.7681 - 81.714 - 3328.5703

-

~3

1.00 1.00 1.00 1.00 5.654

- 1.00 - 1.00 -9.522

Stability

Parameters

- 0.6198 0.8977 - 4.981 x 10- 2 5.0646 × 10-2 - 1.00

3-D 1-D 3-D 1-D 3-D

stable a 1 = 33961.56 unstable a 2 = 1 6 1 7 2 1 7 . 3 1 8 stable fll = 1.0; 132= 0.5 unstable Y1 = 72 = 20 p = 1.0 stable q = 1.0

- 0.7505 1.2019 -1.00

a 1 = 2.688 x 1041 3-D stable 1-D unstable a 2 = 1.344 x 1042 3-D stable /31 = 0 , 0 4 , /32 = 0 . 0 4 71 = Y2 = 100.0 p = 1.0 q = 1.0 2-D unstable a 1 = 7.20049 X 1011 a 2 = 1.98 × 1011 131 = 10; /32 = 3.725 Y1 = "/2 = 25 p = 252.49 q = 250 3-D stable

a 1 = 7.20049 × 10 n a 2 = 1 . 9 8 × 1011 /31 = 1 0 . 0 ;

/32 = 3.725 Y1 = Y2 = 25 p = 254.98; q = 250.

f u r t h e r f o r e n d o t h e r m i c r e a c t i o n s , as w e s h a l l s h o w

For the remainder consider problems

in the next proposition.

of this paper, we shall only

w h e r e iR c o n t a i n s e x a c t l y o n e

unstable critical point.

Proposition

3.2. F o r

endothermic

reactions

(or

s l i g h t l y e x o t h e r m i c r e a c t i o n s ) , t h e c r i t i c a l p o i n t is unique and stable.

Theorem 3.1. A s i n g l e u n s t a b l e c r i t i c a l p o i n t m u s t have

a

two-dimensional

dimensional

unstable

and

stable manifold. Furthermore,

a

one-

in this

c a s e t h e r e is a n a t t r a c t o r i n iR ( d i s t i n c t f r o m t h e

Proof. I t c a n b e s h o w n t h a t A 1 < 0, A 2 > 0, A 3 < 0 f o r f l l < 0,

r 2 < 0. T h e r e f o r e ,

critical point).

b y (3.5), t h e r e is

o n l y a s i n g l e c r i t i c a l p o i n t . T h i s is a p r e v i o u s l y

Proof. B y p r o p o s i t i o n 3.1, t h e n u m b e r o f p o s i t i v e

k n o w n r e s u l t [6]. N o w , b y R o u t h ' s s t a b i l i t y c r i t e r i a ,

real eigenvalues for a single critical point must be

s t a b i l i t y is a s s u r e d if A1A 2 - A 3 < 0. A c u m b e r -

e v e n . T h e p o s s i b i l i t i e s a r e : (i) T h e r e a r e t w o p o s i -

s o m e c a l c u l a t i o n s h o w s t h a t t h i s is, i n fact, v a l i d .

t i v e r e a l e i g e n v a l u e s . T h i s is c a l l e d a n u n s t a b l e

Furthermore,

s i n c e A1, A2, A 3 a r e c o n t i n u o u s i n

n o d e . (ii) T h e r e a r e n o p o s i t i v e r e a l e i g e n v a l u e s ,

fix a n d f12, t h e r e m u s t e x i s t 81,82 so t h a t t h e a b o v e c o n c l u s i o n s a r e t r u e f o r 0 < / 3 1 < 81, 0 < 132

which means

_< 82 a s well.

the critical point

[]

that a complex conjugate pair with

positive real parts must exist since by assumption is u n s t a b l e . T h i s is c a l l e d a n

142

D.G. Retzloff et a L / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

unstable focus. In both cases, the remaining eigenvalue must be real and negative. Therefore, we have a two-dimensional unstable and a onedimensional stable manifold. To prove the second part of the theorem, we note that all orbits starting on the two-dimensional unstable manifold near the critical point cannot return to the critical point along the stable manifold direction due to the uniqueness of solutions of (2.1). Hence, there must exist orbits leaving a neighborhood of the critical point along the unstable manifold which do not return to the critical point. Furthermore, these orbits cannot leave R by proposition (2.2). Therefore, there must exist an attractor in R which is distinct from the critical point. [] It is natural to divide the region R into four subregions by the critical point plane and the surface given by V3 = 0; these surfaces intersect in the curve of proposition 2.4. Also, the critical point plane separates S into a bounded region which we say is below the critical point plane and an unbounded region which we say is above the critical point plane. Region 1 is the set above the critical point plane which has boundary points in the xlx2-plane; region 2 is the other region above the critical point plane; region 4 is the region below the critical point plane with boundary points in the XlX2-plane; region 3 is the remaining region. For a point P in the interior of region 1 we define the symbol sequence for P to be the sequence in the symbols 1, 2, 3, 4 determined by the order in which the four regions are transversed by the trajectory through P. A circulating orbit is an orbit such that for any t > 0 there is a T > t such that the symbol sequence ~r(P) begins with the form 1, 2, 3, 4 . . . . . Of course, if an orbit is circulating the ordering 1,2,3,4 occurs infinitely many times in the symbol sequence. An invariant set, in particular an attractor, is called circulating if all of its points lie on a circulating orbit.

Theorem 3.2. If R contains exactly one critical point which is a hyperbolic saddle with onedimensional stable manifold then all orbits except

the critical point and those on the stable manifold are circulating orbits. In particular, the attractor of theorem 3.1 is a circulating attractor.

Proof. Consider the foliation of R by the planes Dram < D b < Dmax, defined by (2.5) and let N b denote the normal vector for the plane ~ b ( D b ) . An easy calculation shows

7rb(Db) ,

+P - (ill + fl2)x1

--

1~2X2

--

qx3.

(3.6)

Hence, in view of (2.4), N b • V < 0 at points above the critical point plane and N b • V > 0 at points below the critical point plane. Assume that P is in the interior of region 1. If P lies on the stable manifold of the critical point then P does not lie on a circulating orbit because the stable manifold is tangent to the stable eigenspace of the linearized equations at the critical point. If P does not lie on the stable manifold we observe that the positive invariance of R and the fact that N b ° V < 0 above the critical point plane imply the orbit through P eventually crosses the critical point plane. There are two possibilities: it crosses leaving region 2 or it crosses at a point of F. To see that only the first possibility occurs, consider a crossing at a point on F and observe that by continuity with respect to initial conditions all points in a small disk on the critical point plane centered at the crossing point must flow below the critical point plane. However, some of these points must lie in the boundary of region 1 where all points cross the critical point plane but in the opposite direction. Hence, the orbit through P must enter region 2, cross the critical point plane and enter region 3. A similar argument shows that after entering region 3 the trajectory must eventually cross the critical point plane from region 4 to region 1. Of course, it is possible that the trajectory through P crosses into region 2 and returns to region 1 so that a finite portion of the symbol sequence has form 1, 2,1, 2,1, 2,... before the orbit eventually crosses into region 3. Similarly, finite sequences of the form 3, 4, 3, 4, 3,4 .... are possible. However, since the sequences 1, 4 and 3, 2 are

D.G. Retzloffet aL/ Chaoticbehaviorin the dynamicalsystem of a continuousstirredtank reactor not possible it is easy to see that the sequence 1, 2, 3,4 occurs infinitely often in the symbol sequence. In particular, after passing into region 1 the only possible sequence is a finite number of repetitions of the sequence 1, 2 followed by 3, 4. Thus after the last repetition of the sequence 1, 2 we have the desired result. [] Remarks. 1) The proof of theorem 3.2 shows there are orbits which contain the sequence 2,1 or the sequence 4, 3. We are not able to completely analyze the situation to determine all possible symbol sequences. 2) This circulating attractor could be a periodic orbit or a more complicated attractor with associated chaotic behavior. Previous choices of the bifurcation parameter for (2.1) or (2.2) are: q for A ~ B ~ C [2], and fl2/q for A ~ C, B ~ D [1]. However, in both cases because of the lumped parameter nature of the description it is not possible to interpret varying the bifurcation parameter as varying a single physical variable. Thus, it would be difficult to replicate any predicted results in a physical experiment. We have chosen B as the bifurcation parameter which can be interpreted physically as varying the temperature of the heat exchanger fluid, an easily controlled physical quantity. For the next theorem, we define the stability boundary to be the points Bc in parameter space at which the real parts of at least one eigenvalue of the critical point is zero. Theorem 3.3. For a region of parameter space where the critical point is unique, a Hopf bifurcation with B chosen as bifurcation parameter arises at the stability boundaries, except where d(AIA 2 A3)//dx3c = 0. At least one periodic orbit exists in the vicinity of the boundary.

143

unstable region, therefore, either A t = 0 or A t A 2 A 3 = 0 must happen. If A 1 = 0 occurs first, then A1A 2 - A 3 > 0, which is a contradiction. Therefore, at the stability boundary, A t < 0, A 3 < 0, A 2 = A 3 / A 1 > 0 is the only possibility. The eigenvalues are h 1 = A1, X2 = A~2i, X 3 -- - A~2i. To satisfy the conditions of the Hopf bifurcation theorem, it remains to show that the pair of complex conjugate eigenvalues cross the imaginary axis transversely at the critical point values of the bifurcation parameter [7]. This is equivalent to showing that -

d-dB"( R e X2) sc

_[

d (A1A2_A3)] :#0. 2(A2-+ A 2 ) dB jso (3.7}

For the present problem, the Jacobian matrix depends only implicitly on B, whereupon the derivative of (3.7) must be obtained by the chain rule, i.e. d

(AIA2_A3)=

d

(AIA 2

dS

A3) dx3e -

dB

•

By (3.1) and (3.4), dB dx3c =

1 d F ~3o q dx3 A3 (1 + al e-Vl/~'c)(1 + a2e-'=/"'c)q '

so that (3.7) becomes d B ( R e x2) s~ (1 + ale-Vl/x'o)(1 + et2e-'=/~'c)q =

-

2A,(A2+A~)

-

Proof. By Routh's stability criteria, a critical point is stable if and only if A t < 0, AIA 2 - A 3 < O, A 3 < 0. For a single critical point, we already know that A 3 < 0 is assured. To move into the

d ~( A I A =_A3 ) x-a-~--£3 By hypothesis, d ( A ~ A 2 - A 3 ) / d x 3 ¢ is nonzero. Also it is evident that the remaining term in the above expression is always positive and finite. Hence a Hopf bifurcation exists. []

144

D.G. Retzloff et a L / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

Table V Representative list of parameters that give rise to Hopf bifurcation l~ ' ' (0)

p

A ---)C, B - o D a 1 = 7.20 × 10 n a 2 = 1.98 × 1011 fll = 10.0; r2 = 3.725

- 1.0858 X 1012

252.69

-6.1308

225.47

71 = 72 = 2 5 0 q = 250 a t = 4.85 × 109

X 1012

a 2 = 2 . 9 4 X 1013

fll = 11; flz = 3.498 y1=20; T2=30 q = 220 al = 4.85 × 109 a 2 = 2.94 × 1013 B1 = 10.0; B2 = 3.18

-3.3389 × 1012

204.23

- 1.1026 × 1012

252.80

-- 3.4419 × 1013

261.03

- - 2 . 1 2 6 1 × 1014

211.96

Y1 = 2 0 . 0 ; Y2 = 3 0

q = 200 A---) B --* C a x = 7.20 × 10 n a 2 = 1.98 × 10 n f ll = 10.0;

f12=

3.725

=Y2 =25 q = 250.0; m 2 = 0.0 "Yl

a 1 = 7.20 × 10 I1 a 2 = 1 . 9 8 × 1011 fll = 10.0; r2 = 3.725 "Y1 = "/2 = 2 5 . 0 q = 2 5 0 . 0 ; m 2 = 1.0 a 1 = 4 . 8 5 × 109

a 2 = 2.94 × 1013 fll = 10.0; /32 = 3 . 1 8 3'1 = 20.0; "/2 30.0 q = 200.0; m z = 1.0 =

a n y r o o t s exist, they m u s t exist in pairs (two, four, etc.), a s s u m i n g that they are nondegenerate. E a c h r o o t then gives a B~ which is a c a n d i d a t e for the s t a b i l i t y b o u n d a r y . F u r t h e r m o r e , at any degenera t e r o o t of A I A 2 - A 3, the derivative is zero which is p r e c i s e l y the exception of t h e o r e m 3.3. A t these p o i n t s the H o p f b i f u r c a t i o n t h e o r e m does n o t a p p l y a n d m o r e exotic b e h a v i o r such as the a p p e a r a n c e o f two or m o r e p e r i o d i c orbits m a y occur. F o r the d y n a m i c a l b e h a v i o r investigated n u m e r i c a l l y this d e g e n e r a c y d i d n o t occur. R e p r e s e n t a t i v e p a r a m e t e r values which give rise to a H o p f b i f u r c a t i o n are presented in table V. F i n a l l y , we p o i n t out that an algorithm exists [7] to a s c e r t a i n as to whether the p e r i o d i c orbits r e s u l t i n g f r o m the H o p f b i f u r c a t i o n are stable or u n s t a b l e . T h e significance of this calculation lies in the fact t h a t if the p e r i o d i c orbit is unstable, t h e n it w o u l d b e unlikely to be observed either n u m e r i c a l l y o r p h y s i c a l l y as n u m e r i c a l noise or p h y s i c a l fluctuations will cause the d y n a m i c s to evolve to a stable orbit. Thus, only the stable p e r i o d i c o r b i t s are directly a m e n a b l e to experim e n t a l o b s e r v a t i o n . Briefly, the p r e s c r i b e d alg o r i t h m involves first t r a n s f o r m i n g to a new c o o r d i n a t e system, a n d then d e t e r m i n i n g the sign o f a c e r t a i n f u n c t i o n 17"' (0). The r e a d e r is referred to [7] for details; it suffices to note that this is an e x t r e m e l y c u m b e r s o m e calculation a n d m u s t be p e r f o r m e d n u m e r i c a l l y in each case. F o r the p a r a m e t e r values listed in table V, I7 ' ' (0) is f o u n d to b e a large negative n u m b e r which indicates that the p e r i o d i c o r b i t is indeed stable.

F o r a given set of p a r a m e t e r values, the calculation of the stability boundary Be - Be

4. Periodic and chaotic behavior

( a l , a2, fla, f12, "/1, "/2, m2, q) is in general a s t r a i g h t f o r w a r d p r o c e d u r e which involves solving F = 0 a n d A I A 2 - A 3 = 0 s i m u l t a n e o u s l y to o b t a i n x3c a n d B c. F o r the present p r o b l e m , we n o t e that A 1 A 2 - A 3 d o e s n o t involve the b i f u r c a t i o n p a r a m eter B a n d will therefore give x3c directly. This m a y t h e n b e s u b s t i t u t e d into (3.1) to calculate B~. It is o f i n t e r e s t to n o t e that the f u n c t i o n A I A 2 - A 3 is n e g a t i v e at x3~ = 0 a n d as x3~---, 00, so that if

F u r t h e r insight into the n a t u r e of the a t t r a c t o r is o b t a i n e d f r o m a numerical d e t e r m i n a t i o n of the o r b i t structure. A s it is a l r e a d y shown that chaos c a n n o t o c c u r for sufficiently large B, we n o w start w i t h a sufficiently large value of B for which the critical p o i n t is unique a n d stable, a n d decrease the v a l u e o f B to observe the H o p f bifurcation, p e r i o d i c p h e n o m e n a a n d chaos. This is consistent

D.G. Retzloff et al. / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

O, 021 0,0l 0.00

, ......... 20

, ..... .... , ......... "21 22

, ,

23

i

24

. . . . . . . . .

1

25

F i g . 3. T i m e - t e m p e r a t u r e plot of fundamental period for A ---, B ---, C, x~ = x3 - 1; B = 1.002988. q = 250.0, a l = 7.20 X 1011, a 2 = 1.98 × 1011,/31 = 10,/32 = 3.725, Vl = 3'2 = 25.0, m 2 =

0.0.

with the approach generally taken in literature. For the reaction A ~ B ~ C, fig. 3 exhibits the graph of x 3 for the stable periodic orbit arising from the Hopf bifurcation. For subsequent analysis this period which is evident in fig. 3 is considered as the fundamental period. Upon decreasing the value of the bifurcation parameter, periodic orbits whose period is an integer multiple of the fundamental period are observed. For the parameter values that are chosen, the observed bifurcation sequence is (stable critical point, 1, 2, 3, 2, 10, 2, 10, 4, 5, chaos, 3, stable critical point}. (Note: In the remainder of this paper we will use the term chaos to denote the dynamics of (2.1) or (2.2) when the attractor is not a stationary point or a periodic orbit.) The corresponding orbits, Poincar6 sections and time-temperature

145

plots for period 10 and chaos are given in figs. 4 and 5. The general evolution of the attractor can be characterized by the addition or subtraction of one or more extra turns in the center hole of the attractor. This phenomenon has also been numerically observed for the Belousov-Zhabotinsky reaction [8]. Furthermore, this periodic to chaotic behavior, which occurs as the bifurcation parameter is varied, has thus far been numerically observed only when there is a single unstable critical point and ceases when the bifurcation parameter has decreased sufficiently that the single unstable critical point becomes stable again. This is consistent with theorem 3.1 and therefore the orbits shown depict the structure of the attractor identified by theorem 3.1. Further evidence of the chaotic nature of the orbit in fig. 5 is given in fig. 6 which shows a blow up of the Poincar6 section for this orbit. An important point to note is that the values of the parameters chosen are physically realizable. Specifically the infinite activation energies approximation is not employed. In addition, while m 2 = 0 is a typical experimental condition, fig. 7 shows that chaotic behavior similar to that depicted in fig. 5 also occurs when rn 2 ~ 0. Thus the choice of m 2 = 0 has no special significance for the existence of chaotic dynamics. The sequence of orbits in figs. 4 and 5 shows that a more pronounced spiraling structure occurs when bifurcation parameter B is decreased. This may then evolve to a homoclinic orbit of the type shown in fig. 8. When the dilation along the unstable directions is greater than the contraction along the returning stable direction Shil'nikov has shown the existence of chaos in the form of infinitely many horseshoes in the appropriate Poincar4 map in a parameter neighborhood of the condition for homoclinicity [9, 10, 11]. This is suggestive of the mechanism for chaos in this reaction sequence. Current efforts are being directed to numerically establishing that the homoclinic orbit suggested above actually exists. The period-doubling phenomena [12] or period-adding phenomena [8] which are well established routes to chaos were not numerically observed. Changes in the period

146

D. G. Retzloff et a L l Chaotic behavior in the dynamical system of a continuous stirred tank reactor 0,39-

I 0,36-

0.330.30-

0.27-

x;

0.24-

0.21-

0.18

0110

0.15

0,~22

~ "

o 0.06

o o

O.OOn,J

a

.

.

.

.

.

.

.

.

.

20

,

.

.

.

.

.

.

.

.

.

21

,

.

.

.

.

.

.

22

.

.

.

,

23 '

.

.

.

.

.

.

.

.

.

,

24

.

.

.

.

.

.

.

.

.

,

25

b

4-4. 4. 4.

4.4-

÷ 4-

O

÷

+ .000

C .004

.008

.012

.016

.020

.024

Yf Fig. 4. Period 10 for A -* B ---, C; (a) 3 D plot; (b) T i m e - t e m p e r a t u r e plot, x~ = x 3 - 1; (c) Poincar6 section, Yi, Y2 are coordinates tangent to ~rc(Dc); B = 1.0017118. Other constants same as fig. 3.

D.G. Retzloff et a l . / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

147 u~

a

oAe

o.~

) 0.10

O.O$

0,06

0.04

2

0.02

b

0.50-

Q. 02

0,0~-

~ ° ~ 0.45-

O.O,S-

0.40-

0.07-

°" " ~ * ~

~'; .~'~"

~"~'~

"" ~'..

0.350,06-

N

./

0.30-

O. 05-

f

#

0,25-

0.O4-

/

0.200.03.

I i

0.15t

I

0.02. O0

0.01.

0.05-

!

0.00

O. 00, . . . . . . . . .

20

, . . . . . . . .

2t

, . . . . . . . . .

22

i . . . . . . . . .

23

i . . . . . . . .

24

,

25

. . . . . . . .

, . . . . . . . . .

26

r . . . . . . . . .

27

i . . . . . . . . .

28

i . . . . . . . . .

29

,

30

0.000

l 0,004

0.008

0.0t2

0.016

0.020

0,024

t

Fig. 5. Chaos for A - * B ~ C; (a) 3D stereo plot; (b) Time-temperature plot, x~ = x 3 - 1 ; coordinates tangent to qrc(Dc); B = 0.9999004. Other constants same as fig. 3,

(c) Poincar6 section, Yl, Y2 are

D.G. Retzloff et a L / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

148

¢,

(JD

,,¢

o~

00802

.00806

.00810

.008|4

.00818

.00822

.00826

Fig. 6. B l o w - u p of fig. 5c. All c o n s t a n t s same as fig. 5. Straight lines are d r a w n in to i n d i c a t e the observed trend.

0

0.I0

C

o.1o

o'.o+

o'.o~

l

o.o{

G'.o2

Fig. 7. 3 D s t e r e o p l o t of chaos for A ~ B --* C; B = 1.0007968, m 2 = 0.01. O t h e r c o n s t a n t s s a m e as fig. 3.

of the orbits were found to occur for changes in the seventh decimal place of the value for the bifurcation parameter. It is possible that period doubling or more complicated behavior occurs between two integer values of our bifurcation sequence. Of course, such bifurcations are probably beyond experimental observation.

F o r all the orbits studied it was found that the time from one intersection with the Poincar~ section to the next is a constant. This is apparent in the graph of x 3 for these orbits. However the physical reason that this should occur is currently unknown. The Lyapunov exponents for figs. 4 and 5 were also calculated using the established

D.G. Retzloff et aL / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

i!;

the orbits shown in figs. 4 and 5 are calculated. For the periodic orbits the fractal dimension was determined to be 1 while for the chaotic behavior of fig. 5 the fractal dimension was calculated to be 2.021. Chaotic dynamics for this reaction were also reported by Jorgensen [2]. However, the results are obtained for the situation where the activation energies 3'1 and 3'2 are considered infinite while their ratio is taken to be finite. In addition q is chosen as the bifurcation parameter. Infinite activation energies are physically unrealistic while the choice of q as a bifurcation parameter implies that several physical variables must be simultaneously changed in order to experimentally realize the predicted behavior. The dynamics investigated by Jorgensen are significantly different from what are shown here and do not exhibit the constant crossing frequency found in our results. The resuits shown in figs. 4 and 5 correspond to realistic choices for the physical variables, and changing the bifurcation parameter B while holding the other parameters fixed corresponds to changing the cooling water temperature of the process. In an experiment, this is easily realized. For the case of the parallel reactions A ~ C and B ~ D, fig. 9 shows the graph of x 3 for the stable periodic orbit arising from the Hopf bifurcation.

r

×2

F i g . 8. H o m o c l i n i c orbit.

numerical methods [13, 14, 15, 16] and are listed in table VI. For periodic orbits these are the usual Floquet multipliers. The Lyapunov exponents exhibit the sign behavior consistent with the type of orbit structure that occurs, i.e., (0, - , - ) for periodic orbits and ( + , 0 , - ) for chaotic flows [16]. From these exponents and the conjecture of Kaplan and Yorke [17] the fractal dimension for Table VI Lyapunov exponents

Exponents Period

Figure

Reaction: A --, B ~ C

Reaction: A -o C, 1 3 16 Chaos

1

2

3

B

0.0 0.0 0.0 1.4508 x 1 0 - 3

- 1.088 X 1 0 - 2 -5.116x10 -3 - 5 . 3 2 7 x 10 - 3 0.0

- 1.0907 X 1 0 - 2 -1.842X10 2 - 2 . 2 2 8 x 10 - 2 - - 6 . 9 1 5 6 × 10-2

1.002988 1.0021787 1.0017118 0.9999004

0.0 0.0 0.0 4.7419 × 10 - 2

- 1 . 1 4 6 4 X 10 - 2 - 2.7609 X 1 0 - 2 - 1.1806 × 1 0 - 2 0.0

- 1 . 6 0 3 1 X 10 - t -- 2 . 8 8 4 x 1 0 - 1 - - 2.796 × 1 0 - 1 - 4 . 0 9 2 x 10 - 1

1.001992 1.0008964 1.0004615 1.0

( m 2 = 0.0)

1 3 10 Chaos

3 4 5

B --, D 9 10 11

149

P a r a m e t e r s : q = 250; a 1 = 7.20 X 10xl; a z = 1.98 x 1 0 H ; fix = 10.0; f12 = 3,725; "/1 = Y2 = 25

D.G. Retzloff et aL/ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

150

x;

20

21

22

23

24

25

t

Fig. 9. T i m e - t e m p e r a t u r e plot of f u n d a m e n t a l period for A~C, B~D, x~=x 3-1; B=1.001992, q=250.0, a 1= 7.20 × 10 ix, a 2 = 1.98 × 10 u , fll = 10.0, r2 = 3.725, "rl = 2'2 = 25.

This figure exhibits the fundamental period. Decreasing the bifurcation parameter B yields periodic orbits with periods which are integer multiples of the fundamental period. For the parameter values studied for this reaction, the observed integer sequence is (1, 2, 3, 4, 8, 16, 32,64,128,256, chaos}. The orbits, Poincar6 sections and graph of x 3 for selected values in the integer sequence are shown in figs. 10 and 11. The orbits are similar to those found for the reaction A ~ B--* C in that the attractor evolves by the addition of extra turns in the center hole. They also share the property that the periodic to chaotic behavior has so far been numerically observed only when there exists a single unstable critical point and ceases when the bifurcation parameter has decreased sufficiently that this critical point

becomes a sink. Fig. 12 shows a blow-up of the Poincar6 section for the chaotic orbit of fig. 11. The self-similarity arising from the folding and stretching normally present in chaotic behavior stemming from the presence of Smale horseshoes [18, 19] is apparent in this figure: The general orbit structure evolution for this reaction is similar to that observed for A ~ B ~ C. Thus the conjecture of the existence of a homoclinic orbit with its associated Smale horseshoes as being responsible for the observed chaotic dynamics is applicable to this reaction. If the period three orbit is omitted, we can define a period doubling sequence for this reaction. The Feigenbaum number can then be calculated by the standard method [20, 21] and determined to be 4.66446. This is to be compared with the established Feigenbaum number which is 4.66420160 . . . . It is currently not understood why the only odd period that is numerically found is three. This may be due to the accuracy with which the various orbits can be calculated since it is not unusual to observe a change in the period when the bifurcation parameter changes in the seventh decimal place. From the work of Lynch et al. [1] it appears that the existence of period doubling may be just a fortunate situation. The Lyapunov exponents for figs. 9-11 were calculated and are listed in table VI. Their sign behavior is exactly what would be expected for the type of orbit structure that occurs. Using these results and the conjecture of Kaplan and Yorke [17] the fractal dimension for the orbits of figs. 9-11 were calculated. The fractal dimension of the periodic orbits was determined to be 1. The chaotic orbit of fig. 11 yielded a fractal dimension of 2.11. Lynch et al. [1] have also studied this reaction numerically and have observed a similar kind of chaotic behavior. They observed that the to-limit set of the dynamics appeared to be a Mobius band. However, it must be noted explicitly that the to-limit set cannot be any manifold that is topologically equivalent to R × S1. Their analysis involves the use of the next amplitude plots, which, although not stated explicitly by these authors, is

D.G. Retzloff et al. / Chaotic behavior in the dynamical system of a continuous stirred tank reactor

151

x;

~

a

10

11

12

13

14

15 '

16

1T

18

~9

20

b

4*

4

0

~Y tO

o.

0 0

0 .005

#

,00?'

.009

.011

.01 3

.015

Fig. 10. Period 16 for A - , C, B--, D; (a) 3D plot; (b) Time-temperature plot, x~ = x 3 --1; (c) Poincar~ section, Yl, Y2 are coordinates tangent to %(Dc); B = 1.0004615. Other constants same as fig. 9.

D.G. Retzloff et a l l Chaotic behaoior in the dynamical system of a continuous stirred tank reactor

152

!

0 O[3

]

T

O L)~ -

-

'O. tO

-

0.08

'

O.OS

0.04

O. Q2

a

0.02

2

0-29] 0,28

9.a75~

0.27

0.26

0.0~0" 0.965 0.0609.055. 0.950

Y2 ×3

/

0.045 0.O4O 0.035 0.030.

!

!

0.025

i

0.020

0.015 o.alo

I

I |

! 20

21

22

23

24

t

25

26

27

28

29

b

30

0.005 0.002

C 0.005

0.008

r

O.OII

t Yl 0.014

i

. . . . .

O.OIZ

i

. . . . .

0,020

¢

0,023

Fig. 11. C h a o s for A ~ C, B ~ D; (a) 3 D stereo plot; (b) T i m e - t e m p e r a t u r e plot, x~ = x 3 --1; (C) Poincar6 section, Yt, Y2 are coordinates tangent to ~rc(Dc); B = 1.0. Other constants same as fig. 9.

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

153

.g

°090°

.oo;o~

.oo~,o

.Od,,5

Y~

.ooGo

.oo~2~

.oo~,O

Fig. 12. Blow-up of fig. 11c. All constants same as fig. 11. Straight lines are drawn in to indicate the observed trend.

equivalent to choosing a subset of the surface V3 = 0 as a Poincar6 section, and then projecting the crossings onto the x3-axis to obtain a onedimensional map. In our opinion, this simplification is ad hoc and should not be expected to work in general. Furthermore, it was demonstrated in the previous sections that ,re(De) as a Poincar6 section exhibits properties that are useful in a detailed analysis, and is therefore preferable to the surface V3 = 0. Furthermore, Lynch et al. [1] choose fl2/q as the bifurcation parameter, so that their results are not amenable to experimental verification. Finally, it should be noted that the bounds for the trapping block given by propositions 2.2 and 2.3 are reasonably sharp as evidenced by figs. 5 and 11 in which the orbits approach the given bounds. Furthermore, we point out that the condition D~i~ < D b < D*~x is always satisfied, and in fact, our numerical calculations show that D b will assume values that are very close to Dn*~ or D*~x.

5. Conclusions In this work, we have shown that the dynamics of the CSTR in which the reaction A ~ B ~ C occurs share many dynamical characteristics with a CSTR in which the reaction A ~ C, B ~ D occurs. Each situation has a finite invariant trapping region in the physically accessible domain and a critical point plane in this domain which contains all the critical points of the problem. In both problems the adiabatic case leads to the critical point plane being positively invariant and hence chaos in these circumstances cannot occur. F o r the nonadiabatic case we have shown that each problem admits regions of parameter space for which there exists a single two dimensionally unstable critical point arising from a Hopf bifurcation. As a result, there is a circulating attractor distinct from the critical point. It appears from our numerical studies that the attractor responsible for the observed chaotic behavior is a strange

154

D.G. Retzloff et al./ Chaotic behavior in the dynamical system of a continuous stirred tank reactor

attractor [22]. Moreover, this attractor may be associated with a homoclinic orbit. The significance of this latter observation lies in the fact that if it can be proved to be true then a complete understanding of the chaotic dynamics and characterization of the chaos structure will be possible. From this it should be possible to explain what role, if any, period doubling plays in these problems, and the extent to which the structures of the attractors are similar, as evidenced by the fact that they both appear to evolve numerically by the addition of turns in the center hole of the attractor. There clearly are differences in the shapes of the attractors for the two problems. A complete understanding of these differences and their origin remains to be attained. Finally, it is known that multiple critical points exist in select parameter domains for both reactions [3, 4, 5]. However, a characterization of the dynamical behavior in these regions of parameter space remains to be accomplished. Thus the current results represent a beginning in the complete characterization of the dynamical behavior present in these problems.

References [1] D.T. Lynch, T.D. Rogers and S.E. Wanke, Mathematical Modelling 3 (1982) 103-116. [2] D.V. Jorgensen, Ph.D. Thesis, Dynamics and Exotic Behavior of a Stirred Tank Reactor, University of Minnesota (1983). [3] D.G. Retzloff and C. Chicone, Chem. Eng. Commun. 12 (1981) 365-376.

[4] C. Chicone and D.G. Retzloff, Nonlinear Analysis, Theory, Methods and Applications 6 (1982) 983-1000. [5] D.G. Retzloff, P.C-H. Chan, R. Mohamed, C. Chicone and D. Oflin, to appear in Journal of Mathematical Analysis and Applications. [6] C.A. Pikios and D. Luss, Chem. Eng. Sci. 34 (1979) 919-927. [7] J.E. Marsden and M. McCracken, The Hopf Bifurcation and its Applications (Springer, New York, 1976). [8] J.C. Roux, P. Richetti, A. Arneodo and F. Argoul, Chaos in a chemical system: toward a global interpretation, preprint. [9] L.P. Silnikov, Soc. Math. Dokl. 6 (1965) 163-166. [10] L.P. Silnikov, Math. U.S.S.R. Sbornik 10 (1970) 91-102. [11] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcation of Vector Fields, (Springer, New York, 1983). [12] P. Eckman, Routes to chaos with special emphasis on period doubling, in: Chaotic Behavior of Deterministic Systems, Les Houches Session XXXVI, G. Iooss, R.H.G. Hellemann and R. Stora, eds. (North-Holland, New York, 1983). [13] G. Benettin, L. Gulgani, A. Giorgilli and J. Strelcyn, Meccanica, March (1980) 9-20. [14] G. Benettin, L. Gulgani, A. Giorgilli and J. Strelcyn, Meccanica, March (1980) 21-30. [15] I. Shimadu and T. Magashima, Prog. Theoretical Phys. 61 (1979) 1605-1616. [16] A.J. Lichtenberg and M.A. Lieberman, Regular and Stochastic Motion (Springer, New York, 1983). [17] J.L. Kaplan and J.A. Yorke, Chaotic behavior of multidimensional difference equations, in: Functional Differential Equations and Approximation of Fixed Points, AI Dold and B. Eckmann, eds., Lecture notes in Math., vol. 730 (Springer, New York, 1979), pp. 204-237. [18] S. Smale, Bull. Amer. Math. Soc. 73 (1967) 747-817. [19] J. Palls Jr. and W. deMelo, Geometric Theory of Dynamical Systems (Springer, New York, 1982). [20] M.J. Feigenbaum, J. Stat. Phys. 19 (1978) 158-185. [21] M.J. Feigenbaum, Los Alamos Science 1 (1980) 4-27. [22] J.C. Roux, R.H. Sinoyi, H.L. Swinney, Physica 8D (1983) 257-266.