## Abstract

The propagation of planar combustion waves in an adiabatic model with two-step chain-branching reaction mechanism is investigated. The travelling combustion wave becomes unstable with respect to pulsating perturbations as the critical parameter values for the Hopf bifurcation are crossed in the parameter space. The Hopf bifurcation is demonstrated to be of a supercritical nature and it gives rise to periodic pulsating combustion waves as the neutral stability boundary is crossed. The increase of the ambient temperature is found to have a stabilizing effect on the propagation of the combustion waves. However, it does not qualitatively change the behaviour of the travelling combustion waves. Further increase of the bifurcation parameter leads to the period-doubling bifurcation cascade and a chaotic regime of combustion wave propagation. The chaotic regime has a transient nature and the combustion wave extinguishes when the bifurcation parameter becomes sufficiently large. For Lewis numbers of fuel close to unity, the parameter regions where pulsating solutions exist become very close to each other and this makes it difficult to experimentally observe the period-doubling. It is shown that the average velocity of pulsating waves is less than the speed of the travelling wave for the same parameter values.

## 1. Introduction

Pulsating regimes of flame propagation in premixed solid combustion have been observed as early as 1950 (Belyaev & Komkova 1950), in experiments with chrome-magnesium thermite. In that work, the influence of pressure variation of the inert dilutant upon the speed of the combustion wave was examined for various thermites. It was found that for chrome-magnesium thermite, the increase of pressure above a certain critical value resulted in pulsating behaviour of the flame propagation with the alternating intervals of fast combustion and long depression. A tendency towards flame extinction for pulsating regimes was also reported. However, no explanation of this phenomenon has been given and it was treated as an unusual type of flame propagation. A regular description of pulsating combustion waves was given by Maksimov (1963) in an experimental study of nitroglycerine propellant combustion. It was determined that pulsations were of a truly auto-oscillatory nature, which was related to the mechanism of combustion only. In that paper, the pressure of the inert component was used as a bifurcation parameter and it was found that the increase of this parameter could lead to complex regimes of flame oscillations with irregular time intervals between successive temperature peaks. In the paper by Shirko & Nersesyan (1978) on solid premixed combustion of tantalum–carbon mixture, experimental observation of a period-doubling bifurcation was reported. The dilution of the combustible mixture with an inert dilutant resulted in the onset of pulsations, which were determined as almost small harmonic oscillations of the flame speed. As the dilution ratio increased, the oscillations became more relaxational and at a certain stage the period two solution emerged. It was also observed that more complex oscillatory regimes occurred as the dilution rate was further increased. Pulsations are often observed in solid combustion in experiments on self-propagating high-temperature synthesis (Merzhanov & Rumanov 1999). By contrast, pulsating waves are rare in experiments with combustible gas mixtures. One example was given by Pearlman & Ronney (1994), where radial pulsating waves with a complex time behaviour were observed for the high Lewis number premixed gas combustion, which were attributed exclusively to thermal-diffusive instabilities.

Analytical description of the emergence of pulsating planar combustion waves from the travelling planar combustion waves was formulated in terms of single-step diffusional-thermal models. In numerical investigations of gasless solid combustion, Shkadinskii *et al.* (1971) demonstrated that the pulsating regime of flame propagation existed in a model with single-step kinetics of the reaction. Based on numerical data, the phenomenological approximation of the neutral stability boundary was found and it was also shown that as the bifurcation parameter increased, the dynamics of pulsations became more and more complicated. In Aldushin *et al.* (1973), the numerical investigation of the solid gasless combustion was carried out for various kinetic laws of Arrhenius type that modelled the effect of heterogeneous condensed media on the flame propagation regimes. The authors clearly demonstrated period-doubling bifurcations in the dynamics of flame pulsations. The pulsating waves of various periods were reported to be followed by more sophisticated regimes of flame oscillations, which could not be traced accurately as the bifurcation parameter was increased. It was also noted in Aldushin *et al.* (1973) that the number of period-doubling bifurcations observed before the onset of complex pulsations depended upon the specific kinetic law: for the case, where there was strong effect of heterogeneity, only period two solutions were observed, whereas for the case of weak heterogeneity, solutions with periods higher than two were found. Further advances in the understanding of flame oscillations were related to activation energy asymptotics, which was one of the cornerstones of combustion theory. Using this approach, the properties of pulsating waves in both solid (Matkowsky & Sivashinsky 1978) and gaseous (Joulin & Clavin 1979; Matkowsky & Olagunju 1980) premixed combustion were studied in detail based on diffusional-thermal models with one-step kinetics and the Arrhenius type temperature-dependence of the rate of the reaction. The linear stability analysis was carried out for solid adiabatic flames (Matkowsky & Sivashinsky 1978) and for both gaseous adiabatic flames (Sivashinsky 1977; Joulin & Clavin 1979; Matkowsky & Olagunju 1980) and gaseous non-adiabatic flames (Joulin & Clavin 1979) under a condition close to equivalence. It was demonstrated that the onset of pulsating instabilities was because of the Hopf bifurcation and the condition for the stability loss of the travelling planar waves as well as the dispersion relation for the pulsating instabilities were established in terms of the control parameters. In Matkowsky & Sivashinsky (1978) and Matkowsky & Olagunju (1980), the nonlinear bifurcation analysis of the pulsating solutions was also carried out and it was shown that the Hopf bifurcation had a supercritical character. The amplitude, frequency and speed of the pulsating combustion wave were estimated near the point of bifurcation in the parameter space. The pulsating wave was demonstrated to have a mean velocity less than the travelling wave speed. These results were also reviewed by Margolis & Matkowsky (1983). Although the activation energy asymptotic approach has proved to be successful for the analysis of the primary bifurcation of pulsating waves from the travelling combustion waves, investigation of period-doubling bifurcations usually can be done only numerically.

