## Abstract

We investigate a diffusional-thermal model with two-step competitive exothermic reactions for premixed combustion wave propagation in one spatial dimension under adiabatic conditions. A criterion based on the crossover temperature notion was used to qualitatively predict the region in the space of parameters where three travelling combustion wave solutions coexist, which are further studied via numerical means. It is demonstrated that under certain conditions the flame speed is an ‘S’-shaped function of parameters. The fast branch is either stable or is partly stable and exhibits the Andronov–Hopf (AH) bifurcation before the turning point is reached. The mid-branch is completely unstable. The slow solution branch is either unstable or partly stable and exhibits a single or a pair of AH bifurcations. The AH bifurcations are shown to be supercritical giving rise to stable pulsating waves. Bistability and hysteresis phenomena are also demonstrated.

## 1. Introduction

Combustion waves with two-stage competing exothermic reactions and have been studied experimentally [1,2], analytically [3–7] and numerically [5]. In [1,2], the experimental observation of bistability of flame propagation is reported for combustion of Me−C−H_{2}, where Me is either Ti or Zr.

Throughout this work, we use the term bistability to mean the coexistence of two stable travelling combustion waves. In this system, the reaction can proceed along two competing paths with different end products, either metal carbo-hydrides or metal hydrides. The first reaction channel is characterized by high combustion velocity and temperature, whereas the latter by low combustion wave speed and burning temperature. Depending on the experimental parameters, hydrogen pressure, mixture composition, etc., there exist conditions that both regimes of flame propagation can be realized by using high- or low-temperature ignition sources. It was also demonstrated that changing the control parameters may also lead to the disappearance of the bistable regime of flame propagation.

The first analytical investigation of combustion waves with two-stage exothermic reactions dates back to 1975 [3]. In this paper by Berman & Riazantsev [3], the diffusional-thermal model for the planar flame propagation is studied by using the activation energy asymptotic (AEA) method. The flame temperature and speed are calculated for cases of gaseous and solid fuel. Although it is not explicitly formulated in this paper and it is not directly acknowledged in the latter publications [5,6,8], the analytical results of Berman & Riazantsev [3] suggest that the dependence of the flame speed and burning temperature on parameters can be multi-valued ‘S’-shaped function of parameters. There are either three solutions travelling with different velocities or a single solution depending on the parameters of the problem.

Non-uniqueness of the travelling wave solutions in this model is discussed in [4–7]. In [5], the analysis based on the AEA is carried out for both solid and gaseous fuel for the cases of the zero and first-order reactions. It is clearly demonstrated that there can exist three solutions travelling with different velocities and burning temperatures. Also in [5], a brief report is given on the results of the numerical analysis of the model, which was done by the authors earlier. However, to the best of our knowledge, these results have not appeared as a regular publication. In particular, it is stated that the propagation regime corresponding to the intermediate flame speed is always unstable, whereas the branches with fast and slow flame speeds can be stable or lose stability in an oscillatory manner. In [7], the model is generalized to include an additional step , so that the resulting scheme is a combination of the two-step sequential and competitive reaction mechanisms. It is argued that the resulting model can be applied to the description of combustion of coal and mixtures of transition metals with boron. Further, the properties of the travelling combustion waves were studied and it was shown that regimes of flame propagation multiplicity are also possible. In [4], various two-step reaction models for combustion wave propagation are considered and in particular, the model with parallel-competing schemes is investigated. The simple mechanistic explanation based on the notion of the crossover temperature and the balance between the rates of the reactions and heat release is presented, which allows one to qualitatively describe the existence of two different types of the response curves: monotonic with the single travelling wave solution and ‘S’-shaped with either one or three solutions. The AEA analysis is also carried out in [4]. In the most interesting case when the rates of both reactions are comparable, the authors expand the solution near the crossover temperature and obtain the leading order equations for the temperature and fuel concentration in the reaction zone. However, the resulting inner problem has not been solved either analytically or numerically, and the authors were not able to give any estimates for flame speed and temperature. In the case when one reaction is dominating, the problem is reduced to the usual one-step approximation which is similar to the conclusions in [3,5–8].

The linear stability analysis of the travelling combustion wave in the model with two parallel competing reactions is carried out in [9]. Considering only zero-order reactions with an infinitely thin reaction zone (delta-function sources are the simplification adopted in [9]) allows the dispersion relation to be obtained explicitly. The intermediate branch is again found to be unstable, whereas the fast and slow solution branches are either stable or exhibit pulsating instabilities, which agrees with the results communicated in [5].

In our earlier papers [10,11], we have studied a similar problem of flame propagation in a model with competitive exothermic–endothermic reactions. Both properties and flame stability have been investigated numerically and analytically by means of AEA. However, the behaviour of combustion waves in this model is significantly different as the second reaction works as a chemical buffer absorbing the heat produced in the course of the first exothermic reaction. The result is that the endothermic–exothermic model exhibits behaviour more like the one-step non-adiabatic model.

