## Abstract

In this paper, we numerically investigate the stability of propagating combustion waves in the competitive exothermic–endothermic reaction model. The analysis is based on the Evans function method and direct numerical integration of the governing partial differential equations. The critical conditions for the onset of instability are found for a broad range of parameter values of the model. It is demonstrated that for the parameter values for which the combustion wave is unstable in the one-step reaction model, the inclusion of the endothermic step can lead to flame stabilization.

## 1. Background and introduction

Models of premixed combustion with two-step competitive reactions are usually designed to understand the phenomenology of decomposition or dissociation reaction of fuel. Significant attention has been paid to combustion systems with competitive exothermic reactions (e.g. [1] and references therein). By contrast, the competitive exothermic–endothermic reaction models have attracted less interest, despite the fact that these models can be appropriate for modelling pyrolysis [2] or decomposition [3] processes.

In this paper, we report results of a stability analysis of combustion waves in a model for the cases when the endothermic step plays a significant role in the reaction mechanism, and is shown to be able to limit or eliminate instabilities which may be inconvenient or even dangerous in practical applications. This is particularly relevant to the combustion of solid fuels, important in such applications as self-propagating high-temperature synthesis (SHS) [4], solid propellants [5], thermopower waves [6], etc. It is well known that flame propagation in solid fuel is prone to pulsating diffusional–thermal instabilities; in SHS the emergence of flame pulsations was demonstrated as early as [7]. Besides flame pulsations, rich dynamical wave patterns are also reported in literature on SHS [8–10]. As a rule, the emergence of pulsations is detrimental for products of combustion synthesis. Therefore, the control of the regime of flame propagation is an important issue; equally, there is certain demand to extend the limits of stationary combustion of solid fuels in propellant and thrust applications. The results of our investigation, reported fully in §3, suggest that it is possible to significantly modify the boundaries for the onset of oscillatory instability in certain combustion mixtures through the introduction of a competitive endothermic reaction, without substantial reduction of the rate of combustion process. We discuss these potential applications more fully in §4. The endothermic stage can act as an effective suppressant of flame oscillations [11]. This may offer a novel, and perhaps even a cost-effective, approach to undertake the SHS synthesis for conditions which are not accessible otherwise. An approach like this may be able to advance this technology significantly.

Besides the issues mentioned above where pulsations may have detrimental effects in product synthesis, any onset of pulsations in flames may also have safety implications. Oscillations of flame speed may result in flashback or blow-off phenomena while pulsations in temperature may lead to the emergence of superadiabatic temperatures and possible local overheating. It is demonstrated in [12] that as the conditions for the onset of pulsating regimes are approached the sequential stages of relatively slow, cool combustion and short bursts of a hot, fast-moving reaction front are observed. These bursts can be considered to be practically ‘unsafe’ regime of combustion even though such behaviour can be transient [13].

In earlier mathematical models for the propagation of deflagration in condensed material with competitive exothermic and endothermic reactions, e.g. [14], the endothermic reaction was treated as a perturbation and large activation energy asymptotic analysis was carried out to obtain the dependence of the burning velocity and temperature on the perturbation parameter and to identify the boundary of stability of combustion waves. It was demonstrated that the boundary of stability shifts towards larger values of the Zel'dovich number with the increase of the heat absorbed by the endothermic reaction. In [15], the case of general Lewis number of the mixture was numerically studied. As in [14] most of the results were obtained for the range of parameters, where the endothermic step is weakly activated and serves as a perturbation of the flame propagation model by the single-step endothermic reaction. It was found that the combustion wave solution is unique for given set of parameters. The stability was also investigated and it was demonstrated that flame pulsations may emerge for high Lewis number mixtures as the Zel'dovich number for the exothermic step is increased. The inclusion of compressible flow in the competitive exothermic–endothermic model has been considered in [16]. A weakly nonlinear approximation was developed for the combustion waves with burned temperature close to the ambient temperature. This allowed the construction of the appropriate travelling wave solution and was also used to demonstrate its stability.

