Chapter 5 Global behaviour in the oxidation of hydrogen, carbon monoxide and simple hydrocarbons

Chapter 5

Global Behaviour in the Oxidation of Hydrogen, Carbon Monoxide and Simple Hydrocarbons S.K. SCOTT

5.1 INTRODUCTION

This chapter describes in some detail the information and approaches to data analysis that can be drawn from the relatively new conceptual structures of non-linear mathematics in the context of combustion chemistry. Chemical reactions that evolve with time are part of the class of dynamical systems: other examples arise throughout physical, biological and engineering science, but chemical kinetics is an area almost tailormade for the new mathematics developed in the past 15 years or so. The two simplest oxidation/combustion reactions, H 2 and CO, have been the most intimately studied via the non-linear mathematics approach, but some features of the combustion of hydrocarbons, the so-called cool-flame phenomena, can also be usefully discussed in this context. The main interest here will not be in the details of the individual elementary steps by which a chemical reaction occurs but in the global behaviour of the reaction under particular experimental operating conditions such as the ambient temperature, mixture composition and total pressure. This "behaviour" may be for the reaction to respond to the operating conditions by establishing a steady-state reaction rate in which only a small fraction of the initial fuel and oxidant reacts or, for slightly different conditions, to exhibit an explosive reaction. The changes in the qualitative nature of the global behaviour, e.g., from slow reaction to ignition, in response to a small change in the experimental operating conditions (known technically as the system parameters) are known as bifurcations of the system. This chapter will aim to show how the concepts and general understanding of non-linear dynamical systems can be used to explain how such bifurcations from one type of reaction behaviour to another arise, and to predict

440

Global behaviour in simple oxidations

Ch. 5

the existence of additional responses not yet observed as the operating conditions are changed. Additionally, some relatively new types of response, such as chaos, will be discussed. If a chemical reaction is evolving chaotically, then traditional methods for predicting the future evolution of the system in time are no longer appropriate and new approaches and statistics are required. Before discussing particular reactions, some basic concepts will be introduced in the context of combustion chemistry. The first is non-linearity, which we will see is actually an "old friend" amongst chemists and is the rule for chemical reactions rather than anything special. The second feature is that of feedback. There are two main mechanisms by which feedback occurs in combustion processes: chemical feedback and thermal feedback and often these occur together in the systems of interest here to produce thermokinetic feedback. Without these feedback processes, none of the phenomena we associate with the science of combustion would exist (and presumably there would not be such an active field as combustion chemistry). Because the coupling of non-linear feedback, and the various other processes such as fluid flow and spatial inhomogeneity that arise in combustion, simplified representations of the chemical feedback processes can be very useful in seeking generalized theories. Some of the more widely exploited models or caricatures of the global chemistry will also be discussed. The previous chapter in this book describes more formal approaches to the representation of complex mechanistic sets of elementary reactions by reduced but accurate and "lumped kinetics".

5.2 NOTATION