The aim of this paper is to undertake both linear and nonlinear stability analysis for the model of flame propagation with two-stage competing parallel reactions of the first order. Our attention is mainly focused on the parametric region, where multiple solutions exist. Under these conditions, interesting dynamical behaviour can occur. Further, we focus only on pulsating instabilities; therefore, the case of Lewis numbers for fuel greater than unity is always assumed. The paper is organized as follows. In §1, the governing model equations are introduced and the methods used for the investigation are described. Section 3 focuses on the travelling wave solutions to the model detailed in §2 and where regions of multiplicity may exist in the parameter space. Section 4 presents the actual results of the search for travel wave solutions and multiplicity while §5 details the solutions stability and bistability. The article concludes with a brief discussion of the results.

## 2. Model and methods

We consider a diffusional-thermal model with two-step kinetics for premixed combustion wave propagation in one spatial dimension under adiabatic conditions (table 1). It is assumed that the reactant undergoes two competitive exothermic reactions: (R1) and (R2) with Arrhenius kinetics. The non-dimensional equations governing this process can be found in [10] and are the following:
2.1where *t* and *x* are non-dimensional time and space coordinates, *u* and *v* are the dimensionless temperature and fuel concentration, *β* is the dimensionless activation energy of the first exothermic reaction (R1), *q* is the ratio of the enthalpies, *r* is the ratio of pre-exponential factors, *f* is the ratio of the activation energies of the second (R2) to the first (R1) exothermic reaction, respectively, and *L* is the Lewis number for the fuel (symbols are summarized in table 1). Here, *q* is positive in contrast to [10,11], where the second reaction is endothermic and thus *q* was negative.

Equations (2.1) are considered subject to the boundary conditions
2.2On the right boundary, we have a cold (*u*=0) and unburned state (*v*=1). The non-dimensionalized ambient temperature is taken to be equal to zero. The assumption of zero ambient temperature is assumed here to circumvent the cold boundary problem, which has been discussed by many authors (e.g. see the discussion in [12]). On the left boundary (), neither the temperature of the mixture nor the concentration of fuel can be specified. We only require that there is no reaction occurring so the solution reaches a steady state of (2.1). Therefore, the derivatives of *u* and *v* are set to zero for .

In this work, we consider travelling wave solutions to (2.1). We make the substitution *ξ*=*x*−*ct*, where *c* is the constant speed of the travelling wave, into (2.1) and produce the system of ordinary differential equations
2.3System (2.3) consists of two second-order ordinary differential equations and in this sense is similar to one-step models [12]. However, in contrast to one-step models, here there is no additional integral, which can be used to reduce the order of the system (2.3) or to determine the downstream flame temperature as it is done in [12]. Therefore, the burning temperature, is undefined and should be considered as an eigenvalue of the problem as well as the flame speed.

The system of ordinary differential equations (2.3) together with boundary conditions (2.2) constitute the two-point boundary value problem, which is solved numerically by using a standard shooting algorithm with a fourth-order Runge–Kutta integration scheme first, and then the results are corrected by employing the relaxation algorithm. Here, the moving frame coordinate *ξ* is discretized using finite differences. As the domain required by the travelling fronts varies greatly from one parameter set to the next, *ξ* is scaled to unity using an unknown parameter *δ*. After providing an initial guess for the solution profile the resulting nonlinear algebraic equations are solved iteratively using the Newton–Kantorovich method until the change from one iteration to the next is less than 10^{−15}. In addition to solving for *u* and *v*, *δ*_{ξ}=0 and *c*_{ξ}=0, for the scaling parameter *δ* and the front speed *c* are added to the system and solved.

The stability of the combustion waves is investigated in a way similar to [13]. We linearize the governing equations (2.1) near the travelling wave solution. We seek solutions of the form *u*(*x*,*t*)=*U*(*ξ*)+*ϵϕ*(*ξ*) e^{λt}, *v*(*x*,*t*)=*V* (*ξ*)+*ϵψ*(*ξ*) e^{λt}, where [*U*(*ξ*),*V* (*ξ*)] represent the travelling combustion wave, i.e. a solution to the problems (2.2) and (2.3). Here, terms proportional to the small parameter *ϵ* are the linear perturbation terms, *λ* is a spectral parameter governing the time evolution of the perturbation. Substituting this expansion into equation (2.1), leaving terms proportional to the first order of *ϵ* only, and introducing the vector function with components **v**(*ξ*)=[*ϕ*, *ψ*, *ϕ*_{ξ}, *ψ*_{ξ}]^{T} we obtain
2.4where
2.5Here, is a Wronskian of source terms in equations (2.1) calculated at *U*(*ξ*), *V* (*ξ*) and is 2×2 identity matrix