In Weber *et al.* (1997), the propagation of combustion waves in solid and gaseous mixtures was studied numerically for an adiabatic model with one-step kinetics. It was demonstrated that for the case of Lewis number equal to unity (gaseous mixture), the pulsating waves did not exist, whereas for solids the pulsating solutions of period 1, 2 and more complex dynamics were observed. In Mercer *et al.* (1998) and McIntosh *et al.* (2004), non-adiabatic combustion waves in a solid fuel were investigated. The emergence of pulsating waves owing to the Hopf bifurcation was found and the period-doubling bifurcations giving rise to pulsating solutions of period 2 and 4 were determined in the parameter space (Mercer *et al.* 1998). The authors speculated that further increase of the bifurcation parameter could lead to a period-doubling route to chaos. However, as the bifurcation parameter became large the flame extinguished. It was not clear if the flame quenching was attained through the chaotic pulsations or otherwise. The fact that the pulsating behaviour was observed for solids and did not emerge in gases at least under the condition of equivalence was explained in Gubernov *et al.* (2003, 2004), where it was shown that the onset of pulsating instabilities was due to the Bogdanov–Takens bifurcation from which the Hopf locus originated in the parameter space. For a given value of the heat loss parameter, the Bogdanov–Takens bifurcation was located at the definite Lewis number greater than unity. Subsequently for Lewis numbers above this critical value, the pulsating instabilities emerged with an increase in the activation energy. On the other hand, for Lewis numbers below the critical value corresponding to Bogdanov–Takens bifurcation no pulsating behaviour occurred. A similar scenario has been observed recently in diffusion flames (Gubernov & Kim 2006). This result correlates with the analysis in Aldushin *et al.* (1973). In Bayliss & Matkowsky (1990) two adiabatic models of solid premixed combustion were studied numerically. One of the models allowed melting of the solid component while the other one did not. In both models, the increase of the activation energy led to the onset of pulsations owing to the Hopf bifurcation. For the model without melting, the sequence of period-doubling bifurcations was shown to occur as the bifurcation parameter became larger. The pulsating solutions of period 2, 4 and 8 were found. A period-doubling route to chaos was suggested with the increasing values of bifurcation parameter. However, for the model with melting of the solid fuel only bifurcation to period 2 solution was observed, the period 1 and period 2 solutions alternated as the bifurcation parameter was increased and a transition to chaos via intermittency was predicted as the bifurcation parameter became sufficiently large. In a number of papers (Brailovsky & Sivashinsky 1993; Frankel *et al.* 1994, 2000), free-boundary models were considered in order to investigate thermal instabilities in solid combustion. Depending on the kinetic functions used in these systems two scenarios of stability loss with respect to pulsating perturbations were observed. The travelling combustion wave either became unstable owing to the Hopf bifurcation followed by a period-doubling cascade leading to a chaotic solution, or the stability loss was due to the infinite period bifurcation of the Shilnikov type. In Frankel *et al.* (2000), both sub- and supercritical Hopf bifurcations and multiple cascades of period-doubling were reported to occur.

Although the literature overview presented above is far from being complete, it demonstrates that significant success has been achieved recently in understanding the reasons for the onset of planar pulsating regimes of combustion waves propagation on the basis of one-step diffusional-thermal models. The aim of our current research is to investigate how the scenarios of the stability loss and the onset of pulsations are projected onto the dynamical behaviour of the two-step model with chain-branching reaction mechanism and first-order recombination reaction introduced recently by Dold (2007). In Dold (2007), the properties and stability of the travelling combustion waves were studied using activation energy asymptotics. In our recent papers (Gubernov *et al.* 2008*a*,*b*, 2009), we numerically investigated the properties of travelling combustion waves and their stability with respect to pulsating perturbations, and compared the results with the predictions of the corresponding one-step models and the results presented by Dold (2007). We summarize our results from Gubernov *et al.* (2008*a*,*b*, 2009) in §3 of the current paper. In §2, the model equations are introduced and in §4 we perform the analysis of pulsating solutions emergence and investigate the bifurcations leading to stability exchange exhibited by pulsating waves. In §5, we study the emergence of a chaotic regime of flame propagation owing to the Feigenbaum’s period-doubling cascade and demonstrate that it has a chaotic transient nature, which results in abrupt extinction of combustion waves. In §6, the results are summarized and discussed, including prospects of future work.

## 2. Model