symbol a, b a{ a0, b0 A cp E F J

quantity concentrations of species A, B relative third body efficiency of species i inflow concentrations of species A, B pre-exponential factor specific heat capacity activation energy steady-state condition Jacobian matrix

units mol dm 3 mol dm 3 cm 3 molecule - 1 s" J K_1 kg-1 J mol-1

Notation

k P P

s tch *N

'res

T T 1

a

^ad TATT

To V Xi

a

P A>

#ad

© K

A A; V

p T



reaction rate coefficient ("rate constant") pressure vector of parameters = - A H , exothermicity of reaction i rate of reaction / surface area = l/k(T0), chemical timescale = CppV/xS, Newtonian cooling timescale mean residence time thermodynamic (absolute) temperature ambient temperature adiabatic (flame) temperature =E/R, Arrhenius temperature inflow temperature volume mole fraction of species / = a/a0, dimensionless concentration of A = b/a0, dimensionless concentration of B = b0/a0, dimensionless inflow concentration of B = (T - T0)E/RTQ9 dimensionless temperature rise = (Tad ~ T0)E/RTQ9 dimensionless adiabatic temperature rise function in Gray-Yang model parameter in Salnikov model primary bifurcation parameter eigenvalue parameter in Salnikov model stoichiometric coefficient extent of reaction density = t/tch, dimensionless time = *isr/*ch> dimensionless Newtonian cooling time = tTes/tch> dimensionless residence time net branching factor

cm3 molecule Pa Jmol-1 molecule cm - 3 ™2

m s s s K K K K K ™3

m _ -

kgm"

442 X
Ch. 5

Global behaviour in simple oxidations

surface heat transfer coefficient = 0adTN> Semenov number

Wm

2

K

subscripts and superscripts a ambient ad adiabatic cr critical i relating to species / 0 inflow values ss steady-state

5.3 NON-LINEARITY AND FEEDBACK IN CHEMICAL KINETICS: STOICHIOMETRY AND ELEMENTARY STEPS

Chemical reactions can be represented by an overall stoichiometric equation which indicates the relative molecular proportions with which the reactants combine to form the products of the particular reaction. Some simple examples from combustion chemistry are (i)

2H 2 + 0 2 ^ 2 H 2 0

(ii)

2CO + 0 2 - ^ 2 C 0 2

(hi)

C 4 H 10 + 6§0 2 -> 4C0 2 + 5H 2 0

In each case, these stoichiometric equations are appropriate to the complete combustion of the fuel and, hence, assume the availability of sufficient oxidant. If there is less than a stoichiometric proportion of 0 2 , then other partial combustion products such as CO or HCHO etc., may be important in the stoichiometric equation. Thus, there is some sense of pragmatism in specifying whether a particular chemical species is best defined as a product and, hence, featured in the stoichiometric equation or as an intermediate formed perhaps in the early stages but effectively completely consumed again in later stages. Such intermediates do not feature in the stoichiometric equation and do not contribute significantly to the overall thermodynamics of the reaction. The intermediate species may, however, be significant (perhaps even dominant) in determining the

Stoichiometry and elementary steps

443

reaction kinetics, i.e., how fast the conversion of reactants to products is proceeding at any given instant. In general, the reaction stoichiometry tells us nothing about the kinetics. The reaction does not proceed in a single concerted chemical process in which all the molecular rearrangements are achieved in one single stroke. Instead, the chemistry occurs through a sequence of elementary steps with which a variety of other chapters in this volume are primarily concerned. Even for the stoichiometrically-simple oxidation of hydrogen, governed by equation (i) above, it is possible to identify something of the order of 90 different elementary steps that contribute to the overall reaction over a sufficiently wide range of pressure and temperature. The phenomena discussed later in this particular chapter arise over a limited range of experimental conditions, but even so involve approximately 30 elementary steps. These elementary steps involve the intermediate species as well as the initial reactants: in some cases intermediate species are produced, in others they are transformed and in yet others they may combine to form the final products. In the H 2 + 0 2 reaction, even under the limited range of conditions of interest in this chapter, the intermediate species H, O, OH, H 0 2 and H 2 0 2 can be invoked along with the final product H 2 0 . (Under some conditions sufficient H 2 0 2 survives for this to be appropriately classified as a minor product.) Examples of elementary steps from this reaction are (iv)

H2 + 0 2 - ^ H 0 2 + H

(v)

OH + H 2 ^ H 2 0 + H

(vi)

H + 02^OH + 0

(vii)

H + OH + M-» H 2 0 + M

In step (iv) two intermediate species are produced from the initial reactants. At least one of these species, the H atom, is a reactive radical intermediate (it possess an unpaired electron) that plays a major role in carrying the chain of the subsequent chemical transformations. Thus, reaction (iv) is classified as an initiation step (see Chapter 1). Step (v) sees one reactive chain-carrying intermediate, the OH radical, transformed into one other, the H atom, and so it is a propagation step. Step (vi) sees the production of two chain carriers from one: such a step, which brings about

444

Global behaviour in simple oxidations

Ch. 5

an increase in the number of chain carriers, is termed a branching step and we will see that such steps underpin the chemical feedback processes in combustion reactions. Step (vii) is a termination step as the number of chain carriers decreases. The species M in such reaction steps is termed a third body. Its role is to extract some of the energy released in the bondformation process. This acts to stabilize the nascent product molecule, in this case H 2 0 , which would otherwise re-dissociate on the first vibration of the newly-formed bond (see Chapter 2). Any molecule present in the reaction mixture can play the role of a third body, although some are more efficient at participating in particular energy-transfer processes than others and, hence, different species have different third-body efficiencies, and the rate of such steps thus depends on the total concentration of molecules present, i.e., on the pressure and composition of the reacting gas mixture. In the case of elementary steps (only), the reaction kinetics can be inferred from the stoichiometric equations such as (iv)-(vii). Thus, the rate at which step (v) proceeds will be proportional to the local instantaneous concentration of OH radicals and to the local instantaneous concentration of H 2 molecules. We can write, therefore, rv = * v [OH][H 2 ],

(5.1)

where rv is used to denote the rate at which OH and H 2 are converted to H and H 2 0 through this particular step. This term will contribute to governing reaction rate equations for each of these species, i.e., to the differential equations specifying the overall (local) rates of change of these particular concentrations. The proportionality coefficient kv is the reaction rate coefficient (sometimes called "rate constant") for this step. This is likely to be a temperature-dependent quantity (and hence not appropriately called a constant) as described elsewhere in Chapters 1-3 of this volume and to which we return later in this section. The form of equation (5.1) suggests that any bimolecular elementary step will naturally give rise to a term in the reaction rate equations that involves the product of two concentrations. Such a quadratic term in a differential equation provides for a non-linearity and so we see that chemical kinetics naturally produces non-linear terms and equations. Steps (iv) and (vi) are also bimolecular (involving two molecules) and, hence, give rise to quadratic terms: step (vii) gives rise to a cubic term as the total concentration [M] is the sum of the instantaneous individual concentrations, al-

Stoichiometry and elementary steps

445

though this is frequently more conveniently expressed in terms of the local pressure and temperature [M] = p/RT. The elementary steps contributing to a chemical mechanism explicitly specify the (local) rates at which individual species concentrations are varying by summing all the terms of the form of equation (5.1) that contribute to the production, or removal, of the individual species of interest. However, the overall reaction rate is a quantity that must first be defined. In this case, it is appropriate to return to the overall stoichiometric equation (i), for the H 2 + 0 2 reaction. A generalized version of a stoichiometric equations is (viiia) aA + bB -* cC + eE Here the numbers a, b, c and e are related to the stoichiometric coefficients for this overall reaction. These important quantities are obtained by rewriting the stoichiometric equation in the form (viiib)

0 = cC + eE - aA - bB = 2 vjJ J

The final form indicates a summation over all species J appropriate to the particular reaction under consideration. Thus, the stoichiometric coefficient for species C in reaction (viii) is given by vc = c while the stoichiometric coefficient for species A is vA = -a. The stoichiometric coefficients for the reactants and products in the overall reaction (i) are similarly, for the reactants, J^H2=-2,

^o2=-l

for the product ^H2O

= +2

The stoichiometric coefficients for reactants are typically negative while those for products are positive. We can use the stoichiometric coefficients to define the overall rate based equivalently on the rate of consumption of reactants or the rate of production of products. In terms of the general

446

Global behaviour in simple oxidations

Ch. 5

stoichiometric equation (viii), the equivalent definitions of the "reaction rate" R have the form reaction

rate,

« - ifiS , I M , - !4U c dt e dt a dt

=

_ Iffil, b dt (5.2)

or, in general

Vj dt

where vj and [J] are the stoichiometric coefficient and concentration of species J. In some cases, it is appropriate to divide such a rate expressed in terms of a product species such as C by the initial concentration of one of the reactants such as A divided, in turn by its stoichiometric coefficient. The term (^u/^c)([C,]/[A]0) represents the extent of reaction £, which is a dimensionless measure of the progress of the reaction varying from 0 at the beginning of the reaction to 1 at the end (complete conversion). The form of equation (5.2) is useful in that it relates the rates of removal or production of the species in the stoichiometric equation. To specify how this rate depends on the concentrations of the various species involved in the reaction, we must return to the individual elementary steps that comprise the mechanism for this reaction. When dealing with a set of elementary steps, we can use stoichiometric coefficients of the form v{J that specify the coefficient for species ; in step i. The overall rate equation will then involve the rates of the individual elementary steps summed, as they involve the stoichiometric coefficients of the species in terms of which the overall rate has been expressed. As an explicit example, we can consider a truncated version of the mechanism for the H 2 + 0 2 reaction that may be of relevance at very low pressures:

(0)

H2 + 0 2 - ^ 2 0 H

rate = A:0[H2][O2]

(1)

OH + H 2 ^ H 2 0 + H

rate = ^![OH][H2]

(2)

H + 02^OH + 0

rate = fc2[H][02]

Stoichiometry and elementary steps

(3)

O + H 2 -> OH + H

(4)

H-*§H2

447

rate = * 3 [0][H 2 ] rate = *4[H]

The last step here is apparently a unimolecular process and depends only on the concentration of H atoms. This is a step that occurs on the surface of the reaction vessel and the above is a simplification of the rate expression for such a process. We can apparently define the reaction rate R in three different ways: R = - \ ^ 2 at

= ^o[H 2 ][0 2 ] + ^ 1 [ O H ] [ H 2 ] 2 2

+ h3[0][H2]--k4[H], 2 4 R=~

^

(5.3a)

= MH 2 ][0 2 ] + k2[U][02],

(5.3b)

at or R=ld[H£l_l

2

At

2

L

JL

J

(5.3c)

It is not apparent that these are the same and, in fact, in general this will not be the case. The reason for this is that the stoichiometric equation does not consider the intermediate species. If any of these are being formed or removed at significant rates, then the above three choices are not equivalent. If, however, the intermediate species manage to establish a quasi-steady-state in which their rates of formation and removal effectively balance, then these three different definitions of the overall rate become equivalent. This point has been discussed as an example in Chapter 4. Even if the quasi-steady-state approximation is valid, all of the forms of the reaction rate equation (5.3) given above involve the concentrations of at least some of the intermediate species. If we are to be concerned with the intermediate species, and particularly when we wish to employ the methods of non-linear mathematics (or even direct numerical computation) then we must find ways of including such quantities. The appropriate rate equations for the concentration of H, O and OH are readily constructed

448

Global behaviour in simple oxidations

Ch. 5

in a similar manner. There are various ways forward. The six rate equations (one for each different chemical species) can be integrated numerically forwards in time from a given set of initial concentrations and for a particular set of the experimental parameters. Such computations could then be repeated for different parameter values or initial concentrations and a map of the global behaviour in the parameter space be constructed. Experimentally, this is equivalent to determining the pressure-temperature ignition diagram for a particular mixture and vessel size. We may note that it is not, in theory, necessary to integrate all six equations: there are two constraints that arise from the conservation of atoms as an H atom is never converted to an O atom in a chemical process. Thus we know that there is a total atom balance of the form 2[H2] + [H] + [OH] + 2[H 2 0] = 2[H 2 ] 0 ,

(5.4a)

2[0 2 ] + [O] + [OH] + [H 2 0] = 2[O 2 ] 0 ,

(5.4b)

relating the instantaneous concentration of various intermediates to the initial concentrations of the reactants. For complex situations, such as for spatially-inhomogeneous systems, this atom balance must be expressed as an integral over all space and may not hold locally and in many cases, even for well stirred systems, it is often not particularly helpful invoking these constraints. The quasi-steady-state approximation works by replacing the differential equations for the rates of change of the intermediate species by algebraic conditions obtained by setting d[H]/df = 0 etc. (see Section 4.8.5). In some cases, the resulting equations can be solved and manipulated algebraically allowing substitution into the overall rate equation to obtain a form that only involves explicitly the concentrations of the reactants and (perhaps) products. Such rate equations can then be compared with the empirical rate equations determined from experiment to test the validity of the assumed mechanism and to obtain quantitative values for the rate coefficients involved. The non-linear mathematics approach begins with the set of governing differential equations but puts great emphasis in the initial stages of analysis on the true steady-states of such systems, i.e., the combinations of concentrations for which all the rate equations, including those for the reactants and products, become zero. This is only appropriate in so-called

Stoichiometry and elementary steps

449

open systems in which there are inflows of fresh reactants. In closed vessels, the system "simply" has a single steady-state corresponding the state of chemical equilibrium: such states are basically "uninteresting" in the present context, although the way in which the system evolves "early on" in its approach to the equilibrium state may be "interesting". Other approximate approaches can be developed to analyze such situations. The simplest type of open system of interest in combustion is the continuous-flow, well-stirred tank reactor or CSTR, which is an idealization of tank reactors used widely in industry. In essence, this is simply a tank into which reactants flow continuously at some known volumetric flow-rate and the reactant-intermediate-product mixture is efficiently stirred so that there are no spatial concentration or temperature gradients. In order to maintain a constant reaction volume, there is a matching volumetric outflow of the mixture from the CSTR so that molecules spend only a finite time in the reactor. This is known as the mean residence time tTes and is determined by the volumetric flow-rates and the reactor volume. In terms of the governing mass balance equations which again determine the rates of change of individual species concentrations then, under certain mild assumptions, there are simply additional flow terms added to the kinetic terms. For the reactant species there will be both inflow and outflow terms while for intermediates and products there are likely to be only additional outflow terms as these species do not have significant concentrations in the inflow mixture. The nature of these extra flow terms can be exemplified with the simplified H 2 + 0 2 mechanism discussed above. For the fuel H 2 , the mass balance equation appropriate to reaction in a CSTR can be written in the form

ffl=ffl±da_yH!l0il at

tres

- ^[OH][H 2 ] - * 3 [0][H 2 ] + ik 4 \H\.

(5.5)

The first term on the righthand side of this equation gives the net inflow rate of H 2 : here [H 2 ] 0 is the "inflow concentration" of the fuel. This must be specified with some care if, as is likely, there are separate inflows of fuel and oxidant (and perhaps also of other diluents) to the reactor. The appropriate definition of the "inflow concentration" in such equations is,

450

Global behaviour in simple oxidations

Ch. 5

in fact, the concentration of the fuel (or any other species) that would be established in the reactor under the prevailing flow conditions but in the absence of chemical reaction. In other words, we must allow for the dilution effects that arise when more than one inflow stream is mixed in the reactor. For intermediate or product species, such as H 2 0 , there is typically no inflow term so the mass balance equation has the form

ffioi.,i[OHIH2]-a°i, ill

(5,6) *res

with an outflow term appended to the kinetics. We can write a general form (cf. equation (4.3)) for the mass balance equation of any species J

afl.SBkzJffi.s^, Qt

tres

(5.7)

i

where [J]0 is the inflow concentration as defined above (which may be zero for intermediates and products) and the final term is the sum of the kinetic terms with vx j being the stoichiometric coefficient for species J in reaction step i which has a rate rt that will depend in some (typically non-linear) way on the concentrations of the various species in the reactor. If the mechanism is not known in detail, the kinetic terms may be replaced by empirically-determined rate laws, i.e., by approximations to the reaction rate term that typically will be some (non-linear) polynomial fit of the observed rate to the concentrations of the major species in the reaction (reactants and products). Such empirical rate laws have limited ranges of validity in terms of the experimental operating conditions over which they are appropriate. Like other polynomial fitting procedures, these representations can rapidly go spectacularly wrong outside their range of validity, so that they must be used with great care. If this care is taken, however, empirical rate equations are of great value. In many respects, it is surprising that complex chemical mechanisms often underlie relatively simple empirical rate expressions. This is a consequence of the way that the individual elementary steps, of which there may be several hundred, frequently appear to group themselves into a much smaller number of "timescale bands". Then, the kinetic tune may be called, to a great extent, by only one or two such groups, each of which

Chemical feedback: Branched-chain ignition

451

introduces only a simple term to the approximate rate equation. This is discussed more formally in Section 4.8. Examples of such empirical rate laws for the H 2 + 0 2 reaction are found in the literature [1-4]: rate = fc[H2]2[02], rate =

fc[H2]12[02f8,

rate = fc[H2][02]1/2, rate =

k[H2][02]2,

each for different pressure, temperature, mixture composition and vessel surface preparations. For our general reaction (viii) we might be able to provide an empirical rate equation of the form ld[C] c At

=

ld[£]= _ld[A]= e At a At

ld[£] b At

= kv^Ar[BY[CY[EY.

(5.8)

Here the coefficients a-e are the individual (empirical) reaction orders with respect to the species A-E, respectively. The overall reaction order would then be given by a + fi + \ + e> but this may not be a particularly significant quantity. Other possible (and indeed quite common) forms for the empirical rate law involve rational polynomial terms such as [A]a/(1 + ^[C]*) in which case the concept of an overall order, or of an order with respect to C, is not appropriate.

5.4 CHEMICAL FEEDBACK: BRANCHED-CHAIN IGNITION

Empirical rate laws are particularly useful for introducing the idea of feedback. For this, we can proceed by plotting a graph of how the overall reaction rate varies with the degree of advancement of the reaction. The latter can be simply represented by the extent of reaction described above, i.e., by the ratio of how much of a selected reactant has been used up with respect to its initial concentration. In terms of species A in our general reaction, which for convenience we can simplify here to

452 (ix)

Global behaviour in simple oxidations

A + B -» C + D

Ch. 5

rate = f([A], [B])

where we have restricted the rate law to be some function of the reactant concentrations only, we could sensibly choose the extent of conversion as

f-W£4

(5.9)

The form of the dependence of R on £ will be determined by the order of the reaction with respect to A and B. Provided the reaction does indeed follow the overall stoichiometry in the above equation (in this particular case we have vA = vB = - 1 ) then, once we know the initial concentrations of A and 5 , the concentration of B at any time can be determined uniquely if the concentration of A is known, i.e., these concentrations are not independent, and so the rate can be formally expressed as a function solely of one concentration or extent of reaction variable. For example, if the reaction is first order in A (a = 1) and zero order in all other species (/3, *, e = 0) then the rate law R = k[A] = k[A]0(l - £) gives rise to a linear relationship as indicated in Fig. 5.1(a). This is a rather special case and, in the context of the previous section, we have a reaction that shows linear dynamics rather than non-linear dynamics. For any other rate law, the plot of R versus £ will give rise to a curve. Figure 5.1(b) shows the appropriate form for a second-order reaction, e.g., with a = j3 = 1, [A]0 = [B]0, or a = 2 /3 = 0, and also shows the curves for a number of other overall reaction orders of the general form R = kix[A]o (1 - £)", with n not necessarily specified as being an integer. The curves in Fig. 5.1(b) with n + 1 are "non-linear" but do share the feature that the rate is highest at the beginning of the reaction (£ = 0) and the rate decreases monotonically in each case as the extent of reaction increases. This feature is characteristic of deceleratory reactions. As the number of elementary steps and intermediate species (particularly if these involve reactive species such as radicals) increases, so the possibility of more "interesting" shapes for the reaction rate curve increases. Two such "interesting shapes" are illustrated in Fig. 5.1(c). These are characteristic of reactions that display an acceleratory phase at low extents of reaction before attaining a maximum prior to a final deceleratory phase at high extents at the end of the reaction as the state of chemical equilibrium is

453

Chemical feedback: Branched-chain ignition

k[A]0 11=1/2

06

3

2

o

-a o

Fig. 5.1. Variation of reaction rate R with extent of reaction £: (a) linear relationship for first-order reaction; (b) non-linear deceleratory reactions of overall order n; (c) reactions showing chemical feedback in the form of autocatalysis; (d) comparison of chemical and thermal feedback curves.

approached. These two curves have approximately "quadratic" (parabolic) and "cubic" character, and many observed rate curves for reactions that have some form of chemical feedback mechanism can be approximated quite closely by one or other, or by some linear combination of these two curves. In order to give these idealized feedback curves a "chemical face", they are frequently represented in terms of autocatalysis [5, 6]. A purely quadratic autocatalytic curve would arise if the A + B reaction does not give rise to products C and Z), but instead, produces two molecules of species B quadratic autocatalysis

A + B —>2B

rate = A:^[A][5]

We can think of this as a reaction in which the conversion of a reactant species A to an intermediate or product B is assisted or "catalyzed" by B

454

Global behaviour in simple oxidations

Ch. 5

itself. An example of how such a rate law arises in the H 2 + 0 2 will be presented below, but it should be emphasized that the above representation is simply a convenient shorthand rather than a suggestion that any species B directly combines with a reactant A to reproduce. If we start with pure A so that [B]0 « 0, then the reaction rate will be very low (R would actually be equal to 0 if [B]0 is actually equal to zero and there is no other route from A to B). As A is converted (we might initiate the reaction by adding a small trace of B), then [A] decreases as [B] increases with [A] oc (l - £) and [B] <* £ so that we obtain quadratic autocatalysis

R = kq[A]l^(l - £)

giving rise to the parabola in Fig. 5.1(c). In the initial phases of the reaction the deceleratory effect of the decreasing concentration of A is more than compensated for by the increase in the rate as [B] increases. These effects balance as the two species concentrations become approximately equal and, in the later stages, the consumption of the initial reactant is more significant in determining the sign of the slope of R. In a similar way, the feedback mechanism giving rise to the cubic curve can be given the stoichiometric representation cubic autocatalysis

A + 2B —>3B

rate = fcc[y4][l?]2

We still have [A] « ( l - £) and [B] « £, so now cubic autocatalysis

R = kc[A]l^2(l — £)

The higher order, with respect to the autocatalyst, skews the rate curve so that the maximum lies at higher extents of conversion and there is a longer induction phase during which the reaction rate is close to zero at low extents of conversion. Cubic autocatalysis is apparently less significant than quadratic which is relatively common as chemical feedback in combustion systems, although cubic-type curves have been reported and exploited in the oxidation of H 2 for which a rate expression of the form d[H20]/d£ = &[H2][H20]2 was observed [7] and also the oxidation of CS2 in heavilydiluted air mixtures [8]. An important example of chemical feedback in combustion reactions is provided by the so-called branching cycle in the H 2 + 0 2 reaction [9] which

Chemical feedback: Branched-chain ignition

455

is comprised of steps (l)-(3) given earlier. Under typical conditions, the slowest or rate determining step in this sequence is step (2) involving the H atom and 0 2 molecule. The net stoichiometry of the branching cycle is obtained by taking step (2) + step (3) + 2 x step (1) which may be written as, 3H2 + 0 2 + H -» 3H + 2H 2 0

rate = fc2[H][02]

The H atom is retained on the lefthand side of this stoichiometric equation to emphasis its role (because its concentration occurs explicitly in the rate law). Thus, the overall reaction can be regarded as the conversion of H 2 to 2H, "catalyzed" by H and with the simultaneous conversion of 2H2 + 0 2 to 2H 2 0. The latter provides the free energy driving-force for the process. In terms of the reaction stoichiometry, the stoichiometric coefficient for H atoms is now +2 (i.e., +3 - 1) so that H is a "product" but also appears with a positive exponent (reaction order) in the rate expression. Thus, it clearly plays a role that parallels B in the simplified representation with A corresponding to 0 2 as the latter appears as the "reactant" in the rate law. In the simple mechanism for the H 2 + 0 2 reaction given earlier, the branching cycle which increases the H-atom concentration competes with the termination step (4) which removes H atoms from the system. This competition appears more clearly if we make steady-state assumptions on O and OH (see Section 4.8 for a justification of this procedure) and then substitute for these into the rate equation for d[H]/d£, to give ^ 1 = 2* 2 [H][0 2 ] - k4[U] = [H], at

(5.12)

where = 2k2[02] - k4 is termed the net branching factor. If > 0, i.e., if the coefficient for branching exceeds that for termination, then d[H]/dr is positive and there is (exponential) growth in the H-atom (and other radical) concentration. The condition = 0 is important in determining the classic p-Ta ignition limits as we discuss later in this chapter. A simplified and more generalized representation of the chain-branching chain-termination competition involving a single reactant A, an intermediate X and a final product P is the three-step scheme

456

Global behaviour in simple oxidations

initiation

A—>X

rate = k{ [A]

branching

A + X-+2X

rate = fc*[i4][Z]

termination

X—>P

rate = kt[X]

Ch. 5

The concentration of the reactant [A] is frequently assumed to be constant (this is known as the "pool chemical approximation"). Only one reaction rate equation, for the chain carrier X, need be considered. This has the form M = r{ + kb[A][X] - kt[X] = r{ + 4{X]9 at

(5.13)

where r{ = k{[A] is the initiation rate and is the net branching factor as before. Provided (/> + 0, this integrates to give [X] = ^(e+-l). 9

(5.14)

If kt > kb[A], then < 0 and the chain carrier concentration evolves to a steady-state value given by [*]ss=-ri/0.

(5.15)

If is not small, this will be a low steady-state concentration, but as cf) tends to zero so [X]ss increases. For 0, [X] grows exponentially with time again until the consumption of A must be taken into account. The condition for an acceptable steady-state then is clearly (/> < 0, with the condition 0 = 0 being the "critical" case separating steady-state from explosive runaway. The "critical" concentration of the reactant obtained from the condition <\> = 0 is then [A]cr = kt/kb.

(5.16)

Although this mathematical manipulation leads to an explicit criterion

Chemical feedback: Branched-chain ignition

457

for the "ignition limit", and has contributed immensely to the understanding of branched chain reactions, this approach has been criticized in one respect. As is increased towards zero from some negative value, so [X]ss increases smoothly towards infinity, as indicated in Fig. 5.2. This seems to contrast the "expected" behaviour in chemical systems that show a "steady-state rate" that generally increases as the system parameters are brought close to the critical values, but then show a sharp, discontinuous change in behaviour as the limit is crossed. We will see examples of such discontinuous bifurcation structures with other models below but should also note that in closed systems there is, in fact, no such discontinuous change, rather the rate simply increases, usually very rapidly, over a typically very narrow range of experimental conditions (so narrow and with such a dramatic increase that it appears to be discontinuous). The interpretation of "ignition" in closed systems is still a subject of academic study. The initiation step is necessary in the above system to maintain a nonzero steady-state chain carrier concentration if < 0. The problem can be re-specified by setting r{ = 0 and examining how [X] varies if we have a non-zero initial chain-carrier concentration [X]0. If c/> < 0, the [X] decreases exponentially to zero; if >0, then [X] increases exponentially. So again, we see that = 0 is a "critical" value separating qualitatively different responses.

|[ IX].

net branching factor,

ty

0

Fig. 5.2. Variation of the quasi-steady-state radical concentration with net branching factor for the simple model of chain branching and chain termination showing [X]ss -» °° as -* 0 (note: the steady-state only exists for (f> < 0).

458

Ch. 5

Global behaviour in simple oxidations

If the same model is examined in a CSTR, then the governing equations are

d[A]_([A]0-[A]) Qf

d\X]

k^A] - kb[A][X],

(5.17a)

*res

=

_[X]

(5.17b)

+ k-M] + kb[A\[X\ - k,[X\.

In this case we can incorporate the consumption of the reactant and still study (true) steady-state behaviour. Setting d[A]/df = d[X]/df = 0 we obtain [A]0 - [A]ss = (i + A^ res )[;r] ss ,

(5.18)

and, hence, the steady-state condition become

ki[A]0 + { kb[A]0 - ki(l + kttrcs)

k, l res

in

- kb(l + kttres)[Xfss = 0.

(5.19)

The first term in this quadratic equation is the initiation reaction rate based on the inflow concentration of the reactant. The coefficient for the term in [X]ss has something of the character of the previous net branching factor. The above equation has a single positive solution for any set of rate constants, residence time and inflow concentration: a typical variation of [X]ss with [A]0 is shown in Fig. 5.3(b) and shows a rapid increase in the vicinity of some "critical" concentration [^4]o,crThe behaviour can be quantified if we make the approximation of ignoring the (probably small) initiation terms, setting k{ = 0. The steady-state condition can then be written in the form fiow is given by

(5.20)

459

Chemical feedback: Branched-chain ignition

X

(a)

(b)

[A]„,r

[A] 0

Fig. 5.3. Variation of steady-state radical concentration as a function of reactant concentration appropriate to chain branching and chain termination reactions in a flow reactor (a) ignoring initiation reaction, there is a root [X]ss = 0 for all [A]0 and a non-zero root that is positive for [A]0 > [A]0,cr; (b) including initiation step provides a non-zero positive radical concentration for all [A]0 but the radical concentration increases rapidly over a relatively small range of [A]0.

0f low = kb[A]0 - kt

1 ^res

= 0

1

.

(5.21)

'res

The flow term adds an extra chain-carrier loss term with a pseudo-firstorder rate constant given by the inverse of the residence time. Equation

460

Global behaviour in simple oxidations

Ch. 5

(5.20) has a root [X]ss = 0 for all parameter values, with a non-zero root whose sign is determined by the sign of f low . For $ f low < 0, the non-zero root is negative and, hence, physically unacceptable. For fiow > 0, the non-zero root is positive and physically acceptable. The character of the two roots in this case is such that [X] will evolve to the non-zero root from any non-zero (positive) initial value. The variation of the two roots [X]ss with [A]0 is shown in Fig. 5.3(a). There is a transcritical bifurcation where the two roots become equal and their loci cross at [A]0 = [A]0,cr which corresponds to the condition for f low = 0 [A]o,a = (kt + —)/kb. \

(5.22)

^ res' '

Provided rt is small, then the critical inflow concentration for this branching-termination model under CSTR conditions differs slightly from the so-called "pool chemical" result which is obtained by assuming [A] = constant. For typical chemical systems the residence time will be such that kt > VtTes, so the two results are not significantly different but the extra influence of the flow is clearly evident in the above forms. In neither of the analyses above, however, is there a discontinuous jump in the steadystate response as the parameters are varied.

5.5 THERMAL FEEDBACK: IGNITION, EXTINCTION AND SINGULARITY THEORY

Very few chemical processes are precisely thermoneutral and combustion reactions are characteristically highly exothermic processes. If the reaction rate is not negligible, therefore, there will be a significant nonzero rate of chemical heat evolution. Unless this energy is transferred equally rapidly from the system, it will give rise to local increases in the temperature of the reacting mixture above that of the surrounding fluid or ambient heat reservoir. This thermal effect is highly significant because of the other characteristic feature of combustion reactions, i.e., the high sensitivity of the reaction rate coefficients to temperature. A common form of representing the temperature dependence of reaction rate "constants" is the Arrhenius equation k = Ac'E/RT.

(5.23)

Thermal feedback

461

Here A is the pre-exponential factor, E is an effective activation energy or temperature coefficient for the overall reaction, R = 8.314 J K _ 1 m o l _ 1 is the Gas Constant and T is the local absolute temperature. The quotient EIR has units of temperature and is sometimes known as the Arrhenius temperature TArr. It is a characteristic feature of combustion processes that if Ta is the ambient temperature, then TArr > Ta so that the group RTJE and, more generally, the scaled or "dimensionless" temperature RTIE are typically (very) small quantities in the systems of interest. The evolution of the local temperature in a reacting system is governed by the heat balance equation. This may become quite complex in unstirred systems, especially if convective heat transfer processes develop as a consequence of local heating. For the simple CSTR described earlier we can proceed with an ordinary differential equation of the form ( c ^ ^ ^ C c ^ ^ ^ ^ +I^^-^Cr-Ta). at tTes t V

(5.24)

Here cpp is the heat capacity per unit volume and is assumed to be independent of temperature in this simplified formulation, T0 is the temperature of the inflowing reactants. The second term on the righthand side is the sum of the rate of chemical heat evolution due to the individual elementary steps: q{ is the exothermicity (-A//: typical units = kJmol - 1 ) of the ith step and plays a similar role to that of the stoichiometric coefficients in the mass balance equations and rx is, as before, the rate of the ith step and, hence, depends on various species concentrations and, now, also on the temperature T of the reacting mixture. The final term is a Newtonian heat transfer term indicating heat transfer from the reacting mixture at temperature T to the surrounding heat bath at the ambient temperature Ta subject to a heat transfer coefficient \ a n d a surface-tovolume ratio S/V. For a well-lagged reactor we may approach the adiabatic case ^ = 0, for which the only heat transfer process is the inflow of reactants that are cool relative to the outflowing fluid. For the simplest chemical case imaginable, that of a single first-order reaction A—>B

rate = ka

exothermicity = q

where we use a for the concentration of A and the rate constant k has the

462

Global behaviour in simple oxidations

Ch. 5

Arrhenius form given above, the two governing equations for reaction in a CSTR have the form ^ =^ ^ - ^ K At tres

(5.25a)

(cpp) ^ = (cpP) H^H at tr^

+ qk(T)a

- & (T - r a ) . V

(5.25b)

This provides a pair of coupled, non-linear (through the Arrhenius temperature dependence) ordinary differential equation for the two variables a and T. If the temperature increases, the reaction rate increases through the increase in k. The consequent increase in T will lead to increases in the heat transfer rates and also to a decrease in the concentration of A, which in turn tends to decrease the reaction rate term ka. To quantify this effect, we can examine the adiabatic case x= 0. In this situation, the temperature rise above the inflow is uniquely linked to the extent of reaction £ = (a0 - a)/a0 through the condition T-T0 ad

a0-a ?0

=

^

(5 26)

00

Here, Tad - T0 = qa0/cpp is the temperature rise that occurs under adiabatic conditions accompanying complete consumption of the reactant A. Because of this relationship, the temperature and the reactant concentration are not independent (this relationship is only true under adiabatic conditions), and so the instantaneous reaction rate can be expressed in terms of one of these quantities alone, e.g., we can express R as R(€) as in the previous section. A typical form for the reaction rate curve is shown in Fig. 5.1(d). This shows the nature of thermal feedback even in this very simple chemical example. At low extents of reaction, the increase in k as T increases dominates and the rate, which is relatively low when £ = 0 so that T = T0, increases as the reaction proceeds. Only at very high extents of reaction, i.e., close to complete consumption of the reactant, does the rate fall, approaching zero as £ —> 1. The curve has some similarities with the cubic autocatalytic curve but thermal feedback tends to have a shallower initial development and its maximum at higher £.

Thermal feedback

463

The governing mass and heat balance equations can be usefully recast in the following dimensionless forms: da dr

1— a .,. a/(0), T

/e. ^ x (5.27a)

res

^ = eadaf(d) " ( — + - ) 0. ar \r r e s TNJ

(5.27b)

Here 6= (T - T0)E/RTl is a dimensionless measure of the rise in the temperature above the inflow value (for simplicity T0 = Ta has been assumed). The unspecified function f(0) then reflects the increase in the rate coefficient k above its value when T= T0 that arises for a particular dimensionless temperature excess 0, i.e., k(T) = k(T0)xf(d).

(5.28)

The Arrhenius temperature dependence can be represented exactly in this way with /(0) = exp — — . U + eflJ

(5.29)

For the purposes of qualitative discussion, however, this form can be usefully approximated by the simple exponential dependence f(0) - e e ,

(5.30)

which holds as the quantity e = RT0/E = T0/TArT is typically very small compared with unity as discussed earlier. The quantity 0ad = (T ad - T0)E/ RTl is the dimensionless adiabatic temperature excess and the dimensionless concentration used here a = a/a0 = 1 - £. The remaining dimensionless quantities Tres and r N are the dimensionless residence time tres/tch and Newtonian cooling time tN/tch, where tN = cppV/xS and the chemical timescale tch = llk(T0) is the inverse of the rate constant evaluated at the inflow temperature T0. Using the relationship between the temperature rise and extent of con-

464

Global behaviour in simple oxidations

Ch. 5

version determined earlier we have for the adiabatic case (i.e., with TN^oo)

^ - = l - « = £.

(5.31)

The reaction rate term R — k{T)a can thus be written as R = i? 0 (l " €) e"ad*,

(5.32)

where R0 = k(T0)a0 and gives rise to the form in Fig. 5.1(d). The behaviour of this model has been widely studied in the chemical reactor engineering literature, where it is a standard form. Under non-adiabatic conditions it gives rise to complex steady-state responses and to sustained oscillations in the concentration and temperature. A detailed description can be found elsewhere (see e.g., Chapters 6 and 7 in reference [5] and references therein). In the "pool chemical" formulation of this model, the consumption of the reactant A is ignored and so there is only one variable, the temperature T. The governing heat balance equation has the form (cPp) ^ = qk(T)a0 ~^(Tat V

Ta),

(5.33)

there being no flow terms and the reactant concentration is now constant and equal to its initial concentration a0. The dimensionless form of this equation can be written as

f = m~e-r

(5-34)

where ifj is known as the Semenov number and is given by 0adTN- This model is the basis for the theory of thermal explosion or thermal runaway because of chemical self-heating [10]. The neglect of reactant consumption is generally justifiable due to the large magnitude of the dimensionless adiabatic temperature excess 0ad, which indicates that only small extents

Thermal feedback

465

of reaction are need to achieve temperature rises that make 0 of order unity. Criticality for ignition in either of the two formulations of this simple model of thermal feedback can be interpreted by means of a thermal diagram. The variation of the rate of heat evolution with temperature for the simple Semenov model above is, effectively, the graph of the function f{6). For relatively low-temperature excesses, the curve is well approximated by the exponential form/(0) = t6 and so has the appearance shown in Fig. 5.4. The heat transfer rate corresponding to the term dlifjin equation (5.34) is simply a straight line on the thermal diagram. For low values of the Semenov number, corresponding to low reaction exothermicities or to high heat transfer rates, the loss line is steep and intersects the heat release line twice on the diagram as shown. The system starting with 0 = 0, i.e., with T= Ta initially, will evolve to the lower of these intersection points at which d0/dr = 0 and, hence, a steady-state is attained. The corresponding steady-state temperature excess 0SS will have a value less than unity, typically indicating a temperature excess of less than 10 K. As if/ is increased, so the slope of the heat loss line decreases and for high if/, there are no

\|/ small

m

o

e

Fig. 5.4. Thermal diagram for Semenov model of thermal explosion: the rate of chemical heat release varies with the dimensionless temperature excess 0 according to/(0) « e e ; the rate of heat transfer is given by the straight line with a gradient of 1/I/J. For small ifj the loss line is steep and makes two intersections corresponding to two steady-states; for large ip the loss line has a low gradient and does not allow steady-state intersection points; the critical case corresponds to tangency of the heat release and heat loss lines.

466

Ch. 5

Global behaviour in simple oxidations

intersections between the two curves on the diagram. In this case, the heat release rate always exceeds the heat loss rate so the extent of self-heating increases continuously until a high temperature rise is attained (at which point the assumption that the reactant concentration remains constant becomes totally invalid). This latter behaviour corresponds to thermal runaway. The "critical" case separating these two different responses (low 6 steady-state for low i/> and thermal runaway for large if/) arises when the heat release and heat loss curves just touch tangentially. The two intersection points at low i/> approach each other as \\f is increased and merge at the tangency condition. These steady-state intersections vanish for larger \\f. If we plot the variation of the steady-state intersection points with the parameter \\f, to give the bifurcation diagram as shown in Fig. 5.5, then we see two branches that meet at a vertical turning point at the point of tangency. (There is an additional intersection point giving rise to a third branch at large 0: for the Semenov model in which reactant consumption is ignored, this branch relies on the "saturation" of the Arrhenius function at large T>TArr and is totally physically unrealistic. For the CSTR equations, consumption is included explicitly and the third branch of intersection points corresponds to a high, but physically acceptable, temperature and an "ignited" steady-state.) Following the steady-state behaviour then, if we begin with some low value for I/J, the system will evolve to the lowest branch in Fig. 5.5. If if/ is increased slowly, the system effectively moves along the lowest branch, with the steady-state temperature excess increasing slowly and smoothly

9s*cx

Vex

V

Fig. 5.5. Variation of the steady-state temperature excess 0SS with the Semenov parameter \ft indicating the turning point at the critical condition.

Thermal feedback

467

with the parameter. At the turning point, however, the system must jump away from the low 6 steady-state solution since this solution vanishes as tangency occurs. Thus, there is a discontinuous jump in the steady-state of this model, corresponding mathematically to a genuine saddle-node bifurcation point. Such a genuine discontinuity does not arise in the simple quadratic chain-branching model, and this absence leads to the earlier (but probably unwarranted) criticism alluded to in the previous section. (In fact, if we allow for reactant consumption, which must clearly be a feature of the real situation, the distinct bifurcation structure vanishes from the Semenov model too as the only steady-state for any i/> is that with all the reactant consumed and the temperature returned to the ambient. Only open systems have true bifurcations between non-(chemical) equilibrium steady-states.) The critical condition can be located in a straightforward manner using the simultaneous conditions for a steady-state and for tangency of the two curves steady-state condition

f{&) = - ,

tangency condition

= —, d0

(5.35a) (5.35b)

\fj

with the latter being the condition for equal slopes for the heat release and heat loss curves on the thermal diagram. If we take f{6) = e 0 , then these give 0ss,cr = l,

iAcr = e _ 1 .

(5.36)

For the CSTR model, the highest intersection point is physically realistic and, indeed, important. It corresponds to an intersection point on the deceleratory part of the reaction rate curve as indicated in Fig. 5.6(a, b). There are now two possible tangencies of the heat loss line with the heat release curve as some parameter such as the residence time is varied, as indicated in the figure. The first has a similar implication as in the Semenov case. It corresponds to the merger of two low-lying steady-states and to an ignition point on the steady-state locus and, in this model, arises typically as the residence time is increased. The system now jumps to a high steady-

468

Global behaviour in simple oxidations

Fig. 5.6. Bistability for an exothermic reaction in a flow reactor: (a) flow diagram for system with large adiabatic temperature excess, 0ad > 4 showing heat release curve and indicating the position of the flow line corresponding to ignition and extinction events; (b) corresponding variation in steady-state temperature excess with mean residence time showing region of bistability: the "jumps" associated with the tangencies in (a) are indicated by vertical arrows at the turning points; (c) and (d) flow diagram and steady-state locus for a system with 0ad < 4 for which multistability has been unfolded.

state temperature excess, corresponding to the uppermost branch in the figure. If the experimental conditions are now changed, e.g., by decreasing the residence time, we move back to the left on the steady-state locus, but remain on the upper branch. Only at some lower residence time do the upper and middle branches coalesce at an extinction turning point. This is another tangency (saddle-node) bifurcation between the heat release and heat loss curves, as indicated, and sees the system jump back to a low temperature steady-state. Thus, this model shows both ignition and extinction and accounts for a range of hysteresis in the parameter space over

Thermal feedback

469

which different steady-states coexist. The system can sit at either steadystate, depending on its previous history. If some other parameter of the system, such as the adiabatic temperature excess 0ad is varied, so the shape of the steady-state locus may deform. For low values of 0ad in this model, corresponding to weakly exothermic processes, then the hysteresis loop is unfolded, as indicated in Fig. 5.6(c, d), and a simple smooth variation of the steady-state temperature excess with the residence time is observed. Thus, systems can lose criticality as other experimental parameters are changed. The recipe of locating the critical (ignition or extinction) conditions through tangency in a thermal or flow diagram of the form shown in the previous figures or equivalently as a vertical turning point in the steadystate locus can be written as a general prescription [5]. The steady-state condition can be written in the form F(xss, A; p) = 0,

(5.37)

where F is some set of functions, x is potentially a vector of the variables in the system, A is the primary bifurcation parameter, i.e. the parameter chosen as that most likely to be varied in a given set of experiments and p is a vector of the remaining unfolding parameters. For criticality, this condition will be satisfied along with the condition for a vertical turning point in the xss-X curve, which occurs when dF(xss, A; p)/dx = 0,

(5.38)

where the form dF/dx indicates the matrix obtained by differentiating the vector of functions F with respect to the variables. If there is more than one independent variable in the model, then equation (5.38) represents a condition on the determinant of the Jacobian matrix, det(J) = 0. To illustrate this, we can apply this prescription to the two variable CSTR model for which x = (a, 6)T. The steady-state condition F is simply obtained by setting da/dr = 0

Fx = — - a / ( 0 ) = O, ''"res

(5.39a)

470

Ch. 5

Global behaviour in simple oxidations

F2 = eadaf(6) - ( — + — ) 0 = 0.

(5.39b)

The Jacobian matrix has the form

3a

30

§

§

da

30,

Tres

L 0ade.«e

.

/ 1 +. J-) i\

0 a d afle « - ( J -

(5.40) taking f(6) = e e , for simplicity. If we also restrict ourselves to the adiabatic case so that l/r N = 0, then the tangency/determinant condition becomes det(J) = I — + te

1

dadace

*

+ 0 a d ae z " = O.

(5.41)

Using the steady-state equations, this can be reduced to a quadratic equation in 0ss?cr with roots 0ss,cr = koad ± V0 ad (0 ad - 4)}.

(5.42)

This indicates that we need 0ad > 4 for there to be real roots. The roots correspond to the critical temperature excess at the points of ignition (lower root) and extinction (upper root) respectively. As 0ad is decreased towards the transitional value 0ad,trans = 4, so the two turning points approach each other and merge as the hysteresis loop unfolds (Fig. 5.6).

5.6 THERMOKINETIC FEEDBACK: OSCILLATIONS AND LOCAL STABILITY ANALYSIS

In many combustion systems there will be the possibility of both thermal and chemical feedback, with the two processes coupled together through the heat and mass balance equations just for added fun. A simple model

Thermokinetic feedback

471

involving both effects would be that of an exothermic, quadratic autocatalytic reaction with an Arrhenius temperature dependence A + B -* IB

rate = k{T)ab

for which the governing equations would have the form da _ (a0 - a) k(T)ab, dt tre$ db = (fro - b)

(5.43a)

+ k(T)ab,

(Cpp) ^ = (cpP) at

(T

T)

°

tres

(5.43b) + qk(T)ab - ^ (T - T a ), V

(5.43c)

or, in dimensionless forms noting that a0 + b0 = a + b from the reaction stoichiometry, ^ = — - - «(1 + A, - * ) / ( * ) , dr r res

(5.44a)

^ = 0 ad a(1 + ft " a)/(fl) - ( — - — ) 0, dr

\Tres

(5.44b)

TN/

where a, 0, etc. are as defined before and j80 = bja0. This system has been investigated using singularity theory [11,12]. The condition for thermal runaway under the assumption a = a0 in a closed system was determined by Frank-Kamenetskii [13] as i/fcr = 4e _ 1 in terms of the Semenov number described earlier. Another model studied under the name of thermokinetic feedback is due to Salnikov [14,15]. This has two first-order reaction steps which, in general, could both be exothermic and both have an Arrhenius temperature dependence: (51)

A^X

rate = M

(52)

X-*P

Tate = k2x

472

Global behaviour in simple oxidations

Ch. 5

The most interesting behaviour can be obtained, however, by taking a slightly simplified version in which only step (S2) is exothermic and has a temperature-dependent rate law, so that the governing equations have the form: dx — = k1a-k2(T)x, dt (cpp) ^ = qk2(T)x - & (T - Ta), at V

(5.45a) (5.45b)

and the concentration of the reactant a is assumed to be constant. These equations can again be recast in dimensionless terms, with the following form being particularly convenient [16a]: ^ = fi-Kyf(d), dr

(5.46a)

An

— =yf(B)-e. dr

(5.46b)

Here, /x is a dimensionless measure of the initial reactant concentration a0, y is the dimensionless concentration of X and K is the ratio of the Newtonian cooling time to the chemical reaction time at ambient temperature. Equations (5.46a, b) have a steady-state solution which, taking f(8) = ee for convenience, is given by 6ss = fi'K,

yss = (fi/K)e-^/K\

(5.47)

The variation of the steady-state with the parameter jx is shown in Fig. 5.7 for a particular choice of K. There is no "criticality" in this case as the steady-state solution varies smoothly with the parameter \x across the whole range, so there are no discontinuous jumps. There is, however, a different type of bifurcation that arises in this model. Provided the second parameter K is sufficiently small, there is a range of the parameter /JL over which the steady-state solution becomes unstable. To understand this statement we can examine how a system sitting at

473

Thermokinetic feedback

0ss

/

s

H1

&

H

Fig. 5.7. Variation of quasi-steady-state temperature excess 0SS with the dimensionless reactant concentration JJL for the Salnikov model. The steady-state is stable at high and low fi but is unstable over a region fxf < \x < /JL* as indicated by the broken section of the locus.

the above steady-state will respond to a (very small) perturbation. Imagine that the perturbation takes y and 6 to the following values: y = y ss + Ay,

(5.48)

6 = 0SS + A0,

where Ay and A0 are small quantities. The equations governing the evolution of the perturbations can be expressed in terms of the original governing equation in which the functional forms of the righthand sides can be expanded about their steady-state values in a Taylor series in the perturbations. Thus if, in general, dy/dr = f(y, 6) and dd/dr = g(y, 6) then we can write

^ =/(7ss, 0SS) + 7^ (7ss, 0ss)Ay + ^(yss, dr dy 30

6ss)A0 + (5.49a)

474

Global behaviour in simple oxidations

Ch. 5

dA0 de dz — = g(7ss, ess) + -*- (y ss , 0ss)Ay + -* (y ss , 0SS)A0 + dr dy ad (5.49b) By the definition of the steady-state condition, the first terms on the righthand side are zero and provided the first-order partial derivative terms do not all vanish, we can ignore the additional terms which are second-order in the perturbations. Thus, we obtain a pair of linear equations for the evolution of these perturbations in the vicinity of the steady-state point V * = ? ( 7 s s , 0ss)Ay + ^ ( y s s , 0SS)A0, dr dy dd

(5.50a)

^

(5.50b)

= — (7ss, 0ss)Ay + ^ ( y s s , 0ss)A0,

dr

dy

36

or, in the most general form, dT d(Ax)

= J(Ax),

(5.51)

where (Ax) is the vector of perturbations and J is, as before, the Jacobian matrix of first partial derivatives which is evaluated with the steady-state values for the variables y and 0. The evolution of the perturbations is then given by the sum of exponential terms in the form Ay = ax eAlT + a2 eA2T,

A0 = bx eAlT + b2 eA2T,

(5.52)

where the coefficients arb2 depend on the initial perturbation. The qualitative nature of the time dependence is determined by the exponents Ala. These are obtained as the eigenvalues of the Jacobian matrix, i.e., from the equation J - A I = 0, where I is the identity matrix.

(5.53)

475

Thermokinetic feedback

For the Salnikov model, the partial derivatives can be evaluated and the steady-state solution substituted to obtain the Jacobian matrix for this twovariable system in the form

dg \dy

dg\ ae/

\ ee

yee-l)

I t^ K

^-1 K

(5.54) The characteristic equation for this matrix has the form of a quadratic equation with the eigenvalues then given by the roots A2 - tr(J)A + det(J) = 0,

(5.55)

where tr(J) is the trace (the sum of the terms on the leading diagonal) and det(J) is the determinant. For the present model, we have, therefore, A2 + (1 + K^/K

- ^) A + KS*K = 0.

(5.56)

Four possibilities exist for the roots depending on the sign and magnitude of the trace and determinant: (i)

tr(J) < 0, tr(J) 2 - 4 det(J) > 0 two negative real roots

(ii)

tr(J) < 0, tr(J) 2 - 4 det(J) < 0 two complex conjugate roots with negative real parts

(hi)

tr(J) > 0, tr(J) 2 - 4 det(J) < 0 two complex conjugate roots with positive real parts

(iv)

tr(J) > 0, tr(J) 2 - 4 det(J) > 0 two positive real roots

In case (i), the two exponential terms in the series for the evolution of the perturbations decay monotonically to zero, so the system decays

476

Global behaviour in simple oxidations

Ch. 5

monotonically back to its original steady-state. In this case the system is stable and termed as stable node. In case (ii) the steady-state is again stable as the perturbations decay but now, through the imaginary parts of the exponential terms the decay has a damped oscillatory character and the steady-state is termed a stable focus. In case (hi) there is again oscillatory behaviour, but now the positive real parts indicate that the perturbations grow. This growth will not continue in an unbounded manner as higherorder terms in the Taylor expansion of the full non-linear equations will become important once the perturbation is not infinitesimally small. Nevertheless, the system departs from its initial steady-state which is termed an unstable focus. In case (iv) there is again divergence from the vicinity of the steady-state, but in a direct manner characteristic of an unstable node. The terms node and focus are most easily understood if instead of plotting the perturbations as a function of time, we plot one variable against the other. This gives rise to a trajectory in the y-d phase plane. The steady-state corresponds to a point on this plane and the trajectory indicates the direction in which the system evolves in the vicinity of this singular point (it is termed singular as the slope of the trajectory which is given by d0/dy = (d0/dT)/(dy/dr) = 0/0 at this point). In the case of a stable steady-state, the local trajectory is directed towards the steady-state point either approaching directly (node) or as an inward spiral (focus): for unstable points the flows are in the opposite direction. The phase portraits associated with the four cases above, and also for a fifth case to be discussed below but not present in the Salnikov model, as shown in Fig. 5.8. Returning to the specific case of the Salnikov model, the major qualitative change in behaviour occurs when damped oscillatory decay of the perturbation gives way to oscillatory growth. The condition for this change from case (ii) to case (hi) which is known as a Hopf bifurcation is, in general terms, Hopf bifurcation

tr(J) = 0.

(5.57)

In this model this becomes l + K e M / K - - = 0.

(5.58)

K

Provided K < e~2 there are two values of ^, which we may denote by /A* and ^ 1 ? at which this occurs, and these mark the ends of a range of steady-

Thermokinetic feedback

477

Fig. 5.8. Phase plane portraits of different possible steady-state singularities: (i) stable node, trajectories approach singular point without overshoot; (ii) stable focus showing damped oscillatory approach; (iii) unstable focus showing divergent oscillatory departure; (iv) unstable node showing direct departure; (v) saddle point x showing insets and outsets and typical trajectory paths.

state instability as indicted by the broken portion of the steady-state loci in Fig. 5.7. If K increases towards the value e~2, the ends of range approach each other and instability is lost from the system. If the steady-state is unstable over the range between the Hopf bifurcation points, what happens to the concentration and temperature excess in the region /i* < \x < ju,*? We can answer this in a straightforward manner by integrating the full equations for a particular value of K < e~2 and with JJL then chosen in the above range and with initial conditions close

478

Global behaviour in simple oxidations

Ch. 5

to, but not exactly equal to, the steady-state values. In this particular case, we will then find that the concentration and temperature settle into a regular, periodic oscillation about their steady-state values as shown in Fig. 5.9(a, b). The amplitude and period of this oscillatory motion vary with the parameter JJL across the region of instability, with the amplitude tending to zero at the two ends as the steady-state regains stability. If we plot the oscillatory behaviour on the y-0 phase plane, it draws out a closed loop or limit cycle around the steady-state singular point, as indicated in Fig. 5.9(c). It is a characteristic feature that a limit cycle is born in the phase plane at a Hopf bifurcation point. In the Salnikov model with the exponential approximation we have a particularly simple scenario. At the lower Hopf point, the steady-state becomes unstable as the parameter /x is increased through /i*. As this occurs, a stable limit cycle is born off the steady-state point at the Hopf point and grows in size around the singular point as /JL is increased further. This is characteristic of a supercritical Hopf bifurcation. What we will see in the reaction is a steady-state give way to oscillations in a smooth manner, a so-called soft excitation of the oscillations that grow in size from zero. If we reduce /JL again, the oscillations shrink back to zero amplitude an the steady-state is regained at exactly the same parameter value without any hysteresis. Such a bifurcation is often represented by Fig. 5.10(a) which shows the amplitude of a stable limit cycle grow from zero as the parameter is varied so that the steadystate becomes unstable. The behaviour at the upper Hopf point is also that of a supercritical Hopf bifurcation although the loss of stability of the steady-state and the smooth growth of the stable limit cycle now occurs as the parameter is reduced. This is sketched in Fig. 5.10(b). We can "join up" the two ends of the limit cycle amplitude curve in the case of this simple Salnikov model to show that the amplitude of the limit cycle varies smoothly across the range of steady-state instability, as indicated in Fig. 5.11(a). The limit cycle born at one Hopf point survives across the whole range and dies at the other. Although this is the simplest possibility, it is not the only one. Under some conditions, even for only very minor elaboration on the Salnikov model [16b], we encounter a subcritical Hopf bifurcation. At such an event, the limit cycle that is born is not stable but is unstable. It still has the form of a closed loop in the phase plane but the trajectories wind away from it, perhaps back in towards the steady-state as indicated in Fig.

Thermokinetic feedback

479

•MM Fig. 5.9. Variation of (a) concentration y and (b) temperature excess 6 in time for a value of fx in region of steady-state instability showing sustained oscillations; (c) typical limit cycle lying around unstable steady-state o produced by plotting y against 0; (d) for some parameter values there are two limit cycles surrounding a stable steady-state, one unstable (broken curve) the other stable (solid curve). In (c) and (d) example trajectories are shown (thin curves) that wind either onto the stable limit cycle or, in (d) onto the stable steady-state point.

5.9(d). Such a limit cycle grows so that it surrounds a stable steady-state, as sketched in Fig. 5.10(c, d) and so would grow out of the upper Hopf point /x* as the parameter /JL increases. This then leaves the question as to what happens if we start outside the unstable limit cycle. The system does not wind onto an unstable cycle, nor can the trajectory cross it so it cannot approach the stable steady-state inside. The system must move further away across the phase plane and either find another steady-state (of which there is not one in the Salnikov model) or perhaps to find another, but

480

Global behaviour in simple oxidations

(a)

sic

sss f _uss_

Ch. 5

(b)

I^ J c Luss_ ji sss

^ (c) - - jilc sss

\ uss

n

(d) ulc.- uss

/

\

sss \

^

Fig. 5.10. The four possible types of Hopf bifurcation: (a) a stable steady-state (sss) becomes unstable (uss) as a parameter /JL is increased through the bifurcation point (/A*) and a stable limit cycle (sic) emerges - the growth of the limit cycle is indicated by plotting the maximum and minimum of the variable as it undergoes the oscillatory motion around the limit cycle; (b) the scenario is reversed, with the steady-state losing stability and a stable limit cycle emerging as the parameter is reduced; (a) and (b) are termed supercritical Hopf bifurcations. In (c) and (d) there is an unstable limit cycle emerging to surround the stable part of the steady-state branch: this is characteristic of a subcritical Hopf bifurcation.

stable, limit cycle. The latter occurs in the Salnikov case and the modified bifurcation diagram is shown in Fig. 5.11(b). The stable limit cycle born at the lower Hopf point "overshoots" the upper Hopf point but is extinguished by colliding with the unstable limit cycle born at the upper Hopf point which also grows in amplitude as /x is increased. Over a, typically narrow range, then there are two limit cycles, one unstable and one stable around the (stable) steady-state point. If we start with the system at some large value of JJL, so we settle onto the steady-state locus, and then decrease the parameter, we will first swap to oscillations at the Hopf point juf- At this point there is a stable limit cycle available as the system departs from the now unstable steady-state, but this stable limit cycle is not born at this point and so already has a relatively large amplitude. We would expect to

Thermokinetic feedback

Hi

m

£

Hi*

H»*

p.

Fig. 5.11. Variation of the oscillatory (limit cycle) solution with /x for the simple Salnikov model showing that the stable limit cycle born at one supercritical Hopf bifurcation exists over the whole range of the unstable steady-state, shrinking to zero amplitude at the other Hopf point; (b) in this case, each Hopf point gives rise to a different limit cycle, with a stable limit cycle born at /if growing as /JL increases and an unstable limit cycle born at /x* also increasing in size as \x increases. At some /JL > /xj the two limit cycles collide and are extinguished.

see the system suddenly jump from steady-state to large amplitude oscillations in a hard excitation process. Furthermore, if we now increase the parameter again, we stay on the stable limit cycle locus until the later bifurcation point where this coalesces with the unstable cycle born at /**• Thus, there is a region of hysteresis between steady-state and oscillatory behaviour characteristic of a subcritical Hopf bifurcation, although in many real systems this is too small to be clearly identifiable in practice. There

482

Global behaviour in simple oxidations

Ch. 5

are other important ways in which limit cycles disappear or lose stability that we will discuss later. Experimental confirmation of the existence of oscillations in a system with Salnikov-type kinetics has been obtained [17] by Griffiths and co-workers. We saw earlier that the simple model for a single exothermic reaction in a CSTR gives rise to multiple steady-state solutions and also that a mathematical recipe for locating the critical turning points in the steadystate locus involved, in general terms, the determinant of the Jacobian matrix becoming equal to zero. Referring back to the eigenvalue equation above, we can see that under such a circumstance one of the eigenvalues A will become equal to zero. If the determinant changes sign from positive to negative, then we will have a situation where one root will be positive and the other negative, irrespective of the sign or magnitude of tr(J). The evolution of any perturbation will then be governed by one exponentially decreasing term but also by one exponentially increasing term in equation (5.50). Eventually, the growth term will dominate and the system will move away from the vicinity of the steady-state even though there may some initial "fast motion" in its direction. The steady-state point is known as a saddle point and its phase portrait was given earlier in Fig. 5.8(e). There is one special path that passes through the steady-state along which the coefficient of the term with the positive eigenvalue are exactly equal to zero. If the system is perturbed exactly onto this pair of saddle insets then the system will actually travel along them to the steady-state. The saddle outsets correspond to the case in which the coefficients of the terms with the negative exponent are zero. In general, a trajectory from an arbitrary point may evolve towards the steady-state parallel to the inset but will then begin to move away, ultimately parallel to the outset. The insets and outsets play important "organizing" roles in the phase plane when there are multiple steady-states or other attractors such as stable limit cycles. Saddle point steady-states are those on the middle branches in regions of hysteresis. We never encounter an isolated saddle point in chemical systems: they always occur with a node or focus partner. At a vertical turning point in a steady-state locus, where the zero eigenvalue arising because det(J) = 0 occurs, the steady-state locus splits into two branches (or two branches merge): one branch of saddles and the other which is initially a node although it may gain focal character an perhaps change stability at a Hopf point later.

Thermokinetic feedback

483

If we consider a typical phase plane for a two-variable model in which three steady-state solutions lie (i.e., we have parameter values so that we are in a region of hysteresis) and if the steady-states corresponding to the uppermost and lowermost branches are stable, then the insets to the saddle lie across the plane so as to divide it in to two distinct regions. The inset is then known as a separatrix. If we examine the trajectories in the phase plane then we will notice that none cross the separatrix: those that start on one side end up approaching one of the stable steady-states and those that start on the other approach the other steady-state. Thus, the separatrix defines the edges of the two basins of attraction of the individual stable steady-state attractors on the phase plane as indicated in Fig. 5.12. In the simplest cases, the separatrix lies fairly simply across the phase plane, but it is not unusual for the insets to weave complex dances across the plane, although they avoid intersecting each other except under very special (but also very interesting) circumstances. If the insets are wound in a complex manner so as to divide the plane into interlocked thin strips, the different initial conditions that lie quite close but in different bands can end up evolving to different steady-states. An important event in the phase plane occurs if the inset to a saddle manages to join up with an outset from the same saddle point. This then gives rise to a closed loop with the saddle point lying on it as a "corner". Such a loop is known as a homoclinic orbit as it forms a path connecting

Fig. 5.12. Typical phase plane arrangement of two co-existing stable-steady-states separated by a saddle point. The insets to the saddle divide the phase plane into two parts: trajectories starting above this separatrix tend to the steady-state point in the upper righthand part of the plane; those starting below tend to the lower steady-state.

484

Global behaviour in simple oxidations

Ch. 5

the steady-state to itself. (There can also arise heteroclinic orbits which are closed paths connecting two different steady-states.) A homoclinic orbit only occurs for a particular value of the parameter being varied. If the parameter is varied a bit further, then the orbit may break up so that it is no longer closed, or it may shed off from the saddle point to form a separate closed loop in the phase plane independent from any steady-state point. The latter, of course, is simply a limit cycle of the type seen earlier, so homoclinic orbit formation is important as another mechanism for limit cycle (and hence oscillation) formation or, in reverse, extinction in systems that exhibit multiple steady-states as sketched in Fig. 5.13. One more model scheme is of interest: the Gray-Yang model for some aspects of the low temperature oxidation of hydrocarbons [18-21]. This involves the features of chemical and thermal feedback described previously with a chain-carrier X coupled to the temperature T. Four reaction steps are required: (GY1)

initiation

A->X

rate = kxa

(GY2)

branching

A + X —> 2X

(GY3)

termination I

X—>P

rate = ktXx

(GY4)

termination II

X->Q

rate = kt2x

rate = kbax

In the original formulation, the reactant concentration a is assumed to be constant, so this is not being used to model full ignition problems. The

(a)

(b)

(c)

1 £r \cr Fig. 5.13. Formation of a stable limit cycle about an unstable steady-state through a homoclinic orbit as the inset and outset of a saddle point merge.

485

Thermokinetic feedback

branching reaction (GY2) and the termination step (GY4) are taken to be exothermic, with the termination step being the more strongly exothermic of the two. The termination step (GY3) is taken to be thermoneutral (in fact, it is only necessary that this be less exothermic than the branching step) as is the initiation step (GY1). Steps (GY1), (GY2) and (GY4) are taken to have rate constants following the Arrhenius temperature dependence, with Et2 > Eb. The activation energy Etl for the termination step (GY3) is taken to be effectively zero, so that the rate constant is independent of temperature, but again the requirement is actually less strict as Etl < Eb. The mass and heat balance equation for this scheme have the form dx — = ki(T)a + cl>(T)x, dt (cPp)^= at

qMT)a

+ © ( 7 > ~^(TV

(5.59a) 7 a ),

(5.59b)

where = kb - kn - kt2 is the net branching factor and @ is a related quantity @ = qjtb + qnKi + qtiK2,

(5.60)

involving the reaction exothermcities. Both of these factors are temperature dependent via the rate coefficients. The inequalities described above for the exothermicities and activation energies mean that at low temperatures there is an overall exothermic reaction involving steps (GY1) and (GY3). As the temperature increases, so the initiation rate and the branching rate increase, leading to an increase in the number of chain carriers and to a greater rate of heat evolution. On further increasing the temperature of the reacting mixture, the second termination step becomes significant and, although this is also exothermic, the consequent decrease in the net branching factor sees a fall in the chain carrier concentration and, hence, an overall reduction in the rate of heat release. At very high temperatures, the rate of heat release increases again as the rate of the initiation step increases. The resulting dependence of

486

Global behaviour in simple oxidations

Ch. 5

gas temperature, T Fig. 5.14. Thermal diagram (compare Fig. 5.4) for the Gray-Yang model showing a maximum and minimum and a region of negative temperature coefficient. Also shown are three heat loss lines: Li intersects R four times, with states I and III being potential stable steady-states; L2 is a critical case with I and II merging at a point of tangency, so the system would have to jump to the ntc region (steady glow or cool-flame); L3 has a different tangency corresponding to the critical condition for ignition.

the total rate of heat release on the reacting mixture temperature, i.e., the thermal diagram for this model, may have the form shown in Fig. 5.14, with a maximum and a minimum superimposed on the general exponential increases seen earlier for the single one-step exothermic reaction. In between the extrema is a region in which the rate of heat evolution decreases with increasing temperature, a phenomenon known as the negative temperature coefficient or ntc. This region is also intimately connected with the limit cycle oscillation behaviour in this model scheme. The basic Gray-Yang scheme has been extended by Wang and Mou [22] who wrote the proper forms of the governing equations appropriate to a CSTR and allowed for the consumption of the reactant A and an additional branching step (GY5)

branching II

A + X-+2X

rate = kb2ax

This creates a model with three independent variables, a, x and T and allows for some more complex responses and for behaviour that might be related to the full ignition of hydrocarbons.

487

The H2 + 0 2 reaction 5.7 THE H 2 + 0 2 REACTION: p-Ta IGNITION LIMITS IN CLOSED VESSELS

The reaction between hydrogen and oxygen has been thoroughly reviewed elsewhere on several occasions [2, 4, 9, 23] and so the account of the classical behaviour in closed vessels will be restricted here to the basic features necessary for setting the background and interpreting the "new" behaviour from studies in open (flow) reactors. Figure 5.15 shows the classic form of the three pressure-temperature (p-Ta) ignition limits for the H 2 + 0 2 reaction in a static, closed reaction vessel. The experimentalist has control over the mixture composition, the total pressure and the ambient temperature (i.e., the temperature of the oven in which the reactor is enclosed) as well as other less obvious factors such as the coating on the inside of the vessel walls. Depending on the values for these various parameters, the reaction at reduced pressure typically has one of two qualitative forms. Either the reaction is very slow, perhaps even undetectable, or it occurs rapidly on a millisecond timescale. In either case, the final product is essentially complete conversion to water vapour, but the route from reactants to products is obviously different to the observer. The locus of pressure and ambient temperature conditions that separate

third limit

°: 400 ignition

GO

00 CL,

5

700 800 ambient temperature, Ta/K Fig. 5.15. Schematic representation of the p-Ta ignition limits for the hydrogen + oxygen reaction in a closed reactor.

488

Global behaviour in simple oxidations

Ch. 5

P2

pressure, p Fig. 5.16. Schematic variation of "reaction rate" as a function of pressure at fixed ambient temperature showing dramatic increase in rate approaching the three limits.

these two types of behaviour for a given mixture composition and vessel preparation forms an approximately Z-shaped curve comprised of the first, second and third "explosion limits" in order of increasing pressure. At any pressure, the system moves from "slow reaction" to "ignition" as the ambient temperature is increased, but the response to changes in pressure is more interesting. Starting at very low pressures, the reaction rate is extremely low, but increases with increasing pressure until the first limit is reached, at which point ignition appears with a high instantaneous rate. Ignition is the response across a range of pressures, but as we cross the second limit, so the "rate" falls again. The first and second limits form a peninsula on the ignition diagram. If the pressure is increased further, the rate increases again and becomes high as the third limit is approached. This description is frequently illustrated with a diagram of the form of Fig. 5.16. Here the term "rate" needs a rather careful interpretation as the system does not really attain a steady-state at any point: rather the maximum rate is generally much lower outside the ignition region than inside where it becomes high (although always remains finite). The reaction in the region "above" the second limit (i.e., at higher pressure in the slow reaction zone) is sufficiently rapid to be measured with conventional methods: the exothermic reaction can support transient temperature excesses (gas temperature being heated above the ambient temperature) of several kelvin [24] as indicated in Fig. 5.17.

489

The H2 + 0 2 reaction

(a) 250 , AT = 3K NAT = 5K

(b) S

\s

s

\

(c) N^AT = 3 K

^

4T = 2 K ^ ^

AL

N A T = 2K

AT =

AT = 3K k s

^

iks

WIKV AT = 0 K V > - ^ ^ignition

780

ignition

860

780

AW = OK

860780

ignition

860

ambient temperature, T./K Fig. 5.17. Self-heating in the slow reaction zone above the second limit in the H 2 + 0 2 system: (a) large vessel, equimolar mixture; (b) small vessel equimolar mixture; (c) small vessel stoichiometric mixture. The numbers indicate the maximum instantaneous temperature excess observed at a given p-Ta location and the "isotherms" connect points at which the same AT is observed. (Reprinted with permission from reference [24], © Royal Society of Chemistry.)

5.7.1 First limit The location of the first limit is sensitive to a number of experimental parameters including the vessel diameter, surface coating, 0 2 concentration and also to packing the vessel with glass rods (another way of varying the surface: volume ratio). These features are explained by invoking the competition between the branching cycle (steps 1-3) and the termination step (4) occurring on the vessel surface. As described earlier, the net branching factor for this system would have the form

$ = kb-

kt = 2k2[02] - k4.

(5.61)

The rate constant k2 shows a typical Arrhenius form, at least over a modest range of temperature, while k4 involves a combination of diffusion to the walls and a subsequent surface-phase reaction.

490

Global behaviour in simple oxidations

Ch. 5

The condition for criticality, i.e., for the first limit, will then be parametrized by cf> = 0, i.e., by the condition 2k2[02] = k4.

(5.62)

In unstirred systems, the effective rate of diffusion will be decreased by increasing the total pressure, allowing inert gases to influence the explosion pressure. The concentration of 0 2 will be directly related to the partial pressure, [O2]=p02/RT.

(5.63)

If we ignore the effect of p on k4, which may be much smaller in wellstirred systems, the ignition condition can be written as Po2,.r = W t c c r = ( ^ ) RT = ( ^ )

e+*'*

(5.64)

where ptot,cr is the limit pressure for a given mole fraction x02 of 0 2 , indicating that the limit pressure decreases as T increases. Different surface coatings affect the value of k4, the surface termination coefficient. For some particularly inert or "reflective" surfaces, k4 ~ 0 and so the first limit disappears to zero pressure. 5.7.2 Second limit The important part of the mechanism applicable to pressures and temperatures along the major part of the second limit involves the competition between the branching cycle (l)-(3) and a gas-phase "termination" step (5)

H + 02 + M ^ H 0 2 + M

This is effectively a termination step if the species H 0 2 does not continue the reaction chain. The latter situation arises if, say, the radical is efficiently removed in a surface reaction of the form (6)

H 0 2 —> products at wall

Provided step (5) is rate determining in this termination step then the termination rate rt = /c5[H][02][M] and the net branching factor will have

The H2 + 0 2 reaction

491

the form cf> = 2k2[02]-k5[02][M],

(5.65)

and the condition = 0 for the second limit will be equivalent to 2k2 = k5[M].

(5.66)

In this expression, the term [M] is a representation of the concentration of "third body" species M. The role of M can be played by any molecule or radical in the vessel and is essentially to stabilize the H 0 2 by energy transfer. In general, only species present in significant concentrations will be important here, but different species have different effectivenesses at facilitating the energy transfer. This is accommodated by assigning different values of third body efficiencies a{ to different species /. The termination rate can then be written as rt = kf2 (xUl + a02Xo2 + 2 a{x{ J -^z [0 2 ],

(5.67)

where kf2 is the rate coefficient for step (5) with H 2 acting as the third body M and the summation is over all species present in addition to H 2 and 0 2 . The relative third-body efficiency ax any species / is the ratio kl5/k™2 and, thus, by definition aU2 = 1. For 0 2 , a0l ~ 0.3, i.e., 0 2 is only 30% as effective as H 2 at stabilizing the H 0 2 species, so as H 2 is replaced by 0 2 in a mixture, the net efficiency of the termination process decreases. The ignition limit criterion $ = 0 can now be written as 2k2 = kf2 (xU2 + a02x02 + E axxx \ - ^ f ,

or 2k2RT kf2 (aU2 + x02x02 + 2 axx{

(5.68)

492

Global behaviour in simple oxidations

Ch. 5

This indicates that the limit pressure will increase as the temperature increases mainly through the Arrhenius temperature dependence of the branching step rate coefficient k2. The termination step relies on collisional energy transfer and its net rate is likely to decreases as the temperature (and, hence, average energy of the reacting species) increases, equivalent to a negative temperature coefficient (activation energy). The influence of the mixture composition on the limit pressure arises from the term in the denominator relating to the relative third body efficiencies. If H 2 is replaced by a less efficient third body, so ptot,cr increases, i.e., ignition arises at a lower temperature for a given fixed pressure. Conversely, if H 2 is replaced by a more efficient third-body, the termination rate is enhanced and the limit pressure decreases or, equivalently, the limit moves to higher temperatures. An important case of this arises if H 2 is replaced by H 2 0 , the product of the reaction. Water vapour is particularly efficient as a third body in step (5), with aU20 typically being regarded as having a value of approximately 6.3 (although a recent direct measurement has suggested that this may need to be reduced by a factor of 3: see Section 3.2.2(iii)). If a mixture of H 2 and 0 2 is maintained at some pressure above the corresponding second limit for a significant period of time, the slow reaction will cause the production of some water. This in turn will cause the limit for the instantaneous mixture to change its location. Baldwin and Walker [25] describe this for a series of experiments in which the location of the limit is determined by premixing the reactants at a pressure well above the limit and then withdrawing gas through a capillary. They noted that the observed limit pressure depends on the rate of withdrawal of gas and that, with a sufficiently low withdrawal rate the ignition, could be completely suppressed by the formation of the "inhibitor" product. The inhibiting role of the product on the reaction is also of great importance in interpreting flow-reactor phenomena and we return to this in a later section. 5.7.3 Reactions involving H02 and the third limit This simple interpretation of the second-limit mechanism is appropriate provided the reaction vessel surface is "efficient" with respect to removal of the species H 0 2 . This situation arises with many salt-coated vessels, with KC1 being widely exploited. Provided this is arranged, the rate determining part of the termination process is the gas-phase step (5) and the

The H2 + 0 2 reaction TABLE 5.1 The Baldwin-Walker mechanism H2 + 0 2 H2 + 0 2 H 2 + OH H + 02 0 + H2 H + 02 + M H O OH H02 + H02 H02 + H H02 + H H02 + H H 0 2 + H2 H202 + M H202 + H H202 + H H 2 0 2 + OH H202 + O H20 + 0 H20 + H OH + O OH + H OH + OH H + OH + M H+H+M OH + OH + M O+O+M H02 + M H2 + M 02 + M

-* -* -> -» -» —> -> -> -> -* -» —» -» -> -> -» —> —> -* -> -» —» —> -* -> -> -> -» -> -> ->

20H H02 + H H20 + H OH + O OH + H H02 + M wall wall wall H202 + 0 2 20H H20 + 0 H2 + 0 2 H202 + H 20H + M H2 + H 0 2 H 2 0 + OH H20 + H02 OH + H 0 2 20H H 2 + OH 02 + H H2 + 0 H20 + 0 H20 + M H2 + M H202 + M 02 + M H + OH + M 2H + M 20+ M

limit is relatively insensitive to vessel size. With other coatings, and in particular with "aged boric acid" which has been extensively exploited by the Baldwin-Walker group [26-28] for the determination of reaction rate constants, the removal of H 0 2 at the walls is not significant. In such systems, the H 0 2 concentration will increase with time and further reactions of this species may become important. Of these, the reaction (7)

H 0 2 + H 0 2 -* H 2 0 2 + 0 2

494

Global behaviour in simple oxidations

Ch. 5

sees the formation of hydrogen peroxide. This reaction will also become important even for vessels with efficient surfaces at sufficiently high pressure. Hydrogen peroxide can dissociate thermally (8)

H202 + M->20H + M

returning active chain carriers to the radical pool. This decreases the effectiveness of the termination process and, eventually, will provide for a less negative net branching factor. The latter effect, along with the increasing rate of heat generation, brings the possibility of a combined chain-thermal ignition which is believed to be the character of the third limit. Additional reactions that need to be considered when modelling the H 2 + 0 2 reaction in the vicinity of the second limit are given in the Baldwin-Walker scheme in Table 5.1.

5.8 FLOW REACTOR STUDIES OF THE H2 + 0 2 REACTION

The H 2 + 0 2 reaction has been studied in the vicinity of second explosion limit pressures in well-stirred, continuous-flow reactors over the last 15 years [29-38]. The typical arrangement of the apparatus is shown in Fig. 5.18. Experiments are commonly performed by setting the inflow rates with electronic flow controllers to achieve the desired mixture composition and residence time at a specified total pressure. The inflow channels of fuel and oxidant are preheated separately before reactor entry which is achieved through jet nozzles to enhance mixing. In some experiments, additional mechanical mixing is provided. During an experiment it is generally most convenient to maintain a constant total pressure (controlled by a needle valve in the outflow line) and to vary the experimental conditions through the oven temperature Ta as required. In between small step changes in Ta, the system will typically be allowed to adjust to a steadystate appropriate to the operating conditions, i.e., the system is left for sufficient time for transient features to die out. The p-Ta ignition diagram obtained in the above manner shows an explosion limit, Fig. 5.19, similar in both form and location to that observed for identical mixture compositions in closed vessels and also to similar slow reaction behaviour for pressures above the limit. Use of a flow reactor, however, also allows the behaviour on the "ignition" side of the

495

Flow reactor studies of the The H2 + 0 2 reaction Oven

Preheating Coils

Reaction Vessel

V

Liquid Nitrogen Trap and Vacuum Pump Chart Recorder

Mass Flow Controllers

\ Thermocouple

Photomultiplier

Jet Nozzle

Personal Computer

Fig. 5.18. Diagrammatic representation of flow-reactor apparatus for studying combustion reactions. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

50 r13 o

40

slow reaction

steady ignited state

3.0

20

10 0 650

700

750

800

ambient temperature Ta/K Fig. 5.19. Thep-T a ignition diagram for an equimolar H 2 + 0 2 mixture with mean residence time tres = 5.2 ± 0.7 s showing region of slow reaction separated by second limit from regions of oscillatory ignition and steady-ignited state. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

496

Global behaviour in simple oxidations

Ch. 5

limit to be investigated experimentally. In a closed system, the ignition encountered on crossing the limit sees the rapid completion of all reaction with no further kinetic developments. In a flow system, the products of the ignition will subsequently be replaced by an inflow of fresh reactants. Thus, we can imagine at least one of two scenarios: either the inflowing reactants will support a continuation of the ignition process giving rise to a "steady flame" in which there is a steady combustion of the fuel + oxidant stream at a rate matching its inflow; or, we might expect a periodic sequence of ignition event separated by periods of relatively little chemistry but during which the mixture composition is changing under the influence of the inflow and outflow. For an equimolar mixture ( H 2 : 0 2 = 1:1) and if the total pressure in the reactor exceeds approximately 30Torr, then crossing the ignition limit by increasing the ambient temperature causes the reaction to jump to a steadystate corresponding to a steady, high consumption of the incoming reactants. That is, the steady-state concentration of H 2 in the reactor and outlet is approximately zero, and the 0 2 concentration is one-half of its inflow concentration. The reaction rate is not particularly high (compared with those associated with a normal burning or flame state) as it is limited by the inflow rates which are relatively low in such systems (fres is of the order of seconds and is long compared to the Newtonian cooling timescale). If the ambient temperature is subsequently decreased, the system may stay in this high reaction state to temperatures well below the ignition limit temperature located on the upward sweep, with an extinction event at some lower Ta. This corresponds to the type of steady-state hysteresis loop described in Section 5.3 and we can sketch the bifurcation diagram for this system as in Fig. 5.20(a) with a branch of unstable steady-states (saddle points) in between. Pursuing the dynamical systems theory description of this response, we see that the ignition limit can now be identified with the turning point or "saddle-node bifurcation point" in the steadystate locus where the branch of low reaction rate (high [H2]ss) states meets the middle, saddle point branch. Similarly, the extinction point, which has no analogy in the closed system, is also a turning point in the steady-state locus. If the reaction is perturbed while the system is in the ignited state, e.g., by a momentary disturbance to the inflow, there is typically a damped oscillatory return indicating stable focal character. If the system is in the region of hysteresis and the perturbation is sufficiently large, it can cause a transition to the low reaction rate (steady slow reaction) state. The slow

Flow reactor studies of the The H2 + 0 2 reaction

T x

497

T a,cr

x

a,Hopf

Fig. 5.20. Schematic representations of bifurcation diagram appropriate to H 2 + 0 2 system in flow reactor, (a) A low conversion (high [H2]) steady-state exists at low ambient temperature but terminates in an ignition turning point at r a , cr at which the system must jump to the lower branch of high conversion steady-state. If Ta is subsequently reduced, the system can stay on the lower branch for Ta < Ta,CT until the extinction turning point. There is experimental evidence for hysteresis at the limit between steady slow reaction and steady ignited states, (b) Variation in experimental conditions allows high reaction (low [H2]ss) branch to lose stability over range of ambient temperature (for Ta < Ta,Hop{). System evolves along high [H2]ss (low reaction) branch until Ta = Ta,CT at which point it evolves to large amplitude limit cycle (sic). As Ta is increased, amplitude of oscillations decreases, and a steady-state emerges for Ta > Ta,Hop{. If the ambient temperature is reduced again, there is no hysteresis, with oscillatory ignition re-appearing at r a , Hop f and the oscillations terminating at (the saddle-node) bifurcation at r a , cr . (c) As in (b) except that the limit cycle locus is now folded: two separate branches of stable limit cycles (sic) corresponding to different oscillatory solution now exist, separated by a branch of unstable limit cycles (ulc). Hysteresis is observed between the large- and small-amplitude oscillations over a (typically narrow) range of ambient temperature.

498

Global behaviour in simple oxidations

Ch. 5

reaction branch typically has stable nodal character with perturbations decaying monotonically. At a lower total pressure, the reaction again exhibits a low reaction steady-state at low Ta, equivalent to the slow reaction in a closed vessel. This lies on a branch of steady-states that terminates with a saddle-node bifurcation point corresponding to the ignition limit, as indicated in Fig. 5.20(b). In this case, however, the steady-state corresponding to the ignited state is not stable and so the system does not jump to a steady-flame state. Instead it oscillates around the ignited state, showing an oscillatory sequence of ignitions. The typical waveform for a set of experimental conditions lying just beyond the ignition limit is shown in Fig. 5.21, in which the H 2 , H 2 0 and a nominal record of the gas temperature (which is different from the ambient or oven temperature) and the emitted visible light intensity are plotted as a function of time. The temperature is recorded in such experiments using fine-wire thermocouple junctions (Pt/Pt + 13%Rh, coated with silica to avoid catalysis), but these have significant response times compared with the short timescale of the sharp ignition phase in these "relaxation-type" oscillations. (The response-time distortion of the signal is even more pronounced if the voltage is output directly to a traditional chart recorder, but this is reduced with an oscilloscope or A/D converters coupled to computer-based data acquisition systems.) A more realistic record of the temperature excursion is obtained indirectly from the rotational level populations from laser absorption studies on the OH radical concentration, Fig. 5.22, indicating that temperature rises approaching the adiabatic temperature excess are achieved. The H atom concentration has been followed using resonance enhanced multiphoton ionization (REMPI) [34-38]. For conditions close to the limit, there is virtually complete H 2 consumption during the ignition stage of the oscillation, followed by a flow-based recovery during which [H2] approaches its inflow concentration. The closer the conditions are held to the limit condition, the closer [H2] approaches [H 2 ] 0 and also the longer the period between ignition events. In a careful measurement of the variation of the oscillatory period with the ambient temperature in this region, Griffiths et al. [32] showed that the period varies approximately according to the relationship

499

Flow reactor studies of the The H2 + 0 2 reaction Light Output ( arbitary units ) AT\ K 200-1

100

100 Time \ seconds

200

100 Time \ seconds

200

[H:0]

50 100 Time \ seconds

10,]

mm 50 100 Time \ seconds

Fig. 5.21. Typical instrumental records for oscillatory ignition; (a) AT as measured by finewire thermocouple; (b) light emission intensity as measured by photomultiplier; (c) [H 2 0] and (d) [0 2 ] as measured by mass spectrometer. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

500

Global behaviour in simple oxidations

Ch. 5

8.0-

-2500

6.0-

h2000

4.0-

-1500 ^

2.0-

-1000

-500

t/ms

Fig. 5.22. Records of [OH], [H] and T from rovibrational spectrum and REMPI studies. (Reprinted with permission from reference [37], © Royal Society of London.)

where Ta,CT is the value of the ambient temperature exactly on the ignition limit, implying an infinite period as Ta —> Ta^CT. Experimentally, a period of 158 s was recorded for a system with tres = 2 s by holding the ambient temperature 0.05 K above r a?cr . The lengthening of the period close to the limit is illustrated in Fig. 5.23. If the ambient temperature is raised so that the system is taken further into the ignition region, the amplitude of the oscillations decreases and their frequency increases. For Ta sufficiently high, the amplitude decreases to zero and a steady-state (steady flame) emerges. The period of the oscillations remains finite as the ambient temperature approaches this oscillatory extinction. Figure 5.24 shows how the waveform changes with Ta close to this upper limit to the oscillatory range. The waveforms become noticeably more sinusoidal as evidenced by the fourier transforms with the waveform for Ta = 800.5 K having only a single frequency component. The limit cycles corresponding to these oscillations can be obtained from a single experimental record (most conveniently, the thermocouple signal) using the "method of delays". In this, the value of the signal at some time t is plotted against the value of the same signal at a later time t + td where

501

Flow reactor studies of the The H2 + 0 2 reaction

T./K

727.80

727.00

724.00

723.70

723.50

158±8

1

T a (t0.05)/K

I

I

I

723.20

I

L

extinction

$723.15

Fig. 5.23. Variation of oscillatory period in vicinity of ignition limit (saddle-node bifurcation point) showing extreme lengthening as limit is approached. (Reprinted with permission from reference [32], © Manchester University Press.)

td is known as the "delay time". This plot is performed for all the points in a time series and gives rise to the limit cycles shown in Fig. 5.25 for the present system. The choice of td is not particularly important but is typically taken to be approximately one-half of the oscillatory period. (The additional "chaining" on the trajectories in this figure corresponds to a frequency of 50 Hz and is due to electrical interference picked up through the thermocouple leads.) As Ta increases, so the limit cycle shrinks to a point (corresponding to the stable steady-flame state. If the ambient temperature is reduced again, oscillations return at the same temperature at which they were lost, indicating a lack of hysteresis. These features are, therefore, characteristic of a supercritical Hopf bifurcation (see Section 5.5.d). We can now sketch the corresponding bifurcation diagram for this pressure and mixture composition, Fig. 5.20(b). This again shows a folded steady-state curve but with the additional feature of the Hopf point along the high reaction rate (low [H2]ss) branch at some temperature r a?Hopf . The high reaction state is unstable for Ta < ra?HoPf a n d stable for higher ambient temperatures. For Ta < Ta^op{, the high reaction state is surrounded by a stable limit cycle whose amplitude is indicated on the bifurcation diagram by plotting the maximum and minimum values of [H2]

502

Global behaviour in simple oxidations

Ch. 5

784 OK

AT \ K 100 AT\ K

T = 799.8K

60H I/VVVVWSA/VN/VVVVNAA/VVVNA/VVVWNA'VVS/WVVW^

Ta = 797.1 K

30

60 10

0

30

AT \ K 0 AT \ K

TQ = 800.5K

™ 60798.5K

30 J

60 A 10

30 J Al

10 60 "

0 AT \ K

\ = 800.6K

T = 799.3K

30 -

60 J 0 -

i

'

•' r

10

30 JWWVV\AAAAAAAAAAAA^^ Time \ seconds 0 10 Time \ seconds

Fig. 5.24. Variation of oscillatory waveform in vicinity of boundary between oscillatory and steady ignition showing characteristic nature of a supercritical Hopf bifurcation. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

attained during the oscillation. This limit cycle increases as Ta is decreased away from r a?Hopf but undergoes an extinction at lower ambient temperatures. Experimentally, it is observed that this extinction typically occurs at approximately the same ambient temperature as that at which the saddle-node bifurcation (ignition limit point) is located as Ta is increased

503

Flow reactor studies of the The H2 + 0 2 reaction

800.5K 3

799.3K

o O

o

o

T3

798.5K

797.1K

digitised thermocouple output Fig. 5.25. Reconstructed limit cycles for oscillations shown in Fig. 5.24. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

from the low reaction steady-state. This, along with the observed scaling law for the oscillatory period given above, suggests that the limit cycle is lost by collision with the turning point that marks the ignition point on the diagram, as illustrated in Fig. 5.20(b). Technically, this is known as the formation of a homoclinic orbit to a saddle-node point (earlier we saw the extinction of a limit cycle by formation of a saddle point for which the inset and outset connected into a loop). The inverse square-root scaling is characteristic of such as structure. This would account for the lack of hysteresis between the slow reaction state and the oscillatory ignition at the limit. Thus, the ambient temperature range between r a c r and ra?Hopf bounds the range of oscillatory ignition in this instance. Both of these points move as we vary the total pressure between experiments: the variation of TacT corresponds to the ignition limit itself while the variation of r a H o p f describes the second boundary on Fig. 5.19 giving the upper temperature limit to the region of oscillatory ignition. The latter curve appears to approach the ignition limit at a definite angle at the upper end of the pressure range. Arguing from the shape of the bifurcation diagrams

504

Global behaviour in simple oxidations

Ch. 5

sketched previously, and noting that there must be smooth transitions from one type to another, leads us to conclude that as ptot increases, the homoclinic orbit formation must move away from the saddle-node (ignition) point and up along the saddle point branch as T^Uopf approached r a c r . In such cases, the oscillatory ignition response will survive to ambient temperatures below that corresponding to the ignition limit, as indicated in Fig. 5.19. We should then see hysteresis between the low reaction state and the oscillatory ignition. Furthermore, at some pressure, the Hopf point may move along the high reaction branch so that it lies at a lower ambient temperature than the ignition limit. The system should then jump to a high reaction steady-state for Ta > r a>cr on the initial upward sweep in ambient temperature, but on the way down will show hysteresis in which the system loses stability of the steady flame to oscillatory ignition. In practice the region of the experimental conditions under which such relatively complex behaviours are exhibited are very small and so these are unlikely to be observed by accident. However, armed with the above argument it has been possible to search successfully for such behaviour in a narrow range of pressure above the crossing point of the two loci on the bifurcation diagram [33,39].

5.9 COMPLEXITY IN THE OSCILLATORY IGNITION REGION

Returning to the variation of the limit cycle amplitude across the region of oscillatory ignition, for the equimolar mixture this is relatively simple, with the amplitude decreasing continuously. Typically, there is a narrow range of Ta over which the amplitude decreases very rapidly as the waveform changes from the relaxation type observed close to the ignition limit, to the sinusoidal character observed near to r a?Hopf . Such a region of rapid change is known rather misleadingly as a "canard" (i.e., false) bifurcation (as there is no qualitative change). Experimentally, it is often difficult to distinguish between a very rapid change and a genuine discontinuity, so care is needed in these instances. With different mixture compositions, however, the range of behaviour supported becomes richer. The variation of the oscillatory period with the ambient temperature for a stoichiometric mixture (2H2 + 0 2 ) at an operating pressure of 16 Torr and a residence time of 4 s is compared with

Complexity in the oscillatory ignition region

505

8

T3

a. 4

0

725

750 775" TJK Fig. 5.26. Variation of oscillatory period with ambient temperature for different H 2 + 0 2 systems: (a) p = 14 Torr, tres = 4 s; (b) p = 14Torr, tres = 2 s, showing region of birhythmicity. (Reprinted with permission from reference [40], © Combustion Institute.)

that of the equimolar mixture in Fig. 5.26. Whereas the period simply decreases across the range for the latter, there is a region of oscillatory hysteresis for the stoichiometric case. Over a range of approximately 730 K < Ta < 745 K, the reaction exhibits either relatively large amplitude ignition events or small amplitude, more sinusoidal oscillations, depending on whether the ambient temperature is being increased from low values or decreased from the steady-flame state at high ambient temperature [40]. The coexistence of two different oscillatory states for the same experimental conditions is analogous to the coexistence of steady-states (multistability) seen earlier, and is known as birhythmicity. We can imagine this arising from a folding of the limit cycle locus in the bifurcation diagram, as indicated in Fig. 5.20(c), just as steady-state multi-stability arises from a folding of the steady-state locus. This also indicates that there will be a third oscillatory state, corresponding to a branch of unstable limit cycle solutions between the two different stable oscillatory states. Thus, there are three limit cycles arranged concentrically around the (unstable) steadystate point in the phase plane, as indicated in Fig. 5.27. The fold in the

506

Global behaviour in simple oxidations

Ch. 5

Fig. 5.27. Schematic representation of three co-existing limit cycles surrounding an unstable steady-state. The middle limit cycle is unstable: trajectories that start outside this evolve to the outer stable cycle; those that start inside the unstable limit cycle evolve to the smallest, stable cycle.

limit cycle locus must arise as the mixture composition, which plays the role of an unfolding parameter in this system, is varied from equimolar to stoichiometric. The unfolding of the birhythmicity feature can also be accomplished by varying other parameters, such as the residence time (with birhythmicity favoured by low fres) or the total pressure. Another type of complexity is a feature for stoichiometric mixtures at certain pressures and residence times [33]. Thep-T a ignition limit diagram for such a system with tTes = 2.0(±0.2)sis shown in Fig. 5.28. Within the region of oscillatory ignition lies a subregion denoted "complex oscillations". The complex oscillations in this system have a mixed mode waveform. These comprise of a single large ignition followed by a number of small amplitude oscillations before the process repeats. A selection of observed mixed-mode states are shown in Fig. 5.29. Also displayed are the associated limit cycles obtained by the method of delays described previously. These show an apparent crossing corresponding to the small amplitude events. In fact, as stated earlier, trajectories in the phase plane cannot cross, except at the singular points of the system. In this case, then, the trajectory must be embedded in a phase space of more than two dimensions, so the 2-D portraits shown are just projections of a trajectory winding around in at least 3 dimensions. The crossings are, therefore, nothing more than one part of the cycle passing above or below another

Complexity in the oscillatory ignition region

507

50 40

slow reaction steady ignited state

20 10

650

region of complex oscillations

700 750 800 ambient temperature, Ta/K

Fig. 5.28. The p-Ta ignition diagram for a stoichiometric 2H2 + 0 2 mixture with mean residence time tTes = 2.0 ± 0.2 s showing additional region of complex oscillatory ignition. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

part but not actually intersecting. (We should remember that there is really one independent variable representing each chemical species involved in the reaction mechanism + one more for the reacting mixture temperature, so the H 2 + 0 2 system is genuinely multidimensional from the point-ofview of dynamical systems analysis. It is striking how for much of its behaviour this, and most other chemical reactions, confine themselves to behaviour that lies in low-dimensional phase spaces, e.g., a 0-D steadystate or a 1-D limit cycle embedded in a 2- or 3-D phase space.) As the ambient temperature is increased across the region of mixedmode states, so the number of small amplitude oscillations in the repeating unit increases. A particular waveform can be denoted by the symbol 1", where n is the number of small events, as n increases with Ta. As we approach the righthand limit (high Ta) so that n effectively tends to infinity, i.e., the large amplitude ignition disappears and we have simply repetitive, small amplitude oscillations which undergo a Hopf bifurcation to yield the stable steady-state at some higher ambient temperature. Mixed-mode oscillations are also favoured by shorter residence times and by mixtures

508

Global behaviour in simple oxidations

T

.•

200-

Ch. 5

734.5K

100-

l

00

jjLiL

I 10

20

200

-600

-1400

-2200

DIGITALISED THERMOCOUPLE OUTPUT Xlt)

30

Time \ seconds

AT \ K 1r a = 7 4 4.1K

200-

K

0.

O

1 j\

0-

u1_J 10

o o o

) THERM

100-

ID

>

1|

-MOO

X

-600

200 •

,

r

J>^~\

///

\

\\1 \

\

(

,

i

1

|

DIGITALIZED THERMOCOUPLE OUTPUT Xlll

20

Time \ seconds

E OUTPUT )

AT \ K

-1400 •

3

o oo

1 '\\

cc X

-600

, / _ _

o

>-

200 J 200

,

-600

, ,

-1400

-2200

DIGITALIZED THERMOCOUPLE OUTPUT X|t) 10

20

30

Time \ seconds

Fig. 5.29. Typical mixed-mode oscillatory ignition waveforms for stoichiometric 2H2 + 0 2 mixtures with corresponding reconstructed limit cycles. (Reprinted with permission from reference [33], © Royal Society of Chemistry.)

Mechanistic modelling of complexity in H2 + 0 2 reaction

509

of composition close to stoichiometric. Some experimental evidence now exists for so-called concatenated states. These are waveforms apparently comprising two different mixed-mode "parents", e.g., a mixing of l 1 and l 2 states to form a repeating unit of the form 1*12 lying in a region of the parameter plane between the l 1 and l 2 states. Between the 1*12 and l 2 , it is also possible to find 1*1212 oscillations in which the repeating unit involves eight events, three of which are of large amplitude and five small ignitions. Some very complicated bifurcation sequences involving hysteresis (birhythmicity) between small amplitude oscillations and mixedmode state have also been observed but not characterized in any quantitative sense to date.

5.10 MECHANISTIC MODELLING OF COMPLEXITY IN H2 + 0 2 REACTION

The first feature modelled for the flow system [29, 30] was thep-T a ignition limit (corresponding to the locus of saddle-node points in terms of the bifurcation diagram). This would correspond to an application of singularity theory as outlined in Section 5.3. In fact, these studies were predated by a singularity theory approach to computing the ignition limits for H 2 + 0 2 mixtures directly by such techniques for closed systems with reactant consumption ignored [41]. Formally, the approach would involve setting up the mass balance equations for the chemistry and flow rates for each species involved in any postulated mechanism and then solving, simultaneously, the steady-state condition and the vanishing of the determinant of the corresponding jacobian matrix. Such computations are now feasible but for the present system represent taking a large sledgehammer to a fairly innocuous nut. We can proceed with sufficient accuracy by building on the earlier analyses involving the net branching factor and how it becomes modified appropriate to flow systems. The ignition limit lies in the region of the p-Ta plane, corresponding to the second limit in classical closed vessels, and so we may surmise that the dominant features of the mechanism will be the competition between the branching cycle (1-3) and the gas-phase termination step producing H 0 2 , step (5). A full steady-state analysis on the intermediates OH and O would introduce (out)flow terms for each species and a fairly complex polynomial in terms of tTes. The full analysis appears in Chapter 4 of this volume. For now, we can note that the typical residence times of interest, 1 to 10 s,

510

Global behaviour in simple oxidations

Ch. 5

mean that the flow terms are relatively insignificant compared with the fast kinetic terms for these reactive species and, hence, can be neglected. Thus, we can use a similar approach to that given earlier in terms of a net branching factor governing the H atom concentration, 4>=(2k2-k5[M])[02l

(5.71)

as before. To allow for the loss of H by the outflow, we can modify this to a form appropriate to a flow reactor as in Section 5.2, to yield fciow = (2fc2 " * 5 [M])[0 2 ] -

fcfiow,

(5.72)

where kf low = l/fres- Implicit in this is the assumption that the steady-state concentrations of the major reactant species do not differ significantly from their inflow concentrations along the branch of low reaction steady-states. This is unlikely to fail by more than 5%. The condition for the limit then reduces to 4>{iow = 0. In fact, for the experimental studies described in the previous sections, the influence of the flow term on the (p-Ta) location of the limit is almost insignificant. The individual values of 2k2 and, hence, A:5[M] are many orders of magnitude greater than the flow term at the typical pressure and ambient temperatures of interest, for instance 2£;2[02] ~ 103 s _ 1 for an equimolar mixture at 16 Torr and 750 K compared with /cflow ~ 0.1 s" 1 . Furthermore, the relatively high activation energy for step (2) ensures that only very small changes in Ta will be needed to allow for the -/cfiow term in equation (5.72). Thus, the problem of predicting the p-Ta limit, and the influence of mixture composition etc., is effectively the same as for equivalent closed vessels. In order to model the behaviour on the "ignition side" of the limit, however, it is clearly necessary to be more sophisticated, at least during the transient ignition stages, as there will be high radical concentrations, temperature excursions and significant temperature rises. There has been an interest in discovering the "minimal" mechanism that will predict at least some form of oscillatory event. A scheme involving the branching cycle (1-3), the termination step (5) and one extra step (9)

H + H02^H20 + 0

gives oscillatory solution for some parameter values provided certain cm-

Mechanistic modelling of complexity in H2 + 0 2 reaction

511

cial features pertaining to reactant consumption and, in particular, product formation were also included [42-44]. For some ambient temperature Ta slightly in excess of the limit condition r a c r , the basic nature of the oscillatory ignition process involves a branched chain ignition that leads to consumption of H 2 and 0 2 and formation of H 2 0 . Remembering that the product species has a significantly higher third-body efficiency for the termination step (5) than the fuel and oxidant mixture it was formed from, the ignition serves to increase the effective value of the term k5[M] in the net branching factor. (During the actual ignition event, the temperature Twill increase dramatically, increasing the instantaneous value of /c2, but this is short-lived as the temperature falls rapidly back to the vicinity of Ta.) Prior to the ignition, fiow will have been marginally positive (as we have just crossed the limit), but the change in mixture composition will now arrange that 0 flow (based on the instantaneous, product-dominated mixture composition) will now be negative. This causes the reaction rate to fall to a low value which we can treat as zero to a first approximation. Thus, the system thus enters a quiescent phase in which changes in composition within the reactor arise solely from the flow processes. The product concentration decreases while those of the original (inflowing) reactants recover. If these are truly flow-controlled processes with no chemical contribution, then the concentrations will vary exponentially with the time elapsed since the ignition event, i.e., they will have the form [H2] = [H 2 ] 0 (l - e - ' H , [0 2 ] = ([O 2 ] 0 - [0 2 ] resid )(l - e - " - ) + [0 2 ] r e s i d ,

(5.73)

[H2O] = [ H 2 ] 0 e - ^ s ,

(5.74)

and

for a fuel-lean mixture, where [0 2 ] r e s i d is the "residual" concentration of 0 2 after the ignition event ([02]resid = 2[O2]0 for an equimolar mixture) and the concentration of H 2 0 immediately after the ignition event is equal to the concentration of the fuel before the ignition. Thus, the concentration of the inhibitor will fall to e _ 1 of its post-ignition value and the reactant

512

Global behaviour in simple oxidations

Ch. 5

concentrations recover to the same factor of their initial concentration within one residence time and e _ n within ^-residence times. During this period, the effective (mixture-composition dependent) value of the net branching factor will be increasing as H 2 0 is replaced by the less efficient reactant mixture. Eventually, flow will pass through zero again and another ignition will develop. From the above argument we can see why the period between successive ignition events increases the closer we are to the limit, r a?cr . The smaller the "degree of supercriticality", the smaller c/>flow based on the initial reactant concentrations will be (it is exactly equal to zero at the limit and slightly positive for Ta marginally in excess of T^cr). This means that the mixture will need to recover virtually to its original composition through the exponential flow processes above, requiring more residence times to elapse before the next ignition, the observation of Griffiths et al. [32] of a period equivalent to 60 x tres indicates such a marginally supercritical system that only once the H 2 0 concentration has decreased to e~60 of its post-ignition value can another ignition develop. As the difference between Ta and r a?cr increases, so the branching term involving the rate coefficient k2 increases and 4>f\ow based on the inflow concentrations is more positive. This means that following an ignition (/>fiow will become zero before the initial composition is regained, so ignition develops "earlier". The reactant concentrations do not recover to their inflow values at any stage, so the amplitude of the ignition will be slightly smaller. The above argument can be used to predict the oscillatory period reasonably accurately, provided Ta is neither too far above T^cr (such that there is significant chemistry occurring in the period before the actual ignition event) nor too close to r a?cr . The latter case must involve additional features as the exponential variation of the mixture composition with time would lead to a logarithmic lengthening of the period as Ta approaches ^a,cr? whereas, the experimentally determined form involves the inverse square root of the difference as given earlier. The latter scaling arises because once the flow has returned the reactant concentrations to their original values, the reaction almost manages to achieve a steady-state and so "hovers" in the vicinity of a state in which mass and heat balance equations almost vanish. The same scaling is observed in the vicinity of any saddle-node point, even in the simple theory of thermal runaway (see Section 5.3). The loss of oscillatory ignition can also be rationalized on the above

Mechanistic modelling of complexity in H2 + 0 2 reaction

513

arguments. If the ambient temperature is increased sufficiently, then the branching term 2k2 will become so large that $fiow remains positive even for the system in which the H 2 0 concentration is at its highest value, i.e., [H 2 0] = [H2]o just after the ignition event. In this case, the "ignition" will continue rather than be quenched and the inflowing reactant streams will simply support a steady burning state. In order to model the oscillatory waveform and to predict the p-Ta locus for the (Hopf) bifurcation from oscillatory ignition to steady flame accurately, it is in fact necessary to include more reaction steps. Johnson et al. [45] examined the 35 reaction Baldwin-Walker scheme and obtained a number of reduced mechanisms from this in order to identify a minimal model capable of semi-quantitative p-Ta limit prediction and also of producing the complex, mixed-mode waveforms observed experimentally. The minimal scheme depends on the rate coefficient data used, with an updated set beyond that used by Chinnick et al. allowing reduction to a 10-step scheme. It is of particular interest, however, that not even the 35 reaction mechanism can predict complex oscillations unless the non-isothermal character of the reaction is included explicitly. (In computer integrations it is easy to examine the "isothermal" system by setting the reaction enthalpies equal to zero: this allows us, in effect, to examine the behaviour supported by the chemical feedback processes in this system in isolation TABLE 5.2 The "minimal complex oscillator" model for the H 2 + 0 2 reaction Rate constant, k = Ae~ E/RT

H2 + 0 2 - > 2 0 H H2 + O H - + H 2 0 + H H + 02^OH + 0 0 + H2^OH + H H + 02 + M->H02 + M H-^wall OH -> wall H02 + H ^ 2 0 H H02 + H-*H2 + 02 OH + 0 - > 0 2 + H

A*

(E/R)/K

g/kJmol - 1

1 x 108 2.2 x 107 2.2 x 108 1.37 x 107 5.1 x 103 12 12 2.5 x 108 2.5 x 107 1.3 x 107

35194 2590 8450 4480 -500 0 0 950 350 0

-77.748 64.0 -70.66 -8.247 199.451 0 0 159.707 238.651 8.247

*Units for A are (m3 mol *)" 1 s step in the mechanism.

1

where n is the order of each elementary

514

Global behaviour in simple oxidations

Ch. 5

from the thermal feedback routes. The isothermal system supports ignition and simple oscillations but not complex oscillations.) Chapter 4 of this volume discusses the formal aspects of such mechanism reduction further. The minimal complex oscillator mechanism of Johnson et al. is given in Table 5.2.

5.11 THE CO + 0 2 REACTION

The reaction between carbon monoxide and oxygen is the other "stoichiometrically simple" combustion process and has been studied for about as long as the H 2 + 0 2 reaction [46,47]. Unlike the latter, however, the progress towards a full mechanistic understanding - or even agreement amongst different workers as exactly how to interpret each other's observations in a consistent manner - has only emerged in the last two decades. The main obstacle to studying this reaction is its extreme sensitivity to traces of any species containing H atoms [48-50]: H 2 0 and CH 4 being the primary source of such contaminants. This was recognized relatively early and there were extensive programmes of research dedicated to "drying" the reactants by storage in liquid 0 2 etc. for long periods prior to mixing the gases [51-54]. In these studies the p-Ta location of the "ignition limit" shifted to higher ambient temperature as the reactants were dried further. The reaction also shows phenomena not exhibited by the H 2 + 0 2 system, including trains of oscillatory light emission that may continue over several hours even in a closed system. The reaction is sensitive to the state of preparation of the vessel surface to a degree beyond that found for the H 2 + 0 2 reaction. Relatively detailed review of the classical ignition limit phenomena can be found elsewhere [4, 5, 23], and so we concentrate here on features related to the oscillatory and steady glow in this system in closed vessels and on the behaviour in flow reactors. 5.11.1 Closed vessel studies In addition to showing branched-chain ignition, with (first and second) p-Ta limits similar to those of the H 2 + 0 2 system, the CO + 0 2 reaction supports a response known as glow. This is a spontaneous chemiluminscent state, with a weak, pale-blue emission arising from electronically-excited

The CO + 0 2 reaction

515

CO* formed from the three-body recombination reaction (10)

CO + O + M ^ C O | + M

with the triplet O atom giving rise to a triplet C 0 2 molecule. This step is a termination process (there have been occasional suggestions that excited C 0 2 can support a process of "energy-branching" with its excess enthalpy being sufficient to allow the reaction (11)

COf + 0 2 - > C 0 2 + 2 0

but there is no actual evidence in support of this type of process in this reaction). The glow can take one of two forms. The earliest reports of steady glow are due to Prettre and Laffitte [55] who observed a long-lived chemiluminescence from a mixture of the reactants in a closed vessel. Ashmore and Norrish [56] were the first to report observations (made in 1939) of oscillatory glow, which they termed the "lighthouse effect". In this mode, "bursts" of chemiluminscent emission are separated by periods of relative "darkness". Depending on the total pressure and oven temperature, short or long trains of emission were observed. Even in a closed system, they were able to record trains consisting of over 100 pulses. The reactor used by Ashmore and Norrish had been treated by exposure to chloropicrin (CC13N02) a known inhibitor of radical chains. If the surface was washed with acid, the oscillatory glow response was lost, but reappeared on further pre-treatment with the inhibitor. Linnett and various co-workers [57, 58] also observed oscillatory glow, somewhat intermittently, in their studies of the reaction with extensivelydried reactants (see above) in vessels with untreated surfaces. These authors unfortunately used the phrase "multiple explosions" to describe the phenomenon which they observed only by eye. Later, meticulous work exploiting instrumental monitoring of the reaction by Bond et al. [59, 60] revealed that the "ignition limit" reported by the Linnett group on the basis of a detectable emission actually corresponds to the limit for steady glow. The work by Bond has provided a firm framework for the interpretation of the isolated reports by Ashmore and Norrish, Linnett and co-workers and those of McCaffery and Berlad [61] and Aleksandrov and Azatyan [62], by paying particular attention to controlling the H 2 content of the

516

Global behaviour in simple oxidations

Ch. 5

mixture and by following light emission intensity (photomultiplier), selfheating (fine-wire thermocouple) and stable species concentrations (mass spectrometry) during the reaction. The extent of reaction accompanying a single oscillatory excursion is typically less than 0.5% if the reaction mixture is "dry" (<0.1% H 2 ), and there is then no detectable accompanying temperature excursion, i.e., the oscillatory glow can be an isothermal process. Even when the oscillatory response dies in a given experiment, there is typically 50% of the original fuel remaining in an initially stoichiometric mixture. The reaction may exhibit a period of steady glow (generally of diminishing intensity) following the end of the oscillatory phase, but this also ceases well before complete reactant consumption. With higher concentration of H 2 , the oscillations are more intense (sufficiently so to be confused with an ignition if only monitored by eye) and accompanied by significant reactant consumption and self-heating. This means that only short trains, perhaps only 2 or 3 pulses, are possible with the oscillatory phase again ceasing before complete consumption. Some exmples of the oscillatory glow behaviour observed in these systems are shown in Fig. 5.30. An important feature of these experiments was that oscillatory and steady glow were observed for the same ranges of pressure and ambient temperature, lying in a p-Ta peninsula with "first" and "second limits", as indicated in Fig 5.31. It appears that a conditioning of the surface arising from repeated experiments in the same vessel is responsible for determining whether a given experiment will exhibit oscillatory or steady luminescence. A newly-cleaned vessel will generally exhibit a steady glow for early experiments at a given pressure and ambient temperature, but eventually the system will support the oscillatory response for, perhaps, 100 consecutive experiments at the same/? and Ta. Eventually, the surface will age too much and the response will return to the steady glow until the vessel is cleaned again by acid washing. At high H 2 concentration (wet systems) the genuine ignition response invades the p-Ta region for glow. However, the inevitable effects of reactant consumption again indicate that these various phenomena are better studied quantitatively and characterized under continuous flow conditions. 5.11.2 Open systems In well-stirred flow reactors, the CO + 0 2 reaction supports five different modes of response [63,64]: steady slow reaction, steady glow, oscil-

517

The CO + 0 2 reaction

(a) 82

17 c

O)

JUu uuu LLULiUUuLJu 1000

i94 J_JL

Time/s

9000

(b)

(c) c £

O

o o o

2:

0

200 400 Time/s

Fig. 5.30. Representative trains of oscillatory glow for CO + 0 2 mixtures in a closed system; (a) chloropicrin-treated reactor with p = 33 Torr, Ta = 841 K giving a total of 94 oscillations over 3 hr; (b) and (c) light emission and reactant consumption in a clean reactor with p = 22 Torr and Ta = 900 K. (Reprinted with permission from reference [60], © Royal Society of London.)

Global behaviour in simple oxidations

518

Ch. 5

-(b)

duration decreases 800

900

1000

1000

ambient temperature, T,/K Fig. 5.31. The co-existence of (a) steady glow and (b) oscillatory glow for the CO + 0 2 reaction (dry mixtures) in an ageing reactor. (Reprinted with permission from reference [60], © Royal Society of London.)

latory glow, oscillatory ignition or a steady high reaction state ("steady flame") - the first and last two of these being the counterparts of the same phenomena in the H 2 + 0 2 system. We will also see that the oscillatory ignition response can show complex periodicities although perhaps of a different nature from those described in Section 5.7 for hydrogen. Figure 5.32 shows the p-Ta diagrams appropriate to five different mixtures of increasing H 2 content, indicating the different regions of the parameter plane for the basic responses. Steady slow reaction is always found at the lowest pressure and ambient temperature conditions. In this mode, the reaction may sustain a steady temperature excess of the reacting gases above the ambient, vessel temperature with corresponding extents of reaction of up to 10% of the inflow concentration of the reactants. If the system is disturbed, perhaps by injecting a cold inert gas or by perturbing the gas flow, the system returns to the same steady-state after the perturbation is removed and does so in a simple, exponentially decaying manner, indicating stable nodal character of this steady-state. Steady and oscillatory glow responses occupy separate regions of the parameter plane under continuous flow conditions. It is possible that the flow of reactants in the period before the reaction is monitored (and this period must be rigorously reproduced in any given set of experiments if reproducible results are to be obtained) provides a consistent degree of surface conditioning in such experiments [65]. Oscillatory glow, which is

519

The CO + 0 2 reaction

750

775

800

TJK

Fig. 5,32. The p-Ta regions for steady dark reaction, steady glow, oscillatory glow and oscillatory ignition, and the influence of added H 2 for the CO + 0 2 reaction in a flow reactor with fres = 8s: (a) no added H 2 ; (b) 150 ppm added H 2 ; (c) 1500 ppm added H 2 ; (d) 7500 ppm added H 2 ; (e) 10% H 2 in final mixture. (Reprinted with permission from reference [64], © Royal Society of London.)

now sustained for as long as the flows are maintained, typically occupies a closed region at lower pressure within a steady glow peninsula. There is a different response to perturbation of the steady glow from that for steady dark reaction, with a damped oscillatory return in the case of glow indicating stable focal character. As the oscillatory glow is also clearly a stable limit cycle response, we can assume that the boundary between steady and oscillatory glow corresponds to the locus of Hopf bifurcation

520

Global behaviour in simple oxidations

Ch. 5

points, Section 5.4, with these generally being of a supercritical character as the oscillations have vanishingly small amplitude at the boundary itself. For systems with low concentrations of H 2 (i.e., if no H 2 is added and for which the natural impurities are probably less than 250 ppm H-containing species, or with 150 ppm H 2 added to these impurity concentrations), the oscillatory glow is effectively isothermal as no temperature excursion can be measured (the level of detection is approximately 4 K) and no reactant consumption can be detected. As the concentration of added H 2 increases, so temperature excursions and reactant consumption accompanying the oscillatory glow events can be recorded. Ignition is also a feature for such systems. Ignitions are readily distinguished from glow by having temperature excursion well in excess of 200 K, effectively complete consumption of the fuel (in equimolar or stoichiometric mixtures) and significant emission at wavelengths other than that corresponding to CO*. With increasing H 2 concentrations, the p-Ta peninsula for ignition invades the accessible p-Ta region and begins to "cover" the regions of oscillatory and steady glow. The boundaries for the latter do not appear to be very sensitive to the H 2 concentration and even with 10% H 2 added, there is a remnant of the steady glow region with the limit separating it from steady dark reaction in much the same position as for the mixture with no added H 2 . In a sequence of experiments, Johnson [66, 67] has investigated the development of the oscillatory ignition waveform for mixtures containing a "significant" addition of H 2 (typically 2% of the fuel and so 1% of the total mixture). The p-Ta diagram for a system with an average mean residence time of 16 s is shown in Fig 5.33. The slow reaction at low Ta gives way to oscillatory ignition as the (second) ignition limit is traversed. As r a is increased further, so the waveform changes quantitatively (in amplitude and period). More interestingly, perhaps, is that for a range of operating pressure from approximately 15 to 45 Torr, there is a qualitative change in the waveform as the system traverses a region of "complex oscillation". A flavour of the evolution of the waveform can be obtained by allowing the ambient temperature to drift upwards through this region at a fixed total pressure. Such a sweep through this region at a relatively high pressure, p = 40 Torr, where the region of complex response is relatively narrow, is shown in Fig 5.34(a). Complexity here evolves through a perioddoubling bifurcation, with the period-1 type oscillation giving way to a

521

The CO + 0 2 reaction

700

750

800

ambient temperature, Ta/K Fig. 5.33. The p-Ta regions of slow reaction, simple and complex oscillatory ignition and steady ignited state for CO + 0 2 reaction with 0.5% H 2 in a flow reactor with fres = 16 s. (Reprinted with permission from reference [66], © Royal Society of Chemistry.)

waveform that has one large + one small peak separated by approximately the same period per repeating unit. Thus, the repeat period is approximately twice that for the original period-1 response. Notice that as this change is encountered, the difference in amplitude between "large" and "small" grows smoothly: it is possible to sustain the system just inside this region such that the two peaks have amplitudes that differ by as small an amount as desired (at least within the experimental control allowed in this experiment). Thus, this is a supercritical period-doubling bifurcation. The experimental record in Fig. 5.34(a) shows that the system has virtually settled into a sustained period-2 responses before the ambient temperature is increased further. The response to this further increase in Ta is that the system undergoes a supercritical period-halving bifurcation and settles onto a small amplitude period-1 oscillation. A further increase in Ta causes the system to traverse a supercritical Hopf bifurcation and the oscillatory nature of the ignition event gives way to a stable steady-ignited state. At a lower total pressure, e.g., p = 20Torr, the region of complex oscillatory behaviour is of wider extent in Ta. The traverse through this region uncovers more exotic responses after the initial period doubling, as indicated in Fig. 5.34(b). To characterize the various responses in more detail, the system must be allowed to settle to a steady ambient temperature (and all other experimental parameters) at various points within this range. Selected time series obtained in this way at different Ta are shown

522

Global behaviour in simple oxidations i

0

1

Ch. 5 1

1 2 time/min

120 H

P 80H < I

806.1K

808.3K

809.6K •

Ta slowly increasing

i

i

i

120 ^80

< 0 777.3K

- - -•-•• 778.0K 781.8K

782.1K

782.6K

782.8K

• T, slowly increasing

Fig. 5.34. Traverse through region of complex ignition in Fig 5.33 for p = 40Torr (upper trace) and p = 20Torr (lower trace) showing evolution of complex oscillatory waveform. (Reprinted with permission from reference [67], © American Institute of Physics.)

in Fig. 5.35(a)-(f). The original period-1, (a), undergoes period doubling to give period-2 (b). The period-2 response also period-doubles as Ta is increased to give a period-4 response, (c). Evidence for period-8, following another period doubling, and even for period-16 has also been reported, but the control over the experimental conditions in the apparatus is not sufficiently high that such high periodicities can be sustained for data capture. Figure 5.35(d) corresponds to a complex response of no apparent periodicity, but for certain very narrow ranges of ambient temperature within this basic chaotic evolution, the system manages to find periodic states which repeat every 5- or 3-ignition events, (e) and (f), respectively. As the region of chaos is exited at high Ta a reverse period-doubling cascade is observed down through period-8, period-4 and period-2 to a small amplitude period-1 oscillation that eventually undergoes a Hopf

The CO + 0 2 reaction

.(d)

— i

1—

Time / mins

Fig. 5.35. Transient-free records of different periodic states from region of complex oscillation in Fig. 5.33: (a) period-1; (b) period-2; (c) period-4; (d) aperiodic trace; (e) period5, (f) period-3. (Reprinted with permission from reference [66], © Royal Society of Chemistry.)

bifurcation as the steady-ignited state about which all this temporal complexity is evolving gains stability. The reconstructed limit cycles for the five periodic states from Fig. 5.35(a-c, e, f) are shown in Fig. 5.36. The aperiodic time series in Fig. 5.35(d) has no repeating unit, but is not a random sequence of ignition. That there is an underlying order can be revealed by plotting the so-called next-maximum map. In this construction, the amplitude of one ignition (as measured by any appropriate and convenient device, so here we use the signal from the internal thermocouple junction which is an approximate measure of the self-heating

524

Global behaviour in simple oxidations

dtgtialfeed ihermocouple record

Ch. 5

digiiatized thermocouple record

Fig. 5.36. Reconstructed attractors corresponding to time series in Fig. 5.35. (Reprinted with permission from references [66 and 67], © Royal Society of Chemistry and American Institute of Physics.)

of the gas) is plotted against the amplitude of the next peak. This process of finding x-y pairs is repeated for every ignition in the sequence. If the response were random, then the points on the graph obtained in this way would be randomly scattered. As shown in Fig. 5.37, however, the CO ignition lie so as to construct a relatively simple, single-humped maximum

The CO + 0 2 reaction

525

+ o

is

amplitude of n-th maximum Fig. 5.37. Next maximum map obtained from chaotic sequence in Fig. 5.35(d). (Reprinted with permission from reference [66], © Royal Society of Chemistry.)

map. This is a characteristic feature of deterministic chaos and helps distinguish complex responses from noise. We can also reconstruct the attractor corresponding to the time series using the method of delays described above. The strange attractor arising for this chaotic sequence is also shown in Fig. 5.37 and its relationship to the simple limit cycles from the periodic behaviour can be recognized. Complex behaviour is favoured in the CO + 0 2 system by long residence times and by low concentrations of H 2 . With residence times of 10 s or less, the region of complex responses includes only period-2: no further period doublings are evidenced. At tres = 16 s, the region of complex ignition moves to higher ambient temperatures as the H2-content decreases. The full period-doubling cascade to chaos described above for 0.5% H 2 is lost as the concentration of H 2 increases. With 1.33% H 2 , period-2 and period-4 are found, while for 2% H 2 we find only period-2 responses. The region of complex responses covers a decreasing range of pressure as [H2] increases. The behaviour of "pure" CO, i.e., the cylinder gas (CP grade) without added H 2 is of interest in that no period-doubling is observed, but a remnant of the mixed-mode type of oscillation seen for the H 2 + 0 2 system arises.

526

Global behaviour in simple oxidations

Ch. 5

5.11.3 Mechanistic interpretation and modelling The reaction steps involving CO and 0 2 are rather limited: (12)

CO + 0 2 - * C 0 2 + 0

(10)

CO + O + M ^ C O f + M

the first being an initiation step and the second a termination. Thus, a "pure" CO + 0 2 system thus has no branched chain processes and could only support a thermal ignition at best. The fact that all systems include some hydrogen-containing species, and the experiments described in detail above actually have H 2 added, indicates that we should also consider the H 2 + 0 2 reactions and any additional reactions involving CO with intermediate species in that mechanism. The most important step is (13)

CO + O H ^ C 0 2 + H

which is the primary source of C 0 2 in such systems. It is a chain propagation step, converting one chain carrier OH into another H and so neither enhances nor diminishes the propensity to chain-branching. Another step (14)

CO + H 0 2 -> C 0 2 + OH

is relatively unimportant, but has the potential significance that it provides a reaction channel for the relatively unreactive radical H 0 2 and thus would help sustain the radical chain. CO and C 0 2 can also play the third-body role in the reaction (5)

H + 02 + M ^ H 0 2 + M

(a) Modelling the ignition limit The ignition limit is a feature of systems containing approximately 0.5% H 2 or more. If we begin the discussion by imagining a pure H 2 + 0 2 system in which the H 2 is then systematically replaced by CO, we find initially that the (second) limit moves to lower ambient temperatures and higher pressures. This feature can be modelled quantitatively simply by including the relative third body effect of CO, which is less than 1 (i.e., CO is a less

527

The CO + 0 2 reaction

55

(ii) (i)

fa o

15

700

ambient temperature, Ta/K

780

Fig. 5.38. Variation of CO + H 2 + 0 2 ignition limit with fuel mixture composition, all mixtures have 50% 0 2 : (i) 50% H 2 ; (ii) 10% H 2 + 40% CO; (iii) 2% H 2 + 48% CO; (iv) 1.33% H 2 ; (v) 0.5% H 2 ; (vi) 0.33% H 2 , (v) 0.167% H 2 . (Reprinted with permission from reference [67], © American Institute of Physics.)

efficient third body than H 2 in the primary termination process in the H 2 + 0 2 second limit mechanism) so replacing H 2 by CO decreases the effective termination rate and favours ignition. The ignition limit reaches its "minimum" position, furthest to the left on the p-Ta diagram as indicated in Fig. 5.38, for a mixture composition with approximately 10% H 2 . A further replacement of H 2 with CO sees the limit moving back to higher ambient temperature, as indicated in the figure. This "inhibitory" effect arises for systems with low [H2] for which the extra termination step (10) involving CO becomes able to compete with the branching process (3). Using the reactions

(i)

OH + H 2 ^ H 2 0 + H

(2)

H + 02^OH + 0

(3)

0 + H2^OH + H

(13)

CO + O H - * C 0 2 + H

(5)

H + 02 + M ^ H 0 2 + M

(10)

CO + O + M -> C 0 2 + M

528

Global behaviour in simple oxidations

Ch. 5

a quasi-steady-state analysis on the radical intermediates yields the following expression for the net branching factor [67]

* = 1 |

ir r r n i r ^ i ~ ^ [ M ] ' fc10[CO][M]

(5 75)

'

k3[H2] where k5 = k™2{xU2 + a02x02 + acoXco + K}, (5.76) with xU2 + xCo = 0.5 for an equimolar fuel + 0 2 mixture. At relatively high H 2 , the influence of CO is felt mainly through the final term in equation (5.75) and with aco < 1, increasing x c o , and, hence, decreasing jtH2> has the effect of decreasing the effective value of k4 and, hence, increasing the value of
Hydrocarbon oxidation

529

mode is reset by the direct inflow of fresh H 2 and, for ignition, of CO [68]. In a closed system, much the same chemistry must operate (there may be some differences arising from the slightly different p-Ta range for the oscillatory glow in these systems compared with flow reactors and this, in turn, may be related to different states of the reactor surface). The major problem in interpreting oscillatory glow in the above terms, i.e., as an H 2 + 0 2 ignition in a system heavily diluted with CO which is acting primarily as an inert gas, is to understand the resetting mechanism by which H 2 is regenerated. Suggestions due to Babushok et al. [69] and to Gray et al. [64] invoke the "water gas shift" process (15)

CO + H 2 0 ^± C 0 2 + H 2

perhaps occurring on the surface, hence, the sensitivity of closed vessels studies to the pre-treatment processes mentioned above. Chinnick and Griffiths [70] have successfully simulated some of the closed vessel behaviour in these terms. The complex periodicities in the CO + 0 2 reaction are only now yielding to simulation, but are again arising from a mechanism based on an H 2 + 0 2 scheme with the few steps involving CO described above.

5.12 HYDROCARBON OXIDATION

"Exotic" oscillatory and other types of non-linear behaviour are also features of most hydrocarbon oxidations [71-74]. The next chapter will provide a detailed mechanistic description of the basis for cool-flames etc., and their relevance in various situations. It is interesting, however, to apply the classification system developed in the previous sections to the global behaviour in these systems. We start with a description of the oxidation of acetaldehyde (ethanal) and again concentrate on modern studies in flow reactors where the effects of reactant consumption (which are much more significant in closed systems for these cases than for CO) are not a feature.

530

Global behaviour in simple oxidations

Ch. 5

5.12.1 Oxidation of acetaldehyde (a) Experimental features The oxidation chemistry of small, partially-oxygenated fuels is of great interest in combustion chemistry as these are important intermediates in the combustion of virtually all commercial hydrocarbon fuels. Fuels with long carbon backbones react in their early stages mainly through a sequence of reactions that cause chain rupture, yielding smaller hydrocarbon fragments such as radicals. These then typically react with 0 2 to produce precursors of aldehydes, ketones etc. Not all of the features of acetaldehyde chemistry are completely representative of hydrocarbon oxidation, but this point is developed in the next chapter. The basic global behaviour of a mixture of acetaldehyde vapour in 0 2 is illustrated by reference to the p-Ta "ignition" diagram, Fig. 5.39. Up to five regions of qualitatively different responses are characteristic [75]. At low ambient temperature and pressure, the system exhibits a steady dark reaction, Region I. This may support a measurable steady-state tem-

150

100

50

450

500

550 Ta/K

600

650

Fig. 5.39. Experimental p-Ta diagram for CH 3 CHO + 0 2 in a flow reactor with tTes = 3 s showing five regions of different qualitative behaviour. (Reprinted with permission from reference [75], © Royal Society of London.)

531

Hydrocarbon oxidation

perature excess, as the reacting gases self-heat through the exothermic reaction beyond the oven temperature, of up to 40 K. The temperature excess ATSS increases with increasing Ta through this region, indicating that the rate of chemical heat release di has a positive temperature coefficient. The oxidation process is incomplete under these conditions, with the major product species detectable including per ace tic acid (CH3CO3H), methanol (CH 3 OH), formaldehyde (CH 2 0) and methane (CH4) in addition to C 0 2 and H 2 0 . If the system is perturbed in this region, the steady-state reveals a stable nodal character. Region V corresponds to the other steady-state form, that of steady glow, and lies at the highest ambient temperatures. As the name indicates, reaction is accompanied by a steady chemiluminescence. In contrast to the CO oxidation, however, the molecular source is not excited C 0 2 but excited formaldehyde HCHO*. Reaction is also accompanied by a steadystate temperature excess. In this region, however, ATSS typically decreases as the reacting gas temperature increases, as illustrated in Fig. 5.40. This phenomenon, termed the negative temperature coefficient (n.t.c), indicates that the rate of chemical heat release decreases as Tgas increases. (It is important to be careful on two counts here: first it is the gas temperature rather than the ambient temperature that is referred to even though it is the ambient temperature that determines whether the system is in region V or not: second, the rate of heat release is not the same as the "rate of

c o

2"

40

20

v

80

X

c o U

3 cj

-Uo

©

cu

£ S c

3 8

500 600 700 © reactant temperature, T/K Fig. 5.40. The steady-state heat release rate from CH 3 CHO + 0 2 in regions I and V as a function of reacting gas temperature showing negative temperature dependence in region V. (Reprinted with permission from reference [75], © Royal Society of London.) 400

532

Global behaviour in simple oxidations

Ch. 5

reaction" - the latter being a somewhat ambiguous quantity anyway.) The major products in this region are CH 3 OH and other partially oxidized, C 0 2 , CO and H 2 0 . The steady-state is a stable focus as perturbations decay in a damped oscillatory manner. The degree of damping decreases as we approach the boundary with region IV. If the ambient temperature is decreased, the system enters region IV in which the reaction exhibits oscillatory cool-flames. Each oscillatory event is accompanied by emission from HCHO* and by a temperature excursion. Cool-flames were first described by Davy [76] who noted that they were not able to ignite paper held in the flame (he stabilized a cool-flame above a bed of sand in which the liquid reactant was absorbed). The temperature rise may be as much as 200 K, but is typically less than this: cool-flames in which the temperature excursion is only 10 K or so are regularly observed. Indeed the oscillatory amplitude seems to tend to zero as we approach the boundary with region V. This latter point indicates that the cool-flame corresponds to a simple limit cycle oscillation with the species concentrations and gas temperature oscillating around an unstable steady-state. The steady-state regains stability as the ambient temperature Ta is increased emerging as the stable focal, steady glow state in region V. Thus, the boundary between regions IV and V is that of a locus of (supercritical) Hopf bifurcation points in the terminology developed in Section 5.5. The products detected from reaction in region IV include CO, H 2 0 , CH 2 0, CH 3 OH and CH 4 and lower quantities of peracetic acid, ethane (C2H6) and hydrogen peroxide (H 2 0 2 ) all of which oscillate with the coolflame period as indicated in Fig. 5.41. On decreasing the ambient temperature further, we enter the regions II and III (the convention of splitting this region into two subregions has some convenience, but there is in fact a coherence to the behaviour across this whole range). At the lowest Ta in this region (i.e., in region II) the behaviour corresponds to a two-stage ignition event (see Chapter 6). Basically, the behaviour is that of an ignition in which a high transient temperature excess develops accompanying consumption of the primary reactants and the production of CO and H 2 0 and other "final" products. There is, however, some "fine structure" to the ignition waveform, with a distinct shoulder preceding the ignition proper. During this initial development phase, the products of the cool-flame response are observed to build up. Thus, we have a cool-flame as a precursor stage, building up intermediate species and the gas temperature before the transition to a

Hydrocarbon oxidation

533

-IJLLLL «

*

0

« 20

Time/s

Fig. 5.41. Typical experimental records in cool-flame region IV. (Reprinted with permission from reference [75], © Royal Society of London.)

fully-fledged ignition. The transition from regions I to II is accompanied by a "hard excitation" similar to that observed at the transition from slow reaction to oscillatory ignition in the H 2 + 0 2 system in the sense that it corresponds to a saddle-node bifurcation of the low reaction steady-state with the system then moving to a large amplitude limit cycle solution existing around the unstable high reaction steady-state (in this case corresponding to a complex oscillation). Region III corresponds to other types of multistage ignition events [77]. In this type of response, a two-stage ignition event is preceded by a number of distinct cool-flame oscillations, with some examples shown in Fig. 5.42. Thus, just across the border from region II, we observe a single cool-flame followed by a two-stage ignition event, giving rise to an overall 3-stage ignition. The number of cool-flames interspersing the two-stage ignition peak increases as we move across region III to higher Ta, giving 4-stage, 5-stage ignitions etc., with the ignition peak then disappearing from the sequence as we enter region IV. A schematic bifurcation diagram summarizing these changes in behaviour as a function of the ambient temperature is sketched in Fig. 5.43. We may also note that at yet higher Ta, beyond region V, a simple ignition limit may be encountered.

534

Global behaviour in simple oxidations

Ch. 5

CH3CHO

nn n 1

AC

HoO

t

1

1 •

1 !

;

•

i

1

•

\ \ CH20 CH3OH

20 Time/s

j

1

^ •M

mm

AT

1

i ;

AT

40

Fig. 5.42. Typical experimental records in two-stage ignition region III. (Reprinted with permission from reference [75], © Royal Society of London.)

(b) Basis of mechanistic interpretation A full mechanistic account follows in the next chapter. Here we simply indicate the important features of the currently-accepted interpretation of the above facts. The key feature in these thermokinetic phenomena is that there are both thermal and chemical feedback processes combining to produce the various "exotic" responses, including the ntc. At the heart of the clockwork is the equilibrium involving the methyl radical CH 3 and molecular oxygen [78] (16)

CH 3 + 0 2 ^ CH 3 0 2

It is evident that A5?6 < 0, indicating that the equilibrium will shift from

Hydrocarbon oxidation

535

Fig. 5.43. Schematic bifurcation diagram for CH 3 CHO + 0 2 reaction showing separate branches corresponding to dark reaction and steady glow, with upper branch losing stability at a Hopf bifurcation as Ta is reduced to give limit cycle (cool-flame) oscillations. The simple limit cycle also loses stability as Ta is reduced further and a complex oscillation corresponding to the multi-stage ignition will emerge but cannot be adequately represented in this 1-D diagram.

the right (CH 3 0 2 ) to the left (CH 3 + 0 2 ) as the temperature of the reacting mixture increases. Thus, at low gas temperatures we can expect the chemistry to be that characteristic of the peroxy species CH 3 0 2 but at high Tgas it will be that of the methyl radical. The reaction steps at the pressures and temperatures of interest that arise from the formation of CH 3 0 2 are the following: (17)

CH 3 0 2 + CH 3 CHO -> CH 3 0 2 H + CH 3 CO

a propagation step in the sense that one radical is consumed and one (CH3CO) is produced. However the "molecular" product CH 3 0 2 H is unstable to decomposition into two more radical species (18)

CH 3 0 2 H -» CH 3 0 + OH

so there is a degenerate branching in which three chain carriers eventually emerge from the CH 3 0 2 species. The increasing radical/chain carrier pool leads to an increasing reaction rate and, hence, to an increasing rate at which heat is released by this exothermic reaction channel.

536

Global behaviour in simple oxidations

Ch. 5

As r g a s increases in response to the increased rate of heat release, so the equilibrium (16) swings over to the lefthand side in favour of the methyl radical. Under the conditions typical of this reaction, the major fate for CH 3 is recombination to form ethane (19)

CH 3 + C H 3 ^ C 2 H 6

This termination step reduces the radical pool and alters the overall reaction exothermicity, both having the effect of reducing the rate of chemical heat release. In turn, the drop in 9t causes the gas temperature to fall, allowing the equilibrium to shift back to the right. For some conditions, this cycling in Tgas and the CH 3 or CH 3 0 2 chemistry continues to give sustained, limit cycle oscillations (cool-flame); for other operating conditions, a balance can be attained so that a steady-state (steady glow) is established. The negative temperature coefficient corresponding to steady glow also emerges naturally as the equilibrium lies further in the direction of the unbranched CH 3 chemistry leading to a lower 9t, the higher Tgas. Full scale numerical computation on relatively large mechanisms [78] have reproduced the major features (see also refs [79-84]). Simplified model interpretations have also been built on the Gray-Yang scheme discussed earlier [18-21]. This is a two-variable model (chain carrier x and gas temperature T) and can account qualitatively for the dark reaction, cool-flame, steady glow and n.t.c. phenomena and provide hints of the ignition-type events. For the complex oscillations, a three-variable model is needed and an elaboration on the Gray-Yang scheme has been provided by Wang and Mou [22] (see also refs [85-87]). Interestingly, although the computations and extended models can produce "mixed mode" multistage ignitions, one feature of the experimental responses is not matched. In the theoretical approaches, the increasing number of cool-flames as we move across region III is achieved by adding two extra cool-flames at each stage, so we see sequences of 2-, 4-, 6-stage ignitions etc., or of 3-, 5-, 7-stage etc., but not the 2-, 3-, 4-stage etc., reported in the experimental studies. 5.12.2 Non-isothermal oxidation of alkanes Cool-flames arise for some ranges of experimental conditions for many other hydrocarbon oxidation reactions. Chapter 6 presents many p-Ta ignition diagrams for different hydrocarbon fuels, all of which have a

Hydrocarbon oxidation

537

richness matched by the acetaldehyde system. An earlier discussion of the detailed chemistry relevant to such oxidations has already been presented in Chapter 1, but it is worth also stressing that the contributions of the heat generation and heat loss (thermal feedback) processes are as important as the kinetic complexity. If we denote a general hydrocarbon as RH, then initiation processes involving 0 2 may give rise to the corresponding alkyl radical R. We may anticipate an analogous equilibrium to that involving CH 3 produced in acetaldehyde oxidation of the form (20)

R + 02 ^ R02

with the same general dependence of the equilibrium on reactant temperature. This is indeed the case, but the subsequent reactions are not simple generalizations of the acetaldehyde schemes [74]. The major fate of R if the equilibrium lies to the left (the high Tgas regime) involves abstraction of an H atom from the C atom adjacent to the radical site (if this has a C—H bond), giving rise to the conjugate alkene and H 0 2 . More significant, however, is the difference in the R 0 2 chemistry in the low Tgas regime. The feasibility of the reaction (21)

R 0 2 + RH -> R 0 2 H + R

which would then be followed by the degenerate branching (22)

R 0 2 H ^ R O + OH

is greatly reduced for a typical hydrocarbon. For an aldehyde, the —CHO group provides a relatively labile H atom for abstraction by R 0 2 , but in the absence of an activating functionality, the basic C—H bond is too strong to allow this to be a sufficiently fast process at the cool-flame temperatures. This point is developed in the next chapter where the route involving intramolecular hydrogen abstraction is described. This mechanism allows us to account for the different cool-flame characteristics of different fuels and, in particular, for the different propensities of different isomers, such as n- and /-butane to show cool-flame behaviour (see Fig. 6.14) on the basis of their molecular structure.

538

Global behaviour in simple oxidations

Ch. 5

5.13 CONCLUSIONS AND FUTURE DIRECTIONS

The application of the methods of non-linear dynamics to chemical reactions is still in its relative infancy, but it's advocates should not rely solely on its novelty. It is time to ask the question "what can this approach really do for combustion chemists"? At the lowest level, it provides a new system with which to label and categorize different forms of global response (stable node, limit cycle etc.). Beyond this, however, is the power of then being able to bring to bear the body of knowledge pertaining to generic features of how, say, a stable node can give way to a limit cycle as some (experimental) parameter is varied. Important in this is the basic recognition of the role played by steady-state solutions of the governing reaction rate equations. First, if we use flow reactors, such steady-state solutions genuinely exist! (If we operate in a batch process, quasi-steadystates may arise in which the gradual background change in reactant concentrations plays the role of a slowly varying parameter rather than a significant reactant: we can then construct steady-states by ignoring reactant consumption as a first approach to modelling, see Chapter 4.) Steadystate solutions are, in principle at least, obtainable from solutions of algebraic equations (in spatially homogeneous systems) rather than requiring solutions of differential equations (in spatially-distributed systems, we may hope to reduce to problem to the solution of ordinary rather than partial differential equation if there are symmetries or restrictions on the geometry). Thus, we can exploit more efficient computational approaches and obtain the steady-state solutions directly. Second, it is possible (and not particularly unusual) to have more than one steady-state solution for the same set of experimental conditions (parameters). If steady-state multiplicity occurs, there are usually 3 (or less frequently 5 or 7) steady-states for the system to choose from. Multiplicity is the underlying cause of the phenomena of ignition and extinction. In the simplest case, with three co-existing steady-states, two may be stable (and, hence, physically attainable for the reaction system) and the other "middle" solution is unstable. The unstable solution is not attainable but does play an important role in separating the basins of attraction of the competing stable states, i.e., in separating initial states (concentrations and gas temperature) that evolve to the slow reaction state from those that evolve to the ignited steady-state. Thus, this separatrix is of potential real interest: methods for constructing the separatrix in the vicinity of the

Conclusions and future directions

539

unstable state and following it away across the concentration "state space" are already well developed in simple situations. As an experimental parameter is varied, so the different steady-states change smoothly in response. However, it is possible that two of the steady-states can approach each other and merge at some "critical" parameter value. Beyond this, the two states vanish and the system must move to the only remaining state: a discontinuous response to a small, smooth change in the parameter characteristic of an ignition. Thus, we can thus compute ignition points in terms of a condition in the steady-state behaviour: a saddle-node bifurcation. Such points can be located directly by adding the condition for a zero eigenvalue or, equivalently, a zero determinant for the corresponding Jacobian matrix as an additional algebraic condition to the steady-state solution. Furthermore, we can then follow this condition as some other parameter is varied, allowing our numerical package (or even our algebraic analysis in the simplest cases) to construct an ignition limit directly. In addition to merging with other steady-states, a stable steady-state can also give rise to a qualitative change in the behaviour of the chemical system by losing its stability. An important example of this is the Hopf bifurcation which often gives rise to the onset of oscillatory behaviour. By adding the condition for this bifurcation (a pair of purely imaginary eigenvalues) to the steady-state condition, such points can again be solved for directly and followed as a second parameter is varied. Additional characterization of a Hopf point will suggest whether the system shows soft or hard excitation to oscillations. At present, computational techniques really only allow us to work in this way with stiff systems of small numbers of intermediate species (<6). As model reduction techniques become further developed and more sophisticated, and as computing power increases, there is the prospect of more combustion systems being tackled in this way. We also now know that complex oscillations evolve as simple limit cycles become unstable, bifurcating to more complex limit cycles. Only a small number of bifurcation sequences account for all known scenarios. We have seen examples of mixed-mode sequences (H 2 + 0 2 ) and period-doubling cascades (CO + 0 2 ). A third route involving quasi-periodic responses is known and arises in some chemical system [88], but has not yet been observed in combustion systems (except in some special studies in which the ambient temperature or some other parameter is "forced" to vary in some sinusoidal or other periodic manner [89]). The important lesson then

540

Global behaviour in simple oxidations

Ch. 5

is simply that systems do not (except in very, very rare special cases) change straight from a simple response to a very complex response in one step. Except in pathological cases there will always be one of the above sequences that effect the transition from simplicity to complexity and a small number of additional observations (made with some guidance form the above theory) should suffice to decide which is actually responsible in a particular instance. Once the bifurcation sequence has been identified, it is usually possible to predict which other types of response will be found "close by" in terms of the experimental conditions. For instance, a concatenated mixed-mode state consisting of alternate l 1 and l 2 waveforms has been observed recently for a narrow, but finite range of ambient temperature between the l 1 and l 2 states in the H 2 + 0 2 system. Another important lesson arises for modelling/mechanism development. In the past, the observations of oscillations could be regarded as unusual and modellers were perhaps tempted to conclude that any hand-waving argument that supported such behaviour must be correct. It is now realized that the requirements for oscillations are not particularly stringent (some moderately strong feedback) nor unusual (chain-branching, self-heating) and that typically several different oscillatory clockwork mechanisms can be produced for chemically-complex systems. It is now imperative that such suggestions are subjected to quantitative tests. Simply reproducing some oscillating concentration records is not enough. However, we do now know in what form we should obtain data for such quantitative tests: we can observe how oscillatory amplitudes grow as we move away from the onset conditions for oscillations; we can see whether the period lengthens logarithmically or as the inverse square-root of the ambient temperature as we approach an extinction limit. Similarly, although not discussed at any great length here, we now have answers to the question "what should I do with a very complex (perhaps chaotic) record"? First, there are well-specified tests to determine whether the signal is actually chaotic or just random noise (or chaos with noise superimposed). One important, but frequently overlooked approach in this regard, is to vary the experimental conditions slightly to see if period4, period-2 or mixed-mode states etc. arise - a good indication that the signal does indeed have a chaotic basis. Next, we can gain from the data some coefficients known as Lyapounov exponents that give as a measure of how quickly the system is "losing information". This is relevant in chaotic systems for making forward predictions of the subsequent evol-

References

541

ution. System with parameter values such that the response is chaotic are extremely sensitive to their initial conditions. (Note that this sensitivity is not about whether the system will be chaotic or not: that depends only on the parameter values, not on the initial conditions.) This sensitivity, means that two virtually identical systems will typically begin to evolve in a similar way, but will sooner or later begin to behave differently, i.e., they will become uncorrelated. Thus, systems in their chaotic mode are essentially unrepeatable and also consequently unpredictable in the long term. The Lyapounov exponent tells us whether these differences develop quickly or more slowly. This is important because it allows us to calculate a timescale on which the behaviour is predictable within a certain tolerance. For instance, we might find that the value of the exponent is such that we will be able to predict the evolution to ±10% accuracy for lOmin, but if we require ± 1 % accuracy, we will only be able to predict for 10 s. It is perhaps too early to expect major modifications to rate data to be made on the basis of trying to match oscillatory or other non-linear behaviour (although this has happened in other areas of chemistry where some "guestimated" rate constants were revised by a factor of 105 on the basis of an analysis to match the observed speed of chemical waves in solution). Nevertheless, the semi-quantitative matching of experiment and numerical predictions for a whole p-Ta ignition diagram is a stringent test both of a proposed mechanism and of the associated rate data (and any assumptions about the heat and mass transfer processes): and the more "interesting" the p-Ta diagram in terms of regions of oscillations, complex oscillations etc., the more stringent the test. We can hope to see chemical kineticists and combustion scientists becoming more familiar with the concepts of nonlinear dynamics to the extent that they become recognized as additional extra tools in our "detective kit" for teasing out our understanding of the behaviour of complex chemical systems.

References [1] H. Richtering, in: Low Temperature Combustion, ed W. Jost, Chapter 2 (Gordon and Breach, 1965) pp. 37-82. [2] B. Lewis and G. von Elbe, Combustion, Flames and Explosions of Gases, 3rd edn. (Academic Press, 1985). [3] N.N. Semenov, Chain Reactions (Clarendon Press, Oxford, 1935).

542

Global behaviour in simple oxidations

Ch. 5

G. Dixon-Lewis and D J . Williams, Comprehensive Chemical Kinetics, eds C H . Bamford and Tipper, Vol. 17, (1977) pp. 1-239. P. Gray and S.K. Scott, Chemical Oscillations and Instabilities (Oxford University Press, 1994) S.K. Scott, Oscillations, Waves, and Chaos in Chemical Kinetics (Oxford University Press, 1994). N. Chirkov, Acta. Phys. Chim. URSS 6 (1937) 915. V.G. Vorokonov and N.N. Semenov, Zh. Fiz. Khim. 13 (1939) 1695. R.R. Baldwin and R.W. Walker, in: Essays in Chemistry, eds J.A. Barnard, R.D. Gillard and R.F. Hudson, Vol. 3, (1972) pp. 1-37. T. Boddington, P. Gray and G.C. Wake, Proc. Roy. Soc. Lond. A 357 (1977) 403. S.R. Kay, S.K. Scott and A.S. Tomlin, Chem. Eng. Sci. 44 (1989) 1129. J.L. Liu and S.K. Scott, Dynam. Stab. Sys. 8 (1993) 273. D.A. Frank-Kamenetskii, Diffusion and Heat Transfer in Chemical Kinetics (Plenum Press, 1969). I.Ye. Salnikov, Dokl. Akad. Nauk SSSR 60 (1948) 405. I.Ye. Salnikov, Zh. Fiz. Khim. 23 (1949) 258. (a) P. Gray, S.R. Kay and S.K. Scott, Proc. Roy. Soc. Lond. A 416 (1988) 321. (b) S.R. Kay and S.K. Scott, Proc. Roy. Soc. Lond. A 416 (1988) 343. (a) J.F. Griffiths, S.R. Kay and S.K. Scott, 22nd Symp. (Int.) Comb. (1988) 1597. (b) P. Gray and J.F. Griffiths, Comb, and Flame 78 (1989) 87. C H . Yang and B.F. Gray, J. Phys. Chem. 73 (1969) 3395. C H . Yang and B.F. Gray, Trans. Faraday Soc. 65 (1969) 1614. B.F. Gray, Trans. Faraday Soc. 65 (1969) 1603. C H . Yang, J. Phys. Chem. 73 (1969) 3407. X.-J. Wang and C.Y. Mou, J. Chem. Phys. 83 (1985) 4554. P. Gray and S.K. Scott, in: Oscillations and Traveling Waves in Chemical Systems, eds R.J. Field and M. Burger, Chapter 14, (John Wiley, 1985) pp. 493-528. J.F. Griffiths, S.K. Scott and R. Vandamme, J. Chem. Soc. Faraday Trans. 1, 77 (1981) 2265. R.R. Baldwin, D. Jackson, R.W. Walker and S.J. Webster, Trans. Faraday Soc. 63 (1967) 1665. R.R. Baldwin, D. Jackson, R.W. Walker and S.J. Webster, Trans. Faraday Soc. 63 (1967) 1676. R.R. Baldwin, D. Jackson, R.W. Walker and S.J. Webster, 10th Symp. (Int.) Comb. (1965) 423. R.R. Baldwin, M.E. Fuller, J.S. Hillman, D. Jackson and R.W. Walker, J. Chem. Soc. Faraday Trans. 1, 70 (1974) 635. P. Gray, J.F. Griffiths and S.K. Scott, Proc. Roy. Soc. Lond. A 394 (1984) 243-58. P. Gray, J.F. Griffiths and S.K. Scott, 20th Symp. (Int.) Comb. (1984) 1809. P. Gray, J.F. Griffiths, A.J. Pappin and S.K. Scott, in: Complex Chemical Reaction Systems, eds J. Warnatz and W. Jager (1987) pp. 150-159. [32] J.F. Griffiths, W. Kordylewski, A.J. Pappin and A.F. Sykes, in: Spatial Inhomogeneities and Transient Behaviour in Chemical Kinetics, eds P. Gray, G. Nicolis, F. Baras, P. Borckmans and S.K. Scott (Manchester University Press, 1987), pp. 237-253.

References

543

[33] J.F. Griffiths, B.R. Johnson and S.K. Scott, J. Chem. Soc. Faraday Trans. 87 (1991) 523. [34] D.L. Baulch, J.F. Griffiths, B.R. Johnson and R. Richter, Proc. Roy. Soc. Lond. A 430 (1990) 151. [35] D.L. Baulch, J.F. Griffiths and R. Richter, Chem. Eng. Sci. 46 (1991) 2315. [36] D.L. Baulch, J.F. Griffiths and R. Richter, Proc. Roy. Soc. Lond. A 434 (1991) 399. [37] D.L. Baulch, J.F. Griffiths and R. Richter, Phil. Trans. Roy. Soc. Lond. A 337 (1991) 199. [38] D.L. Baulch, J.F. Griffiths and R. Richter, Comb, and Flame 85 (1991) 271. [39] B.R. Johnson, PhD thesis (University of Leeds, 1991). [40] D.L. Baulch, J.F. Griffiths, A.J. Pappin and A.F. Sykes, Comb, and Flame 73 (1988) 163. [41] W. Kordylewski and S.K. Scott, Comb, and Flame 57 (1984) 127. [42] D.L. Baulch, J.F. Griffiths, A.J. Pappin and A.F. Sykes, J. Chem. Soc. Faraday Trans. 1, 84 (1988) 1575. [43] K. Chinnick, C. Gibson, J.F. Griffiths, and W. Kordylewski, Proc. Roy. Soc. Lond. A 405 (1986) 117. [44] K. Chinnick, C. Gibson and J.F. Griffiths, Proc. Roy. Soc. Lond. A 405 (1986) 129. [45] B.R. Johnson, S.K. Scott and A.S. Tomlin, J. Chem. Soc. Faraday Trans. 87 (1991) 2539. [46] H. Kuhl, Z. Phys. Chem. 44 (1903) 385. [47] M. Bodenstein and F. Ohlmer, Z. Phys. Chem. 53 (1905) 166. [48] B. Topley, Nature 125 (1930) 560. [49] V.E. Cosslett and W.E. Garner, Trans. Faraday Soc. 26 (1930) 190. [50] E.J. Buckler and R.G.W. Norrish, Proc. Roy. Soc. Lond. A 167 (1938) 292. [51] G. Hadman, H.W. Thompson and C.N. Hinshelwood, Proc. Roy. Soc. Lond. A 138 (1932) 297. [52] G. Dixon-Lewis and J.W. Linnett, Trans. Faraday Soc. 49 (1953) 756. [53] D.E. Hoare and A.D. Walsh, Trans. Faraday Soc. 50 (1954) 37. [54] A.S. Gordon and R.J. Knipe, J. Phys. Chem. 59 (1955) 1160. [55] M. Prettre and P. Laffitte, C.R. Hebd. Seances Acad. Sci. 189 (1929) 177. [56] P.G. Ashmore and R.G.W. Norrish, Nature 167 (1951) 390. [57] P.G. Dickens, J.E. Dove and J.W. Linnett, Trans. Faraday Soc. 60 (1964) 539. [58] J.W. Linnett, B.G. Reuben and T.F. Wheatley, Comb, and Flame 12 (1968) 325. [59] J.R. Bond, P. Gray and J.F. Griffiths, Proc. Roy. Soc. Lond. A 375 (1981) 43. [60] J.R. Bond, P. Gray, J.F. Griffiths and S.K. Scott, Proc. Roy. Soc. Lond. A 381 (1982) 293. [61] B.J. McCaffery and A.L. Berlad, Comb, and Flame 26 (1976) 77. [62] E.N. Aleksandrov and V. Azatyan, Combust. Explos. Shock Waves 12 (1977) 407. [63] P. Gray, J.F. Griffiths and S.K. Scott, Proc. Roy. Soc. Lond. A 397 (1985) 21. [64] P. Gray, J.F. Griffiths and S.K. Scott, Proc. Roy. Soc. Lond. A 402 (1985) 187. [65] P. Gray, J.F. Griffiths and S.K. Scott, J. Chim. Phys. 84 (1987) 49. [66] B.R. Johnson and S.K. Scott, J. Chem. Soc. Faraday Trans. 86 (1990) 3701. [67] B.R. Johnson, J.F. Griffiths, and S.K. Scott, Chaos 1 (1991) 387.

544

Global behaviour in simple oxidations

Ch. 5

[68] J.F. Griffiths and A.F. Sykes, J. Chem. Soc. Faraday Trans. 85 (1989) 3059. [69] V.I. Babushok, E.A. Novikov and V.S. Babkin, Dokl. Akad. Nauk SSSR 271 (1982) 878. [70] K. Chinnick and J.F. Griffiths, J. Chem. Soc. Faraday Trans. 2, 82 (1986) 881. [71] J.F. Griffiths, in: Oscillations and Traveling Waves in Chemical Systems, eds R.J. Field and M. Burger, Chapter 15 (John Wiley, 1985) pp. 529-564. [72] J.F. Griffiths, Ann. Rev. Phys. Chem. 37 (1985) 77. [73] J.F. Griffiths, Adv. Chem. Phys. 64 (1986) 203. [74] J.F. Griffiths and S.K. Scott, Prog. Energy Comb. Sci. 13 (1987) 161. [75] P. Gray, J.F. Griffiths, S.M. Hasko and P.G. Lignola, Proc. Roy. Soc. Lond. A 374 (1981) 313. [76] H. Davy, Phil. Trans. Roy. Soc. Lond. A 107 (1817) 77. [77] P. Gray, J.F. Griffiths, S.M. Hasko and P.G. Lignola, Comb, and Flame, 43 (1981) 175. [78] J.F. Griffiths and A.F. Sykes, Proc. Roy. Soc. Lond. A 422 (1989) 289. [79] M.P. Halstead, A. Prothero, and C.P. Quinn, Proc. Roy. Soc. Lond. A 322 (1971) 377. [80] M.P. Halstead, A. Prothero and C.P. Quinn, Comb, and Flame, 20 (1971) 211. [81] M.P. Halstead, L.J. Kirsch, A. Prothero and C.P. Quinn, Proc. Roy. Soc. Lond. A 346 (1975) 515. [82] L.J. Kirsch and C.P. Quinn, 16th Symp. (Int.) Comb. (1976) 233. [83] A.H. Harrison and L.R. Cairnie, Comb, and Flame 71 (1988) 1. [84] A.J. Harrison, R.M. Furzeland, R. Summers and L.R. Cairnie, Comb, and Flame 72 (1988) 119. [85] X.-J. Wang, Comb, and Flame 75 (1989) 107. [86] X.-J. Wang and P. Gaspard, in: Spatial Inhomogeneities and Transient Behaviour in Chemical Kinetics, eds P. Gray, G. Nicolis, F. Baras, P. Borckmans and S.K. Scott (Manchester University Press, 1987) pp. 687-690. [87] P. Gaspard and X.-J. Wang, J. Stat. Phys. 48 (1987) 151. [88] S.K. Scott, Chemical Chaos (Oxford University Press, 1993). [89] S.A. Pugh, M.A. Schell and J. Ross, J. Chem. Phys. 85 (1986) 868.