In [17], the properties of the travelling combustion waves were studied both numerically and analytically when the assumption of weak endothermic reaction is relaxed. It was demonstrated that the flame speed, as a function of parameters, is either a single-valued monotonic function or a double-valued ‘C’-shaped function with the turning point type behaviour depending on the parameter values. It is expected that interesting instability characteristics can occur for these different regimes. In [15], we have analysed the stability of combustion waves in the regime where the exothermic reaction is dominating over the endothermic step; it is the purpose of the current investigation to carry out a similar systematic stability analysis of combustion waves for the cases when the endothermic step plays a more significant role in the reaction mechanism.

## 2. Model

Diffusional–thermal model for premixed combustion wave propagation in one spatial dimension under adiabatic conditions is considered. The kinetics is described by the two-step competitive reaction mechanism. Fuel, *A*, is assumed to decompose along two competitive paths: *t* and *x* are non-dimensional time and space coordinates, respectively; *u* and *v* are the dimensionless temperature and fuel concentration, respectively; *β* is the dimensionless activation energy of the exothermic reaction; *q* is the ratio of the enthalpies of the endothermic to exothermic reaction; *r* is the ratio of pre-exponential factors of the endothermic to exothermic reaction; *f* is the ratio of the activation energy of endothermic to exothermic reaction; *L* is the Lewis number for the fuel. The non-dimensional time and space coordinates used here, *t* and *x*, are related to the non-dimensional variables, *t*′ and *x*′ from [15] via the equations *t*′=*t* e^{β} and *x*′=*x* e^{β/2}. This additional scaling of time and coordinate is more convenient for the numerical analysis. As shown in [17], the adiabatic flame temperature for the model including the exothermic step only is equal to *β*^{−1}. In the non-dimensionalization adopted in [15,17], the rate of the exothermic reaction and flame speed are scaled with variation of *β* as e^{−β} and e^{−β/2}, respectively. The updated parametrization taken in the current work eliminates this strong exponential dependence on the activation energy, which is one of the key parameters in the model.