We will call a set, *Σ*, of all *λ* values for which there exists a solution to equation (2.4) bounded for both a spectrum of linear perturbations. In the general case, *Σ* is a set on the complex plane and it consists of the essential and the discrete spectrum. If there exists at least one *λ*∈*Σ* such that Re *λ*>0, then the travelling wave solution is linearly unstable, otherwise, if for all *λ*∈*Σ* the real parts are not positive, then the travelling wave solution is linearly stable. Therefore, in order to investigate linear stability of the travelling wave solutions to equation (2.1), the spectrum *Σ* of the problem (2.4) has to be found. It can be shown (see [14] for details) that the essential spectrum consists of parabolic curves in the complex plane with Re*λ*≤0. This implies that it is the discrete spectrum of the problem (2.4) that is responsible for the emergence of instabilities. The linear stability problem is solved by finding the location of the discrete spectrum on the complex plane using the Evans function method [15], as described in [13].

The nonlinear stability analysis as well as investigation of properties of the solutions bifurcating from the travelling waves are performed by direct integration of equations (2.1). In order to simulate equations (2.1), we employ the so-called method of lines [16]. The spatial dimension is discretized using the standard second-order centre finite difference scheme. The resulting system of ordinary differential equations are stiff. The stiffness index of the discretized system behaves like ∼2/(Δ*x*)^{2}, where Δ*x* is the distance between spatial points. Thus, increasing the resolution of the spatial dimension greatly increases the stiffness of the system to be integrated. Further, the associated Jacobian of the system is sparse. Therefore, we used the specialized stiff-solver for sparse problems mebdfso detailed in [17]. In all simulations presented here, the relative and absolute tolerance was fixed to 10^{−6} and controlled by the mebdfso.

## 3. Travelling waves

In this section, we focus on the investigation of the travelling wave solutions, which are found by numerically solving equations (2.3) subjected to boundary conditions (2.2) as described in §2. Firstly, we discuss the conditions for the existence of multiple solutions. Following [4], the crossover temperature is introduced as a temperature, *u**, at which the rates of the first and second reactions are equal. By means of simple algebraic manipulations with the second equation in (2.3), it can be shown that . This definition has physical significance to the present competitive reaction scheme only if *u** is positive. Otherwise, if *u** is negative, one of the reactions is always dominating, the model behaves effectively as a one-step reaction model and no multiplicity occurs. For *u**>0, there are two parameter regions: {*f*>1,*r*>1} or {*f*<1,*r*<1}, which can be mapped to each other by the change of notation of the first and second reaction. Hence, in what follows, we restrict the consideration to the case of {*f*>1,*r*>1}.

If the flame temperature is substantially below the crossover temperature then the first reaction is dominating and the flame speed and burning temperature are mainly defined by (R1). Thus, the reaction terms proportional to *r* can be omitted in equations (2.3), and we derive a standard one-step model [12,13] for which the flame speed and temperature can be written by using the AEA as
3.1

Similarly, if the flame temperature is much higher than *u**, then (R2) is dominating, and the flame speed and burning temperature are mainly defined by the second reaction step. In this case, the first reaction terms in equations (2.3) can be neglected and we once again derive the one step model, which after introducing new coordinate, , time, *t*′=*trq*/*f* and activation energy, *β*′=*βf*/*q*, can be written exactly in the same form as in [12,13]. Thus once again, we can use AEA formula (3.1) and obtain the flame speed and temperature, which in the original variables can be expressed as
3.2

Comparing the equations for the flame speed (3.1) and (3.2) and definition of the crossover temperature, we see that there are two conditions for the activation energies,
when the crossover temperature is equal to the burning temperatures for the one-step models governed by reactions (R1) and (R2), respectively. It is seen that depending on the value of *q*, the temperatures, *u*^{(1,2)}_{b}, and activation energies, *β*^{(1,2)}, can be ordered in different ways.

If *q*<1, then and *β*^{(2)}<*β*^{(1)}. In this case, as we increase *β* from the values when *u*^{(2)}≫*u** (or *β*≪*β*^{(1)}), the reaction (R2) is dominating, and the flame temperature follows the dependence in equation (3.2). As *β* approaches *β*^{(2)} the reaction (R1) becomes important. The flame temperature switches from the colder branch corresponding to (R2) to the hotter branch corresponding to (R1) as *β* is increased from *β*^{(2)} to *β*^{(1)} in such a way that *u*_{b} is staying close to *u**. For the values of the activation energy higher than *β*^{(1)}, the dependence *u*_{b}(*β*) approaches the limiting behaviour (3.1) corresponding to (R1). In other words, in the case when *q*<1, there is no multiplicity, and the flame temperature is monotonic function of activation energy with a characteristic plateau behaviour in the region *β*^{(2)}<*β*<*β*^{(1)} when both reactions are significant. This situation is also discussed in [4].