We consider an adiabatic model for premixed flame propagating in one spatial dimension that includes two steps: autocatalytic chain branching and recombination . Here *A* is the fuel, *B* is the radical, *C* is the product and *M* is a third body. It is assumed that all the heat of the reaction is released during the recombination stage and the chain-branching stage does not produce or consume any heat. Following Gubernov *et al.* (2008*a*), the governing equations for the non-dimensional temperature *u*, concentration of fuel *v* and radicals *w*, can be written in a non-dimensional form as
2.1
where *x* and *t* are the dimensionless spatial coordinate and time, respectively, *L*_{A} and *L*_{B} are the Lewis numbers for fuel and radicals, respectively, *β* is the dimensionless activation energy of the chain-branching step (which corresponds to the definition for the one-step model (Gubernov *et al.* 2004), *r* is the ratio of the characteristic time of the recombination and branching steps (which cannot be reproduced in one-step approximations of the flame kinetics).

Equations (2.1) are considered subject to the boundary conditions
2.2
On the right boundary (fresh mixture) we have cold (*u*=*u*_{a}) and unburned state (*v*=1), where the fuel has not been consumed yet and no radicals have been produced (*w*=0). In contrast to one-step models, the ‘cold boundary problem’ does not exist for the current model. The ambient temperature is allowed to take any value as long as the branching reaction is not activated in the fresh mixture, i.e. the reaction term in the third equation in equation (2.1) is negative. More formally this implies that *r*>*e*^{−1/ua} or the ambient temperature is below the temperature at which the fresh mixture becomes self-ignitable. 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 equation (2.1). Therefore, the derivatives of *u*, *v* are set to zero and *w*=0 for .

## 3. Properties of the travelling wave solution

The solution to the problem (2.1) and (2.2) is sought in the form of a travelling wave *u*(*x*,*t*)=*u*(*ξ*), *v*(*x*,*t*)=*v*(*ξ*) and *w*(*x*,*t*)=*w*(*ξ*), where a coordinate in the moving frame, *ξ*=*x*−*ct*, is introduced and *c* is the speed of the travelling wave. Substituting the solution of this form into the governing equations, we obtain
3.1
The boundary conditions (2.2) can be modified if we multiply the first equation in equation (3.1) by *β*, add it to the second and third parts of equation (3.1) and integrate it once over *ξ* from minus to plus infinity. This yields a condition: , where *S*=*βu*+*v*+*w*. Combining this condition with equation (2.2) results in
3.2
where is the adiabatic burned mixture temperature and *σ* denotes the residual amount of fuel left behind the wave and is unknown until a solution is obtained.

In order for the travelling wave solution of equation (3.1) to exist, several conditions have to be satisfied. The reactions in front of the wave, in the preheat region and behind the wave, in the product region, should be in chemical equilibrium. The latter restriction implies that the radicals cannot be produced in these regions effectively. In other words, the production of radicals is overwhelmed by their consumption of radicals, which are governed by the second and the third terms in the right-hand side of the last equation in equation (3.1). In front of the flame, the temperature of the fresh mixture is limited from above by *u*_{i} as discussed earlier. In the product region, the reaction must be completed; therefore, the branching term is less than the recombination term in the third equation of equation (2.1), that is . This condition was first derived in Gubernov *et al.* (2008*a*) for the adiabatic case and defines the region in the parameter space, where the travelling combustion waves exist. Once this condition is violated the adiabatic combustion waves exhibit extinction.

The properties of the travelling combustion waves were investigated analytically by Dold (2007), Gubernov *et al.* (2008*a*) and numerically by solving the set of equations (3.1) subject to boundary conditions (3.2) by Gubernov *et al.* (2008*a*, 2009). It was found that varying the Lewis number for fuel, *L*_{A}, has a substantial effect on the properties of the premixed flames, whereas the variation of Lewis number corresponding to the radicals, *L*_{B}, affects only quantitative behaviour of the combustion waves. For the case of Lewis number for fuel less than unity (*L*_{A}<1), the dependence *c*(*β*) is a monotonically decaying function exhibiting extinction as the flame speed reaches zero at a certain value of the activation energy, *β*_{e}, corresponding to extinction. For parameter values sufficiently away from the extinction condition, the residual amount of fuel can be neglected. Almost all fuel is converted to radicals and no fuel leakage is observed (). As we increase *β* and approach the extinction point, the value of *σ* becomes significant. At the extinction condition the residual amount of fuel reaches its maximum value. The dependence of the flame speed, *c*, on *β* becomes C-shaped for *L*_{A}>1, i.e. *c*(*β*) is a double-valued function. There are either two solutions travelling with different flame speeds or the solutions cease to exist owing to the extinction when the fast solution branch meets the slow solution branch at the turning point of the *c*(*β*) curve. We denote the coordinates of the turning point with subscript ‘tp’, i.e. *β*_{tp} and *c*_{tp}=*c*(*β*_{tp}). For small values of *β*, the fast solution branch is characterized by a negligible amount of fuel. As the activation energy is increased, the residual amount of fuel grows and it becomes significant as the turning point of *c*(*β*) is approached. The slow solution branch is characterized by a considerable fuel leakage. As we move along the slow solution branch by decreasing *β* from the turning point value, the flame speed decreases and at a certain value of *β*_{e} it becomes equal to zero. As shown in Gubernov *et al.* (2008*b*), variation of *L*_{B} over two orders of magnitude does not affect the qualitative behaviour of the solution branches in the parameter space. The behaviour of the travelling combustion waves described above is also illustrated in figure 1, where the dependence of the wave speed, *c*, is plotted against the dimensionless activation energy, *β*, for two values of the Lewis number for fuel *L*_{A}=1 and *L*_{A}=10 as shown in figure 1. The value of the Lewis number for radicals is fixed as *L*_{B}=1 throughout this paper, whereas *r* is taken to be equal to 10^{−4}. For the case *L*_{A}=1, the dependence of *c* on *β* is a monotonically decaying function. Once *L*_{A} becomes greater than one, the flame speed *c* turns into a double-valued function of *β*. The effect of the ambient temperature variation is also investigated. In figure 1, the dashed lines represent the solutions obtained for *u*_{a}=0 and the solid lines show *c*(*β*) for *u*_{a}=0.01. As in the case of the one-step kinetics studied in Gubernov *et al.* (2005), as the fresh mixture is preheated to higher initial temperatures, the domain of existence of the travelling combustion waves becomes larger. In other words, the increase of the ambient temperature results in the growth of both the flame speed and the *β*_{tp} values, expanding the domain of *β* values for which the travelling combustion waves exist. At the same time, the *u*_{a} variation does not affect the *c* on *β* dependence qualitatively.

The linear stability of the travelling combustion waves is investigated in Gubernov *et al.* (2009) in detail for the case *u*_{a}=0. Here, we extend the previously obtained results to include the effect of non-zero ambient temperature on the flame stability. The governing equations (2.1) are linearized near the travelling wave solution in order to obtain the linear stability problem, which describes the evolution of the infinitely small perturbations of the travelling wave. We seek the solution of the form , where ** u**=[

*u*,

*v*,

*w*]

^{T}is the solution in vector form,

**(**

*U**ξ*)=[

*U*(

*ξ*),

*V*(

*ξ*),

*W*(

*ξ*)]

^{T}represents the travelling combustion wave and terms proportional to the small parameter

*ϵ*are the linear perturbation terms. Substituting this expansion into equation (2.1), leaving terms proportional to the first order of

*ϵ*only, we obtain 3.3 where is the linear differential operator and we have introduced , a 3×3 identity matrix, 3.4 the reaction terms are defined as . In equation (3.4), the derivative denotes the Jacobi matrix evaluated at

**(**

*u**ξ*,

*t*)=

**(**

*U**ξ*). The linear stability problem, equation (3.3) subject to boundary conditions as , is solved by finding the location of the discrete spectrum of in a complex plane using the Evans function method (Evans 1972) implemented using a compound matrix approach as discussed in Gubernov

*et al.*(2009). Since the technique of the stability analysis that is used in this paper is essentially identical to the analysis in Gubernov

*et al.*(2009), we skip it here to avoid the repetition of our earlier paper and proceed to the description of the new results.

For the case *L*_{A}>1, the travelling combustion wave is either stable or exhibits the onset of pulsations as a result of the Hopf bifurcation. The results of this study are summarized in figure 2 where the critical parameter values for the extinction and Hopf bifurcation are depicted with the solid and the dotted lines, respectively, on the *L*_{A} versus *β* plane for *L*_{B}=1 in figure 2*a* and on *β* versus *u*_{a} plane for *L*_{A}=10 in figure 2*b*. In figure 2*a*, the Lewis number for fuel is plotted in decimal logarithmic scale in order to map the range of its variation uniformly on the graph. The term extinction here means the case when the turning point of the *c*(*β*) curve is reached for *L*_{A}>1. In figure 2*a*, the stability diagram is plotted for two values of the ambient temperature *u*_{a}=0 and *u*_{a}=0.01. For each value of *u*_{a} the region to the right from the solid line corresponds to the parameter values, where the travelling wave solutions do not exist. The parameter values to the left from the extinction curve correspond to the region where two travelling wave solution branches may coexist *L*_{A}>1. The stability analysis shows that in this parameter region, the slow solution branch is always unstable and the fast solution branch is either stable or exhibits the onset of pulsating instabilities via the Hopf bifurcation. The dotted line corresponds to the Hopf bifurcation locus. Between the dotted (Hopf) and the solid line (extinction) lies a region where the fast travelling solution branch is unstable. Once the dotted line is crossed in the parameter space, for example, by increasing *β* for fixed *L*_{A}, the fast branch of the travelling wave solutions becomes unstable with respect to pulsating instabilities. The point where the Hopf bifurcation curve meets the extinction curve is a bifurcation of co-dimension two and is known as the Bogdanov–Takens bifurcation. There are two zero eigenvalues of the linearized stability problem (3.3), and a Hopf bifurcation and a saddle-node bifurcation meet at this point. It is also a point from which the Hopf bifurcation locus in the parameter space originates and therefore this bifurcation is directly responsible for the onset of pulsations in the model. As seen in figure 2*a*, the increase of *u*_{a} shifts both the boundary of the existence of travelling solutions and the Hopf loci towards higher values of the dimensionless activation energy. In other words, the ambient temperature rise has a stabilizing effect on the propagation of the combustion waves, which correlates with the results in Gubernov *et al.* (2005). This issue is further clarified in figure 2*b* where the loci for extinction (solid line) and Hopf bifurcation (dashed line) are plotted in *β* versus *u*_{a} plane for *L*_{A}=10 and *L*_{B}=1. The region above the solid curve is marked as ‘no solutions’ and represents the parameter values, where there are no travelling combustion waves. For parameter values located below the extinction locus, there are two travelling wave solutions. The slow solution is unstable. The fast solution is stable below the Hopf locus and exhibits pulsating instabilities in the parameter region between the solid and dashed curves. Since the variation of *u*_{a} does not change the qualitative behaviour of the combustion waves, in the rest of the paper we fix the value of the ambient temperature to zero.

Gubernov *et al.* (2009) demonstrated that the variation of *L*_{B} does not affect the qualitative properties of the stability diagram, although the critical parameter values for the Hopf bifurcation substantially shift towards the larger values of *β* with the increase of *L*_{B}. The extinction curve is only slightly influenced by variations in *L*_{B}, namely, the extinction curve rotates clockwise around the Bogdanov–Takens bifurcation point with the increase of *L*_{B}. It should also be noted that the location of the Bogdanov–Takens bifurcation is not affected by variation of *L*_{B}.

## 4. Period-doubling bifurcations of pulsating solutions

We investigate the properties of pulsating combustion wave solutions emerging as a result of the Hopf bifurcation when the parameters reach critical values. The governing equation (2.1) are solved in a sufficiently large coordinate domain with the boundary conditions (2.2) imposed at the edge points of the space grid. For our numerical algorithm we use the method of splitting with respect to physical processes. Initially we solve the set of ordinary differential equations that describe the temperature and the species concentration variations owing to the branching and recombination reactions by using the fourth-order Runge–Kutta algorithm. As a next step, equations of mass transfer for fuel and radicals are solved with the Crank–Nicholson method of the second-order approximation in space and time. The initial conditions for the numerical scheme are taken in the form of the travelling wave solution (or autowave) of equation (3.1).

As the Hopf curve is crossed in the parameter space, the emerging oscillatory instabilities lead to the formation of pulsating waves, which on average travel with definite speed, *c*_{drift}, however distribution of *u*(*x*,*t*), *v*(*x*,*t*) and *w*(*x*,*t*) are periodic functions of time. In order to describe these solutions it is convenient to use coordinates *ξ*=*x*−*c*_{drift}*t* in the frame travelling with the speed *c*_{drift} and trace the location of the maximum value of the radical concentration, *ξ*_{max} and its value, *w*_{max}. An example of a periodic pulsating wave is presented in Gubernov *et al.* (2009). Here, we investigate the properties of the Hopf bifurcation in greater detail. We use a standard perturbation analysis that can be found in a number of textbooks (e.g. Volpert *et al.* 1991). Therefore, below we only briefly outline the derivation of the main results and refer the reader to literature on this subject for more details. The pulsating solution bifurcating from the travelling wave solution owing to the Hopf bifurcation is sought in the form of the series
4.1
where the small parameter *A* is the amplitude of oscillations, *c* is a travelling wave speed, *β*_{h} is the activation energy, index ‘*h*’ here and in what follows implies that the parameter value is taken at the Hopf bifurcation point. The other parameter values are fixed, namely, *u*_{a}=0, *r*=0.001, *L*_{B}=1 and the value of *L*_{A} is specified below. As shown in Volpert *et al.* (1991), in the first order of the asymptotic expansion *v*_{1}(*ξ*,*t*) can be found in the form , where *v*_{1}(*ξ*) satisfies the linear stability problem (3.3). In the case of the Hopf bifurcation considered here it has three solutions: , and corresponding to *λ*=0, *iω*_{h} and −*iω*_{h}, respectively. Together with the problem (3.3) an adjoint problem is considered, where . The adjoint problem also has three bounded solutions *z*^{(0,1,2)} corresponding to *λ*=0, *iω*_{h} and −*iω*_{h}, respectively. In the first order of perturbation analysis it can be shown that *c*_{1}=0 and *β*_{1}=0 and in the second order, the following solvability condition has to be satisfied
4.2
where ∂^{2}** N**/∂

*u*^{2}denotes the vector function that can be obtained by expanding

**(**

*N***) subject to equation (4.1) into series over a small parameter**

*u**A*and taking terms of the order

*O*(

*A*

^{2}) containing

*v*_{1}. The standard inner product is used in equation (4.2), i.e. and [·,·]

_{T}means time average over a period

*T*of pulsations.

We calculate the vector functions involved in equation (4.2) numerically using the technique described in Gubernov *et al.* (2004) and then evaluate *c*_{2} according to formula (4.2). In figure 3*a*, the dependence of Δ*c*=*c*_{drift}−*c*, the difference between the pulsating and travelling wave speed, is plotted as function of the amplitude *A* of *w*_{max} oscillations for *L*_{A}=3.0, *L*_{B}=1, *r*=0.001 and *u*_{a}=0. The activation energy is varied from *β*_{h}≈4.0703 to *β*=4.075. The solid line shows the parabola Δ*c*=*c*_{2}*A*^{2}, where *c*_{2} is calculated using equation (4.2) as outlined above. The results obtained by direct numerical integration of equation (2.1) are shown with crosses. As *β* is increased from the value corresponding to Hopf bifurcation the pulsations emerge. Initially, the amplitude *A* is small and correspondence between the results derived from both approaches is good. However, as *β* is further modified to approximately 4.075 the discrepancy between the data obtained by direct integration of equation (2.1) and by using formula (4.2) becomes visible. It should be noted however that for *β*≈4.075, the amplitude *A* of oscillations becomes of the order of 50 per cent of the *w*_{max} value for the travelling combustion wave for the same parameter values. Therefore equation (4.2) describes the properties of pulsating waves reasonably well. In figure 3*b*, the dependence of *c*_{2} on *L*_{A} is plotted for other parameters being the same as in figure 3*a*. It is seen that for large values of the Lewis number for fuel, *c*_{2} tends to a constant value. In the opposite limit, as , the value of *c*_{2} sharply increases. This implies that for *L*_{A} approaching 1, the system becomes very sensitive to small variations in the parameters.

In Gubernov *et al.* (2009) it is demonstrated that as we move away from the Hopf locus in the parameter space, by increasing *β* for fixed *L*_{A}, the shape of the limit cycle deforms and becomes triangular indicating that *w*_{max}(*t*) and *ξ*_{max}(*t*) contain higher harmonics in a Fourier series expansion. It is also shown that further increase of *β* results in the formation of period 2 and period 4 solutions. In our current paper, the sequence of period-doubling bifurcations and the emerging solutions are investigated in greater detail. In figure 4, the Fourier time series of *w*_{max}(*t*) and Poincaré maps are plotted for *L*_{A}=3, *L*_{B}=1, *r*=0.001 and several values of the activation energy as is shown in the figure caption. The Fourier amplitudes are plotted in decimal logarithmic scale and the frequency is counted in inverse periods of oscillations, *T*≈3.05×10^{5}, so that *ω*=1 in figure 4 is equivalent to *T*^{−1} in normal units. In figure 4*a*(i),*b*(i) *β*=4.08 which is just above the critical value for the Hopf bifurcation *β*_{h}=4.0703…. The oscillations of *w*_{max}(*t*) are clearly periodic as seen in figure 4*a*(i). The spectrum consists of discrete lines located at equally spaced frequency intervals so that *ω*_{i}=*i* for each line, where index *i* is a natural number used to enumerate the spectral lines. The first line of the spectrum with *i*=1 is marked as *T* in order to signify that it contributes to oscillations of period *T*. The amplitudes of the higher order harmonics with *i*>2 of the Fourier spectrum of *w*_{max}(*t*) are also significant and exponentially decay with increasing *i*. This results in the triangular form of the limit cycle for *β*=4.08 as shown in fig. 6 of Gubernov *et al.* (2009). Figure 4*b*(i) depicts the Poincaré map on the *w*_{max} versus *ξ*_{max} plane. The points for the map were sampled at the moments of time *t** when *w*_{max}(*t*) reaches a local maximum. The corresponding values of *ξ*_{max}(*t**) and *w*_{max}(*t**) are then plotted. In figure 4*b*(i), all data are located in a single spot near the point with coordinates (0.0, 0.117) indicating that the dynamics of the system is of a limit cycle nature. In figure 4*a*(ii),*b*(ii), the parameters are further shifted towards the unstable region so that *β*=4.0818. In this case, the Fourier spectrum contains both the harmonics *ω*=*i* of frequency 1 and of the frequency 0.5, i.e. *ω*=0.5*i*. The spectral lines corresponding to *ω*=0.5 and to *ω*=1.0 are marked as 2*T* and *T*, respectively. The emergence of the spectral lines with frequencies multiple to 0.5 is related to the birth of limit cycle with period 2*T*. In the Poincaré map 4*b*(ii), this type of solution corresponds to a diagram with two distinct points or to a two loop curve on the (*ξ*_{max}, *w*_{max}) plane. A slight variation in the activation energy to *β*=4.0826 and to *β*=4.08274 results in sequential bifurcations to solutions with period 4*T* shown in figure 4*a*(iii),*b*(iii) and period 8*T* shown in figure 4*a*(iv),*b*(iv), respectively. The Fourier spectra of these solutions contain additional sets of lines with frequencies multiple of 0.25 marked as 4*T* and 0.125 marked as 8*T* in figure 4*a*(iii),*a*(iv). The Poincaré maps of 4*T* and 8*T* solutions are represented in figure 4*b*(iii),(iv), respectively, where it is seen that the map of period 4*T* cycle contains four points and the map of period 8*T* solution is comprised of eight distinct points. The most remarkable change in the system dynamics from regular to chaotic is observed in figure 4*a*(v),*b*(v) for *β*=4.0829. As seen from figure 5*a* the Fourier spectrum of the *w*_{max}(*t*) is becoming continuous up to spectral resolution that was 0.003 in the units used in figure 4. On a visibly irregular background, a set of spectral lines of higher amplitudes and frequencies *ω*=*i* can be observed indicating that the underlying dynamics of the system should bear the pattern of fundamental oscillations of period *T*. The Poincaré map in figure 4*b*(v) possesses the behaviour of chaotic dynamics. The successive points of the map randomly occupy a certain region on the *w*_{max} versus *ξ*_{max} plane. The distribution of the map points does not exhibit any visible pattern and appears to be uniform, which gives evidence of the irregular character of the dynamics.

In figure 5, the critical parameter values for the first three successive period-doubling bifurcations are plotted for *L*_{B}=1, *r*=0.001 on the plane with crosses connected with the solid line. The crosses represent the results obtained from the numerical integration of the governing equations (2.1) for *L*_{A}=1.5, 2, 3 and 10, whereas the solid curves are the results of the numerical data interpolation. Using the critical parameter values for period-doubling bifurcations, *β*_{i}, we estimated the border for the onset of the chaotic regime, , from the Feigenbaum’s universal period-doubling cascade (Feigenbaum 1978). Then, the predicted results were verified using the numerical integration. The data obtained in such a way is also presented in figure 5 as crosses connected with the solid line. The dotted and dashed curves show the locus for the Hopf and fold bifurcations, respectively. The region of the parameter plane located to the left from the dotted curve corresponds to the parameter values for which the stable fast autowave solution branch exist and is marked as deflagration. In the region of parameters situated between the dotted and the dashed curves, the travelling waves are linearly unstable and various pulsating regimes are observed. Finally, for parameter values below the dashed line, no travelling wave solutions are observed. As is seen in figure 5, the pulsating solutions of various periods are located in a narrow region along the Hopf bifurcation locus. Among the various periodic solutions, the strip of parameter values where the pulsating solutions of period *T* exist has the largest width and is marked as ‘*T*’. The parameter region where the period 2*T* solution is observed is just visible on a full-scale graph. The critical parameter values, *β*_{1}(*L*_{A}), *β*_{2}(*L*_{A}) and *β*_{3}(*L*_{A}) corresponding to the bifurcations of period 2*T*, 4*T* and 8*T* solutions, respectively, are located too close to be distinguished between each other; therefore, the magnification of the region of the figure near the point *L*_{A}=10 and *β*=4.125 is also plotted as an inlet to figure 5. In figure 6, loci for the Hopf and the period-doubling bifurcations, *β*_{1,2,3}(*L*_{A}), are plotted on the (, ) plane with crosses connected with interpolation lines that are marked as 1, 2, 4 and 8, respectively. The coordinates of the graph are scaled in such a way as to resolve the fine structure of the bifurcation cascade for various *L*_{A} values. The width of the parameter regions where the period *T*, 2*T*, 4*T* and 8*T* solutions exist (marked accordingly in figure 6) decrease in an exponential manner along the *β* direction. The distance between successive bifurcations *β*_{i+1}−*β*_{i} decays according to the Feigenbaums cascade (Feigenbaum 1978) with the Feigenbaum constant estimated from the numerical data as *δ*=4.7±0.1. It is important to note that the difference between *β*_{i+1} and *β*_{i} is a monotonic function of *L*_{A}: for high Lewis number for fuel, the width of the regions where the pulsating solutions exist are wider and *β*_{i+1}−*β*_{i} vanishes as *L*_{A} tends to 1.

After the border between periodic and chaotic region has been crossed, the strange attractor is formed. Figure 7*a* illustrates the projection of the strange attractor onto the plane (*ξ*_{max},*w*_{max}). In this figure, the first 200 loops of the trajectory approaching the strange attractor are plotted for *L*_{A}=3, *L*_{B}=1, *r*=0.001 and *β*=4.0829. As seen from the figure, the time oscillations of the maximum value of concentration of radicals and its location grow as *β* is varied from 4.08 for figure 4 to 4.0829 for figure 7. These oscillations are accompanied by the descent in average flame speed, and growth of the oscillations of instant velocity of the wave in the travelling coordinate frame, . For instance, for the parameter values corresponding to chaotic behaviour depicted in figure 7*a*, the mean flame velocity is , whereas the corresponding travelling wave speed is *c*=0.01733. This is seen in figure 7*c* and is zoomed in figure 7*d*, where the dependence of *c* on *β* is plotted. In figure 7*c*, the stable solution branch is shown with the solid line and the dashed line represent the unstable travelling wave solution. In figure 7*d*, the solid lines represent the numerical data and the dashed lines show the imaginary continuation of the solution branches which become unstable owing to period-doubling bifurcations. In the enlarged figure 7*d*, the sequence of period-doubling bifurcations can be clearly observed. The instant velocity along the direction of flame wave propagation is reaching maximum values of and in the backward direction the speed minimum is . It should be noted that such a significant change in the propagation velocity is attained as the control parameter *β* is altered from 4.0703 for Hopf bifurcation to 4.0829, i.e. for about 0.3 per cent of its value.

As the activation energy is further increased, the abrupt change in the dynamics of the system can be observed as illustrated in figure 7*b* for *β*=4.0831. After encircling several loops in *w*_{max} versus *ξ*_{max} plane in a region where the strange attractor is expected to be situated, the trajectory departs away to the point located at *w*_{max}=0. In other words, the solution of equation (2.1) reaches the stationary state *u*(*x*,*t*)=0, *v*(*x*,*t*)=1, *w*(*x*,*t*)=0 corresponding to extinction. The same scenario is observed for other values of *L*_{A}∈(1,10]. For all *L*_{A} values, the extinction has been found close to the border of chaotic region in parameter space and its width along the *β* direction is estimated to be of the order of .

## 5. Chaotic transient

As shown in §4, when the border of chaotic region is crossed and *β* becomes greater than , the strange attractor emerges in the phase space of the system. In order to study the properties of chaotic attractor, we calculate the largest Lyapunov exponent, *σ* (Temam 1993). In figure 8, the dependence of *σ* on *β* is plotted for *L*_{A}=10, *L*_{B}=1, *r*=0.001 and *u*_{a}=0. The numerical error corresponding to the standard deviation is shown with bars. The value of the Lewis number is fixed to *L*_{A}=10 for illustrative purposes, since the distances between the successive bifurcations of period-doubling, chaos and extinction borders in the parameter space are the largest in this case. The estimated critical value of the activation energy calculated by using the Feigenbaum scenario as discussed in §4 is . This number is very close to the border of chaotic region in figure 8 where it is found between *β*=4.144371 for which *σ*=0 and *β*=4.144372 for which *σ* is positive. Therefore, this implies that the estimate of the critical parameters for the onset of chaos using the Feigenbaum criterion is very accurate. For *β* greater than the critical value, *σ* becomes positive and grows with further increase of *β* confirming the chaotic nature of the system dynamics.

Further increase of *β* up to the values about 4.1444 results in the occurrence of the abrupt changes in the dynamics of the system, i.e. sudden extinction of the flame propagation as discussed in §4. In the phase space of the system, there always exists an asymptotic regular attractor related to fresh mixture state, where the reaction is frozen and the temperature is equal to its ambient value. This stationary state has its own basin of attraction and at certain parameter values it becomes connected to the basin of the strange attractor. As a result, the spatio-temporal chaos collapses after some time to the regular state of fresh mixture. This type of dynamics is called the chaotic transient and is known to manifest itself in various dynamical systems including the reaction–diffusion systems (Tel & Lai 2008). The collapse of chaotic oscillations occurs in an abrupt and random manner. We have studied the statistical properties of this process. For a given fixed parameter values, the solutions are calculated by direct numerical integration of equation (2.1) for various initial conditions taken on the basin of the strange attractor and the lifetime of pulsating solutions before the extinction is estimated. In figure 9, the dependence of the number of extinction events versus the time of integration is plotted as a histogram, i.e. time is split into equal intervals of length, Δ*T*=1.26×10^{5} and the number of evens happening during this interval is counted. The parameters are taken as *L*_{A}=10, *L*_{B}=1, *β*=4.14435, *r*=0.001 and *u*_{a}=0. The overall number of tries exceeds 500 for figure 9. The *N*(*t*) dependence shows characteristic exponential distribution behaviour that describes the times between events in a Poisson random process. The linear regression allows to obtain the average lifetime, which is found to be *τ*=4.75×10^{5}, the corresponding Pearson correlation coefficient 0.94. In figure 9, the resulting fit of the data to the exponential distribution is shown with the dashed line. As *β* is increased, the average lifetime of propagating solution decays very rapidly. Even the slight variations of the activation energy cause substantial change in the values of *τ*: for *β*=4.14435 the value of *τ* is 4.75×10^{5}, for *β*=4.144375 *τ* is 3.1×10^{5}, and for *β*=4.1444 it drops to 2.35×10^{5}. Such a high sensitivity to changes of control parameters is typical for supertransients (Tel & Lai 2008). However, the clarification of the detailed characteristics of the chaotic transient is the subject of ongoing work.

## 6. Conclusions and discussion

In this paper, we have investigated the properties of travelling and pulsating combustion waves in a model with a two-step chain-branching reaction mechanism in the adiabatic limit. The effect of the ambient temperature variation on the properties of the travelling waves is investigated. We demonstrate that the increase of the ambient temperature results in the growth of both the flame speed and expands the domain of *β* values for which the travelling combustion waves exist. At the same time the *u*_{a} variation does not qualitatively affect the behaviour of *c* on *β* dependence. It is shown that the pulsating instabilities can emerge in the model for *L*_{A}>1 owing to the Hopf bifurcation. The critical parameter values for the Hopf bifurcation are found in the parameter space and it is shown to originate from the Bogdanov–Takens bifurcation at *L*_{A}=1, when the locus for the Hopf and fold bifurcations meet. The increase of *u*_{a} shifts both the boundary of the existence of travelling solutions and the Hopf loci towards higher values of the dimensionless activation energy, i.e. the ambient temperature rise has a stabilizing effect on the propagation of the combustion waves, which correlates with the results in Gubernov et al. (2005) for the one-step model. The variation of *u*_{a} does not change the qualitative behaviour of the travelling combustion waves, however, it alters the quantitative properties.

In the current paper, it is found that as the locus for the Hopf bifurcation is crossed in parameter space the pulsating waves emerge as a result of the supercritical Hopf bifurcation. The properties of the Hopf bifurcation are investigated in detail by means of the perturbation analysis and direct numerical integration. The pulsating waves are propagating with certain average speed, and in the coordinate frame moving with this speed the temperature, concentration of fuel and radicals are periodic functions of time with period *T*. It is demonstrated that the pulsating solution propagates with the average velocity smaller than the corresponding speed of the travelling wave solution. Also it is shown that for , the system becomes very sensible to small variations of parameters, i.e. small changes of *β* alters substantially. It is convenient to describe the pulsating waves in terms of the maximum values of radicals concentration, *w*_{max} and the coordinate of this maximum in the co-moving frame, *ξ*_{max}, which become the periodic functions of time for period *T* solutions. In the plane (*ξ*_{max}, *w*_{max}), the period *T* solutions correspond to the limit cycle. For the parameter values close to the Hopf locus the limit cycle has an elliptic form and *ξ*_{max}(*t*) and *w*_{max}(*t*) are harmonic functions with the period of oscillations governed by the frequency of oscillations of the unstable modes of the linear stability problem (Gubernov *et al.* 2008*b*).

As we move away from the Hopf locus in the parameter space by increasing the activation energy, the limit cycle expands in the (*ξ*_{max}, *w*_{max}) plane and deforms to a triangular form. At the same time in the Fourier time series of *w*_{max} the high-order frequencies, which are multiples of the basic frequency *T*^{−1}, become significant. Further increase of *β*, while keeping the other parameters fixed, leads to the emergence of the period 2*T* solutions as a result of the period-doubling bifurcation at *β*_{1}. In the (*ξ*_{max}, *w*_{max}) plane these solutions correspond to two-loop trajectories, the Fourier spectrum is shown to include the lines with both the multiples of *T*^{−1} and (2*T*)^{−1}, and the Poincaré map consists of two points. Similarly, for higher values of activation energy, we find the second and third period-doubling bifurcations at *β*_{2} giving birth to period 4*T* solutions and *β*_{3} leading to period 8*T* solutions, which correspond to four and eight points map on Poincaré map, respectively. The distance between successive bifurcations *β*_{i+1}−*β*_{i} decays according to the Feigenbaum universal cascade (Feigenbaum 1978) with the Feigenbaum constant estimated from the numerical data as *δ*=4.7±0.1. Based on the knowledge of the critical parameter values *β*_{1,2,3}, and by assuming that the sequence of period-doubling bifurcations follow the Feigenbaum scenario, the boundary of chaotic region is estimated. If the parameter values are taken above this boundary, the dynamics of the system becomes chaotic and is characterized by the continuous Fourier spectrum of *w*_{max}(*t*), irregular distribution of the points of the Poincaré map and positive maximum Lyapunov exponent. The boundary of the chaos onset obtained from the Lyapunov exponents calculations agrees well with the estimate based on the Feigenbaum scenario.

The period *iT* solutions (*i*=1,2,4,…) as well as chaotic solution propagate with certain average speed and can also be described by an instant speed, , of the pulsations of the wavefront location in the co-moving coordinate frame travelling with . Here, we denote the index *i* to each solution branch, i.e. *i*=0 to travelling wave solution, *i*=1 to solution of period *T*, etc. On the *c* versus *β* plane, the solution of period *T* emerges from the travelling solution branch at the Hopf bifurcation point and each sequential solution branch of period *iT* originates from the previous branch at the point of corresponding period-doubling bifurcation, so that forms the tree-like structure. All solution branches are monotonically decaying functions of activation energy; moreover, the pulsating solutions on average always travel with the speed smaller than the flame speed corresponding to the travelling wave solution. We demonstrate that such a decay of the flame speed is very sensitive to the increase of *β*, i.e. small variations of the activation energy results in significant change in values.

The width of chaotic region along *β* is limited by the extinction of the propagating flame. As *β* is increased over the strange attractor is formed, however at certain value of activation energy, *β*_{e}, the propagating chaotic solution disappears and the flame extinguishes to a trivial stationary state *u*=0, *v*=1, and *w*=0. This phenomenon is known in the literature as chaotic transient. As the activation energy is increased above , the basin of attraction of the strange attractor changes and at a certain value *β*_{e} it becomes connected owing to the unstable–unstable pair bifurcation with the basin of attraction of the regular attractor—the trivial solution. In this case, the initial conditions taken in the form of the propagating solution evolves similarly to chaotic waves; however, at a certain moment of time it enters the basin of attraction of the trivial state and the combustion wave extinguishes. The collapse of the chaotic oscillations occurs in an abrupt and random manner. Our studies of the statistical properties of this process indicate that the lifetime of the transient is described by the exponential distribution. The average lifetime *τ* is very sensitive to the variations of control parameters, which is rather common to supertransients. However, the quantitative description of the extinction limit and chaotic transient characteristics are the subjects of our future work.

The location of the regions of parameters with various types of pulsating solutions are found in the parameter space. It is demonstrated that in the (*L*_{A}, *β*) plane they are situated along the Hopf bifurcation locus in the form of narrow strips that expand with increasing *L*_{A} and shrink to a point corresponding to the Bogdanov–Takens bifurcation of the travelling wave solution as . More formally this can be formulated as *β*_{i+1}−*β*_{i} vanishes for and monotonically grows as *L*_{A} increases. This conclusion correlates with the experimental observations that pulsating waves of periods 1*T*, 2*T* and 4*T* are readily observed in solid phase combustion, where *L*_{A}≫1, and the width of the regions in parameter space corresponding to various pulsating solutions is relatively high. By contrast, for *L*_{A}∼1 even the period *T* solutions exist for such a narrow region in the parameter space that the experimental observation of one-dimensional pulsating waves is not possible in practice. The distance between the successive bifurcations, *β*_{i+1}−*β*_{i}, vanishes in an exponential manner as the index *i* is increased. The structure of the parameter space in the regions where the pulsating solutions with periods higher than 2*T* exist is complex. A small variation of parameters here can lead to a change in the type of solution and even to extinction. Our preliminary studies of the extinction phenomenon show that the width of the chaotic region is small in comparison with the width of the region, where the unstable travelling wave exists, i.e. the region between the Hopf and fold bifurcations. Therefore, we expect that the pulsating solutions extinguish before the travelling wave solution branch ceases to exist in the space of parameters. In future, we plan to investigate in detail the route to extinction via chaotic oscillations.

## Acknowledgements

V. V. Gubernov, A. V. Kolobov and A. A. Polezhaev would like to acknowledge the financial support from the Russian Foundation for Basic Research: grants 08-02-00806, 08-02-00682, 08-01-00131; Homeland Science Support Fund: program ‘Young Candidates and Doctors of Science’. H. S. Sidhu and G. N. Mercer would also like to acknowledge the support of the Australian Research Council Grant DP0878146. V. V. Gubernov and H. S. Sidhu are grateful to G. A. Gottwald for fruitful discussions.

## Footnotes

- Received December 24, 2009.
- Accepted March 5, 2010.

- © 2010 The Royal Society