Equations (2.1) are considered subject to the boundary conditions
*u*=0) and unburned state (*v*=1). The non-dimensionalized ambient temperature is taken to be equal to zero. On the left boundary (

In order to obtain the travelling wave solutions, equations (2.1) are reduced to the system of ordinary differential equations using the ansatz *ξ*=*x*−*ct*, which introduces the coordinate in the frame co-moving with the wave. The resulting system of ordinary differential equations is solved numerically using the standard shooting-relaxation algorithms. The linear stability of the travelling combustion waves is investigated by using the Evans function method as described in [19]. The non-stationary solutions are obtained by direct integration of the partial differential equations (2.1), which are solved by the variant of the predictor–corrector scheme, namely, the method of splitting with respect to the physical processes. Initially (predictor) by assuming spatial uniformity, we solve the set of ordinary differential equations which describe temperature and fuel concentration variations by using the fourth-order Runge–Kutta algorithm. Secondly, as the corrector step, the equations of heat and mass transfer are solved with the Crank–Nicolson method of the second-order approximation in space and time using the predictor as a first iterate.

## 3. Stability analysis

The properties of the travelling wave solutions are analysed in [17] both numerically and asymptotically. To put our current work into context, we briefly summarize the main results from Gubernov *et al*. [17]. In this earlier work, it is found that the dependence of the flame speed, *c*, on parameters is either double-valued C-shaped function with a single turning point, or monotonic single-valued function. The regions of existence of different flame regimes in the *rq* versus *f* parameter plane are schematically illustrated in figure 1. As follows from the first equation in (2.1), the combination of parameters *rq* shows the characteristic ratio of the heat consumed and released through the endothermic and exothermic reactions, respectively. It is further referred to as the ‘ratio of heat effects’. In the case when *f*<1, the activation energy of the endothermic reaction is lower than the activation energy of exothermic reaction. The travelling wave solutions do not exist if *rq*>1 since the rate of heat consumption by the endothermic step exceeds the rate of heat release by the exothermic step for all values of *u*. When *rq*<1 and *f*<1, the propagation of deflagration is possible even though the activation temperature is lower for the endothermic reaction when compared with the exothermic reaction. This is possible, since the balance of rates of heat release and consumption is always in favour of the former. In this case, the dependence of the flame speed on the parameters is C-shaped, and there are either two solutions travelling with different velocities, or no travelling wave solutions exist. If *f*>1, there is always a single travelling wave solution and *c*, as function of parameters, is single-valued.

In figure 2, the results of the stability analysis for the case when *rq*<1 and *f*<1 is presented in the *L* versus *β* parameter plane for *r*=1, *q*=0.01 and several values of *f*=0.93, 0.95, 0.97 and 0.99. From figure 1, this choice of parameters represents the case when the dependence of *c* on parameters is double-valued. The dashed lines show the critical parameter values for the fold bifurcation. For given values *r*, *q* and *f*, there are two solution branches corresponding to different velocities for *L* and *β* located to the left from the dashed curve. The fast and slow solution branches merge as the parameter values reach the critical values for the fold bifurcation [17]. This is schematically shown in the left-bottom corner of figure 1. The travelling wave solutions cease to exist for values of *L* and *β* located to the right of the dashed curve. The solid lines represent the neutral stability boundaries for the emergence of flame pulsations due to the Andronov–Hopf bifurcation. This bifurcation causes the loss of stability of the travelling wave solution and emergence of pulsating waves. It is associated with crossing of the imaginary axis by a pair of complex conjugate eigenvalues of the linear stability problem, which is obtained through the linearization of system (2.1) around the travelling wave solution [20]. The slow solution branch is always unstable, this is due to the presence of the positive real eigenvalue of the linear stability problem which emerges at the fold bifurcation point. The fast solution branch is either stable for the parameter values below the solid curves in figure 2 or loses stability for *L* and *β* above the neutral stability boundaries for each value of *f*. The loci of parameter values for the Andronov–Hopf and fold bifurcations intersect as a result of the Bogdanov–Takens bifurcation. The Bogdanov–Takens bifurcation [20] is characterized here by the situation when a pair of complex conjugate eigenvalues of the linearized problem approach the origin along the imaginary axis, merge with the origin at the point of bifurcation and diverge along the real axis after that. This scenario is similar to the case of the one-step non-adiabatic model [21], where the detailed picture of the location of the eigenvalues of the linearized problem is given and their transformations associated with fold, Andronov–Hopf, and Bogdanov–Takens bifurcations are described. As the value of *f* is decreased, the amount of heat lost by the system due to the endothermic reaction grows and the critical parameter values for the fold, Andronov–Hopf and the Bogdanov–Takens bifurcation shift to the lower values of *β*. On the other hand, increasing the ratio of activation energies to *f*=0.99 results in the critical values of *β* for the fold and Bogdanov–Takens bifurcations to occur for values greater than *β*=12 and thus not shown in figure 2. As *f* tends to 1 the fold bifurcation disappears and *c* becomes a single-valued function of parameters (figures 1 and 2).

The numerical results of the stability analysis for the case of *f*>1 are presented in figure 3 in the (*β*, *L*) parameter plane for *r*=1. The critical parameter values for the emergence of Andronov–Hopf bifurcation are plotted for *f*=1.1, 1.5, 2 with the solid line for *q*=0.95, the dotted line for *q*=2.5, the dashed line for *q*=5. The dependence of the critical Lewis number on the dimensionless activation energy of the exothermic step is a monotonic function, which qualitatively agrees with the neutral stability boundary for the one-step adiabatic model [22]. It is shown in figure 3 that the neutral stability boundary strongly depends on the ratio of activation energies, *f*. Decreasing *f* to unity, i.e. to the limit of existence of travelling wave solutions has a strong destabilizing effect. The neutral stability boundary shifts towards smaller values of *β*. As *f* is increased from 1.1 to 1.5 the Andronov–Hopf locus moves to larger values of activation energy, *β*. Further variation of *f* from 1.5 to 2.0 has almost no effect on the stability boundary as illustrated in the inset of figure 3. For *f*=2.0, the curves corresponding to different values of *q* cannot be distinguished. It is worthwhile to note that for *f*≥2 the competitive exothermic–endothermic reaction reduces to the degenerate version of the one-step adiabatic model [17] and the Andronov–Hopf locus approaches the classical boundary of stability [19,22].

The location of the stability boundary in the (*β*, *L*) parameter space non-trivially depends on *q*, which is a ratio of the heat effects of the reactions. For *f*=1.1, the variation of *q* from 0.95 (the solid line) to 2.5 (the dotted line) and then to *q*=5 (the dashed line) shift the critical curve for the Andronov–Hopf bifurcation to smaller and then back to larger values of *β*. Thus, the increase in parameters *q* and *f*, which are specific for the competitive exothermic–endothermic model, can both stabilize and de-stabilize the combustion wave. We now focus on these parameters and how they affect the flame stability in greater detail.

The non-trivial behaviour of the system on the ratio of the heats of the endothermic to exothermic reactions is illustrated in figure 4, where the stability diagram is shown in the (*f*,*rq*) parameter plane for *L*=5, *β*=7 and *r*=1.0. The neutral stability boundary is plotted with the solid curve and is obtained via the Evans function method. The stable travelling waves are found for parameter values located to the right of the solid curve, whereas the combustion waves lose stability for the parameter values located to the left of the curve. It is seen in figure 4 that for certain parameter values such as *f*=1.17 shown with the dashed line, an increase in *q* leads to the transition of the combustion wave from stable to unstable and back to a stable state. For *f*≈1.155, the lower branch of the stability boundary crosses the line *rq*=0. As a result, when *f*≤1.155 there is a single transition from the unstable to stable behaviour as *q* is increased. If *f* exceeds the fold point value for the stability boundary *f*≈1.197, then any increase in *q* has no effect on the stability of the travelling combustion waves.

The results of direct numerical integration of (2.1) are also plotted in the inset of figure 4, where the dependence of the amplitude of oscillations *A* on *q* is displayed for *f*=1.17. Using the instant solution profile *v*(*ξ*) in a travelling frame co-moving with the combustion wave, we determine *ξ*_{1/2} as a coordinate such that *v*(*ξ*_{1/2})=1/2 and define *u*_{1/2}=*u*(*ξ*_{1/2}). Thus, *u*_{1/2} is the function of time only. The amplitude is defined as *q* is varied near the lower point on the stability boundary for *f*=1.17. As *q* is increased, the amplitude *A* grows according to the root law, which is typical for the Andronov–Hopf bifurcation. The curve *A*(*q*) can be used to obtain the critical points for the loss of stability from the results of direct integration of (2.1). They are shown in figure 4 with diamonds for *f*=1.17. It is seen that the correspondence of the data calculated with the direct integration of the governing PDEs (2.1) and with the Evans function method is reasonably good.

In figure 5, the results of stability analysis are presented in the plane of parameters *rq* versus *f* for *L*=5. The plane is split into three regions (as in figure 1): (i) {*rq*>1,*f*<1} is the shaded region with no travelling wave solutions; (ii) {*rq*<1,*f*<1} corresponds to parameter values for which two solution branches coexist; (iii) { *f*>1} is the region of parameters with single travelling wave solution. The critical parameter values for the Andronov–Hopf bifurcation are plotted with the solid lines numbered 1,2,…,6 for *β*=2, 5, 7, 8, 8.15, 8.3, respectively. For ease of interpretation of the data in figure 5, we also show the isocontours of critical flame speed at the Andronov–Hopf bifurcation, i.e. these curves satisfy *c*=const. and the condition for the Andronov–Hopf bifurcation simultaneously. Symbols ‘plus’, ‘diamond’, ‘asterisk’ and ‘cross’ connected with the dashed-dotted line represent the critical parameter values corresponding to the flame speed *c*=0.5, 0.1, 10^{−2} and 10^{−3}, respectively.

In [17], it is demonstrated that for *rq*>1 the velocity of combustion wave vanishes as *rq* versus *f* plane, *c* decreases significantly. This situation is observed in figure 5 for the stability boundaries shown with the curves 1, 2 and 3. We note that the curves are discontinued when the value of the flame speed, *c*, drops below 10^{−4}. The stable combustion waves correspond to parameter values located to the right of these curves for each value of *β*=2, 5 and 7. The neutral stability boundaries 1, 2 and 3 show similar behaviour. They all have a turning point where the maximum value of *f* is reached.

The curve 4 in figure 5 is plotted for *β*=8, which is close to the critical value *β*≈8.082 which corresponds to the loss of stability in the one-step adiabatic model for *L*=5. In this case, the neutral stability boundary is almost a straight line in the (*f*, *rq*) parameter plane which separates the parameters corresponding to the stable and unstable combustion wave solutions (to the right and the left sides of the curve 4, respectively). The parabolic-shaped curves 5 and 6 represent the Andronov–Hopf loci for *β*=8.15 and 8.3, respectively, for which the travelling wave solutions are unstable in the one-step adiabatic model. The stable combustion waves correspond to the parameter values located above the curves 5 and 6. Thus, it is possible to stabilize the flame propagation for the parameter values where pure exothermic reaction leads to the emergence of pulsations by introducing the competitive endothermic reaction step. To the right from the local minima of the parabolic curves 5 and 6, the value of *f* becomes sufficiently large. It is shown in [17] that with the increase of *f*, system (2.1) tends to the one-step adiabatic reaction scheme, therefore the propagation of flame should become unstable at certain value of *f* and this explains the emergence of the right part of the stability boundary. On the other hand, in the left part of the curves 5 and 6 the flame speed and temperature decrease which then results in the emergence of pulsating instabilities. A key result of the analysis shown in figure 5 is that a combustion wave can be stabilized without a drastic reduction of the flame speed. This may have important implications for physical processes and applications.

In figure 6, the stability diagram in the (*f*, *rq*) parameter plane is plotted for *L*=5, *β*=8.15 and different values of *r*=0.5, 1 and 2 with the dash-dotted, solid and dashed lines, respectively. The stabilization of flame exists for different values of *r*, and it is not a special case for the specific parameter value *r*=1.

It can be seen that the Andronov–Hopf bifurcation loci for different values of *r* are close to each other in the (*f*,*rq*) parameter plane (as evident from the three curves shown in figure 6). Therefore, the critical parameter values for the loss of stability are mostly governed by the combination of parameters *rq*, which represents the characteristic scale for the rate of heat consumption in the endothermic step. The results of the direct numerical integration of the governing equations (2.1) are also shown in figure 6 with the bold circles and ‘x’ symbols for *r*=1 and *q*=4. The system of equations (2.1) is integrated for a sufficiently long time to remove the transient behaviour until either stable travelling wave (bold circles) or periodic in time pulsating solutions (denoted by ‘x’) are obtained, which are presented in figure 6 for different *f*. It is seen that the results obtained using the Evans function and the direct numerical integration of the governing equations (2.1) are in good agreement.

The travelling and pulsating solutions are illustrated in figure 7 for *L*=5, *β*=8.15, *r*=1, *q*=4 and *f*=1.25, *f*=1.35, *f*=1.4. The time history *u*_{1/2}(*t*) is shown in figure 7*a*, whereas the dependence of *u*_{1/2}(*t*) on *ξ*_{1/2}(*t*) is plotted in figure 7*b*. The constant value of *u*_{1/2} in figure 7*a* represents the travelling wave solution for *q*=1.4. This choice of parameters is marked as ‘1’ in figure 6 and in figure 7*b* it appears as a dot. For pulsating solutions at *f*=1.35 and 1.25, the dependences *u*_{1/2}(*t*) are periodic oscillating functions as seen in figure 7*a*. In figure 6, these parameter values are numbered as ‘2’ and ‘3’, respectively. The travelling combustion waves become unstable for these choices of parameters and pulsating waves emerge due to the Andronov–Hopf bifurcation. As *f* is decreased from the boundary of stability, the amplitude and period of pulsations grow. In figure 7*b*, the pulsating solutions are represented with trajectories forming closed loops or limit cycles in the *u*_{1/2} versus *ξ*_{1/2} plane of parameters.

## 4. Conclusion and discussion

In this paper, the stability of combustion waves has been investigated in the model with competitive exothermic–endothermic reactions, and neutral stability boundaries have been found for a wide range of parameter values. Previously, in [17], we demonstrated that if the activation energy of the endothermic reaction is less than the activation energy of the exothermic reaction the dependence of the flame speed on the parameters is a double-valued function. In this paper, for the first time, the stability of deflagration in this regime has been investigated. We have shown that the slow solution branch is always unstable, whereas the fast branch is either stable or loses stability through an Andronov–Hopf bifurcation. The switching between these two regimes is due to the Bogdanov–Takens bifurcation. This scenario for the onset of instability is qualitatively similar to the scenario for the one-step *nonadiabatic* model.

For the case when the activation energy of the endothermic step is greater than the activation energy of the exothermic step, the dependence of the flame speed on the parameters becomes a monotonic single-valued function. In terms of parameters of the one-step reaction model, i.e. the Lewis and Zel'dovich numbers, the neutral stability boundary is found to be a monotonic curve such that the Lewis number tends to infinity as the Zel'dovich number tends to one, and the Lewis number tends to one as the Zel'dovich number becomes infinitely large. The transition to instability is again found to be caused by an Andronov–Hopf bifurcation. This behaviour is now qualitatively similar to the scenario of the onset of instability in the one-step *adiabatic* model.

Despite the qualitative similarities, it is important to note that the usual analogy with the one-step adiabatic and non-adiabatic models does not always work in the case of the competitive endothermic–exothermic model; the critical conditions for the transition to instability can depend nontrivially on the specific parameters for the exothermic–endothermic model, i.e. to the ratio of activation energies, the heats of the endothermic and exothermic reactions. The endothermic reaction absorbs a certain amount of heat, which is released in the exothermic step. However, this occurs locally in the thin reaction zone. By contrast, the conductive or radiative heat loss mechanism has a definite heat flux along the whole product zone, which cools down to the ambient temperature. As a result, increasing the rate of heat consumption of the endothermic reaction does not necessarily lead to the loss of stability. In fact, as is shown here, this can even cause the stabilization of the propagation of combustion waves. This property may have important applications such as flame inhibition by endothermic reactions, or more general frontal combustion in heterogeneous systems [11]. Our results are also of considerable relevance to recent developments in thermochemically coupled SHS reactions [23,24].

To conclude, we summarize the main results obtained in this investigation. For parameter values for which the combustion wave is unstable in the one-step reaction model, the flame propagation can be stabilized by the inclusion of the competitive endothermic reaction. This is achievable without a significant reduction in the flame speed, and, by allowing the steady propagation of deflagration in mixtures which would otherwise exhibit flame pulsations or quenching, has potential importance for applications such as SHS or safe and efficient use of solid-state fuels.

## Data accessibility

This work does not have any experimental data. All figures are available as Electronic supplementary material.

## Authors' contributions

All authors participated in formulation of the mathematical model and wrote the paper. V.V.G. and A.V.K. did most of the numerical simulations. All authors gave final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

V.V.G. acknowledges the financial support from the Russian Foundation for Basic Research grant no. 13-03-00282. H.S.S. acknowledges the support of the Australian Research Council grant no. DP0878146.

- Received May 5, 2015.
- Accepted July 1, 2015.

- © 2015 The Author(s)

Published by the Royal Society. All rights reserved.