In the case *q*>1, the burning temperature, *u*^{(2)}, associated with (R2) is higher than *u*^{(1)}, and activation energies are also ordered differently, *β*^{(1)}<*β*^{(2)}. This situation is depicted in figure 1 for *L*=1, *f*=1.5, *q*=5 and *r*=25. The dependencies *u*^{(1,2)}(*β*) are plotted with the dashed lines, where index ‘1’ corresponds to the curve marked as ‘1-step slow branch’ and index ‘2’ to the curve labelled ‘1-step fast branch’, respectively. The dotted-dashed line corresponds to the crossover temperature, *u**, which crosses the curve *u*^{(1)}(*β*) at *β*^{(1)}≈6.3, while the intersection with *u*^{(2)} is for *β*^{(2)}≈32 and is not shown in figure 1. Numerical calculations of the temperature, *u*, for parameter values *L*=1, *f*=1.5, *q*=5, and *r*=25 are also plotted in figure 1 with the solid line. It is seen that as the activation energy is increased from the values below both *β*^{(1,2)} and the mixture is becoming less exothermic, the flame temperature follows formulae (3.2) for the one-step model governed by (R2). As the temperature decreases and approaches *u**, the first reaction comes into play, and the solution cannot follow the behaviour given by equations (3.2). This causes the appearance of the first turning point close to the crossover temperature. Further decrease in the burning temperature causes the solution to turn backward and to approach the slower and colder branch governed by reaction (R1), which results in folding and occurrence of the intermediate lag of the *u*(*β*) dependence. As *u*_{b} gets colder, the second reaction freezes and the combustion wave is controlled by the reaction (R1). This causes the second turning point after which the solution rapidly approaches the limiting behaviour defined by equations (3.1).

Thus, we see that the condition *q*>1 is necessary for the onset of multiplicity. However, this condition is not sufficient as it is demonstrated in figure 2, where the location of the turning points is plotted in the *β* versus *r* plane for *L*=1, *f*=1.5 and *q*=5. The solid lines correspond to numerical solutions, whereas the dashed lines represent *β*^{(1,2)}. The region of multiple solutions is enclosed by the solid curve. As seen in figure 2, the simple analytic approximation of the conditions for emergence of multiplicity substantially overestimates the region where triple solutions exist; however, it qualitatively correctly predicts the location of the turning points.

## 4. Results

The insights gained from the crossover temperature concept does away with the need to search a five-dimensional parameter space for regions of multiplicity. Provided *f* and *q* are greater than unity, we need only to find the value of *r* for which multiplicity exists. This is illustrated in figure 3. For a fixed value of *r*, system (2.3) is solved for different *β*, using the numerical methods mentioned earlier. For low activation energies, the flame front speed tends asymptotically to the values predicted by (3.2) indicated by the upper dashed line in figure 3. Likewise, for high activation energies, the flame front speeds are well described by (3.1), indicated by the lower dashed line. For *r*=3 (thin solid line) and *r*=7 (dotted-dashed line), the curve is monotonic, and, therefore, there is only one travelling wave solution for these parameters. At *r*=9, the curve (dotted) is ‘S’-shaped and there are in fact three solutions of system (2.3) for a thin range of *β*: the fast, slow and intermediate solutions. The crossover temperature analysis indicates that increasing *r* expands the range of *β* for which multiplicity exists. The remaining curve, for *r*=25, in figure 3 confirms this analysis and the width of the multiplicity region is approximately unity.

Increasing *f* (requiring more energy to ignite R2 relative to R1) has the effect of increasing the crossover temperature. This, in turn, means that any multiplicity regions will appear at smaller values of *β* as the fast branch reaction will extinguish earlier. Larger *f* will also have the effect of compressing the width of the multiplicity region. The crossover temperature does not depend on *L*, so theoretically multiplicity should exist for competing parallel reactions regardless of whether gases or solids are involved. The weak dependence of the location of the turning points on the Lewis number is also shown later in this paper in figures 5 and 7.

Although it is not shown here for brevity, the effect of variation of *q* on the existence of multiplicity was also investigated. It follows from equations (3.1) and (3.2) that decreasing *q* makes the width of the region of multiplicity smaller in terms of *β* and at a certain threshold value the turning points coalesce implying the disappearance of the multiplicity. These results were also confirmed by means of numerical calculations.

