Estabilizacion de Reactores CSTR
Estabilizacion de Reactores CSTR
Estabilizacion de Reactores CSTR
Abstract
The steady states of the dynamical Endex CSTR problem are analysed in terms
of degenerate bifurcation (singularity) theory. Steady-state analysis is used in two
ways: (1) as a device to assist in understanding and predicting complex dynamical
behaviour, and (2) as a practical design and operational tool, by applying the concept of quasi-static parameter variation. In the bifurcation analysis, consideration
is given to the effects of thermal mismatching, kinetic mismatching, and variations in the thermal exchange, mass flux, and heat loss rate parameters, on the
structure of the parameter space. The qualitative equivalence of the bifurcation
structure of the kinetically matched, perfectly coupled adiabatic Endex scheme
and the single-reaction adiabatic CSTR is demonstrated, with reference to the
role of the reaction enthalpy effects. Numerical analysis shows that either kinetic
mismatching of the reactions or imperfect heat exchange may introduce Hopf bifurcations into the adiabatic system. This is a result that is both philosophically and
practically important because it shows that limit cycles are not restricted to nonadiabatic thermokinetic systems. The coefficients of thermal exchange, mass flux,
and heat loss are found to induce distortions of the surface of saddle-node bifurcations (the limit-point shell) through codimension 2 bifurcations. The steady-state
and oscillatory-state degeneracies are discussed with reference to the design and
operational implications for a working system.
Keywords: Endex; bifurcation theory; singularity theory; thermokinetics; reactive
thermal coupling; limit-point shell
Present address: School of Chemistry, Macquarie University, NSW 2109 Australia. E-mail:
bgray@laurel.ocs.mq.edu.au
Introduction
In a previous paper (Gray & Ball (1999), hereafter referred to as Part I), we presented an
exposition of the Endex CSTR as a dynamical system. It was shown how experimentally
feasible physical constraints could be imposed to reduce the original four-dimensional
system to one with the relative familiarity of planar dynamics. The concept of perfect
thermokinetic matching was applied to this reduced system, which gave an explicit solution that was shown to define bounds for the state variables over a certain parameter
range. For the unconstrained Endex system, however, global dynamical information is
almost impossible to obtain. This is generally true of most complex dynamical systems
dependent on many parameters. The classical technique for studying the time-evolution
of dynamical systems is that of phase-space analysis, in which trajectories, stable and
unstable manifolds, and separatrices in the space of the state variables (the phase space)
for various initial conditions are deduced or plotted numerically. However, phase space
analysis, and its extension, sensitivity analysis, have severe limitations with regard to
the global study of thermokinetic systems. This is because initial conditions are usually
prescribed, either by practical expedience or by virtue of being built into the model; and
also because thermokinetic systems are typically constrained by at least several variable
parameters, so that it is nearly impossible to obtain a qualitative portrait of the system
dynamics over the physical, or even practicable, parameter realm by means of phase
space analysis alone.
In this paper we turn to the steady-state analysis of the Endex system. Bifurcation
theory, and its extension singularity theory, will play a central role in this study. The
theory and techniques of degenerate bifurcation or singularity theory were developed
by Guckenheimer & Holmes (1983) and Golubitsky & Schaeffer (1985), and applied to
thermokinetic systems by Gray & Roberts (1988) and Ball (1999), among others. A
major aim of this work is to show how the steady-state analyses required in the implementation of bifurcation theory provide qualitative insight into the causes of thermal
instability that is not obtainable by other methods.
The scope of bifurcation theory: dynamical information from steady-state
analysis.
The behaviour of solutions to the dynamical system 1 is the general problem of interest.
dv
(1)
= G v,
dt
v (t=0) = v (0) .
This is a concise representation of a coupled system of first-order, autonomous, ordinary
differential equations, where G is a vector of real-valued functions of the components of
v, the vector of state variables, and , the vector of parameters . In the thermokinetic
context of this work the physical quantities represented by the components of the state
vector v are species concentrations and temperatures, while the vector consists of
time-invariant parameters such as the ambient conditions and thermal properties.
A dynamical process that is subject to mass or energy flux through the boundaries
may attain a nonequilibrium steady state. Setting the left-hand side of equation 1 to zero
2
allows solution in terms of a single state variable v, and the steady states of equation 1
may be expressed as a scalar algebraic equation:
(2)
G v, = 0.
Bifurcation theory is concerned with the study of the qualitative properties of solutions
to equation 2. Specifically, the provenance of bifurcation theory is the multiplicity,
stability, and parameter dependence of solutions to the problem (2). Multiplicity refers
to the number of solutions to (2) at any given . To define criteria for stability we
need to refer back to the dynamical system (1), since the notion of stability involves
assessing the effects of arbitrarily small, transient perturbations to the steady state. A
steady-state solution is said to be asymptotically stable if a perturbation decays to zero
as t . Hirsch & Smale (1974) give precise criteria for stability: a solution to (2) is
locally stable if the eigenvalues of the Jacobian matrix of (1) evaluated at the solution
have nonpositive real parts, and is unstable otherwise.
One of the great strengths of the CSTR, both as a mathematical model and as a
chemical conversion device, is the existence in it of non-trivial steady states. Much of
what we know about time-dependent behaviour in thermokinetic systems actually comes
from steady-state bifurcation analysis. Two instances of how steady-state analysis has
enhanced our understanding of dynamical phenomena were described in Ball & Gray
(1995): the critical damping condition, which marked the appearance of spiral trajectories in the phase plane; and the thermal runaway curve which, although based on a
system-specific and strictly dynamical criterion, records the consequence of perturbations to the steady state. Another example comes from combustion theory, where the
dynamical phenomenon of ignition is most successfully explained when the capacity for
multiplicity is built into the model (Yang & Gray (1969)).
The use of steady-state analysis to extract information about dynamical processes is
closely tied to the concept of quasi-static variation of parameters. Bifurcation diagrams,
in which a system variable in the steady state is plotted as a function of a parameter
(and which sometimes include other useful information such as indications of stability
and the amplitude maxima of limit cycles), are often interpreted as portraits of the states
of a system that is subjected either to deliberate slow parameter tuning or affected by
incidental parameter drift.
Figure 1 illustrates one possible start-up procedure for coaxing the simulated example Endex CSTR from Part I, in which the reactions 2,3-epoxy-1-propanol 1,2,3propanetriol and 2,3-dichloro-1-propanol 2,3-epoxy-1-chloropropane are thermally
coupled, to a desired steady state of operation. For reference the dimensionless adiabatic Endex equations are:
dx1
= x1 e1/u1 + f (1 x1 )
(3)
d
du1
(4)
= x1 e1/u1 + f 1 (uf u1 ) lex (u1 u2 )
1
d
dx2
(5)
= x2 e/u2 + f (1 x2 )
d
du2
(6)
= x2 e/u2 + 2 f (uf u2 ) + lex (u1 u2 )
2
d
3
(Note: The symbols and notation are defined at the end of the paper.)
[Figure 1 about here.]
The substance that reacts exothermally (2,3-epoxy-1-propanol) is titrated into the feed
stream to one compartment of the reactor in contiguous doses, while the coupled compartment houses the endothermal reaction. Directly after each dose there is a transient
spike which decays to a steady state, whereupon the next dose is added. The relatively
large transients in the middle range of the dosage steps warn of possible rogue thermal
behaviour, if the doses are too large. Conversely, it can be seen that, were the doses
made smaller, the transients would become less of a problem. In the theoretical limit of
quasi-static addition of the reactant the system is always at a steady state and the figure
becomes a bifurcation diagram, with the reactant feed concentration replacing time on
the horizontal axis. (Note: The doses of exothermal reactant represent variations in the
specific reaction enthalpy, (H1 ) c1,f . This quantity is obtained as an independently
variable parameter very simply; by defining = 1/, multiplying equations 4 and 6
through by , and redefining 1 , 2 , and lex accordingly.)
Curves and surfaces of limit points or other codimension 0 bifurcations also may be
interpreted as plots of quasi-static variations of two or three parameters. The importance
of bifurcation diagrams to steady-state analysis lies in their value as devices for understanding and predicting dynamical behaviour, rather than their curiosity value as objects
to be collected, sorted, and named. This, then, is the rationale for the steady-state analysis of the Endex problem. It is not intended here to construct a catalogue of all of the
bifurcation possibilities. We shall use steady-state bifurcation analysis to reinforce the
concept of quasi-static parameter variation and to provide a qualitative appreciation of
the relationship of the Endex to the single-reaction CSTR. This is still a formidable task,
because the system is, in general, analytically intractable. The bifurcation analyses of
the single CSTR in Ball & Gray (1995) and Ball (1999) and the constrained reductions
of the current system in Part I provide the formal framework within which the largely
numerical steady-state analyses in this paper are carried out.
At this stage we shall return to the constrained planar system that was derived in Part
I, with reaction exponents (m and n) equal to 1:
dx1
= x1 e1/u + f (1 x1 )
d
du
(1 + 2 )
= x1 e1/u x2 e/u + f (1 + 2 ) (uf u)
d
x2 = (1 x1 ) (1 + 2 ) (uf u)
u(0) = uf , x1 (0) = 1, x2 (0) = 1.
(3)
(7)
(8)
(9)
This system in two dynamical state variables may be derived from equations 36 by
taking the limit 1/lex 0, and making use of the fact that it is integrable to obtain an
4
(10)
G = Gu = Guu = 0, Guuu 6= 0,
yield the coordinates of the cusp of the limit point curves:
r
1
1
us =
2 +1
+
1
+
1
f = e2 1+/(1)
+1 1
1
.
uf =
p
2 1 + / (1 ) + 1
(13)
(14)
All of these reduce to the corresponding equations of the single-reaction CSTR when
= 0 (no endothermal effect). It would seem that this version of the Endex CSTR
is qualitatively identical to the single reaction adiabatic CSTR, because no new degeneracies of multiplicity are introduced. However, in two-variable systems there is the
5
possibility of oscillatory instabilities arising from the Hopf bifurcation. To what extent
is there an inbuilt potential for thermal oscillations in the kinetically matched planar
Endex scheme? It is easily shown that this system cannot show Hopf bifurcations. The
Jacobian matrix can be written
e1/us f
J =
xs e1/us /u2s
1/us
( 1) (uf us )
/u2s
1/us
(15)
The Hopf bifurcation condition tr J = 0 gives an expression for uf that can be positive:
uf = u s +
1
+ 2u2s 1 + f e1/us ,
the rate of a reaction with E = 90kJ/mol, all other things being equal. Reactions which
differ in activation energy by only a few percent may not see each other because they
occur on quite different time-scales, and the efficaciousness of endothermal stabilization
then is much diminished or non-existent. Where the activation energies are of comparable magnitude, but not identical, it is reasonable to expect to find bifurcations that may
introduce thermal instabilities. A previous analysis of a thermokinetic system containing
two Arrhenius terms (Gray & Forbes (1994)) showed a great diversity of phase portraits
and bifurcation diagrams, due to both a dissipative term and unmatched reaction rates.
By numerical analysis we can find no evidence of degeneration with respect to multiplicity in the current system. However, Hopf bifurcations can be introduced where 6= 1,
which contrasts strikingly with their absence where 6= 1. Figure 3 shows a numerically
computed curve of saddle-node bifurcations, for = 1.06, with an associated curve of
Hopf bifurcations that terminates at a double zero eigenvalue (DZE) point.
[Figure 3 about here.]
The short branches of limit cycles emanating from these Hopf bifurcations are unstable,
of low amplitude, and terminate at homoclinic orbits.
The numerical calculations indicate that the origin of the Hopf curve may be a DZE
DZE bifurcation at a point on the saddle-node curve. Computationally it is impossible
to follow the fate of the Hopf curve accurately enough find such a point because the curve
is extremely sensitive to parameter variations. Note that the DZEDZE bifurcation is of
codimension 2, whereas the highest order singularity of multiplicity in the steady states
of equations 38 evidently is the hysteresis variety, of codimension 1. If the Hopf curve
does originate at a codimension 2 singularity then it would seem that parameters which
may not introduce degeneracies of multiplicity, because they are qualitatively already
accounted for, may nevertheless introduce degeneracies of the Hopf bifurcation. We
shall have the opportunity to observe further consequences of Hopf bifurcations in the
Endex CSTR in the next sections. It is interesting to note that Hopf bifurcations and
hence periodic solutions which occur here, do so in an adiabatic system.
In this section, the parameter space of the unconstrained adiabatic Endex CSTR, equations 36, is investigated numerically. (Note: Since it is an integrable system, equations
36 may be regarded as having three dynamical variables and one algebraic variable,
hence the phase space is effectively three-dimensional.) We are interested in how the
combination of finite heat exchange between reactor compartments and imperfect reaction matching may affect patterns of multiplicity and Hopf bifurcations.
The efficiency of thermal exchange between endothermal and exothermal reactions is
one of the more pressing questions arising from Part I and the previous section, where a
conservative system in two dynamical variables was obtained by assuming a technological
fait accompli of perfect thermal exchange. However, imperfect heat exchange could
have serious design and operational implications. Even the best heat exchange surfaces
7
(16)
(Note: Here G is a bifurcation function in the state variable v, which may denote either
u1 or u2 , and the subscripts specify nth order differentiation with respect to v or f as
appropriate.) The asymmetric cusp has not been found analytically or numerically but
its effect can be seen in B: the lower hysteresis point of the upper island becomes tilted
in the direction of the lower region of multiplicity. In C the hysteresis point of this lower
region and the lower hysteresis point of the upper island have converged via a third
codimension 2 bifurcation.
In B, the Hopf bifurcations associated with the lower region of multiplicity are due
to kinetic mismatching; we encountered them earlier in figure 3. By contrast, the Hopf
bifurcations to stable limit cycles associated with the upper region of multiplicity are
due solely to a finite rate of heat exchange between the reactions. These separate loci
begin to interact in C and unite in D. At some point there is also a switch in stability
of the unstable limit cycles emanating from the lower Hopf bifurcations, because in D
the entire locus is stable. The nature of the degeneration of the Hopf loci is therefore
unclear, although a DZE-DZE bifurcation is one obvious possibility.
It is important to remember that the reactor described by equations 36 has no
dissipative losses. The assumption that a dissipative term is necessary for oscillatory
behaviour to occur, implicit in much of the CSTR literature, is incorrect. Limit cycle
8
behaviour is not due to a loss term per se, but may occur in an adiabatic system of two or
more dimensions whenever there is competition between different rate processes. As with
the single-reaction dissipative CSTR, limit cycles in the adiabatic Endex CSTR are associated with Hopf bifurcations and with multiplicity via homoclinic terminii. We should
also be aware of the possible occurence in the Endex system of dynamical instabilities
such as spiral trajectories and wrong-way behaviour.
From this and the previous discussions we may conclude that the cast of multiplicity
has an intrinsic thermal properties component, associated with the thermal and energy
properties of the reaction mixtures, and an exchange component, that moves into the
parameter space from infinity. Oscillatory behaviour may be due to kinetic mismatching
or imperfect thermal exchange, or interaction between the two.
Thus far in the steady state analysis we have retained the conservation conditions
f1 = f2 f and ld = 0
(17)
on the Endex CSTR, firstly because this achieves a reduction of the state space dimensionality, and secondly because adiabatic operation allows the maximum recovery of heat,
and thus possibly optimum conversion. However, quite apart from the fact that the conditions (17) represent idealizations, there may be legitimate control reasons for relaxing
one or both of them. For example, it was shown in Ball (1999) that heat dissipation
allows multiplicity to be avoided in the single diabatic CSTR by increasing the residence
time and/or lowering the feed temperature, in a way that is impossible under adiabatic
conditions. From the above discussion of multiplicity in the adiabatic Endex CSTR it
might be fairly predicted that the cast of multiplicity in the truly dissipative system will
be a distortion of that found in the single dissipative CSTR. A sneak preview of figure 7
confirms this impression, but we may wish to know more about particular causes and
effects of dissipation. To what extent is the adiabatic ideal a workable (or unworkable)
approximation? What is the role of the two feed flowrates in dissipative behaviour? How
and where are oscillatory or other thermal instabilities introduced? The analysis of this
section is based on the following dissipative system:
dx1
d
du1
1
d
dx2
d
du2
2
d
= x1 e1/u1 + f1 (1 x1 )
(18)
(19)
= x2 e/u2 + f2 (1 x2 )
(20)
(21)
flowrate or residence time are very scanty. Li (1994) demonstrated that a winged cusp
singularity is introduced into a model for an isothermal autocatalytic reaction in a CSTR
with two inflows of reactants; but there was no intimation that this singularity might
be allowed through the breaking of a conservation condition. Of the two conservation
conditions in (17) it is the first over which the experimenter has the greater degree
of operational control. It is therefore of great practical interest to observe how two
independent flowrates may affect the steady states.
There is also a theoretical reason why we are interested in uncoupling the residence
times. In Ball (1999), evidence was presented that indicated a correspondence between
the number of physically different quantities in a system and the highest codimension of
degeneracy in the mathematical model. For the dissipative single CSTR it was found that
it is not the separate feed and coolant temperatures per se that increased the codimension of system degeneracy, but rather the (implied) existence of finite coupling between
the two temperature parameters. For the adiabatic Endex CSTR it was found that the
thermal and energy properties of the second reaction are already accounted for in the
energy or thermal properties dimension of the parameter space, in the sense that they
do not introduce higher degeneracies of multiplicity. The flowrates (inverse residence
times) f1 and f2 like the feed and coolant temperaturesare physically identical and
independently variable quantities, yet we saw in Part I that they make or break a conservation condition. How, then, do the residence times affect the structure of the parameter
space?
For simplicity, we shall use the assumption of perfect thermal exchange and retain the
second conservation condition in (17), so that the effect of f2 may be observed without
interference from lex and ld . The steady-state bifurcation equation can be written as:
G=
f1 e1/us
f2 e/us
+ (f1 1 + f2 2 ) (uf us ) = 0.
e1/us + f1 e/us + f2
(22)
A suitable strategy in principle might be to search for limit, hysteresis, pitchfork bifurcations, and so on, by writing down the appropriate derivatives, then compare the
codimension of the highest order singularity with that of the fully conservative Endex
system. However, by inspection we can see that the higher order derivatives with respect
to us will be sufficiently complex to thwart such attempts at analytic solution. Instead,
we shall impose the additional condition = 0. (Admittedly, this is a rather absurd
condition, for it means that no reaction or no thermal effect at all is occurring in the
endothermal compartment. This situation is contrary to the philosophy of the Endex
concept; but it has pedagogical value.) Then, with the bifurcation equation
G=
f1 e1/us
+ (f1 + f2 ) (uf us )
e1/us + f1
(23)
(24)
(25)
2. We may expect the winged cusp to occur in the physical region for zero- or finitelycoupled feed temperatures, uf1 and uf2 . For zero-coupled feed temperatures the
bifurcation equation becomes
G=
f1 e1/us
+ f1 (uf1 us ) + f2 (uf2 us ) .
e1/us + f1
The algebra yields the following equations for the coordinates of the winged cusp,
at fixed uf2 :
2uf,2
us =
1 4uf,2
us (1 + us )
uf 1 =
1 + 2us
1 + 2us
1/us
f1 = e
(26)
1 2us
(1 4u2s ) (1 2us )
=
4u2s
(1 + 2us )
f2 = 2us e1/us
.
(1 2us )2
3. The pitchfork singularity, of codimension 2, is the highest order degeneracy in the
system where uf1 and uf2 are infinitely coupled, as in equation 23. The defining
conditions
G = Gu = Guu = Gf = 0; Guuu Gf u 6= 0
(27)
yield the following parametric equations for the pitchfork locus:
us
1 + 2us
1 + 2us
1/us
f1 = e
1 2us
2
(1 4us ) (1 2us )
=
8u2s
2
1 + 2us
1/us
.
f2 = e
1 2us
uf =
(28)
11
5. Thus we may suspect that there is a rather elusive connection between conservation
conditions, state space dimensionality, and the codimension of the highest order
singularity in a bifurcation problem. Although formally the parameters f1 and
f2 describe the same physical property (mass flux), their uncoupling effectively
introduces a physical property that is qualitatively equivalent to heat loss. Figure 5
displays a limit point shell from equation 23 for a fixed value of f2 . It is formally
identical to that of the single-reaction dissipative CSTR in Ball (1999).
[Figure 5 about here.]
Returning now to the bifurcation equation (22) it can be seen that in general we may
expect to find a variety of codimension 1, 2, and 3 singularities in this system. Can we
hope (or hope not) to find higher order singularities? This is an open question because,
as mentioned above, the higher order derivatives with respect to us become very complex.
Since we have just demonstrated the qualitative equivalence of heat exchange, heat loss,
and a second residence time it appears unlikely that higher order bifurcations are hidden
in the system, because there are no more physically independent properties to vary.
These results also indicate that we should expect to find the usual assortment of
Hopf bifurcations, degeneracies, and other oscillatory instabilities in the Endex CSTR
with variable residence times. This expectation is borne out in figure 6, which is a
bifurcation diagram for an Endex system under the conditions indicated by equation 23.
It shows a branch of limit cycles emanating from the upper Hopf bifurcation and ending
at a homoclinic terminus. (There is also a branch of limit cycles from the lower Hopf
bifurcation that ends at a homoclinic terminus. It is so short that it does not show
on the diagram.) The Hopf bifurcations have occurred solely as a result of uncoupling
the residence times. Engineers should note that under these circumstances attempts to
control reactant conversion by manipulating the feed flow rates separately may induce
thermal oscillations.
[Figure 6 about here.]
This concludes the discussion of the enthalpy flux conservation condition.
Heat losses
The effects of relaxing the second conservation condition in (17) should be simpler to
assess, because ld appears linearly in the bifurcation equation. In this section we shall
look at how heat losses affect multiplicity in equations 1821 with f1 = f2 f . The
structure of the parameter space may be deduced from the sequence of limit-point shell
transections in figure 7, which have been constructed at fixed values of lex and ld .
[Figure 7 about here.]
They should be compared to those in figure 4. Again, as in figure 4, in passing from
A to D we are moving backwards along the -axis. The unexpected appearance of two
islands of multiplicity in A is mundanely resolved as we move back to lower values of
: they are nothing more than transections of knobbles on a surface of multiplicity that
12
begins to resemble the limit-point shell of the single diabatic CSTR. Some of the possible
features of this surface, inferred from Figure 7 and from the previous discussion of the
unconstrained adiabatic system, may be enumerated:
1. There may be two pitchfork points. The upper pitchfork is associated with finite
heat exchange. It moves in from infinity as lex is decreased. The lower pitchfork is
associated with finite heat loss, and it moves out from zero as ld is increased.
2. The upper and the lower pitchforks may each degenerate to a winged cusp on variation of a fourth parameterbut there is a physical condition for this. It was shown
in Ball (1999) that the occurence of the winged cusp in the single CSTR requires
a looser parameter space, in the sense that feed and coolant temperatures must
be finitely coupled. In the dissipative Endex CSTR the association of the lower
pitchfork with heat losses suggests that the winged cusp of the lower pitchfork may
occur for zero- or finitely-coupled feed and coolant temperatures. The association
of the upper pitchfork with heat exchange suggests that an upper winged cusp requires two zero- or finitely-coupled feed temperatures. (Note: This does not imply
that the two hypothesized winged cusps must coexist.)
3. The possibility of bifurcations of codimension > 3 must be countenanced, although
the evidence indicates that lex , ld , and f2 are different incarnations of an overall
loss quantity, which forms the fourth dimension of the parameter space, just as
the combined intrinsic energy properties were found to form the third dimension
of the parameter space. If this is true, then we may expect these and any other
qualitatively similar parameters to introduce smooth bumps and hollows on the
limit point surface due to codimension 2 bifurcations, but not necessarily cause
bifurcations of higher order.
Example
Finally, we turn to a specific example. We shall observe the steady states of the example
Endex reactor described in Part I, which was last seen in the time-series/bifurcation
diagram of figure 1 of this paper. Recall also that the phase portrait of figure 7A of
Part I showed that dynamical stabilization of the exothermal reaction could be achieved
under idealized 2-dimensional conditions at specified operating parameter values. In
general we cannot expect to see stabilization under all operating conditions, so we must
aim for a range of steady states that are thermally resilient yet maximise conversion of
both reactants. A parameter plane diagram of saddle-node and Hopf bifurcations for the
example pair of reactions housed in an adiabatic Endex reactor is shown in figure 8A,
while in figure 8B the thermokinetic parameters have been manipulated to eliminate
multiplicity and Hopf bifurcations in the lower range of f .
[Figure 8 about here.]
The bifurcation diagrams in Figure 9 have been reconstructed from figure 8B at a selected value of f by calculating the steady states and the Hopf branches from the Hopf
13
bifurcation point to as close to the homoclinic terminus as possible. The diagrams indicate the possibilities for optimization, with respect to conversion and temperatures.
(For easy reference we have provided dimensional temperature units on the axes.) It
can be seen by comparison of figures 8 and 9 that multiplicity at low flowrates is built
in to the system by virtue of the mismatch in ratesbut it should be noted that it is
low-temperature multiplicity. Clearly there is potential here for using the hysteresis loop
as a form of control or switch. However, it may be a one-way switch as indicated by the
arrows: we may jump the system to a safe high-conversion state by adjusting the feed
temperature upwards, but a branch of limit cycles may be an impediment to the reverse
procedure of quenching the reactions. Paradoxically, it is the branch of limit cycles in
the thermokinetically matched system in figure 9D, E, and F that is more dangerous because the stable attractor can occur above the boiling point, whereas in the unmatched
sytem of A, B, and C the branch is unstable and contained well below the boiling point.
We see that in this case a certain degree of kinetic mismatching actually tempers the
oscillations indicated by figures 8B and 9D, E, and F that are purely exchange in origin.
[Figure 9 about here.]
Conclusion
The Endex reactor is a promising new development in chemical reactor design and engineering technology. In Part I we showed that unequivocal dynamical stabilization
could be achieved in a two-reaction Endex CSTR over a wide range of thermokinetic
parameters. In this paper we applied a combination of degenerate bifurcation theory
and numerical analysis to the steady states of this system and achieved a comprehensive survey of the parameter space. What seemed at first to be a hopelessly large and
intractable system, in terms of the number of parameters and variables, became somewhat simpler when the bifurcation analysis showed that the parameters fall into only
four physically distinct groups, representing ambient temperature, mass flux, intrinsic
thermal and energy properties, and heat losses.
The bifurcation structure of the Endex parameter space was found to be a smooth
perturbation of the single-reaction CSTR parameter space, with distortions of the saddlenode bifurcation surface, or limit-point shell, occurring via additional codimension 2 bifurcations. Mismatched reaction rates and dissipative terms introduce Hopf bifurcations,
which may be a potential source of thermal instability in reactors with poorly matched
reaction pairs, different flowrates through the compartments, inadequate heat exchange
surfaces, or badly chosen heat dissipation rates. We used an example reaction couple,
hydration of 2,3-epoxy-1-propanol and dehydrochlorination of 2,3-dichloro-1-propanol,
to demonstrate the value of steady-state analysis in reactor design and operation. It was
shown how the concept of quasi-static parameter variation could be usefully applied to
assess a dynamical procedure such as start-up or shut-down. Bifurcation diagrams were
used to compare the behaviour of the selected reaction pair with the matched ideal.
This study is restricted to reactors of the CSTR configuration, although preliminary
work indicates no particular difficulties in extending the Endex concept to other configurations such as tubular, batch, or semi-batch reactors. In future work we shall report
14
on the results of these studies, as well as applications of the concept to reactor and
storage-tank situations where there may be large spatial gradients.
= V2 /V1
= CE1 /c1,f (H1 ) R
= E2 /E1
m1
n1
= A2 c2,f
/A1 c1,f
n1
= tA1 c1,f
Other symbols
q = / (1 )
v =a state variable
W = winged cusp bifurcation
G = a bifurcation function
= 1/
= dimensionless scaling factor (0.0338)
Subscripts
1 = of the exothermal reaction, of the reactor compartment housing the exothermal
reaction
2 = of the endothermal reaction, of the reactor compartment housing the endothermal
reaction
a = ambient, of the surroundings
d = dissipative, of heat transfer to the surroundings
ex = exchange, of heat transfer between reactions
f = of the feed
s = of the steady state
Abbreviations
DZE = double zero eigenvalue
trJ = trace of the Jacobian matrix of differential coefficients of the linearization of
coupled ODEs.
det J = determinant of the Jacobian matrix of differential coefficients of the linearization
of coupled ODEs.
16
Endnote
For computing the the data in figures 1, 3, 4, 6, 7, 8, and 9, the equations were scaled
by the factor e1/ , where is a reference temperature selected only for computational
convenience. Thus the exponential term is written as e1/ e1/ e1/u1 = e1/ e1/1/u1 ,
then dividing through by e1/ also scales the dimensionless time, flowrates, and loss
and exchange coefficients. Scaling the equations in this way forestalls possible numerical
problems, and avoids having very small numbers ( 1010 ) on the axes of the diagrams.
Acknowledgments
The authors thank the Australian Research Council for support of this research.
The success of this work owes much to the computer programs that were used or
adapted to produce the numerical and graphical output. For most of the bifurcation
and parameter plane diagrams the Fortran program Auto (Doedel (1986)) was used
in its (1994) version. This program has the wonderful ability to compute branches
of periodic solutions and to compute second-parameter continuations at codimension
zero bifurcations. It cannot, however, compute its own starting solutions. These were
obtained from another excellent piece of public domain scientific software, the C program
Dstool (Guckenheimer et al. (1993)), which also was used to compute trajectories for
phase portraits (in Part I) and time series. The graphical capabilities of these packages
are limited however, so that the diagrams using the the numerical output were rendered
by the program Gnuplot (Williams & Kelley (1993)), with subsequent (and sometimes
extensive) hand editing of the generated output of Postscript code. For some diagrams,
the fast function-plotting capability of Gnuplot was also used. The bifurcation surfaces,
and some two-dimensional plots of Part 1, were obtained using the commercial package
Mathematica (Wolfram (1991)). Again, the realization of these plots required extensive
post-translational modification of the generated Postscript code.
References
Ball, R. (1999). The origins and limits of thermal steady-state multiplicity in the continuous stirred tank reactor. Proc. R. Soc. Lond. A 455, 141161.
Ball, R & Gray, B.F. (1995). Transient thermal behaviour of the hydration of 2,3epoxy-1-propanol in a continuously stirred tank reactor. Ind. Eng. Chem. Res. 34,
37263736.
Doedel, E. J. (1986). AUTO: Software for Continuation and Bifurcation Problems in
Ordinary Differential Equations. Technical report, Californian Institute of Technology.
Golubitsky, M. & Schaeffer, D.G. (1985). Singularities and Groups in Bifurcation Theory,
volume 1. SpringerVerlag, New York.
Gray, B. F. & Roberts, M. J. (1988). A method for the complete qualitative analysis
of two coupled O.D.E.s dependent on three parameters. Proc. R. Soc. Lond. A 416,
351389.
17
Gray, B.F. & Ball, R. (1999). Thermal stabilization of chemical reactors. I. The mathematical description of the Endex reactor. Proc. R. Soc. Lond. A 455, 163182.
Gray, B.F. & Forbes, L.K. (1994). Analysis of chemical kinetic systems over the entire
parameter space IV. The Salnikov oscillator with two temperature-dependent reaction
rates. Proc. R. Soc. Lond. A 443, 621642.
Guckenheimer, J. & Holmes, P. (1983). Non-linear Oscillations, Dynamical Systems and
Bifurcations of Vector Fields. In Appl. Math. Sciences, volume 42, New York. Springer
Verlag.
Guckenheimer, J., Myers, M.R, Wicklin, F.J. & Worfolk, P.A. (1993). Dstoola dynamical systems analysis and exploration toolkit. Center for Applied Mathematics, Cornell
University.
Hirsch, M.W. & Smale, S. (1974). Differential Equations, Dynamical Systems, and Linear
Algebra. Academic Press, New York.
Li, Ru-Sheng (1994). Continuous flow stirred tank reactor with two inflows of reactants:
A versatile tool for study of bifurcation in chemical systems. Chem. Eng. Sci. 49,
202931.
Williams, T. & Kelley, C. (1993). Gnuplotan interactive plotting program. Copyright
(C) 19861993 Thomas Williams, Colin Kelley.
Wolfram, S. (1991). Mathematica: a system for doing mathematics by computer.
Addison-Wesley, Redwood City, California, second edition.
Yang, C. H. & Gray, B. F. (1969). Unified theory of explosions, cool flames, and two-stage
ignitions. Trans. Faraday Soc. 65, 16141622.
18
List of Figures
1
2
3
5
6
7
8
19
20
21
22
23
24
25
26
27
28
A
390
360
T1/K
330
300
time/10 3s
B
390
360
T1/K
330
stable nodes
stable foci
300
0.4
0.8
Figure 1:
20
1.2
1.6
1.5
=0
f
= 0.3
= 0.6
0.5
-0.1
uf
Figure 2:
21
0.1
0.2
10
1
DZE
f
0.1
0.01
0.022
0.026
uf
Figure 3:
22
0.03
0.034
A. = 80
B. = 50
1E+4
100
f
1
0.01
0.03
0.04 uf
0.05
0.03
C. = 40
uf
0.04
D. = 30
1000
f 10
0.1
0.03
uf
0.04
0.03 uf
Figure 4:
23
0.04
f 2= 0.01
0.3
f1
0.1
uf
0.2
Figure 5:
24
0.044
0.04
us
0.036
HB
0.032
0.03
HB
0.031
Figure 6:
25
uf
0.032
100
100
f
1
0.01
0.01
0.03
0.035
0.03
0.035
D
100
100
f
1
0.01
0.01
0.03
0.035
0.025
F
1000
1000
0.035
10
10
0.1
0.1
0.001
0.001
0.02
uf
0.03
0.04
Figure 7:
26
1e-05
0.02
uf
0.03
0.04
10
f
4
0.025
uf
0.03
0.035
0.03
0.035
10
f
4
0.025
uf
Figure 8:
27
400
x1
355
T1
310
0.5
265
0
C
x2
x2
0.5
0.5
0
E
400
--- boiling point
x1
355
T1
0.5
310
265
265
Tf
310
355
Figure 9:
28
265
Tf
310
355