- Email: [email protected]
Structure of the steady-state solutions of lumped-parameter chemically reacting systems
STRUCTURE OF THE STEADY-STATE SOLUTIONS LUMPED-PARAMETER CHEMICALLY REACTING SYSTEMS
OF
VEMURI BALAKOTAIAH and DAN LUSS* Department of Chemical Engineering, University of Houston, Houston, TX 77004, U.S.A (Received
9 October
1981; accepted
26 May
1982)
Abstract-A new, powerful mathematical technique enables a systematic determmation of the maximal number of steady-state solutions of lumped parameter systems in which several chemical reactions occur simultaneously. The method can predict also the different types of diagrams describing the dependence of a state variable of the reactor on a design or operating variable. The technique is applied to several reaction networks giving new results and insight. For example, it is proven that when N independent, parallel exothermic reactions with equal and high activation energies occur in a CSTR there exist N ! distinct regions of parameters in each of which 2N + 1 steady-state solutions exist.
INTRODUCTION
reactions
in a large number of chemically reacting systems. The multiplicity features of systems in which a single chemical reaction occurs have been studied extensively and are well understood by now[l-51. However, most practical control and start-up problems due to steady-state multiplicity are encountered in systems in which several chemical reactions occur simultaneously. Mathematical difficulties have prevented so far the analysis of the qualitative features of these practical cases. to know the In many applications it is important behavioral features of the reactor, such as the number of possible steady-states and the inIluence of a change in one or more operating variables on these states. Unfortunately, the techniques used so far to answer these questions for the single chemical reaction case cannot be extended directly to multireaction systems. We present here some new, powerful mathematical techniques which enable a systematic determination of the above mentioned features of lumped-parameter systems in which several chemical reactions occur simultaneously. We apply these techniques to several types of chemical reaction networks in order to determine: (a) The maximal number of possible steady-state solutions and the parameters for which they exist; (b) The different types of bifurcation diagrams which describe the dependence of a state variable of the system (say the temperature) on a design or operation variable. The analysis is of lumped-parameter systems such as a continuously stirred tank reactor (CSTR), a non-porous catalytic pellet or wire and a porous catalytic pellet with no intra-particle gradients. To simplify the treatment of multireaction systems we assume that all the reactions are of first-order, so that the steady-state behavior can be described by a single nonlinear algebraic equation. The full power of the new technique is illustrated by its ability to predict all the multiplicity features of a CSTR in which an arbitrary number of independent parallel Steady-state
multiplicity
is encountered
*Author to whom correspondence CBS “a, 37, No. ,,--8
should be addressed.
occur,
a feat which could not be accomplished used previously in the chemical engineering literature to analyze similar problems. Due to its power and ease of application, we expect this technique to become a standard tool for determining the multiplicity features of chemically reacting systems and of other systems described by a single algebraic equation. The analysis is based mainly on the singularity theoryIf%141 with which most chemical engineers are not familiar. The possibility of applying these techniques to the analysis of chemically reacting systems was brought to our attention by [IS]. A formal, rigorous survey of the subject is rather lengthy and abstract. Thus, we present in the next section a simple, heuristic description of many of the mathematical terms and concepts. The reader interested in a more formal approach is referred by any of the techniques
to 1133. MATHEMATICAL BACKGROUND
We consider a nonlinear lem) of the type
equation
F(x, h, p) = 0 x is temperature and p is a theory can vectors we the moment
where
(bifurcation
prob-
(1)
a state variable (such as the steady-state or conversion), A is a bifurcation parameter constant vector of parameters. While the handle cases in which both x and F are do not consider here this general case. For we regard p as fixed and rewrite eqn (1) as F(x, A) = 0.
(2)
We select the origin such that
F(O,O)=0
(3)
and define F to be smooth if it has continuous derivatives of all orders. The graph of x vs ,+ which satisfies equation (1) for a fixed p is defined as a bifurcation 161 I
1612
V. BALAKOTAIAH and D. Luss
diagram. A local bifurcation diagram describes the graph of x vs A in a small neighborhood of some point. A global bifurcation diagram describes this graph for all values of x and A within the domain of interest. Two bifurcation diagrams F(x, A) = 0 and G(y, CL)= 0 are cont(lct equiualent if one can be transformed into the other without changing any of the qualitative features in the neighborhood of the origin by a smooth and invertible change of coordinates (4)
(Y>P)-+ MY. CL).A(P)) and a scaling
factor
T(y, p). Thus,
G(Y, ccl = T(Y, ~r)F(x(y,
P), A(F)).
(5)
When this smooth change of coordinates is used the number and order of solutions of x for each fixed A in the neighborhood of the origin is the same as the number and order of solutions of y for the corresponding p if[13]
is unstable
since the bifurcation
diagram
F,(x,h)ex’-
of (\rl.=%1)
assumes different forms depending on whether l >O or +
=0
(13)
is an unfolding of Fz, but it is not versa1 (since it cannot describe the diagram corresponding to E < 0 in Fig. 1) while FI=x”-~,~2-~L~-A=0
x(0,0) = 0,
h(O)=O;
(6)
~(O,OPO,e (0) > 0. the bifurcation
problems
F(x,h)bx*-A=0
(7)
and G(y,p):l-CL-cosy=O. are contact
equivalent,
(8)
since next to the origin
G(Y, Y) = 2F (sin r/2,@) and F(x, A) = OSG(2 sin-’ x, 2A).
(I4 4 1)
F(O, 0, P? = 0
$
CO,%P9 = 0,
(15)
i= 1,2 ,...,r
(16)
and A 5 $$
(O,O, p‘? $
then in the neighborhood
(0, 0, P’? + 0,
of p”, F(x, A) is contact
(10)
the qualitative features of the bifurcation diagram arc unchanged. On the other hand, the bifurcation diagram of F&h):?-A=0
is versa1 but not universal. Only F, is universal, and the codimensiou of F, is one. The bifurcation diagram of eqn (1) is stable for most values of the parameters vector p. For example the bifurcation diagram of F, is stable for all values of l except zero. A set of p values exists at which the function defined by eqn (1) has the most degenerate singularity. Near these points several qualitatively different stable bifurcation diagrams exist. These singular points are characterized by the vanishing of a finite number of the first derivatives of F with respect to x and A. There exists always a polynomial equation which is contact equivalent to F next to these singular points. We shall use in this work the following result which is a special case of a more general theorem proved in [131. Suppose that the function F given by (1) has a singular point p” at which
(9)
Thus, the local qualitative features of a complex nonlinear equation can be determined by the analysis of the features of a much simpler contact equivalent polynomial equation. The next important concept in singularity theory is that of stability. The bifurcation diagram of F(x, A) is defined to be stable if F(x, A)+ g(x, A, c) is contact equivalent to F(x, A) in the neighborhood’of the origin for all sufficiently small smooth functions g(x, A, s). For example, the bifurcation diagram of F(x, A) = 0 given by eqn (7) is stable in the sense that by perturbing F(x, A) to F&C, A) ” x2- a - A = 0,
(14)
T(0, 0) f 0
and
For example,
(12)
(11)
C&
C-0
c>o
Fig. 1. The bifurcation diagram of eqn (12) (E = 0) and the two corresponding perturbed diagrams (et 0).
Structure of the steady-state solutions of lumped-parameter chemically reacting systems equivalent
1613
to G(x,‘)=~‘+‘i’
whose universal
(
unfolding
- sign if A < 0 +si~“if~>o)
(17)
is
G(x, A, n) =x’+’ - a,_,~~-’ - LY,_-~x~-*--. .
- (Y,XT A. (1%
This theorem is most useful as it proves that if such a singular point exists, the maximum number of solutions of eqn (1) next to the singular point is r + 1 and it also enables finding all the possible local bifurcation diagrams. We present now some results useful for a global analysis of the features of the steady-state equation. We define a solution (x. A, p) of eqn (1) to be a singular point if z
(x, A, P) = 0.
(19)
The set of all singular points in the (x, A, p) space is called the singular set. The projection of the singular set onto the (A, p) space is called the catastrophe sot. A point (x, A) on the bifurcation diagram of eqn (1) is called a bifurcation point if the number of solutions of F(x, A) = 0 near x is not constant as A varies in an arbitrarily small neighborhood. A limit point is the least degenerate singular point, i.e. F = F, = 0 while F, and F, do not vanish. Using the implicit function theorem, it can be shown that a bifurcation point is a singular point. However, the converse is not true in general. Thus, the number of solutions of eqn (1) can change only when the parameters cross the catastrophe set in the (A, p) space. It is proven in [13] that the qualitative features of a bifurcation diagram of (1) change when the parameters p cross one of three surfaces. The first, called the hysteresis variety (If) is the set of all points in the parameter space p satisfying 2 F(x,A,p)=~(x,h,p)=~(x,A,p)=O. (20) Elimination of x and A from these three equations gives a single algebraic equation in p, defining a hypersurface. Similarly, the set of all points in the parameter space p satisfying F(x, A, p) = g
(x. A\,p) =
$
(x. A. p) = 0
(21)
is called the isola variety (I), while the set of points satisfying the four equations
p
F(x,, A, p) = F(xz, A,p) = 0 ~(x,,A,P)=$(x,,A,P)=O;
xlzx2
(22)
is called the double limit variety (DL). When p crosses the hysteresis variety, typically two limit points appear or disappear and the local nature of the bifurcation
Fig. 2. Possibli transitions of a local unstable bifurcation diagram (middle column) at the Hysteresis (A), Isola (B, C) and the Double Limit varieties (D, E) into one of two stable forms (left or right column).
diagram changes as shown in Fig. 2(A). When p crosses the isola variety, typically two limit points appear or disappear, so that either the bifurcation diagram separates locally into two isolated curves (Fig. 2C), or one isolated curve disappears (Fig. 2B). When p crosses the double-limit variety, the number of limit points does not change but the relative position of the limit points changes as illustrated by Figs. 2(D) and 2(E). When for no feasible p value a limit point can exist on the boundaries (extreme values) of the feasible x, A values, the features of the bifurcation diagrams can change only by crossing one of the three varieties. In these cases the three varieties divide the global parameter space into several regions, in each of which a different type of bifurcation diagram exists. Note that when eqn (1) satisfies the condition dF/dA f Ofor all feasible x and y the isola variety does not exist. APPLICATION TO
LLIMFSLB PAMMETER CAEMlCALLYREACTINGSYgTlMs
We apply in this section the mathematical techniques described in the previous section to predict the multiplicity features of a lumped parameter system such as a continuously stirred tank reactor (CSTR) in which several types of chemical reactions occur. To familiarize the reader further with the mathematical tools we start with two simple cases of a single exothermic first and zeroth-order reaction. We then illustrate the full power of the new technique by considering the more complex cases of two consecutive exothermic reactions and of N independent, parallel reactions. (1) A single first-order reaction in a CSTR Consider a cooled CSTR in which a single, exothermic, first-order reaction A-B occurs. The corresponding species and energy balances are q(C,qpc,(T,-
C,) - Vi(T)&
T)+ V(-AH)&IJ&--
= 0
Ua(T-
(23) Tc)= 0
(24)
V. BALAKOTMMI and D. Luss
1614
The corresponding
where f&T) = ff(T,,,)exp
[i
(+-J$)].
(26)
where H = Ualqpc,,
a
(2%
Equations (23) and (24) can be combined to give a single steady-state equation for the dimensionless temperature F(Y, D, B) 5 Y - D(B - Y) exp (Y) = 0
hysteresis
variety
defined by (20) is
=o.
(35)
The isola and the double-limit varieties do not exist. Thus, Q = 0 divides the one dimensional (Yspace into two regions. For all (Y~0, the bifurcation diagrams of x vs A are monotonic while for any (Y>O the bifurcation diagram has an S shape. Thus, next to the singular point ( IID, Do, Bo) eqn (26) has only two types of bifurcation diagrams (Y vs D or Y vs B). The catastrophe set of the steady-state eqn (26) is given by
T, = (To + HTc)/(l + If)
F=E=, aY
which upon elimination
Y =&-(l-+).
(27)
We wish to find the maximum number of solutions and the various types of bifurcation diagrams of Y vs D or Y vs B for eqn (26). We note that the set of equations d2F F=n=ay’=O
B5 + ~(B(B D = B v’(B(B
of Y gives
- B T - 4)) exp 1 2 t/(B(B
4)) 1, (37)
If we consider D as the bifurcation parameter, teresis variety defined by (28) reduces to B = 4.
the hys-
(38)
aF
has a unique
solution Y”=2,
B”=4,
D’=e-‘,
(29)
At this point a’F F=_E=’
aB (31)
It follows that a unique singular point exists in the space (Y, D, B) at which the bifurcation problem (26) is contact equivalent to G(x, A) k x3 - A = 0, whose universal unfolding is G(x,h,a):+CCC-h=O.
(32)
The isola and double-limit varieties do not exist in this case. Since no limit point exists on the boundaries of the feasible Y, D region, the nature of a bifurcation diagram can change only at an intersection with one of the varieties. Thus, in this case, Y is a monotonic function of D for all B <4, while the bifurcation diagram (Y vs D) has an S shape for all B > 4. We conclude that in this case the two local bifurcation diagrams close to the singular point are the only possible diagrams in the global parameter space. If we consider B as the bifurcation parameter, the H variety is defined by D = e?. Thus for ah D > em*, Y is a single valued function of B and for D < em2the bifurcation diagram (Y vs B) has an S shape. Figure 3 shows the surfaces represented by eqns (26) and (32) and the corresponding catastrophe sets.
The singular point of F at (Y”,Do, By or of G at the origin is called a “cusp singularity”. Since G = 0 has either one or three real solutions the same is true for F = 0. The nature of the local bifurcation diagrams of (26) can be obtained by studying the bifurcation diagrams of (32). The catastrophe set of (32) is obtained by a simultaneous solution of
G=~=() ax
giving 4a3 = 271=.
(34)
Fig. 3. The surfaces represented by eqns . (32) and (26) and the corresponding catastrophe Sets.
Structure of the steady-state solutions of lumped-parameter
chemically reacting systems
1615
(2) A single zero&order remlion in a CSTR The steady-state conversion in a cooled CSTR in which a zeroth-order reaction occurs is determined by the equation exp
F(5, Da, Y. B) ” E-Da
(&
) =O
(39)
where 5 = (CArI- C.4)/CAll
DO
Da = Vk(T,)/qC,
(W
y = E/RT,,,
p = (- AH)C~/pcpTm(l
+ H)
and T,,,, H are defined by (27). It should be noted that the zeroth-order rate expression always breaks down at high conversions, and for any set of y and p some critical Da exists above which the steady-state equation has no feasible solution (Oi 5 < 0. In order to determine the maximum number of steadystate solutions and the different possible bifurcation diagrams (t vs Da) we find first the singular points of (39). This is accomplished by noting that the set of equations F=aF=aZF=O a[at has a unique
solution 5” = l/,9.
Fig. 4. Different possible global bifurcation diagrams for a zeroth order reaction. Thus, for any fi > 1, .$ is a monotonic function of Da if y < 4 and a cusp exists if y > 4. For parameters far away from the cusp point, the nature of the bifurcation diagrams depends on the intersection between the surface represented by (39) and the boundary of the feasible region (t = 1). For example, for any p > 1 the multiplicity pattern for y slightly larger than four is 1-3-1-O. However, as y is increased the following multiplicity patterns are obtained: l-3-2-0, 1-2-O. The corresponding bifurcation diagrams are shown in Fig. 4 as cases b, c and d. For /3 < I, the cusp singularity defined by (42) is outside the feasible region. However, the set of equations
(41)
for any p, i.e. Da’ = e-‘JP
yD = 4,
(42)
has a feasible
solution
+r-2-t/trtr-4) 33
and that
’
Do’=Pexp(&)
(48) $
(5”, Da’, y”, P) g
(5”. Duo, y”, P) (0.
(43) if 5” < 1. This occurs
Thus, in the (6 Da, y. p) space, the surface represented by (39) consists of a one- parameter family of cusp singularities. However, p is in the feasible region only if p > 1. Therefore, for any p > 1, a unique singular point exists in the neighborhood of which the bifurcation problem (39) is contact equivalent to (32). It follows that eqn (39) can have at most three steady-state solutions in the feasible region provided p > 1. The local bifurcation diagrams for this case are the same as that for eqn (32). We find now the global bifurcation diagrams. The catastrophe set is given by
if -r ’ (I+ P)‘/P.
(49)
It can be verified easily that
t$ th”,Da’,
Y.
8) g
(5”. ho, y. PI > 0.
(50)
Thus, for every set of g < 1 and y > (1 + p)‘/p a unique singular point defined by (48) exists. At that point the bifurcation problem (39) is contact equivalent to G(x, A) : x2 + A = 0.
(51)
Do=.$*exp(a) of G at the origin is known as the “fold” theory[l2] and as the limit point in Therefore, when p < 1, the maximum number of steady-states is two. Next to the singular point eqn (39) has two or no solutions. Figure 4(d) shows the corresponding bifurcation diagram. It is seen that for all The singularity
where
in the catastrophe bifurcation theory.
5_=~-2rd(~(~-4)) 2P The corresponding
hysteresis
variety
y = 4.
. is
(52)
V. BALAKOTAIAH and D. Luss
1616
a unique solution exists. For all Da in (z, Da@) two solutions exist, while for Da > Da’ no feasible solution exists. When p < 1 and (49) is not satisfied a unique feasible solution exists for 0 < Da < 5 and no feasible solution exists for Da > z (Fig. 4a). When this analysis is repeated taking y as the bifurcation parameter, bifurcation diagrams similar to those shown in Fig. 4 are obtained. (3) Two consecutive exothermic reactions Consider a CSTR in which two consecutive
reactions
. ^ AfiB::C
We rewrite
eqn (57) in the standard
F(B,Da,,u,B,,B3=Da,X(B,-e-be)
t vDuI*X’(B, + Bz- 0) - 6 = 0. taking Da, as the bifurcation parameter. Appendix A that the set of equations
has a unique
The species
and energy
Da,‘=
q(CAo-
CA) = V&(T)Ca
q(C*-
Ce) = vUz3
qpc,(T-
balances
(2 tv'3)e-',
are (53) (54)
+(-AH,)V&T)CB
(*)I.
(56)
With these assumptions eqns (53)-(55) can be combined to give a single equation for the dimensionless steadystate temperature
+ (1
B,” = 6d3(2
vBzDa12X2 t vDa,X)
t Dn,X)(l
(57)
which is in the physically point
feasible
=
Vlf~(Tm) 4
e=&(T-
’
Bz = B I(- A&)/(v = &Tm)I~dTm)~
(61)
region. At this singular
(62)
It follows that a unique singular point exists in the feasible region at which the bifurcation problem (59) is contact equivalent to G(x, A) A x’- A = 0, whose universal unfolding is C(r, A, a) : x5 - ax?-
(1*x2- a,x -A
RTm
= 0.
(63)
The singularity of G at the origin is known as the “butterfly”. The surface represented by (63) in the (x, A, a) space is called the butterfly catastrophe and its properties are reported in detail in the mathematical literature [S, 121. Equation (63) has at most five solutions. This implies that the steady-state equation, which is contact equivalent to G, has five solutions for some p. Now, instead of finding directly all the local bifurcation diagrams of the steady-state equation (59) next to its singular point (defined by (61)), we tackle the simpler equivalent problem of finding all the bifurcation diagrams of G. This is accomplished by finding the various regions bounded by the hysteresis and double-limit varieties in the ((I,, (Y?,eyj) space: The hysteresis variety is obtained by eliminating x and A from (63) and ~=5*‘-3a,*‘-2a2x-a,=0
where Da,
- t/3)
$&=-24(2+3)e'
A detailed analysis of the corresponding multiplicity and multiplicity criteria was presented patterns elsewhere[l6]. The mathematical techniques presented in this work can be used to determine directly the maximum number of the steady-state solutions, and all the possible bifurcation diagrams. For illustration purposes we consider only the case of two exotbermic reactions. To simplify somewhat the algebraic manipulations we assume that the activation energies of the two reactions are identical and that no B is present in the feed. Moreover, we assume that the temperature dependence of each rate constant can be described adequately by the positive exponential approximation [ 171, i.e.
B,Da,X e = I+DalX
Y0 = (2 - d3)2
(9.5)
- Ua(T - T,).
/&(T) = f&T,,,) exp.[&
in
#=3
- V&(T)CA
T,)=(-AH,)Vk^,U-)C.
It is shown
(59)
solution
B,“ = 12(2 - g3), occur.
form
(64)
Tm) 2
T.v.
The elimination of x and A from eqns (63)-(65) is equivalent to eliminating x from (64) and (65) and the hysteresis variety can be expressed in the parametric form
AHI) X = exp (0).
= 20x3 - 6as.x - 2~~2= 0.
(58)
IY, = 3cQx*- 15x’ (12= 10*3-3u3x.
(-m
(66)
Structure of the steady-state solutions of lumped-parameter chemically reacting systems
aI
1617
aI
I
-H ------CL
VARIETY WMETY
I
x c
b
b
0
CY,CO
LY3=-0
Fig. 5. Classification of the x vs A bifurcation diagrams of the butterfly singularity defined by eqn (63).
Similarly, the double-limit can be expressed as
variety,
defined
by eqns (22),
t(4t2+7t +4) (x1 = (3tZ+ 4t + 3)2 a32P &4(1+
o’u +3r+ t3’
a33
(3P + 4t + 3)3 with - m < f i m being a parameter. These two varieties are shown schematically in Fig. 5. The hysteresis variety divides the a space into three regions. In the first region no limit points exist and x is a single valued function of A shown as case n in Fig. 5. In the second region two limit points exist and the bifurcation diagram has the S shape shown as case b in Fig. S. For all a in the third region, four limit points exist. Here, five types of bifurcation diagrams may exist, shown as cases c-g in Fig. 5. The double-limit variety separates between the five regions corresponding to each of these diagrams and in applications it is important to know these transitions. One can construct all the possible bifurcation diagrams in the following way[13]. Let (xi, Ai), i = 1, 2, 3, 4 be the four limit points corresponding to an c~ value in region 3. Without loss of generality assume that x, < xz < xa i x+ Since dx/dh is positive for h + + m the Ai must satisfy the restrictions
(9
Al>&
(ii)
A,>h+
(iii)
A,>A,
Thus, the A,‘S can be arranged only in one of the following five possible sequences (shown in Fig. 5). (c) (d) (4 0-I (Ed
Bl Fig. 6. Global classification of the 0 vs DOI bifurcation diagrams for two consecutive exothermic reactions with no B in feed. Close to its unique singular point (eqn (61)), the steadystate equation F = 0 has seven bifurcation diagrams similar to those of G. A detailed procedure for constructing the hysteresis and double-limit varieties in the global parameters space is reported elsewhere[lS]. Figure 6 summarizes the results of the analysis showing that four different cross sections of these varieties with the v-B, plane exist depending on the values of Bz. The hysteresis variety intersects the B, axis always at B, = 4. A monotonic 8 vs Da, bifurcation diagram, similar to case a in Fig. 5, exists in the dotted region to the left of the hysteresis variety. When B+ 14 this region is outside the positive quadrant. An S type bifurcation diagram similar to case b in Fig. 5 exists in the hatched region to the right of the hysteresis variety. The pocket shaped region. which exists for all B, > B,” is divided by the double limit variety into five regions. The corresponding bifurcation diagrams have four limit points and are similar to cases c-g in Fig. 5. For all B, > 27/8 this region is stretched out to unbounded B, values. We report in Appendix B one set of parameter values corresponding to each of the bifurcation diagrams as well as the corresponding limit points.
Consider a lumped-parameter system such as a nonporous catalytic wire or pellet or a cooled CSTR in which N independent, parallel reactions AiL
Pi,
i=l,2
,...,
N
occur. We shall determine the maximum number of possible steady-state solutions, the number of singular points or the number of s’eparate parameter regions in which this number of solutions exists, and the number of all the possible local and global bifurcation diagrams. To simplify the algebraic manipulations we assume that the activation energies of all the independent reactions are
1618
V. BALAKOTN.~~ and D. Luss
identical and sufficiently large so that the positive exponential approximation is valid. The impact of these assumptions on the multiplicity features is examined in Appendix C. With these assumptions the dimensionless steady-state temperature is given by[18]
paring
eqns (74) and (76) we note that zi = xi
17i= wi,
(77)
(68)
is a solution of (74). Since the weights w, and the zeros xI are uniquely determined for any fixed N, it follows that (77) is the unique solution of (74). Substituting (77) in (75) we get
(69)
Thus, the solution
where Da_
=
vm”J
B_
=
jS_-
AfWL,
’ RT,,, ~c,Tm
’
4
and the other parameters
are as defined before. y = e-e
eqn (68) can be written
of eqns (71) and (74) is
Setting (70)
BO = WJXi(1 -xi)
as
y’=exp[-22
PJ Da.B. F(Y, Dai, Bi) k In Y + z, (y
= O.
(71)
We consider one of the Damkijhler numbers as a bifurcation parameter. (Note that the change of variables defined by eqn (70) transforms an S diagram of 0 vs Da; into an inverse S diagram of y vs Da,. However, the coordinates of the singular point are unaffected by the substitution.) Setting the first 2N derivatives of F with respect to y to zero gives
(79) l/j]
The Gauss-Legendre quadrature integrates exactly any polynomial of degree (2N - I), but not one of degree 2N. Thus
a ‘,+:F ay~~
Moreover,
I
= (2N)!
j-zm [ I- (79 + 0 ,z, wiZN] i 0.
at the singular
point defined by eqns (79).
&>O.
N
,z, +“;f$_+L Yrn-; =0,
(80)
(81)
(y
m=1,2
,...,
2N.
(72)
Defining
Thus, at the singular point, the steady-state is contact equivalent to xzN+‘* A whose folding is G(x, A, a) k xZN+’ ~ oa,._,xzN
zi = __ Y
y+Dai
eqns (72) and (7 1) can be expressed
(73)
as
m=l,2,...;2N
(74) (75)
According
to the Gauss-Legendre
quadrature
-.
. .-a,x+A
equation universal
(68) un-
-’ - (I~+,-~x~~~*
=o.
(82)
(The sign of h in (82) is opposite to that of the derivative in (80)). We conclude that there exists always a set of B, and Do, close to the singular point for which 2N+ I solutions exist. Setting dfldy to zero we obtain a polynomial of degree 2N in y. Thus, the maximum number of solutions in the global space is also 2N + 1. We note that if (y. Bi, hi)
= [Y'.
$?q
+
Y’]
i=1,2,...,N (76) where xi are the zeros of the Nth degree Legendre polynomial and wi are the corresponding weights. Com-
tThe equality of the rhs of (78) was proven by _I.Villadsen and S. Wedel.
is a solution of eqns (74) and (75) then because of the symmetry of the problem every permutation of the (St, Dai) i = 1,2,. . , N is also a solution, and there exist exactly N! distinct singular points. We conclude that there exist N! distinct parameter regions in each of which the sfeady-state equation (68) has (2N + 1) solutions. To illustrate this result consider a system with three independent parallel reactions. Here, the following six singular points exist:
Structure of the steady-state solutions of lumped-parameter chemically reacting systems Da,’ = y”
Da,0 = (4k $q15))y0
Da,’ = (4 5 q(lS))y’
Da,“=(4-c~/(15))y0
Do,” = y”
Daz” = (4 T q/(lS))y’
Du,~=(~s~/(~S))~~
Da>’ = (4 ri:d(l5))yO
Da3’ = y”
1619
(83) B,’ = 1619
B,’ = 25/9
Bz” = 2519
BzO = 16/9
Bz” = 2519
B,‘=
Bs’ = 25/9
B,O = 16/9
2519
where y” = exp (-11/3). The nature of the local bifurcation diagrams of the steady-state equation can be determined by finding the bifurcation diagrams of the polynomial equation (82). Unfortunately, because of the large number of parameters it is not possible to obtain a clear picture of the number of regions formed by the hysteresis and double-limit varieties in the (2N - 1) dimensional a space. The hysteresis variety separates regions of the parameter space which differ by two limit points. Next to every singular point (defined by (79)) there exists a region in which 2N + 1 solutions exist, i.e. one which has 2N limit points. Thus, the hysteresis variety divides the a space into no more than N + I regions. The bifurcation diagrams in region j(j = 1,2,. . , N + 1) have 20’ - 1) limit points. For j 2 3, each region is divided into several subregions by the double-limit variety. One can determine all the possible local bifurcation diagrams using the constructive scheme described in the previous example. For example, if no limit points exist 0 is a single valued function of Da,. When two limit points exist an S pattern is obtained. When four limit points exist the five local diagrams shown as cases c-g in Fig. 5 are possible. Similarly, it can be shown that for six limit points 61 different local diagrams exist, etc. Thus, for two parallel reactions, the total number of possible local bifurcation diagrams is 1 + 1 + 5 = 7 (Fig. 5). while for three parallel reactions, there exist 1 + 1+ 5 + 61 = 68 different possible local bifurcation diagrams. The above example illustrates the complexity of mapping all the bifurcation diagrams for a system in which many reactions occur simultaneously. The global bifurcation diagrams of f3 vs Dai can be found by constructing the hysteresis and double limit varieties in the (2N - 1) dimensional space of Bi (j = 1,2 ,.__, N) and Da! (j= 1,2 ,__. ,N, j#i). Different cross sections of the bifurcation set in the Dal-Da2 plane for the case of two parallel reactions were shown in [ 161. These figures show that a 1-3-5-3-5-3-1 multiplicity pattern can be obtained when one of the Damkahler numbers is varied while the other is kept constant. This indicates that the local and global bifurcation diagrams need not be the same when one of the Damkiihler numbers is taken as the bifurcation parameter. A common way of changing the Damkahler number is by adjusting the residence time in the reactor. Consider an adiabatic reactor in which the flow rate is changed in a continuous fashion so that Dai = viDa,,
i=1,2
,.__,
N.
(84)
B la = 25/9
The steady-state
equation
(68) can be rewritten
as
The local bifurcation diagrams of B vs Da, of eqn (85) are the same as those of (82) and there exist again N! singular points. Surprisingly, the local and global biiurcation diagrams are identical in this case. The reason is that eqn (85) can be written as an Nth degree polynomial in Da,e’ whose coefficients are linear functions of 8. Differentiation of the polynomial with respect to 0 gives another polynomial of Nth degree whose coefficients are also linear functions of 0. Elimination of Dales from these two equations gives a 2Nth degree polynomial in 8. The zeros of this polynomial are the limit points in the global parameter space. It follows that any bifurcation diagram in the global parameter space can have at most 2N limit points. However, all the bifurcation diagrams containing up to 2N limit points exists next to each of the N! singular points. This point is illustrated in Fig. 7 which shows the four different types of cross sections the hysteresis and double limit varieties have with the vrBl plane. This figure is symmetric about the line YZ= 1 in the B,-log v2 plane. Bifurcation diagrams with zero and two bifurcation uoints exist in the dotted and hatched regions. resneckvely. Five different types of bifurcation diagrams with
Bl Fig. 7. Global classification of the 0 vs Da, bifurcation diagrams for eqn (85) with N = 2.
1620
V. BALAK~TAIAH and D. Luss
four points exist in the two pocket shaped regions. Thus, there exist seven different bifurcation diagrams in the global space, which are identical to those existing next to each of the two singular points.
DISCUS!3lON
This work illustrates the power of the singularity theory to predict the number of local steady state solutions. The theory is a local one and additional analysis is needed to determine the maximal number of solutions in the global space. It is proven in [ 181 thatfor all the examples presented here the maximal number of solutions next to the highest order singularity is equal to that in the entire global space. It is not yet clear what properties of the steady-state equations are responsible for this result. In all the examples presented here we analyzed the influence of a single dimensionless group, usually the Damkiihler number, on the steady-state temperature or conversion. The value of this dimensionless group can be varied by changing the value of the rate constant at the reference temperature. In practice, one may want to examine the influence of other parameters, such as the reactants concentration, feed or coolant temperature and residence time (flow rate) on the steady-state temperature. To determine the bifurcation diagrams for any of these cases, the dimensionless steady-state equation should be rewritten so that changes in the specific variables affect only one of the dimensionless numbers. This reparameterization does not change the number of solutions but may increase the number of dimensionless groups appearing in the steady-state equation, and is essential for an efficient determination of the singular point and the corresponding contact equivalent bifurcation problem. For example, in order to determine the bifurcation diagrams of the steady-state temperature of a CSTR in which a single first-order reaction occurs vs the coolant or feed temperature, we should introduce the following parameters
versa1 unfolding is given by (32). Thus, the local bifurcation diagrams in this case are identical to those obtained by changing the parameter D in eqn (26). In all the examples considered so far the derivative of F with respect to the bifurcation parameter did not vanish. Consequently, the hysteresis and double-limit varieties describe all the transitions between the various bifurcation diagrams, which are either single-valued or have one or more hysteresis loops of an S shape. Other types of bifurcation diagrams are possible when the derivatives of F with respect to the bifurcation parameter vanishes at the singular point so that an isola variety exists. This has been illustrated in [15, N-211 by analyzing the influence of changes in the residence time on a nonadiabatic CSTR in which a first-order reaction occurred. A steady-state equation such as (1) describes the dependence of a state variable x on a bifurcation parameter A and a finite set of n other parameters p. Thus, it describes a surface in an n +2 dimensional space. While this surface is continuous and smooth its projection into the (~\.p) space may have certain singularities. The examples discussed in this paper show that the number of local solufions to (1) as well as the di$erent possible bifurcation diagrams can be determined if the highest order singular point on the surface is known. All the examples discussed above show that the highest possible order singular point is characterized by the vanishing of a maximum of (n + 1) derivatives. Consider a steady-state equation F, whose derivative with respect to the bifurcation parameter does not vanish, and for which the highest order singularity is characterized by the vanishing of n + 1 - i (i 2 0) derivatives so that
+E=a*F_ 2 ax
and at the feasible
solution a”“-%‘ ax”
o = RTJE,
6 = RTIE,
s = &(-AHKb E pc,(l +W
flu = VkO/q,
(86)
k(T) = k. exp (- EIRT) and use instead tion
of (26) the following
steady-state
equa-
F(&, (T,fie, p, ” B - v - &(jZ + o - ,_)e_“i = 0
(87)
where (r is the bifurcation parameter. This equation for every b > 0 a unique singular point at
has
i”=
B+t/@(l+B)) 2llf
a + X4&1 + PM
2 = wm1+ Da”=(8’-a4exp(ll~9/(B+ and is contact
equivalent
s,, - 81/2 o”-
(88)
- ---
i-7-i
=
a““_‘F ax”+,_i
uni-
0
(89)
of (89) aF x + 0.
In order to determine the value of i we attempt to solve the set of eqns (89) starting with i = 0. If a physically feasible solution exists then the singular point is of order n + 2 and F = 0 has at most n + 2 local solutions. All the examples in the previous section with the exception of the zeroth order reaction belong to this case. If no feasible solution exists for i = 0 one attempts to find a solution for i = 1 and so on. One should note that when i t 1 the solution of (89) consists of n + 2 - i varidetermined as function of ables Cc, A, pi, pZ, , . . , p,-,) the remaining i parameters @.+,_i,. _ , p,). If, for some physically feasible values of these i parameters, the solution of (89) is feasible, then the maximal number of local steady-state solutions in this case is (n + 2 - i). The example of a zeroth-order reaction in a CSTR, discussed in the previous section, illustrates cases for which n = 2, i = 1 and n = 2, i = 2. In fact, by defining
t?‘)
to G(x, A) = x3 - h whose
=
DaP =
d
rSI=g 1+ PS
Structure of the steady-state solutions of lumped-parameter eqn (43) can be written
chemically reacting systems
as
y = L5(y - P) exp (P)
1621
NOTATION
(92)
which is identical in form to (26). However we have chosen the form of (39) as it enables a more direct apphcation of the feasible boundaries. In certain situations it may not be possible to identify the singular point because of algebraic difficulties. One can still find all the bifurcation diagrams by constructing the hysteresis, double-limit and isola varieties in the parameter space determining one diagram in each region bounded by these surfaces. Finally, we wish to point out that the analysis in the previous section can be extended to systems for which the steady-state is described by a vector of m state variables which arc defined by a set of m equations. This situation is likely to be encountered in the description of the steady-state of a lumped-parameter system in which several reactions with complex rate expressions occur. CONCLUslONS This work illustrates that singularity theory is the most efficient and systematic tool for ascertaining the maximal number of steady-state solutions of lumped parameter systems in which several chemical reactions occur simultaneously, and for a direct determination of the regions in the parameter space corresponding to these solutions. The combined use of the singularity and bifurcation theory enables us to determine some important local as well as global properties of the system, specifically to find out all the possible types of bifurcation diagrams and the boundaries of the region in the parameter space corresponding to each fo the diagrams. The construction of the hysteresis, double-limit and isola varieties is the most efficient direct method of determining the parameter regions corresponding to each type of bifurcation diagram. Knowledge of the bifurcation diagram enables the prediction of the influence of a change in a design or operating variable on the structure and number of steady-state solutions. In practice, one needs to know in addition to the number of solutions the type of corresponding bifurcation diagram. For example, consider a case of two consecutive reactions in a cobled CSTR for which three steady-state solutions exist for some set of parameters. This situation must correspond to one of the six diagrams shown as cases b-g in Fig. 5. Obviously, it is important to know to which of these six diagrams this case corresponds. We presented in this work only a few examples which illustrate the technique and its power. Additional examples are given in 1181. The technique can be used to ascertain the structure of many other reaction networks and it is expected to become the standard tool for a rational analysis of the structure of complex reaction networks, whose steady-state equations contain a large number of parameters. Acknowlcdgcments~We are indebted to professors B. L. Keyfitz and M. Golubitsky for helpful discussions and comments. This work was supported by a grant from the National Science Foundation.
heat transfer area dimensionless parameters defined by (27) and (69), respectively C.4 concentration G heat capacity D parameter defined by (27) Da Damkijhler number E activation energy H dimensionless heat transfer coefficient AH heat of reaction preexponential factor rc$ rate constant P vector of parameters flow rate ; temperature reference temperature, defined by (27) T, u heat transfer coefficient V volume of reactor quadrature Wi ith weight in Gauss-Legendre x a state variable xi ith zero of Legendre polynomial dimensionless temperature defined by (70) : dimensionless temperature defined by (27) .G quantity delined by (73) Greek symbols parameter vector dimensionless temperature rise, defined by (40) dimensionless temperature rise, defined by (86) dimensionless activation energy, defined by (40) quantity defined by (16) quantity defined by (73) dimensionless temperature, defined by (58) dimensionless temperature, defined by (86) bifurcation parameter conversion density dimensionless parameter defined by (86) ratio of Damkijhler numbers defined by (&I) Subscripts 0 inlet conditions i ith reaction or ith element c coolant Superscripts 0 singular
point coordinate
REFERMCES 111 Van den Bosch B. and LUSS D., Chem. Engng Sci. 1977 32 203. (21 Tsotsis T. T. and Schmitz R. A., Chem. Engng Sci. 1979 34 13.5. I31 Chang H. C. and Calo J. M., Chem. Engng Sci. 1979 34 285. [4] Leib T. M. and Luss D., Chem. J&gng Sci. 1981 36 210. [51 Luss D., in Dynamics and Modelfing of Reoctioe Systems (Edited by Stewart W. E., Ray W. H. and Conley C. C.), pp. 131-159. Academic Press, New York (1980). 161 Golubitsky hf. and Guillemin W., Stable Mappings and their Singularities. Springer-Verlag, New York (1973). [71 Lu Y. C.. Singularity Theory and an Introduction to Cotastrophe Thheo~y.Springer-Verlag, New York (I 976).
V. BALAKOTAIAH and D. Luss
1622
[8] Brocker Th. and Lander L., Diferentioble Germs and Cafastrophes. Cambridge University Press (1975). [91 Poston T. and Stewart I. N.. Taylor Expansions and Catastrophes. Pitman, London (1976). [IO] Gibson C. G., Singular Points of Smooth Mappings. Pitman, London (1979). [Ill Poston T. and Stewart I., Catestrophe Theory and its Applications. Pitman, London (1978). [12] Zeeman E. C., Cntosrrophe Theory-Selected Papers, 19721977. Addison-Wesley, Reading, Mass. 1977. [13] Golubitsky M. and Schaeffer D., Commu. Pure Appl. Math. 1919 32 21. [14] Gilmore R., Catastrophe Theory for Scientists and Engineers. Wiley, New York (1981). [ISI Golubitsky M. and Keyfttz B. L., SIAMJ. Math. Anal. 1980 11 216. [16] Balakotaiah V. and Luss D., Chem. Engng Sci. 1982 37 433. [ 171 Frank-Kamenetskii D. A., Diflusion and Heat Transfer in Chemical Kinetics. Plenum Press, New York (1969). [18] Balakotaiah V., Ph.D. Thesis, University of Houston (1982). [19] Zeldovich Ya. B. and Zysin Y. A., I. Tech. Phys. 1941 11 cI\-
APPWIX B Two consecutive, exofhennic reactions with no B present in the feed Examples of (Bz, BI, v) values having a bifurcation diagram which is denoted in Fig. 5 as:
1201 Uppal A., Ray W. H. and Poore A. B., Chem. Engng Sci. 1976 31205. [211 Balakotaiah V. and Luss D., Chem. Engng Commun 1981 13
i olie3-i = 0
When both reactions are exothermic, only the first solution is in the feasible region. Equations (AS) shows that a singular point, next to which a parameter region with five steady-state solutions can be found, exists even when the first reaction is exothermic (BI > 0) while the second reaction is endothermic (Bz
(a) (2,3,0.1)
(b) (6,3,0.1)
(c) (6,5,0.001) (e) (6,6.6,0.002)
(d) (6,6,O.ooZ) (0 (6,4,0.012)
(g) (6,5,0.01). The corresponding limit points satisfy the fourth degree equation (Bl)
111. where
APPENDIX A
(I, = 4v2- u(l+ Y)2
Determiytion of the singular point defined by eqns (59) Carrying out the differentiations we get Da,X[B,-O(l+v)l+
~,=v(l+v)B1-8v’(B1+Bz)+v(l+v)(2B, +Bz+vB,+vB$
~Da,~X’(B,+Bz-0)-0=0
Do,X[B,-(8+1)(1+v)]+vDn,zX2[?.(B,+B~-B)-1]-1
=0
aa = 4v*(B1+ Bz)(B, + BZ + 1) +(u + B,)(B, t Bz))-
Da1X[B1-(8+2)(1tv~ltvDa,~X~~4(B1tB~-~)-4I=0 (Al) Da,X[B,-(B+3)(1tv)l+vDa,‘X*[S(B,+B~-8)-12]=0
~(1 +
v)(Bz
IJB,(B, t (1t v)(Bt
Do,X[B, -(e + 4)(1 t v)] t vDa,‘XZ[16(B, + B2- fI)- 321 = 0.
+vBI(Bz+(Y+BI)(BI+Bz)) a~= u*(B, t B& vB,(Bz+(v+B,)(B,+B,))
Defining
and x, = Dn,X[B - S(l+ v)]
Da) =
x2 = Da,X(l + v) x1 = vDa,*X*(B, + Bz - tJ)
(A21
x., = vDa,*X* x5 = 0
+B&-
v(l+ v)Bt
a~=~B,(B,+(l+v)(B,+B~))~4u*~B,tB3* (B2)
B(BI+B2-B-I)-(I-e)(B,+&-8) (I + ~)(BI +Bz- e)+ (BI+ Bz- 0 - l)(B, - @- ve) eme’
For example, the 0 and Dal values at the four limit points for case (g) are (1.365,0.0953), (4.240,0.0609), (5.505,0.0664), (9.746,0.0227).
eqns (Al) can be written as a linear set of equations
x,-x~+2x,-x*=
1
x,-2~~+4~3-4xq=O
(A31
x, - 3x2 +8x, - 12x4= 0 x, -4x2+
16x,-32x& = 0
APPENDlxc In order to determine the impact of the positive exponential approximation and the assumption of equal activation energies on the case of N parallel reactions we examine first the case of N=2. When the activation energies of both reactions are finite but equal there exist two buttertly singular points at e = 3/(1- 31y)
whose solution is x, = 0, x* = 4, x3=3,
xr=l,
(A41
x3 = 3. Substituting (A4) in (AZ) and solving these equations, we solutions B, = 12(2T d/3)
Bz=-1821243
Dal = (2 + d3)e-’
Y = (2 ? q3)*
tJ = 3.
gettwo
(As)
Thus, for every y 26, the qualitative features (such as the number of steady-states and bifurcation diagrams) of the steady state equation are the same as those obtained by taking y-m.
Stmcture of the steady-state solutions of lumped-parameter chemically reacting systems For y<6 there exist only cusp singularities in the feasible region. We conjecture that for all N>2 the qualitative multiplicity features are not affected by the positive exponential approximation provided y is greater than some critical value. When the activation energies of both reactions are unequal and the positive exponential approximation is used the steady state equation is (C2) where
For every p >O, eqo (C2) has two butterfly singularities in the feasible region at
,,qD~,,~=(3~2+5)~~/(w~*+1))
(I+ 3P2) _&A &*e’B = (5,~’ + 3) 7 ~(24,&* (8‘ + 3) g,(&)+B~(&)=n B +Bz=12(a+1)(114-~‘+4~2-p+1) I P0+3PLhP’+3)
+ 1)) fC3)
1623
Two additional butter@ singularities exist at
@=6
P+t
e=12fw*+P+1) IL& + 0
x,=x2=-1
x,=x,=
12 g1=7 1-p
Bt=- 4P2 P -1 Bz = - B,//L’.
BI=-Btp
1
The first of these is outside the feasible region and the second is feasible only when one reaction is exothermic while the second is endothermic. Thus, the assumption of equal activation energies neither increases nor decreases the order of the most degenerate singular point when the positive exponential approximation is used. We conjecture that the same is true for N > 2. Numerical calculations show that in the case of finite and unequal activation energies the (vl. y2) plane is divided into three regions, In one butterfly singularities exist, in one only cusp singularities exist and in one no singular point exists, i.e. only single valued bifurcation diagrams exist. We conjecture that for any N > 2 there exists a region of y values for which N! singular points of order 2N + 1 exist.
OF
VEMURI BALAKOTAIAH and DAN LUSS* Department of Chemical Engineering, University of Houston, Houston, TX 77004, U.S.A (Received
9 October
1981; accepted
26 May
1982)
Abstract-A new, powerful mathematical technique enables a systematic determmation of the maximal number of steady-state solutions of lumped parameter systems in which several chemical reactions occur simultaneously. The method can predict also the different types of diagrams describing the dependence of a state variable of the reactor on a design or operating variable. The technique is applied to several reaction networks giving new results and insight. For example, it is proven that when N independent, parallel exothermic reactions with equal and high activation energies occur in a CSTR there exist N ! distinct regions of parameters in each of which 2N + 1 steady-state solutions exist.
INTRODUCTION
reactions
in a large number of chemically reacting systems. The multiplicity features of systems in which a single chemical reaction occurs have been studied extensively and are well understood by now[l-51. However, most practical control and start-up problems due to steady-state multiplicity are encountered in systems in which several chemical reactions occur simultaneously. Mathematical difficulties have prevented so far the analysis of the qualitative features of these practical cases. to know the In many applications it is important behavioral features of the reactor, such as the number of possible steady-states and the inIluence of a change in one or more operating variables on these states. Unfortunately, the techniques used so far to answer these questions for the single chemical reaction case cannot be extended directly to multireaction systems. We present here some new, powerful mathematical techniques which enable a systematic determination of the above mentioned features of lumped-parameter systems in which several chemical reactions occur simultaneously. We apply these techniques to several types of chemical reaction networks in order to determine: (a) The maximal number of possible steady-state solutions and the parameters for which they exist; (b) The different types of bifurcation diagrams which describe the dependence of a state variable of the system (say the temperature) on a design or operation variable. The analysis is of lumped-parameter systems such as a continuously stirred tank reactor (CSTR), a non-porous catalytic pellet or wire and a porous catalytic pellet with no intra-particle gradients. To simplify the treatment of multireaction systems we assume that all the reactions are of first-order, so that the steady-state behavior can be described by a single nonlinear algebraic equation. The full power of the new technique is illustrated by its ability to predict all the multiplicity features of a CSTR in which an arbitrary number of independent parallel Steady-state
multiplicity
is encountered
*Author to whom correspondence CBS “a, 37, No. ,,--8
should be addressed.
occur,
a feat which could not be accomplished used previously in the chemical engineering literature to analyze similar problems. Due to its power and ease of application, we expect this technique to become a standard tool for determining the multiplicity features of chemically reacting systems and of other systems described by a single algebraic equation. The analysis is based mainly on the singularity theoryIf%141 with which most chemical engineers are not familiar. The possibility of applying these techniques to the analysis of chemically reacting systems was brought to our attention by [IS]. A formal, rigorous survey of the subject is rather lengthy and abstract. Thus, we present in the next section a simple, heuristic description of many of the mathematical terms and concepts. The reader interested in a more formal approach is referred by any of the techniques
to 1133. MATHEMATICAL BACKGROUND
We consider a nonlinear lem) of the type
equation
F(x, h, p) = 0 x is temperature and p is a theory can vectors we the moment
where
(bifurcation
prob-
(1)
a state variable (such as the steady-state or conversion), A is a bifurcation parameter constant vector of parameters. While the handle cases in which both x and F are do not consider here this general case. For we regard p as fixed and rewrite eqn (1) as F(x, A) = 0.
(2)
We select the origin such that
F(O,O)=0
(3)
and define F to be smooth if it has continuous derivatives of all orders. The graph of x vs ,+ which satisfies equation (1) for a fixed p is defined as a bifurcation 161 I
1612
V. BALAKOTAIAH and D. Luss
diagram. A local bifurcation diagram describes the graph of x vs A in a small neighborhood of some point. A global bifurcation diagram describes this graph for all values of x and A within the domain of interest. Two bifurcation diagrams F(x, A) = 0 and G(y, CL)= 0 are cont(lct equiualent if one can be transformed into the other without changing any of the qualitative features in the neighborhood of the origin by a smooth and invertible change of coordinates (4)
(Y>P)-+ MY. CL).A(P)) and a scaling
factor
T(y, p). Thus,
G(Y, ccl = T(Y, ~r)F(x(y,
P), A(F)).
(5)
When this smooth change of coordinates is used the number and order of solutions of x for each fixed A in the neighborhood of the origin is the same as the number and order of solutions of y for the corresponding p if[13]
is unstable
since the bifurcation
diagram
F,(x,h)ex’-
of (\rl.=%1)
assumes different forms depending on whether l >O or +
=0
(13)
is an unfolding of Fz, but it is not versa1 (since it cannot describe the diagram corresponding to E < 0 in Fig. 1) while FI=x”-~,~2-~L~-A=0
x(0,0) = 0,
h(O)=O;
(6)
~(O,OPO,e (0) > 0. the bifurcation
problems
F(x,h)bx*-A=0
(7)
and G(y,p):l-CL-cosy=O. are contact
equivalent,
(8)
since next to the origin
G(Y, Y) = 2F (sin r/2,@) and F(x, A) = OSG(2 sin-’ x, 2A).
(I4 4 1)
F(O, 0, P? = 0
$
CO,%P9 = 0,
(15)
i= 1,2 ,...,r
(16)
and A 5 $$
(O,O, p‘? $
then in the neighborhood
(0, 0, P’? + 0,
of p”, F(x, A) is contact
(10)
the qualitative features of the bifurcation diagram arc unchanged. On the other hand, the bifurcation diagram of F&h):?-A=0
is versa1 but not universal. Only F, is universal, and the codimensiou of F, is one. The bifurcation diagram of eqn (1) is stable for most values of the parameters vector p. For example the bifurcation diagram of F, is stable for all values of l except zero. A set of p values exists at which the function defined by eqn (1) has the most degenerate singularity. Near these points several qualitatively different stable bifurcation diagrams exist. These singular points are characterized by the vanishing of a finite number of the first derivatives of F with respect to x and A. There exists always a polynomial equation which is contact equivalent to F next to these singular points. We shall use in this work the following result which is a special case of a more general theorem proved in [131. Suppose that the function F given by (1) has a singular point p” at which
(9)
Thus, the local qualitative features of a complex nonlinear equation can be determined by the analysis of the features of a much simpler contact equivalent polynomial equation. The next important concept in singularity theory is that of stability. The bifurcation diagram of F(x, A) is defined to be stable if F(x, A)+ g(x, A, c) is contact equivalent to F(x, A) in the neighborhood’of the origin for all sufficiently small smooth functions g(x, A, s). For example, the bifurcation diagram of F(x, A) = 0 given by eqn (7) is stable in the sense that by perturbing F(x, A) to F&C, A) ” x2- a - A = 0,
(14)
T(0, 0) f 0
and
For example,
(12)
(11)
C&
C-0
c>o
Fig. 1. The bifurcation diagram of eqn (12) (E = 0) and the two corresponding perturbed diagrams (et 0).
Structure of the steady-state solutions of lumped-parameter chemically reacting systems equivalent
1613
to G(x,‘)=~‘+‘i’
whose universal
(
unfolding
- sign if A < 0 +si~“if~>o)
(17)
is
G(x, A, n) =x’+’ - a,_,~~-’ - LY,_-~x~-*--. .
- (Y,XT A. (1%
This theorem is most useful as it proves that if such a singular point exists, the maximum number of solutions of eqn (1) next to the singular point is r + 1 and it also enables finding all the possible local bifurcation diagrams. We present now some results useful for a global analysis of the features of the steady-state equation. We define a solution (x. A, p) of eqn (1) to be a singular point if z
(x, A, P) = 0.
(19)
The set of all singular points in the (x, A, p) space is called the singular set. The projection of the singular set onto the (A, p) space is called the catastrophe sot. A point (x, A) on the bifurcation diagram of eqn (1) is called a bifurcation point if the number of solutions of F(x, A) = 0 near x is not constant as A varies in an arbitrarily small neighborhood. A limit point is the least degenerate singular point, i.e. F = F, = 0 while F, and F, do not vanish. Using the implicit function theorem, it can be shown that a bifurcation point is a singular point. However, the converse is not true in general. Thus, the number of solutions of eqn (1) can change only when the parameters cross the catastrophe set in the (A, p) space. It is proven in [13] that the qualitative features of a bifurcation diagram of (1) change when the parameters p cross one of three surfaces. The first, called the hysteresis variety (If) is the set of all points in the parameter space p satisfying 2 F(x,A,p)=~(x,h,p)=~(x,A,p)=O. (20) Elimination of x and A from these three equations gives a single algebraic equation in p, defining a hypersurface. Similarly, the set of all points in the parameter space p satisfying F(x, A, p) = g
(x. A\,p) =
$
(x. A. p) = 0
(21)
is called the isola variety (I), while the set of points satisfying the four equations
p
F(x,, A, p) = F(xz, A,p) = 0 ~(x,,A,P)=$(x,,A,P)=O;
xlzx2
(22)
is called the double limit variety (DL). When p crosses the hysteresis variety, typically two limit points appear or disappear and the local nature of the bifurcation
Fig. 2. Possibli transitions of a local unstable bifurcation diagram (middle column) at the Hysteresis (A), Isola (B, C) and the Double Limit varieties (D, E) into one of two stable forms (left or right column).
diagram changes as shown in Fig. 2(A). When p crosses the isola variety, typically two limit points appear or disappear, so that either the bifurcation diagram separates locally into two isolated curves (Fig. 2C), or one isolated curve disappears (Fig. 2B). When p crosses the double-limit variety, the number of limit points does not change but the relative position of the limit points changes as illustrated by Figs. 2(D) and 2(E). When for no feasible p value a limit point can exist on the boundaries (extreme values) of the feasible x, A values, the features of the bifurcation diagrams can change only by crossing one of the three varieties. In these cases the three varieties divide the global parameter space into several regions, in each of which a different type of bifurcation diagram exists. Note that when eqn (1) satisfies the condition dF/dA f Ofor all feasible x and y the isola variety does not exist. APPLICATION TO
LLIMFSLB PAMMETER CAEMlCALLYREACTINGSYgTlMs
We apply in this section the mathematical techniques described in the previous section to predict the multiplicity features of a lumped parameter system such as a continuously stirred tank reactor (CSTR) in which several types of chemical reactions occur. To familiarize the reader further with the mathematical tools we start with two simple cases of a single exothermic first and zeroth-order reaction. We then illustrate the full power of the new technique by considering the more complex cases of two consecutive exothermic reactions and of N independent, parallel reactions. (1) A single first-order reaction in a CSTR Consider a cooled CSTR in which a single, exothermic, first-order reaction A-B occurs. The corresponding species and energy balances are q(C,qpc,(T,-
C,) - Vi(T)&
T)+ V(-AH)&IJ&--
= 0
Ua(T-
(23) Tc)= 0
(24)
V. BALAKOTMMI and D. Luss
1614
The corresponding
where f&T) = ff(T,,,)exp
[i
(+-J$)].
(26)
where H = Ualqpc,,
a
(2%
Equations (23) and (24) can be combined to give a single steady-state equation for the dimensionless temperature F(Y, D, B) 5 Y - D(B - Y) exp (Y) = 0
hysteresis
variety
defined by (20) is
=o.
(35)
The isola and the double-limit varieties do not exist. Thus, Q = 0 divides the one dimensional (Yspace into two regions. For all (Y~0, the bifurcation diagrams of x vs A are monotonic while for any (Y>O the bifurcation diagram has an S shape. Thus, next to the singular point ( IID, Do, Bo) eqn (26) has only two types of bifurcation diagrams (Y vs D or Y vs B). The catastrophe set of the steady-state eqn (26) is given by
T, = (To + HTc)/(l + If)
F=E=, aY
which upon elimination
Y =&-(l-+).
(27)
We wish to find the maximum number of solutions and the various types of bifurcation diagrams of Y vs D or Y vs B for eqn (26). We note that the set of equations d2F F=n=ay’=O
B5 + ~(B(B D = B v’(B(B
of Y gives
- B T - 4)) exp 1 2 t/(B(B
4)) 1, (37)
If we consider D as the bifurcation parameter, teresis variety defined by (28) reduces to B = 4.
the hys-
(38)
aF
has a unique
solution Y”=2,
B”=4,
D’=e-‘,
(29)
At this point a’F F=_E=’
aB (31)
It follows that a unique singular point exists in the space (Y, D, B) at which the bifurcation problem (26) is contact equivalent to G(x, A) k x3 - A = 0, whose universal unfolding is G(x,h,a):+CCC-h=O.
(32)
The isola and double-limit varieties do not exist in this case. Since no limit point exists on the boundaries of the feasible Y, D region, the nature of a bifurcation diagram can change only at an intersection with one of the varieties. Thus, in this case, Y is a monotonic function of D for all B <4, while the bifurcation diagram (Y vs D) has an S shape for all B > 4. We conclude that in this case the two local bifurcation diagrams close to the singular point are the only possible diagrams in the global parameter space. If we consider B as the bifurcation parameter, the H variety is defined by D = e?. Thus for ah D > em*, Y is a single valued function of B and for D < em2the bifurcation diagram (Y vs B) has an S shape. Figure 3 shows the surfaces represented by eqns (26) and (32) and the corresponding catastrophe sets.
The singular point of F at (Y”,Do, By or of G at the origin is called a “cusp singularity”. Since G = 0 has either one or three real solutions the same is true for F = 0. The nature of the local bifurcation diagrams of (26) can be obtained by studying the bifurcation diagrams of (32). The catastrophe set of (32) is obtained by a simultaneous solution of
G=~=() ax
giving 4a3 = 271=.
(34)
Fig. 3. The surfaces represented by eqns . (32) and (26) and the corresponding catastrophe Sets.
Structure of the steady-state solutions of lumped-parameter
chemically reacting systems
1615
(2) A single zero&order remlion in a CSTR The steady-state conversion in a cooled CSTR in which a zeroth-order reaction occurs is determined by the equation exp
F(5, Da, Y. B) ” E-Da
(&
) =O
(39)
where 5 = (CArI- C.4)/CAll
DO
Da = Vk(T,)/qC,
(W
y = E/RT,,,
p = (- AH)C~/pcpTm(l
+ H)
and T,,,, H are defined by (27). It should be noted that the zeroth-order rate expression always breaks down at high conversions, and for any set of y and p some critical Da exists above which the steady-state equation has no feasible solution (Oi 5 < 0. In order to determine the maximum number of steadystate solutions and the different possible bifurcation diagrams (t vs Da) we find first the singular points of (39). This is accomplished by noting that the set of equations F=aF=aZF=O a[at has a unique
solution 5” = l/,9.
Fig. 4. Different possible global bifurcation diagrams for a zeroth order reaction. Thus, for any fi > 1, .$ is a monotonic function of Da if y < 4 and a cusp exists if y > 4. For parameters far away from the cusp point, the nature of the bifurcation diagrams depends on the intersection between the surface represented by (39) and the boundary of the feasible region (t = 1). For example, for any p > 1 the multiplicity pattern for y slightly larger than four is 1-3-1-O. However, as y is increased the following multiplicity patterns are obtained: l-3-2-0, 1-2-O. The corresponding bifurcation diagrams are shown in Fig. 4 as cases b, c and d. For /3 < I, the cusp singularity defined by (42) is outside the feasible region. However, the set of equations
(41)
for any p, i.e. Da’ = e-‘JP
yD = 4,
(42)
has a feasible
solution
+r-2-t/trtr-4) 33
and that
’
Do’=Pexp(&)
(48) $
(5”, Da’, y”, P) g
(5”. Duo, y”, P) (0.
(43) if 5” < 1. This occurs
Thus, in the (6 Da, y. p) space, the surface represented by (39) consists of a one- parameter family of cusp singularities. However, p is in the feasible region only if p > 1. Therefore, for any p > 1, a unique singular point exists in the neighborhood of which the bifurcation problem (39) is contact equivalent to (32). It follows that eqn (39) can have at most three steady-state solutions in the feasible region provided p > 1. The local bifurcation diagrams for this case are the same as that for eqn (32). We find now the global bifurcation diagrams. The catastrophe set is given by
if -r ’ (I+ P)‘/P.
(49)
It can be verified easily that
t$ th”,Da’,
Y.
8) g
(5”. ho, y. PI > 0.
(50)
Thus, for every set of g < 1 and y > (1 + p)‘/p a unique singular point defined by (48) exists. At that point the bifurcation problem (39) is contact equivalent to G(x, A) : x2 + A = 0.
(51)
Do=.$*exp(a) of G at the origin is known as the “fold” theory[l2] and as the limit point in Therefore, when p < 1, the maximum number of steady-states is two. Next to the singular point eqn (39) has two or no solutions. Figure 4(d) shows the corresponding bifurcation diagram. It is seen that for all The singularity
where
in the catastrophe bifurcation theory.
5_=~-2rd(~(~-4)) 2P The corresponding
hysteresis
variety
y = 4.
. is
(52)
V. BALAKOTAIAH and D. Luss
1616
a unique solution exists. For all Da in (z, Da@) two solutions exist, while for Da > Da’ no feasible solution exists. When p < 1 and (49) is not satisfied a unique feasible solution exists for 0 < Da < 5 and no feasible solution exists for Da > z (Fig. 4a). When this analysis is repeated taking y as the bifurcation parameter, bifurcation diagrams similar to those shown in Fig. 4 are obtained. (3) Two consecutive exothermic reactions Consider a CSTR in which two consecutive
reactions
. ^ AfiB::C
We rewrite
eqn (57) in the standard
F(B,Da,,u,B,,B3=Da,X(B,-e-be)
t vDuI*X’(B, + Bz- 0) - 6 = 0. taking Da, as the bifurcation parameter. Appendix A that the set of equations
has a unique
The species
and energy
Da,‘=
q(CAo-
CA) = V&(T)Ca
q(C*-
Ce) = vUz3
qpc,(T-
balances
(2 tv'3)e-',
are (53) (54)
+(-AH,)V&T)CB
(*)I.
(56)
With these assumptions eqns (53)-(55) can be combined to give a single equation for the dimensionless steadystate temperature
+ (1
B,” = 6d3(2
vBzDa12X2 t vDa,X)
t Dn,X)(l
(57)
which is in the physically point
feasible
=
Vlf~(Tm) 4
e=&(T-
’
Bz = B I(- A&)/(v = &Tm)I~dTm)~
(61)
region. At this singular
(62)
It follows that a unique singular point exists in the feasible region at which the bifurcation problem (59) is contact equivalent to G(x, A) A x’- A = 0, whose universal unfolding is C(r, A, a) : x5 - ax?-
(1*x2- a,x -A
RTm
= 0.
(63)
The singularity of G at the origin is known as the “butterfly”. The surface represented by (63) in the (x, A, a) space is called the butterfly catastrophe and its properties are reported in detail in the mathematical literature [S, 121. Equation (63) has at most five solutions. This implies that the steady-state equation, which is contact equivalent to G, has five solutions for some p. Now, instead of finding directly all the local bifurcation diagrams of the steady-state equation (59) next to its singular point (defined by (61)), we tackle the simpler equivalent problem of finding all the bifurcation diagrams of G. This is accomplished by finding the various regions bounded by the hysteresis and double-limit varieties in the ((I,, (Y?,eyj) space: The hysteresis variety is obtained by eliminating x and A from (63) and ~=5*‘-3a,*‘-2a2x-a,=0
where Da,
- t/3)
$&=-24(2+3)e'
A detailed analysis of the corresponding multiplicity and multiplicity criteria was presented patterns elsewhere[l6]. The mathematical techniques presented in this work can be used to determine directly the maximum number of the steady-state solutions, and all the possible bifurcation diagrams. For illustration purposes we consider only the case of two exotbermic reactions. To simplify somewhat the algebraic manipulations we assume that the activation energies of the two reactions are identical and that no B is present in the feed. Moreover, we assume that the temperature dependence of each rate constant can be described adequately by the positive exponential approximation [ 171, i.e.
B,Da,X e = I+DalX
Y0 = (2 - d3)2
(9.5)
- Ua(T - T,).
/&(T) = f&T,,,) exp.[&
in
#=3
- V&(T)CA
T,)=(-AH,)Vk^,U-)C.
It is shown
(59)
solution
B,“ = 12(2 - g3), occur.
form
(64)
Tm) 2
T.v.
The elimination of x and A from eqns (63)-(65) is equivalent to eliminating x from (64) and (65) and the hysteresis variety can be expressed in the parametric form
AHI) X = exp (0).
= 20x3 - 6as.x - 2~~2= 0.
(58)
IY, = 3cQx*- 15x’ (12= 10*3-3u3x.
(-m
(66)
Structure of the steady-state solutions of lumped-parameter chemically reacting systems
aI
1617
aI
I
-H ------CL
VARIETY WMETY
I
x c
b
b
0
CY,CO
LY3=-0
Fig. 5. Classification of the x vs A bifurcation diagrams of the butterfly singularity defined by eqn (63).
Similarly, the double-limit can be expressed as
variety,
defined
by eqns (22),
t(4t2+7t +4) (x1 = (3tZ+ 4t + 3)2 a32P &4(1+
o’u +3r+ t3’
a33
(3P + 4t + 3)3 with - m < f i m being a parameter. These two varieties are shown schematically in Fig. 5. The hysteresis variety divides the a space into three regions. In the first region no limit points exist and x is a single valued function of A shown as case n in Fig. 5. In the second region two limit points exist and the bifurcation diagram has the S shape shown as case b in Fig. S. For all a in the third region, four limit points exist. Here, five types of bifurcation diagrams may exist, shown as cases c-g in Fig. 5. The double-limit variety separates between the five regions corresponding to each of these diagrams and in applications it is important to know these transitions. One can construct all the possible bifurcation diagrams in the following way[13]. Let (xi, Ai), i = 1, 2, 3, 4 be the four limit points corresponding to an c~ value in region 3. Without loss of generality assume that x, < xz < xa i x+ Since dx/dh is positive for h + + m the Ai must satisfy the restrictions
(9
Al>&
(ii)
A,>h+
(iii)
A,>A,
Thus, the A,‘S can be arranged only in one of the following five possible sequences (shown in Fig. 5). (c) (d) (4 0-I (Ed
Bl Fig. 6. Global classification of the 0 vs DOI bifurcation diagrams for two consecutive exothermic reactions with no B in feed. Close to its unique singular point (eqn (61)), the steadystate equation F = 0 has seven bifurcation diagrams similar to those of G. A detailed procedure for constructing the hysteresis and double-limit varieties in the global parameters space is reported elsewhere[lS]. Figure 6 summarizes the results of the analysis showing that four different cross sections of these varieties with the v-B, plane exist depending on the values of Bz. The hysteresis variety intersects the B, axis always at B, = 4. A monotonic 8 vs Da, bifurcation diagram, similar to case a in Fig. 5, exists in the dotted region to the left of the hysteresis variety. When B+ 14 this region is outside the positive quadrant. An S type bifurcation diagram similar to case b in Fig. 5 exists in the hatched region to the right of the hysteresis variety. The pocket shaped region. which exists for all B, > B,” is divided by the double limit variety into five regions. The corresponding bifurcation diagrams have four limit points and are similar to cases c-g in Fig. 5. For all B, > 27/8 this region is stretched out to unbounded B, values. We report in Appendix B one set of parameter values corresponding to each of the bifurcation diagrams as well as the corresponding limit points.
Consider a lumped-parameter system such as a nonporous catalytic wire or pellet or a cooled CSTR in which N independent, parallel reactions AiL
Pi,
i=l,2
,...,
N
occur. We shall determine the maximum number of possible steady-state solutions, the number of singular points or the number of s’eparate parameter regions in which this number of solutions exists, and the number of all the possible local and global bifurcation diagrams. To simplify the algebraic manipulations we assume that the activation energies of all the independent reactions are
1618
V. BALAKOTN.~~ and D. Luss
identical and sufficiently large so that the positive exponential approximation is valid. The impact of these assumptions on the multiplicity features is examined in Appendix C. With these assumptions the dimensionless steady-state temperature is given by[18]
paring
eqns (74) and (76) we note that zi = xi
17i= wi,
(77)
(68)
is a solution of (74). Since the weights w, and the zeros xI are uniquely determined for any fixed N, it follows that (77) is the unique solution of (74). Substituting (77) in (75) we get
(69)
Thus, the solution
where Da_
=
vm”J
B_
=
jS_-
AfWL,
’ RT,,, ~c,Tm
’
4
and the other parameters
are as defined before. y = e-e
eqn (68) can be written
of eqns (71) and (74) is
Setting (70)
BO = WJXi(1 -xi)
as
y’=exp[-22
PJ Da.B. F(Y, Dai, Bi) k In Y + z, (y
= O.
(71)
We consider one of the Damkijhler numbers as a bifurcation parameter. (Note that the change of variables defined by eqn (70) transforms an S diagram of 0 vs Da; into an inverse S diagram of y vs Da,. However, the coordinates of the singular point are unaffected by the substitution.) Setting the first 2N derivatives of F with respect to y to zero gives
(79) l/j]
The Gauss-Legendre quadrature integrates exactly any polynomial of degree (2N - I), but not one of degree 2N. Thus
a ‘,+:F ay~~
Moreover,
I
= (2N)!
j-zm [ I- (79 + 0 ,z, wiZN] i 0.
at the singular
point defined by eqns (79).
&>O.
N
,z, +“;f$_+L Yrn-; =0,
(80)
(81)
(y
m=1,2
,...,
2N.
(72)
Defining
Thus, at the singular point, the steady-state is contact equivalent to xzN+‘* A whose folding is G(x, A, a) k xZN+’ ~ oa,._,xzN
zi = __ Y
y+Dai
eqns (72) and (7 1) can be expressed
(73)
as
m=l,2,...;2N
(74) (75)
According
to the Gauss-Legendre
quadrature
-.
. .-a,x+A
equation universal
(68) un-
-’ - (I~+,-~x~~~*
=o.
(82)
(The sign of h in (82) is opposite to that of the derivative in (80)). We conclude that there exists always a set of B, and Do, close to the singular point for which 2N+ I solutions exist. Setting dfldy to zero we obtain a polynomial of degree 2N in y. Thus, the maximum number of solutions in the global space is also 2N + 1. We note that if (y. Bi, hi)
= [Y'.
$?q
+
Y’]
i=1,2,...,N (76) where xi are the zeros of the Nth degree Legendre polynomial and wi are the corresponding weights. Com-
tThe equality of the rhs of (78) was proven by _I.Villadsen and S. Wedel.
is a solution of eqns (74) and (75) then because of the symmetry of the problem every permutation of the (St, Dai) i = 1,2,. . , N is also a solution, and there exist exactly N! distinct singular points. We conclude that there exist N! distinct parameter regions in each of which the sfeady-state equation (68) has (2N + 1) solutions. To illustrate this result consider a system with three independent parallel reactions. Here, the following six singular points exist:
Structure of the steady-state solutions of lumped-parameter chemically reacting systems Da,’ = y”
Da,0 = (4k $q15))y0
Da,’ = (4 5 q(lS))y’
Da,“=(4-c~/(15))y0
Do,” = y”
Daz” = (4 T q/(lS))y’
Du,~=(~s~/(~S))~~
Da>’ = (4 ri:d(l5))yO
Da3’ = y”
1619
(83) B,’ = 1619
B,’ = 25/9
Bz” = 2519
BzO = 16/9
Bz” = 2519
B,‘=
Bs’ = 25/9
B,O = 16/9
2519
where y” = exp (-11/3). The nature of the local bifurcation diagrams of the steady-state equation can be determined by finding the bifurcation diagrams of the polynomial equation (82). Unfortunately, because of the large number of parameters it is not possible to obtain a clear picture of the number of regions formed by the hysteresis and double-limit varieties in the (2N - 1) dimensional a space. The hysteresis variety separates regions of the parameter space which differ by two limit points. Next to every singular point (defined by (79)) there exists a region in which 2N + 1 solutions exist, i.e. one which has 2N limit points. Thus, the hysteresis variety divides the a space into no more than N + I regions. The bifurcation diagrams in region j(j = 1,2,. . , N + 1) have 20’ - 1) limit points. For j 2 3, each region is divided into several subregions by the double-limit variety. One can determine all the possible local bifurcation diagrams using the constructive scheme described in the previous example. For example, if no limit points exist 0 is a single valued function of Da,. When two limit points exist an S pattern is obtained. When four limit points exist the five local diagrams shown as cases c-g in Fig. 5 are possible. Similarly, it can be shown that for six limit points 61 different local diagrams exist, etc. Thus, for two parallel reactions, the total number of possible local bifurcation diagrams is 1 + 1 + 5 = 7 (Fig. 5). while for three parallel reactions, there exist 1 + 1+ 5 + 61 = 68 different possible local bifurcation diagrams. The above example illustrates the complexity of mapping all the bifurcation diagrams for a system in which many reactions occur simultaneously. The global bifurcation diagrams of f3 vs Dai can be found by constructing the hysteresis and double limit varieties in the (2N - 1) dimensional space of Bi (j = 1,2 ,.__, N) and Da! (j= 1,2 ,__. ,N, j#i). Different cross sections of the bifurcation set in the Dal-Da2 plane for the case of two parallel reactions were shown in [ 161. These figures show that a 1-3-5-3-5-3-1 multiplicity pattern can be obtained when one of the Damkahler numbers is varied while the other is kept constant. This indicates that the local and global bifurcation diagrams need not be the same when one of the Damkiihler numbers is taken as the bifurcation parameter. A common way of changing the Damkahler number is by adjusting the residence time in the reactor. Consider an adiabatic reactor in which the flow rate is changed in a continuous fashion so that Dai = viDa,,
i=1,2
,.__,
N.
(84)
B la = 25/9
The steady-state
equation
(68) can be rewritten
as
The local bifurcation diagrams of B vs Da, of eqn (85) are the same as those of (82) and there exist again N! singular points. Surprisingly, the local and global biiurcation diagrams are identical in this case. The reason is that eqn (85) can be written as an Nth degree polynomial in Da,e’ whose coefficients are linear functions of 8. Differentiation of the polynomial with respect to 0 gives another polynomial of Nth degree whose coefficients are also linear functions of 0. Elimination of Dales from these two equations gives a 2Nth degree polynomial in 8. The zeros of this polynomial are the limit points in the global parameter space. It follows that any bifurcation diagram in the global parameter space can have at most 2N limit points. However, all the bifurcation diagrams containing up to 2N limit points exists next to each of the N! singular points. This point is illustrated in Fig. 7 which shows the four different types of cross sections the hysteresis and double limit varieties have with the vrBl plane. This figure is symmetric about the line YZ= 1 in the B,-log v2 plane. Bifurcation diagrams with zero and two bifurcation uoints exist in the dotted and hatched regions. resneckvely. Five different types of bifurcation diagrams with
Bl Fig. 7. Global classification of the 0 vs Da, bifurcation diagrams for eqn (85) with N = 2.
1620
V. BALAK~TAIAH and D. Luss
four points exist in the two pocket shaped regions. Thus, there exist seven different bifurcation diagrams in the global space, which are identical to those existing next to each of the two singular points.
DISCUS!3lON
This work illustrates the power of the singularity theory to predict the number of local steady state solutions. The theory is a local one and additional analysis is needed to determine the maximal number of solutions in the global space. It is proven in [ 181 thatfor all the examples presented here the maximal number of solutions next to the highest order singularity is equal to that in the entire global space. It is not yet clear what properties of the steady-state equations are responsible for this result. In all the examples presented here we analyzed the influence of a single dimensionless group, usually the Damkiihler number, on the steady-state temperature or conversion. The value of this dimensionless group can be varied by changing the value of the rate constant at the reference temperature. In practice, one may want to examine the influence of other parameters, such as the reactants concentration, feed or coolant temperature and residence time (flow rate) on the steady-state temperature. To determine the bifurcation diagrams for any of these cases, the dimensionless steady-state equation should be rewritten so that changes in the specific variables affect only one of the dimensionless numbers. This reparameterization does not change the number of solutions but may increase the number of dimensionless groups appearing in the steady-state equation, and is essential for an efficient determination of the singular point and the corresponding contact equivalent bifurcation problem. For example, in order to determine the bifurcation diagrams of the steady-state temperature of a CSTR in which a single first-order reaction occurs vs the coolant or feed temperature, we should introduce the following parameters
versa1 unfolding is given by (32). Thus, the local bifurcation diagrams in this case are identical to those obtained by changing the parameter D in eqn (26). In all the examples considered so far the derivative of F with respect to the bifurcation parameter did not vanish. Consequently, the hysteresis and double-limit varieties describe all the transitions between the various bifurcation diagrams, which are either single-valued or have one or more hysteresis loops of an S shape. Other types of bifurcation diagrams are possible when the derivatives of F with respect to the bifurcation parameter vanishes at the singular point so that an isola variety exists. This has been illustrated in [15, N-211 by analyzing the influence of changes in the residence time on a nonadiabatic CSTR in which a first-order reaction occurred. A steady-state equation such as (1) describes the dependence of a state variable x on a bifurcation parameter A and a finite set of n other parameters p. Thus, it describes a surface in an n +2 dimensional space. While this surface is continuous and smooth its projection into the (~\.p) space may have certain singularities. The examples discussed in this paper show that the number of local solufions to (1) as well as the di$erent possible bifurcation diagrams can be determined if the highest order singular point on the surface is known. All the examples discussed above show that the highest possible order singular point is characterized by the vanishing of a maximum of (n + 1) derivatives. Consider a steady-state equation F, whose derivative with respect to the bifurcation parameter does not vanish, and for which the highest order singularity is characterized by the vanishing of n + 1 - i (i 2 0) derivatives so that
+E=a*F_ 2 ax
and at the feasible
solution a”“-%‘ ax”
o = RTJE,
6 = RTIE,
s = &(-AHKb E pc,(l +W
flu = VkO/q,
(86)
k(T) = k. exp (- EIRT) and use instead tion
of (26) the following
steady-state
equa-
F(&, (T,fie, p, ” B - v - &(jZ + o - ,_)e_“i = 0
(87)
where (r is the bifurcation parameter. This equation for every b > 0 a unique singular point at
has
i”=
B+t/@(l+B)) 2llf
a + X4&1 + PM
2 = wm1+ Da”=(8’-a4exp(ll~9/(B+ and is contact
equivalent
s,, - 81/2 o”-
(88)
- ---
i-7-i
=
a““_‘F ax”+,_i
uni-
0
(89)
of (89) aF x + 0.
In order to determine the value of i we attempt to solve the set of eqns (89) starting with i = 0. If a physically feasible solution exists then the singular point is of order n + 2 and F = 0 has at most n + 2 local solutions. All the examples in the previous section with the exception of the zeroth order reaction belong to this case. If no feasible solution exists for i = 0 one attempts to find a solution for i = 1 and so on. One should note that when i t 1 the solution of (89) consists of n + 2 - i varidetermined as function of ables Cc, A, pi, pZ, , . . , p,-,) the remaining i parameters @.+,_i,. _ , p,). If, for some physically feasible values of these i parameters, the solution of (89) is feasible, then the maximal number of local steady-state solutions in this case is (n + 2 - i). The example of a zeroth-order reaction in a CSTR, discussed in the previous section, illustrates cases for which n = 2, i = 1 and n = 2, i = 2. In fact, by defining
t?‘)
to G(x, A) = x3 - h whose
=
DaP =
d
rSI=g 1+ PS
Structure of the steady-state solutions of lumped-parameter eqn (43) can be written
chemically reacting systems
as
y = L5(y - P) exp (P)
1621
NOTATION
(92)
which is identical in form to (26). However we have chosen the form of (39) as it enables a more direct apphcation of the feasible boundaries. In certain situations it may not be possible to identify the singular point because of algebraic difficulties. One can still find all the bifurcation diagrams by constructing the hysteresis, double-limit and isola varieties in the parameter space determining one diagram in each region bounded by these surfaces. Finally, we wish to point out that the analysis in the previous section can be extended to systems for which the steady-state is described by a vector of m state variables which arc defined by a set of m equations. This situation is likely to be encountered in the description of the steady-state of a lumped-parameter system in which several reactions with complex rate expressions occur. CONCLUslONS This work illustrates that singularity theory is the most efficient and systematic tool for ascertaining the maximal number of steady-state solutions of lumped parameter systems in which several chemical reactions occur simultaneously, and for a direct determination of the regions in the parameter space corresponding to these solutions. The combined use of the singularity and bifurcation theory enables us to determine some important local as well as global properties of the system, specifically to find out all the possible types of bifurcation diagrams and the boundaries of the region in the parameter space corresponding to each fo the diagrams. The construction of the hysteresis, double-limit and isola varieties is the most efficient direct method of determining the parameter regions corresponding to each type of bifurcation diagram. Knowledge of the bifurcation diagram enables the prediction of the influence of a change in a design or operating variable on the structure and number of steady-state solutions. In practice, one needs to know in addition to the number of solutions the type of corresponding bifurcation diagram. For example, consider a case of two consecutive reactions in a cobled CSTR for which three steady-state solutions exist for some set of parameters. This situation must correspond to one of the six diagrams shown as cases b-g in Fig. 5. Obviously, it is important to know to which of these six diagrams this case corresponds. We presented in this work only a few examples which illustrate the technique and its power. Additional examples are given in 1181. The technique can be used to ascertain the structure of many other reaction networks and it is expected to become the standard tool for a rational analysis of the structure of complex reaction networks, whose steady-state equations contain a large number of parameters. Acknowlcdgcments~We are indebted to professors B. L. Keyfitz and M. Golubitsky for helpful discussions and comments. This work was supported by a grant from the National Science Foundation.
heat transfer area dimensionless parameters defined by (27) and (69), respectively C.4 concentration G heat capacity D parameter defined by (27) Da Damkijhler number E activation energy H dimensionless heat transfer coefficient AH heat of reaction preexponential factor rc$ rate constant P vector of parameters flow rate ; temperature reference temperature, defined by (27) T, u heat transfer coefficient V volume of reactor quadrature Wi ith weight in Gauss-Legendre x a state variable xi ith zero of Legendre polynomial dimensionless temperature defined by (70) : dimensionless temperature defined by (27) .G quantity delined by (73) Greek symbols parameter vector dimensionless temperature rise, defined by (40) dimensionless temperature rise, defined by (86) dimensionless activation energy, defined by (40) quantity defined by (16) quantity defined by (73) dimensionless temperature, defined by (58) dimensionless temperature, defined by (86) bifurcation parameter conversion density dimensionless parameter defined by (86) ratio of Damkijhler numbers defined by (&I) Subscripts 0 inlet conditions i ith reaction or ith element c coolant Superscripts 0 singular
point coordinate
REFERMCES 111 Van den Bosch B. and LUSS D., Chem. Engng Sci. 1977 32 203. (21 Tsotsis T. T. and Schmitz R. A., Chem. Engng Sci. 1979 34 13.5. I31 Chang H. C. and Calo J. M., Chem. Engng Sci. 1979 34 285. [4] Leib T. M. and Luss D., Chem. J&gng Sci. 1981 36 210. [51 Luss D., in Dynamics and Modelfing of Reoctioe Systems (Edited by Stewart W. E., Ray W. H. and Conley C. C.), pp. 131-159. Academic Press, New York (1980). 161 Golubitsky hf. and Guillemin W., Stable Mappings and their Singularities. Springer-Verlag, New York (1973). [71 Lu Y. C.. Singularity Theory and an Introduction to Cotastrophe Thheo~y.Springer-Verlag, New York (I 976).
V. BALAKOTAIAH and D. Luss
1622
[8] Brocker Th. and Lander L., Diferentioble Germs and Cafastrophes. Cambridge University Press (1975). [91 Poston T. and Stewart I. N.. Taylor Expansions and Catastrophes. Pitman, London (1976). [IO] Gibson C. G., Singular Points of Smooth Mappings. Pitman, London (1979). [Ill Poston T. and Stewart I., Catestrophe Theory and its Applications. Pitman, London (1978). [12] Zeeman E. C., Cntosrrophe Theory-Selected Papers, 19721977. Addison-Wesley, Reading, Mass. 1977. [13] Golubitsky M. and Schaeffer D., Commu. Pure Appl. Math. 1919 32 21. [14] Gilmore R., Catastrophe Theory for Scientists and Engineers. Wiley, New York (1981). [ISI Golubitsky M. and Keyfttz B. L., SIAMJ. Math. Anal. 1980 11 216. [16] Balakotaiah V. and Luss D., Chem. Engng Sci. 1982 37 433. [ 171 Frank-Kamenetskii D. A., Diflusion and Heat Transfer in Chemical Kinetics. Plenum Press, New York (1969). [18] Balakotaiah V., Ph.D. Thesis, University of Houston (1982). [19] Zeldovich Ya. B. and Zysin Y. A., I. Tech. Phys. 1941 11 cI\-
APPWIX B Two consecutive, exofhennic reactions with no B present in the feed Examples of (Bz, BI, v) values having a bifurcation diagram which is denoted in Fig. 5 as:
1201 Uppal A., Ray W. H. and Poore A. B., Chem. Engng Sci. 1976 31205. [211 Balakotaiah V. and Luss D., Chem. Engng Commun 1981 13
i olie3-i = 0
When both reactions are exothermic, only the first solution is in the feasible region. Equations (AS) shows that a singular point, next to which a parameter region with five steady-state solutions can be found, exists even when the first reaction is exothermic (BI > 0) while the second reaction is endothermic (Bz
(a) (2,3,0.1)
(b) (6,3,0.1)
(c) (6,5,0.001) (e) (6,6.6,0.002)
(d) (6,6,O.ooZ) (0 (6,4,0.012)
(g) (6,5,0.01). The corresponding limit points satisfy the fourth degree equation (Bl)
111. where
APPENDIX A
(I, = 4v2- u(l+ Y)2
Determiytion of the singular point defined by eqns (59) Carrying out the differentiations we get Da,X[B,-O(l+v)l+
~,=v(l+v)B1-8v’(B1+Bz)+v(l+v)(2B, +Bz+vB,+vB$
~Da,~X’(B,+Bz-0)-0=0
Do,X[B,-(8+1)(1+v)]+vDn,zX2[?.(B,+B~-B)-1]-1
=0
aa = 4v*(B1+ Bz)(B, + BZ + 1) +(u + B,)(B, t Bz))-
Da1X[B1-(8+2)(1tv~ltvDa,~X~~4(B1tB~-~)-4I=0 (Al) Da,X[B,-(B+3)(1tv)l+vDa,‘X*[S(B,+B~-8)-12]=0
~(1 +
v)(Bz
IJB,(B, t (1t v)(Bt
Do,X[B, -(e + 4)(1 t v)] t vDa,‘XZ[16(B, + B2- fI)- 321 = 0.
+vBI(Bz+(Y+BI)(BI+Bz)) a~= u*(B, t B& vB,(Bz+(v+B,)(B,+B,))
Defining
and x, = Dn,X[B - S(l+ v)]
Da) =
x2 = Da,X(l + v) x1 = vDa,*X*(B, + Bz - tJ)
(A21
x., = vDa,*X* x5 = 0
+B&-
v(l+ v)Bt
a~=~B,(B,+(l+v)(B,+B~))~4u*~B,tB3* (B2)
B(BI+B2-B-I)-(I-e)(B,+&-8) (I + ~)(BI +Bz- e)+ (BI+ Bz- 0 - l)(B, - @- ve) eme’
For example, the 0 and Dal values at the four limit points for case (g) are (1.365,0.0953), (4.240,0.0609), (5.505,0.0664), (9.746,0.0227).
eqns (Al) can be written as a linear set of equations
x,-x~+2x,-x*=
1
x,-2~~+4~3-4xq=O
(A31
x, - 3x2 +8x, - 12x4= 0 x, -4x2+
16x,-32x& = 0
APPENDlxc In order to determine the impact of the positive exponential approximation and the assumption of equal activation energies on the case of N parallel reactions we examine first the case of N=2. When the activation energies of both reactions are finite but equal there exist two buttertly singular points at e = 3/(1- 31y)
whose solution is x, = 0, x* = 4, x3=3,
xr=l,
(A41
x3 = 3. Substituting (A4) in (AZ) and solving these equations, we solutions B, = 12(2T d/3)
Bz=-1821243
Dal = (2 + d3)e-’
Y = (2 ? q3)*
tJ = 3.
gettwo
(As)
Thus, for every y 26, the qualitative features (such as the number of steady-states and bifurcation diagrams) of the steady state equation are the same as those obtained by taking y-m.
Stmcture of the steady-state solutions of lumped-parameter chemically reacting systems For y<6 there exist only cusp singularities in the feasible region. We conjecture that for all N>2 the qualitative multiplicity features are not affected by the positive exponential approximation provided y is greater than some critical value. When the activation energies of both reactions are unequal and the positive exponential approximation is used the steady state equation is (C2) where
For every p >O, eqo (C2) has two butterfly singularities in the feasible region at
,,qD~,,~=(3~2+5)~~/(w~*+1))
(I+ 3P2) _&A &*e’B = (5,~’ + 3) 7 ~(24,&* (8‘ + 3) g,(&)+B~(&)=n B +Bz=12(a+1)(114-~‘+4~2-p+1) I P0+3PLhP’+3)
+ 1)) fC3)
1623
Two additional butter@ singularities exist at
@=6
P+t
e=12fw*+P+1) IL& + 0
x,=x2=-1
x,=x,=
12 g1=7 1-p
Bt=- 4P2 P -1 Bz = - B,//L’.
BI=-Btp
1
The first of these is outside the feasible region and the second is feasible only when one reaction is exothermic while the second is endothermic. Thus, the assumption of equal activation energies neither increases nor decreases the order of the most degenerate singular point when the positive exponential approximation is used. We conjecture that the same is true for N > 2. Numerical calculations show that in the case of finite and unequal activation energies the (vl. y2) plane is divided into three regions, In one butterfly singularities exist, in one only cusp singularities exist and in one no singular point exists, i.e. only single valued bifurcation diagrams exist. We conjecture that for any N > 2 there exists a region of y values for which N! singular points of order 2N + 1 exist.