In figure 4, typical examples of flame front profiles are presented. The general shape of travelling waves in this model does not change greatly from one parameter set to another. What does change, by orders of magnitude, is the width of the reaction zone and the wavefront speed. Figure 4 presents solutions from a region of multiplicity with parameters *L*=2, *q*=5, *r*=25, *f*=1.6 and *β*=12. The upper plot is a profile from the slow branch with a speed *c*≈0.0018 and the scaling parameter, *δ*≈4000. By contrast, the lower plot has a profile from the fast branch and the values of *c* and *δ* differ from the slow branch by 2 orders of magnitude: *δ*≈25 and *c*≈0.46.

## 5. Stability analysis

In this section, we present the results of our investigation of flame stability. Having found where in the parameter space the multiple solution branches exist, it is important to know the stability of solutions within the region. Linear stability analysis of the travelling wave solutions of system (2.3) was conducted via the Evans function method. It is feasible to consider a large range of parameters with high resolution as only ODEs are involved. Full nonlinear stability analyses was then conducted for a few select parameter values and comparisons between the two stability approaches were made.

### (a) Linear stability

First, we consider the stability of the fast branch for parameter sets exhibiting multiplicity. The stability analysis indicates that the flame front is stable for small and moderate values of *β*. As *β* is increased closer to the turning point, the solution may lose stability. However, to investigate whether stability is lost before the turning point, we consider the existence of Andronov–Hopf (AH) bifurcations. Using the Evans function, the eigenvalues of the associated linear stability problem for a given parameter can be found. Combining the Evans function and the standard Newton–Raphson method the condition for an AH bifurcation, Re[*λ*(*β*)]=0, can be solved for *β*. The existence of an AH bifurcation at some *β*=*β*_{AH} means that the flame fronts have an oscillatory instability for *β*≥*β*_{AH}.

Two possible scenarios for the fast branch are illustrated in figure 5. For all the curves presented, *r*=25 and *q*=5. For three values of *f*, the position of the AH bifurcation is drawn for a given value of *β* and *L* and is indicated by the solid lines. The dotted lines to which the solid lines approach represent the position of the fast branch turning point for the respective value of *f*. We note here that *L* has only a weak influence on the value of *β* at the turning point *β*_{TP}. For *L*=2 and *f*=1.5, no AH point exists. This means that the fast branch is stable right up until *β*_{TP}. If, however, the Lewis number is increased, to say, *L*=10, then an AH point exists at *β*_{AH}=15.7 and all flame fronts from the AH points until *β*_{TP}=16.8 will have an oscillatory instability. The switching between the two scenarios occurs at the point where the AH and the turning point curves intersect. The point from which the AH locus originates and the fast solution branch changes from the completely stable to partly oscillatory unstable behaviour is the point of the Bogdanov–Takens bifurcation. For *f*=1.5, it is located at *L*=2.265.

The value of *f* has a strong influence on width of the range (*β*_{AH}, *β*_{TP}). Larger values of *f* drive the AH locus towards the turning point. For *f*=1.8 and *L*>8, the AH and the turning point locus are virtually identical. Further, increasing *f* has the effect of lifting the value of *L* for which the AH curve begins, i.e. the Bogdanov–Takens bifurcation point. Thus, we see the fast branch is stabilized completely for mixtures with low Lewis numbers if *f* is large enough. Here, we would like to note that in the single-step adiabatic models the pulsating instabilities also emerge via the Hopf bifurcation if the Lewis number for fuel is greater then unity and is taken to be sufficiently large [18]. However, in this case, the turning point and Bogdanov–Takens bifurcations do not exist. The instability behaviour for the upper branch is qualitatively closer to the single-step non-adiabatic model [19], i.e. there exists a turning point type flame quenching condition and the Hopf bifurcation critical curve originates from the Bogdanov–Takens point located on the turning point loci. In some sense, we can say that for the upper branch the second reaction is playing a role similar to the heat loss in the one-step model, i.e. it serves as a chemical sink for the fuel molecules providing the reaction path to burn fuel in less exothermic way.

The spectra of the associated linear stability problem for *f*=1.5, *L*=3 and 10 are shown against the flame front speed in figure 6. The values of Re(*λ*) are scaled over the flame speed in order to accommodate the range of eigenvalues variation on the same graph for flame speed changing two orders of magnitude. The spectrum is plotted against *c* so as to distinguish the fast to the mid branches. Starting from *c*=1 and moving leftwards to lower speeds, there is no curve drawn as Re(*λ*)≤0 for the fast branch indicating linear stability. An AH bifurcation (squares) occurs at *c*≈0.3. The real part of the spectrum becomes positive beyond the AH bifurcation and the eigenvalue is complex (dashed line) until the imaginary parts of the eigenvalues collide with each other and induce two, positive, purely real eigenvalues (solid lines). One real eigenvalue tends to be zero as solutions approach the turning point (diamonds), at *c*≈0.1, where the fast branch ends and the middle branch begins. The remainder of the curve for *L*=10 is purely real and positive indicating the entirety of the middle branch is unstable. For *L*=3, there is no AH bifurcation so the flame front is linearly stable for the entirety of the fast branch and stability is lost once the middle branch begins immediately after the turning point (*c*≈0.07).

A similar analysis was performed for the lower branch and the results are presented in figure 7. In addition to the AH loci (thick solid and dashed lines) and the slow branch turning point curves (vertical light dotted lines), there is a thick dotted line indicating the existence of AH bifurcations in the one-step model. The latter curve was obtained for the model which contains just the reaction R1. The slow solution branch is located to the right from the turning point curves. The stable solutions correspond to the parameter values below the AH curves. For a given *f*, if *β* is decreased to a certain critical value, the AH locus meets the turning point curve at the Bogdanov–Takens bifurcation point.

All stability curves considered here approach the one-step AH locus in the limit of large *β*, as in this limit the reaction R2 becomes frozen owing to the difference in activation energies. Increasing *f* has a similar effect, thus for *f*≥2.5 the AH curve is indistinguishable from the one-step one (as R2 is not activated).

If *f*≥2.5, the flame front is linearly stable for *β* values in the whole range between the turning point and AH locus. However, for lower *f* values (curves for *f*≤2.1 in figure 7), the AH curve folds over indicating that for a given *f* and *L* there can be several situations. If *L* is above the folding point of the *L*(*β*) curve for the AH bifurcation, then the slow solution is completely unstable right from the turning point. If *L* is between the maximum value of *L*(*β*) curve and the Bogdanov–Takens point, there can be two AH bifurcations. If *L* is less than the Bogdanov–Takens bifurcation point, the slow solution is stable in the range of *β* values between the turning point and AH curves.

Consider the curve for *f*=2.05 (thick solid curve in figure 7), if *L*=8, then no AH bifurcation exists and the entirety of the slow branch is unstable. However, if 4<*L*<7.5, then flame fronts are first unstable for *β*_{TP}≤*β*≤*β*_{AH1}, regain linear stability for *β*_{AH}≤*β*≤*β*_{AH2}, and become unstable again if *β*>*β*_{AH2}. Here, *β*_{AH1} and *β*_{AH2} are the first and second crossing of the AH curve with the line *L*=*constant* (for 4<*L*<7.5). Finally, if *L* is small enough, which for the case of *f*=2.05 is *L*<4, then the AH point exists for large *β* and the slow branch is linearly stable in a wide range of *β* values. In these last two scenarios, there exists the possibility of bistability.

The behaviour of the spectrum for the three scenarios of slow branch stability are presented in figure 8. Once again the eigenvalue is scaled over squared flame speed and plotted against *c* to show the mid and slow branches on the same graph. Three values of *L* are taken: *L*=1.5, 3.0 and 5.0. First, consider the curve for *L*=1.5 which is labelled ‘1’ in the figure 8. At *c*≈0.2, we see the turning point (diamonds) indicating the start of the middle branch. The eigenvalue is positive and purely real (solid curve) and reaches its maximum at *c*≈0.06 returning to zero at the second turning point denoting the beginning of the slow branch (*c*≈0.02). Thus, the entirety of the middle branch of solutions is unstable as expected. For *c*<0.02, there is no curve labelled ‘1’ in the figure so the real value of *λ* is negative, and the slow branch is linearly stable. Curve ‘2’, for which *L*=3.0, behaves in a similar fashion as curve ‘1’ except once we approach the second turning point near the end of the middle branch. At *c*≈0.025, the middle branch ends and a second purely real eigenvalue emerges from the turning point. The two eigenvalues collide soon after the start of the slow branch and produce a complex eigenvalue (dashed line), the real part of which goes to zero at the AH bifurcation (squares) near *c*=0.015. There is now a region for which curve ‘2’ is below the axis Re *λ*=0 and here the flame fronts are linearly stable. A second AH point occurs at *c*≈0.0045 and slow flame fronts have an oscillatory instability. Lastly, curve ‘3’ (*L*=5) shows that the middle branch is unstable, and once again a second real eigenvalue emerges from the turning point between the mid and slow branches. After the collision of the two real eigenvalues, the real part of the resulting complex eigenvalue grows. Thus, the slow branch is also unstable for *L*=5.

The linear stability analysis provided by the Evans function technique predicts the fast branch is either completely stable or loses its stability close to the turning point through an oscillatory instability. How close to the turning point the stability is lost strongly depends on the value of *f*. In cases of interest (i.e. moderate *f* values), the slow branch can be either completely stable, completely unstable, or has a region between two AH bifurcations that is stable. Thus, we see that the instability behaviour for the slow branch is far more complex than it is encountered in the single-step models.

### (b) Pulsating solutions

To confirm the stability predicted in the previous section, we performed a number of numerical simulations of the governing partial differential equations (PDEs) (2.1). In all cases, a narrow Gaussian was used as the initial condition for the temperature profile to mimic an igniting spark or ‘hot spot’. The fuel was initially set to unity for all *x*. During the simulation, the flame front speed was continuously calculated as
5.1where *L*^{′} is the length of the numerical calculation domain and the simulation is continued until the speed has settled to a clear constant value.

Comparisons were made between the results obtained by solving the system (2.3) and associated stability problem and the integration of the full PDE model (2.1). Figure 9 shows the flame front speed determined from equations (2.3), where the solid line corresponds to stable and the dashed line to unstable regimes according to the Evans function analysis. Only the fast branch is present and the turning point for the branch is at *β*≈16.8. The linear stability analysis predicts the fast branch to be stable until the AH bifurcation at *β*≈15.7. The solid dots give the flame front speed calculated during the simulations of the PDEs once transitory dynamics had dissipated. The PDE results are in good agreement with the ODE results up until the AH point. Pulsations in the wavefront speed are evident for simulations with *β* values beyond the AH bifurcation for which the instantaneous flame speed becomes a periodic function of time. In figure 9, this is indicated by the presence of two dots for a given value of *β* marking the local maximum and minimum of *c*. Period doubling of front speed occurs at *β*=16.2. For greater values of *β*, it was not clear what was happening and more simulations will be performed in future work to ascertain the dynamics. Period doublings were also observed earlier for single-step models [12] and led to the formation of the pulsating waves with complex spatio-temporal behaviour.

The flame front speed for *β*=16.2 from figure 9 as a time series is detailed in figure 10. The speed cycles through four distinct extrema values which are highlighted by the solid dots on the time series. The temperature profiles for exactly the times labelled with dots in figure 10 are shown in figure 11. Here, the focus is only on the reaction zone. Instead of a smooth propagation of the front with a flat temperature profile behind the reaction zone, the front proceeds in a ‘stuttering’ fashion, and the temperature profile has a peak near the reaction zone from which it decays and oscillates about some average temperature.

Similar simulations were performed for the slow branch with *L*=3, *f*=1.8, *r*=25 and *q*=5. Figure 12 details the comparison of results from ODE and PDE analysis. There are two AH bifurcations (squares) which occur at *β*=8.48 and *β*=9.96. For *β* values between the AH bifurcations, the flame front is propagating with a constant speed that is confirmed by simulations. The turning point for the fast branch is at *β*=9.88 (dotted line), which is stable everywhere (figure 5), so bistability exists for these parameter values. Pulsations in flame front speed are evident in simulations but decay slowly for stable fronts and grow in unstable fronts. Pulsations are most pronounced for *β* values close to the slow branch turning point at *β*≈8.25.

### (c) Bistability and switching

The ODE analysis demonstrated the existence of not only multiplicity of flame front speeds in the model but also bistability. This was confirmed by simulation of the full PDE model. Bistability allows the possibility of switching rapidly between the fast and slow reactions and vice versa. Figure 13 demonstrates the concept. Starting at a flame front with a speed such that it is near the first turning point (solid dot) the activation energy, *β*, is perturbed to a higher value. The fast reaction will now be extinguished leaving only a slow reaction flame front. The switching is indicated by the arrow in figure 13, and the final flame front speed is shown as the solid dot at the arrowhead. It may be possible to cause switching in a regime of combustion by a global change of the system such as a variation in pressure or ‘doping’ of the combustive material, so that it is no longer homogeneous and the inhomogeneity is enough to perturb the flame front.

We present two numerical simulations of system (2.1) which demonstrate that if such a perturbation could be made then the switching is stable in the sense that after the perturbation the flame front is able to settle to a constant speed. For the parameter values *L*=2, *q*=5, *r*=25 and *f*=3, the temperature and fuel profiles were found for *β*=3.047 (slow branch) and *β*=3.42 (fast branch). Both solutions are linearly stable. These profiles were then used as the initial conditions in the simulations—the only exception made to our usual Gaussian initial condition. Using the exact profile as the initial condition to a simulation avoids the sometimes lengthy process of the Gaussian spike reforming to the flame front solution for the given parameters. The profiles are propagated for 1000 time units and the value of *β* was changed: from *β*=3.047 to *β*=3.03 in the case of the switch from slow to fast branches; and from *β*=3.42 to *β*=3.43 in the case of the switch from fast branch to slow. What can be seen in figures 14 and 15 is the propagation of the initial condition for 0≤*t*≤1000, where the flame front speed is constant and agrees with the values calculated by ODE analysis in figure 13. At *t*=1000 (vertical dashed lines in both figures), the value of *β* is changed, and the flame front speeds go through a transitory phase of 500–1000 time units. The profiles have then fully reformed for the new conditions, and the speeds reach constant values. The new speed values are significantly different to the initial speeds and match the values of the opposite branch, as found by the ODE analysis.

## 6. Conclusion and discussion

In this paper, we numerically investigated a diffusional-thermal model with two-step competitive exothermic reactions for premixed combustion wave propagation in one spatial dimension under adiabatic conditions. The analysis was undertaken by using the shooting–relaxation technique and Evans function method applied to the ODE (2.3) and direct PDE integration methods. We studied the conditions for multiplicity of the travelling combustion waves and examined the stability of the solutions under the conditions of multiplicity.

A simple mechanistic explanation of the multiplicity phenomena was given based on the notion of the crossover temperature in the spirit of the analysis in [4]. This allows to qualitatively predict the regions in the parameter space where triple travelling combustion waves exist, which are then investigated numerically. The numerical analysis confirms that for certain parameter values the dependence of the flame on activation energy is triple-valued, ‘S’-shaped function, i.e. for given external conditions, there are three solutions travelling with different velocities: fast, mid and slow solutions, which are separated by two-turning point (or fold) bifurcations. The width of the multiplicity region in the parameter space is investigated and conditions are found for which the turning points merge and the triple solution is replaced with a single one.

The linear stability analysis was carried out and it was shown that the fast solution is either completely stable up to the turning point or exhibit oscillatory instability owing to the AH bifurcation prior to the turning point has been reached. In this case, two complex conjugate points of the discrete spectrum of the linear stability problem move from the left half to the right half of the complex plane causing the emergence of oscillatory instabilities. The neutral stability boundary and the turning point loci were found in the space of parameters. It was demonstrated that the oscillatory instability originates from the Bogdanov–Takens bifurcation—a point in the space of parameters where the turning point and AH loci merge. The mid branch was found to be completely unstable from one turning point to the other, which is characterized by the existence of purely real eigenvalue in the right half of the complex plane. As for the slow solution branch, the situation is more complex. The analysis of the location of the turning point and AH loci in the parameter space indicates that there are three scenarios resulting in the loss of stability as the activation energy is increased for other parameters being fixed: (i) the slow branch can be completely unstable from the turning point condition, (ii) there can be a single AH bifurcation so that the slow branch is stable for the range of activation energies between the turning point and the AH bifurcation conditions and loses stability in the oscillatory manner for activation energies higher than the critical values for the AH bifurcation, and (iii) there are two AH bifurcations for the slow branch so that the solution is unstable from the turning point to the first AH bifurcation, stable for activation energies between the two AH bifurcations and unstable again for the activation energies higher than the second AH bifurcation point. The critical conditions for the AH and turning point bifurcations was studied in the space of parameters and it was demonstrated that there exist conditions of bistability, i.e. a situation when both the fast and slow solutions are stable.

The nonlinear stability analysis was also carried out. It is demonstrated that the AH bifurcations observed were supercritical. A stable limit cycle emerges as the neutral stability boundary is crossed. In terms of the original problem, this implies that the stable pulsating solutions bifurcate from the travelling solution branches owing to the AH bifurcation. As the activation energy is varied away from the neutral stability boundary, the amplitude of oscillations of the instantaneous flame speed increase and at certain stage the period-doubling bifurcation occurs. It is not clear at the moment what kind of complex dynamical behaviour can be encountered as the bifurcation parameter is further increased and is the focus of future work. It was also demonstrated that the dependence of the flame speed on the activation energy can show the hysteresis type of behaviour. In particular decreasing (or increasing), the activation energy below (or above) the turning point conditions may result in sudden, very abrupt, acceleration (or deceleration) of the flame velocity so that the speed can rapidly change in orders of magnitude in a very short time interval. Such issues may be of practical interest, for example, in fire safety and will be studied in our future work.

As already mentioned in the introduction, the phenomenological model which is the subject of this paper is directly applicable to the modelling of a practical physical process like self-propagating high-temperature synthesis [20]. Therefore, it would be valuable if experimental activity similar to [1,2] could be revived in order to directly observe the dynamical patterns reported in here.

## Funding statements

V.V.G., A.V.K. and A.A.P. acknowledge the financial support from the Russian Foundation for Basic Research grant no. 11-01-00392 and 13-03-00282, and the Dynasty Foundation. H.S.S. acknowledges the support from the Australian Research Council grant no. DP0878146.

## Acknowledgements

I.N.T. acknowledges the hospitality and assistance of the P.N.L. Physical Institute, where much of this work was undertaken.

- Received May 14, 2013.
- Accepted July 24, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.