ISRN Applied MathematicsVolume 2013 (2013), Article ID 580561, 13 pageshttp://dx.doi.org/10.1155/2013/580561

Research Article

Weakly Nonlinear and Numerical Analysis of Auto-Oscillatory Dynamics in a Solid Propellant Combustion Model

Michael Reardon and Jun Yu

Department of Mathematics and Statistics, University of Vermont, 16 Colchester Avenue, Burlington, VT 05405, USA

Received 19 March 2013; Accepted 9 April 2013

Academic Editors: I. Doltsinis, J. R. Fernandez, J. Shen, and H. C. So

Copyright © 2013 Michael Reardon and Jun Yu. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In solid combustion, a planar traveling flame wave may lose stability to oscillatory modes in certain parameter regimes. The loss of stability is commonly achieved through a Hopf bifurcation in which linearly unstable modes grow to fixed amplitudes. This situation has been well studied in models of strictly condensed combustion and has been observed experimentally. In this paper, we describe the onset and evolution of oscillatory modes arising in a free-interface model of solid propellant combustion. After a brief introduction, the structure and stability of small perturbations to a planar flame wave are then investigated. Multiple-scale expansions are utilized to describe changes in the amplitude and the phase of a single, weakly unstable mode, and the results are compared to numerical predictions. It is found the over much of the parameter space that the analytical and numerical solutions predict a supercritical Hopf bifurcation similar to those found in models of strictly condensed combustion.

1. Introduction and Problem Formulation

The departure from steadyplanar combustion results from the amplification of small disturbances to the basic traveling wave structure due to instability. When instability is present in the absence of forcing, it is referred to as intrinsic instability which results from the inability of the solid reagent to adjust instantaneously when the temperature profile of the planar traveling wave is perturbed. If the reaction rate is particularly sensitive to temperature, these perturbations will lead to either increased or decreased reaction zone temperature and propagation speed. This situation is only temporary, however. An increase in the temperature and propagation speed of the reaction zone is eventually counteracted by a decrease in the thermal energy content of the reagent immediately ahead of the front as well as increased abstraction of thermal energy from the reaction zone. Similarly, a decrease in temperature and propagation speed will be followed by a relative buildup of thermal energy ahead of the front which is eventually reabsorbed into the reaction zone. Such a system is therefore possessed of a mechanism for correcting small deviations from the steady, planar configuration. However, the relatively large thermal inertia of the solid reagent may lead to overcorrection as the system attempts to return to the basic configuration. Thus, sustained oscillations of the front temperature and velocity about their steady, planar values can occur. The presence of self-excited oscillations about the basic planar wave is referred to as auto-oscillatory combustion. Oscillations may be planar (pulsating), nonplanar, or a combination of the two.

Intrinsic instability may be present in any reaction-diffusion system where reagents and products are separated by a quasiplanar boundary. However, due to their tendency to produce more violent oscillations, nonintrinsic instability types including acoustic and instability have received the bulk of the research attention in the area of solid propellant combustion and motor design. The study of auto-oscillatory combustion due to intrinsic instability has been explored more thoroughly in models of strictly condensed combustion, though we note the work of Zel’dovich and later Novozhilov in the Russian literature and Denison and Baum in laying the groundwork for similar studies in propellant combustion. The first numerical evidence of auto-oscillatory combustion appeared in the oft-cited work by Shkadinsky et al. [1] wheresustained oscillations were observed in a numerical simulation of a condensed phase combustion model. Soon after, experimental results confirming the existence of auto-oscillatory combustion were produced by Merzhanov et al. in [2]. The increased popularity of self-propagating high-temperature synthesis (SHS) as means for fabricating materials motivated several studies which sought to describe the transition from the steady planar deflagration to an oscillatory regime using both numerical and analytical methods. A good review of both experiment and theory as it pertains to oscillatory combustion in SHS processes can be found in [3]. These studies have predicted a wide variety of nonsteady burning regimes where the temporal form of the oscillations is usually dependent on physical parameters, while their spatial form depends on the model geometry. Rectangular geometries have been found to produce sinusoidal variations in front temperature, location, and velocity [4–6], while cylindrical geometries have been found to produce more complex spatial patterns including “multiple point” and “spinning” combustion where spots or bands of relatively higher or lower temperature appear and disappear or rotate around the combustion zone, respectively [7–9]. In either geometry, a planar, pulsating mode is usually possible as are oscillations containing multiple spatial modes. The temporal evolution of condensed combustion waves has proven equally complex and may include growth to a fixed-amplitude, periodic or quasiperiodic limit cycle, or a transition to relaxational [10, 11] and chaotic [12–14] regimes. The latter is often accompanied by period doubling cascades [15–18] in the parameter space.

Much of the analytical work in auto-oscillatory combustion has focused on describing the evolution of small perturbations containing linearly unstable modes beyond an initial period of exponential growth. Linear stability analyses typically allow for growing or decaying sinusoidal oscillations of the temperature and velocity profile about the basic solution see, for example, [8, 16, 19–21] for stability analyses in various models and geometries. In the case that a particular mode is linearly unstable, its amplitude is usually constrained by nonlinear effects and eventually approaches a finite value. In this case, the basic solution loses stability via a supercritical Hopf bifurcation which cannot be predicted by the linear analysis alone and requires more advanced analytical methods or numerical simulation. The most common analytical technique is the method of multiple scales which has been employed in [8, 22, 23] among others. In this study, we apply similar methods to a mathematical model of solid propellant combustion which was developed by Margolis and Armstrong in [20]. It will be shown that the evolution of the amplitude of a nonplanar perturbation to the laminar flame wave evolves in superslow time according to a nonlinear ODE of Stuart-Landau type.

1.1. The Physical Model

The solid propellant combustion model is initially set in a 3D rectangular spatial geometry in which there exists a region of homogeneous solid reagent and a region of gaseous byproducts separated by a quasiplanar reaction zone containing both sublimated reagent and product gases. The reaction is limited to a band-oriented transverse to its direction of travel. The combustion zone is assumed to move principally in the negative direction with -velocity given by , where denotes the location of the sublimation surface. The physical picture is summarized in Figure 1(a). In the gas phase, the transverse velocities are assumed to be small compared with the component in the positive- direction, denoted , which is assumed to vary with , , and . The authors of [20] assume a small Mach number and that the pressure is large and approximately constant with the majority of the heat rise occurring in the solid phase preheat zone. These assumptions effectively suppress acoustically driven forcing and instability found in other models of propellant combustion. The resulting instability is therefore intrinsic and similar to that found in strictly condensed combustion. With these assumptions the gas velocity in the direction can then be related to the regression rate of the sublimation surface via mass continuity yielding which is then used to remove from the model equations. Here, and are the (constant) densities of the solid reagent and gaseous products, respectively. The authors of [20] also note that this model with the given assumptions is particularly applicable to the burning of ammonium perchlorate as a monopropellant.

(a) The dimensional model geometry

(b) The nondimensional asymptotic model geometry

Figure 1: Comparison of dimensional and nondimensional asymptotic model geometries.

A closed system of equations can then be derived using thermodynamic, hydrodynamic, and pyrolysis laws which relate the unknown reagent concentration function , the sublimation surface location, and the temperature distribution and is given by equations (2.23)–(2.31) in [20]. The temperature of the reagent for is given by , while the temperature of the products approaches the adiabatic flame temperature for . The parameters and represent the gas and solid phase heat capacities, respectively while is the heat of reaction. Margolis and Armstrong proceed by nondimensionalizing the model equations and applying a change of coordinates where the laboratory frame is replaced by a front-attached nondimensional coordinate frame. The physical coordinates and the dimensional functions and are then replaced by their nondimensional counterparts: In the above, is equal to the yet undetermined value in the case of steady planar combustion. The nondimensional parameters , , , and are defined as the gas-to-solid ratios of the density, heat capacity, thermal diffusivity, and activation energy, respectively.

1.2. The Asymptotic Model

Margolis and Armstrong [20] then consider the nondimensional model in the case of a large modified nondimensional activation energy , where is the solid-phase activation energy and is the gas constant. In this well-known limit, the combustion zone collapses to an asymptotically thin sheet located at . A stretched coordinate is introduced and the nondimensional model is divided into an outer nonreacting region for and an inner combustion region for . The functions , , and as well as the nondimensional model equations are then expanded in powers of and the methods of singular perturbation theory are employed. In the outer problem, the thickness of the reaction zone goes to 0 and it is therefore coincident with the sublimation surface which we now denote simply by . The function is replaced by a Heaviside function. Asymptotic matching produces a closed outer problem where the inner combustion zone dynamics are removed in exchange for jump conditions at . We refer to the resulting system of equations as the “asymptotic model” which is given by the piecewise PDE the jump conditions and the -boundary and continuity conditions The dynamics outside of the asymptotically thin reaction zone are essentially captured by (2a) and (2b), while after matching and simplification, the reaction zone dynamics are captured by jump conditions (2c) and (2d) which determine the relationship between front speed, temperature, and reaction rate as well as the heat flux across the interface. In this study we further assume the lateral boundary conditionsThe operator appearing in (2a) and (2b) is the Laplacian in the moving coordinate system and is given by The function appearing in (2c) is determined implicitly from where is defined as and is a constant. The nondimensional parameters , , , and in (2a)–(2f) are defined in terms of the gas-to-solid ratios of density, thermal diffusivity, heat capacity, and activation energy, respectively. Note that the dimensional propagation velocity is dependent on the “burn rate eigenvalue” which must be determined by numerically solving the implicit equation (4). A discussion on steady planar propagation and results on the propagation rate for various parameter values are given in Section  6 of [20].

1.3. Linear Stability of the Basic Solution

An immediate advantage over the dimensional model is that (2a)–(2f) admits a basic solutionwhich describes a laminar combustion wave propagating into the reagent at a constant velocity . We now wish to study the form of small perturbations about the basic solutions (6a) and (6b). We first note that a linear stability analysis was carried out in [20] which we summarize here. First, we introduce a small, nondimensional parameter whose physical significance will be determined later. Let our perturbed solution be given by where and are given by the basic solutions (6a) and (6b). We then substitute these expansions into the asymptotic model equations (2a)–(2f) and obtain a set of equations to be satisfied at each order of . The expansion of in powers of is discussed in Appendix A where the functions at are expressed in terms of where . In particular, we note that and where the constant is given by . The equations are satisfied by the basic solution, and we obtain the following set of equations at :where the differential operator denotes the traditional 3D Laplacian: .

Separable solutions to (8a)–(8c) and (8e)–(8f) are given bywhere , , and . The arbitrary constants and determine the (complex) amplitude of the linear mode and a shift in the position of combustion front, respectively. The complex exponents and are given by with . The linear solutions lead to decaying oscillations of the temperature profile as one moves away from the combustion front in either direction. After inserting these forms into (8d) and canceling common factors, we obtain the dispersion relation For , we find that corresponds to exponential growth of the linear solutions, while for , the linear solutions decay.

The requirement along with the two dispersion relation equations can be used to obtain the values of and corresponding to neutral stability for a given . We denote the neutrally stable value of corresponding to a particular value of as , while the corresponding value of is denoted by . A typical neutral stability curve is shown in Figure 2. With respect to the neutral stability boundary, the analysis predicts that linear modes are stable for . We note that each combination of the nondimensional physical parameters will correspond to a different neutral stability curve, though the general shape of these curves is such that there exists a global minimum for some followed by a sharp increase in for large . This corresponds to the planar mode being generally more unstable than the nonplanar modes with stability increasing for large . However, certain geometries may allow for small , nonplanar modes to lose stability before the planar mode. Of these the modes are the most linearly unstable.

Figure 2: A typical neutral stability curve.

2. Weakly Nonlinear Analysis

We next wish to describe the evolution of a single, weakly unstable perturbation to the basic solution. Let be perturbed by a small amount above the neutral stability boundary, that is, let , . The model geometry is chosen so that is an allowable wave number of the linear solution while all other modes remain stable. Upon expanding the dispersion relation about one obtains , . We expect the behavior of a weakly unstable linear mode in this parameter regime to be similar to its neutrally stable counterpart, oscillating in time with frequency approximately equal to and possessing a similar spatial form while its amplitude slowly grows. For larger amplitudes, however, nonlinear effects may act to reinforce or retard growth and alter the temporal oscillation frequency. The true solution should then exhibit behavior on at least two time scales: sinusoidal oscillation in “fast” time with frequency approximately equal to , and amplitude/phase changes occurring in “slow” time. Changes in the geometry of the combustion front and the temporal pulsations about the basic solution may also result in a departure from the basic propagation rate which may occur on a separate slow time scale. In subsequent analysis, we employ the method of multiple scales which will allow us to naturally incorporate these separate time scales and more accurately describe the dynamics for large times as the solution transitions from a linear to a nonlinear regime.

2.1. The Multiple-Scale Equations

To ease the analytical and numerical burden, we first restrict our analysis to the case of two spatial dimensions, say the and dimensions, which is appropriate in the case of small as compared to . Next, we introduce the time scales , , and . We consider and replace the differential operator in the asymptotic model equations (2a)–(2f) with . We also expand and asymptotically in powers of : where and are given by our basic solution with replaced by . The model equations (2a)–(2f) are then expanded and separated at orders of . The equations at are satisfied while the equations at for are given byThe linear operators , , and in (13a) are defined similarly to those in (8a)–(8f) with the -derivatives removed and replaced by .

2.2. The Problem

The equations at are homogeneous and nearly identical to (8a)–(8f). We propose a solution form which has been modified slightly from (9a) and (9b). We first remove the -dependence of the linear solutions, and therefore . As we have chosen the parameters to ensure all modes with decay exponentially in fast time, we will neglect these modes in our weakly nonlinear solution form and include only the mode. We also replace the “true” (computed with ) with (computed with ) which removes the slow exponential growth from the linear solution. Instead, we allow the arbitrary amplitude constant to become a function of the slow time scales. The constant from the linear solution is also allowed to become a function of the slow time scales in order to allow for an adjustment to the average front location and speed from that of the laminar wave. The weakly nonlinear solution is then given to by for and Note that the constants and are defined using instead of as in (9a) and (9b). Also, the roots and in (14a) and (14b) defined in (10) are computed using the complex as we expect the spatial form of the weakly nonlinear solution to more closely resemble the perturbed linear solution, though the difference from their corresponding neutrally stable values is .

2.3. Solvability Conditions

In subsequent sections, we employ the Fredholm alternative theorem (FAT) to obtain solvability conditions on the inhomogeneous multiscale equations at , . These solvability conditions result in differential equations which describe the evolution of and and preclude the appearance of secular terms in the solution at each order. The FAT implies that a solution to the linear operator equation exists if and only if the RHS is orthogonal to the kernel of the adjoint operator . In the current context, is a differential operator defined on a function space whose elements satisfy appropriate smoothness conditions as well as homogeneous boundary and jump conditions. A full discussion on the adjoint problem including the definition of as well as the derivation of and can be found in [24]. The inner product in the range space of (which contains the domain of ) is defined by where we have assumed . Separable solutions to are of the form where and . There also exists a planar, time-independent solution with given by Only the planar solution and the nonplanar solution whose wave number equals to that of the linear solution will produce nontrivial solvability conditions in the present problem.

2.4. The Problem

Upon inserting our solution forms into the equations, we obtain the differential equations and jump conditions for and as follows: with or .

As the domain of definition of our differential operator contains functions which satisfy homogeneous boundary and jump conditions, we are unable to directly apply the FAT. Therefore, we seek a change of variables such that . We define our new function by where and are two real valued twice-differentiable functions in satisfying the following boundary and jump conditions: After inserting into (19a)–(19d) we arrive at a new set of equations:

We are now in a position to apply the FAT, and upon inserting the RHS of (22a) and (22b) into the inner product (16) along with the time-dependent adjoint solution (16) and simplifying, we obtain the solvability condition Requiring the RHS of (22a) and (22b) to be orthogonal to (16) leads to the solvability condition and therefore The constant is defined in Appendix B.

2.5. The Problem

After obtaining the solutions for and as described in Appendix C we can insert these into the problem at , where we find a system of inhomogeneous equations of the form:where or . Proceeding as in the problem, we again change variables and apply the FAT to the transformed equations. After simplifying the orthogonality condition with the time-dependent adjoint solution, we obtain the Stuart-Landau differential equation for The constant is defined in Appendix D. Applying the FAT with the time-independent adjoint solution yields the solvability condition As , we obtain We note that the constants and that multiply the homogeneous solution to the problem must be determined at higher order which is not pursued here.

3. A Comparison of Numerical and Analytical Results

In this section, we use (27) to make analytical predictions for the evolution of a single unstable mode which are then compared with a direct numerical simulation of the asymptotic model equations (2a)–(2f). Due to its desirable stability properties, the Crank-Nicolson method was employed to carry out the simulation, and its implementation is described in detail in [24]. In all simulations we set , , and .

3.1. Limit Cycles for

From the form of (27), it is evident that there is a fixed point at . The linearization at this fixed point determines that it is unstable for which we have in fact chosen our parameters to ensure. Therefore, we recover the exponential growth for small predicted by the linear stability analysis. This growth is modeled well by the linear solution until the cubic term in the S-L equation becomes important. To examine this scenario further, we make the change of variables where , and we let and . This form is then inserted into (27). After equating real and imaginary parts, the complex differential equation for can be transformed into a pair of partially decoupled real differential equations for and :Note that (31a) does not depend on , while (31b) depends on alone. From the form of (31a), we see that there is a fixed point at corresponding to the aforementioned fixed point at . There is also a possibility of a second fixed point at . As we are given , this fixed point exists whenever in which case there is then a corresponding limit cycle with radius which is a global attractor for .

Note that (31a) is real and separable for for which an analytical solution is given by The constant is real and positive and is defined as This function can then be substituted into (31b) which can be subsequently integrated to obtain the following: where is given in terms of the initial phase by From the form of these equations, it is clear that as , . The limit cycle frequency is then shifted by from the neutrally stable frequency: .

3.2. Structure and Evolution of a Nonplanar Mode

In the linear stability analysis, we found that the oscillations consist of modes which decay as one moves away from the interface in the -direction. Comparing the values of and , one finds the dynamics are stretched in the positive -direction due to the additional velocity of the gaseous exhaust moving away from the front (i.e., the decay is faster in the direction as ) and oscillations about this value caused by the imaginary parts of and die away. The decay rate has important implications for the numerical solution as it is necessary to greatly extend the numerical -domain, especially for small in which case numerical simulation becomes computationally prohibitive for 2 or 3 physical dimensions.

In the case where the lowest nonzero wavenumber mode is the only unstable mode, the temperature on one side of the front oscillates about the basic solution value out of phase with the opposite side with these oscillations propagating away from the reaction front with time. These variations in the temperature profile can be seen both along and propagating away from the flame sheet in the snapshot of the temperature distribution shown in Figure 3. This sinusoidal spatial variation is also found in the front position and propagation velocity. One side of the front will propagate with greater speed during periods of higher temperature and enhanced combustion, advancing further into the unreacted medium only to slow and eventually be overtaken by the opposite side as evidenced by Figure 4.

Figure 3: Numerical solution for .

Figure 4: and : numerical interface location and temperature for . The dashed interface location corresponds to the analytical prediction.

3.3. Varying the Parameters

The analytical predictions for the radius and frequency of the limit cycle for various , , and are shown in Figure 5. We only include the region of parameter space for which the weakly nonlinear analysis also predicts a finite limit cycle radius. Beyond these regions, the analysis predicts a subcritical Hopf bifurcation at the neutral stability boundary. For a single nonplanar mode, only positive limit cycle frequencies are predicted which correspond to positive shifts from the neutrally stable frequency. The qualitative dependence of the solutions on is such that smaller values yield small sinusoidal oscillations which become distorted for larger , eventually leading to relaxation oscillations. The analysis predicts that increased leads to an increase in the limit cycle radii, while the region of the physical parameter space where a supercritical Hopf bifurcation is predicted shrinks. We also find increased contributions from higher harmonics and a subsequent loss of agreement between the analytical and numerical solutions as seen in Figures 6 and 7 where the oscillations in front temperature and speed about their steady planar values are plotted. Again, it can be seen that the oscillations in front speed are in phase with those of front temperature and also of greater (nondimensional) amplitude

(a)

(b)

(c)

Figure 5: Limit cycle frequency and radius for varying , and with .

Figure 6: Numerical and analytical values of (smaller amplitude oscillations) and (larger amplitude oscillations) along with their frequency power spectra for . In this parameter region, the numerical and analytical solutions are indistinguishable even for large times.

Figure 7: Numerical and analytical values of and along with their frequency power spectra for . Higher harmonics are quite evident in the numerical solution, while the analytical solution contains only first and second harmonics as it is determined only to .

To explore the dynamics in a region where the analysis predicts a subcritical Hopf bifurcation, we choose to vary the nondimensional heat capacity ratio . From Figure 5, we expect that as we decrease from for fixed , the analytical and numerical solutions will experience increased amplitudes and frequency shifts. This is what we see in Figures 7, 8(a), and 8(b) which correspond to , and 0.85, respectively. The trend is clearly towards relaxation oscillations characterized by rapid increases in front temperature and velocity followed by extended periods of depressed combustion. We note that in the final figure, the analytical solution is not shown as the limit cycle has disappeared and the analytical solution is unbounded. The numerical solution remains bounded in this parameter region as it appears to be attracted to a higher-order limit cycle. We note that the tendency towards instability and relaxation oscillations as is decreased might be expected physically as a decrease in the gas-to-solid heat capacity ratio would correspond to an increase in the amount of heat required to raise the temperature of the reagents, effectively increasing the thermal inertia of the solid reagents (relative to the products) and thus leading to larger overcorrection and increased instability as discussed in the introduction. Increasing was found to produce similar results as a relatively smaller solid-state thermal diffusivity would also tend to increase the thermal inertia of the solid reagent relative to the products.

(a)

(b)

Figure 8: and for decreasing with and . For , the analytical solution has become unbounded and is not plotted.

4. Concluding Remarks

We have applied a nonlinear analysis to a model of solid propellant combustion and encountered dynamical behavior similar to that found in models of strictly condensed combustion. The analytical results were compared to numerical simulations and found to be in excellent agreement provided that a nondimensional parameter related to the activation energy was chosen to lie a small distance from the neutral stability boundary. As is increased, oscillations about the planar traveling wave take on a more relaxational character. This transition was also observed for fixed as the physical parameters were varied to increase the relative thermal inertia of the reagent with respect to the product. We conjecture that as the parameters are pushed further into these regions increasingly, the dynamics will become increasingly relaxational, possibly leading to more complex behavior including period doubling, cascades, and chaotic combustion similar to those observed in condensed combustion. Unfortunately, exploration of such regions of the parameter space is numerically prohibitive, although we have observed such transitions for a single pulsating mode when the present model is restricted to a single spatial dimension.

Appendices

A. The Expansion of (4)

The jump condition (2d) involves the unknown function whose values are given implicitly by (4). In order to obtain the multiple-scale equations, it is necessary to expand in powers of about the basic solution value and obtain equivalent expressions in terms of , . After substituting these expansions into (4) and expanding the functions and integrals in powers of , we obtain at an equation at each order of . At , we find where the parameter grouping is defined as Note that constant is proportional to the nondimensional activation energy . In the special case when , we have, . At , we have where At we find that the function is given by where Note that for , .

B. Terms from the Problem

The constant appearing in (25) is defined with

C. The Solutions

The solutions consist of a homogeneous and a particular solution. To solve for and , we employ the method of undetermined coefficients and propose a particular solution with the same form as the inhomogeneities: where and The functions lie in the null space of with or given by As , we must choose . Furthermore, which is the reason for the introduction of in (C.1).

After replacing and with their respective values 0 and in the RHS of the equations and substituting our solution forms into the LHS, we obtain 8 equations from equation, 4 equations from equation, and 4 equations each from and equations. We also obtain 4 equations from the continuity condition at . Upon inspection, we find that does not enter into any of these equations, though as the functional form is identical to that of the homogeneous solution, we may choose . Equating terms in the equations yields 9 solutions and 3 equations: where the constants without tildes are then defined similarly to their counterparts in – as: We also obtain 3 equations: As is known, we find for , otherwise, we obtain three equations from the continuity condition: The and equations yield six more equations: for . Equations (C.5)–(C.8) form a linear system of 12 equations for the twelve remaining unknowns which can then be solved numerically. Lastly, we note the existence of a homogeneous solution given by of and where and .

D. Terms from the Problem

The constant appearing in (27) is defined as where In these summations, .

References

K. G. Shkadinsky, B. I. Khaikin, and A. G. Merzhanov, “Propagation of a pulsating exothermic reaction in the condensed phase,” Combustion, Explosion, and Shock Waves, vol. 7, pp. 15–22, 1971.

A. G. Merzhanov, A. K. Filonenko, and I. P. Borovinskaia, “New phenomena during the burning of condensed systems,” Doklady Physical Chemistry, vol. 208, pp. 122–125, 1973.

B. V. Novozhilov, “Nonlinear SHS phenomena: experiment, theory, numerical modeling,” Pure and Applied Chemistry, vol. 64, no. 7, pp. 955–964, 1992.

J. Yu and L. K. Gross, “The onset of linear instabilities in a solid combustion model,” Studies in Applied Mathematics, vol. 107, no. 1, pp. 81–101, 2001.

M. Frankel, L. K. Gross, and V. Roytburd, “Thermo-kinetically controlled pattern selection,” Interfaces and Free Boundaries, vol. 2, no. 3, pp. 313–330, 2000.

S. B. Margolis and B. J. Matkowsky, “Flame propagation in channels: secondary bifurcation to quasiperiodic pulsations,” SIAM Journal on Applied Mathematics, vol. 45, no. 1, pp. 93–129, 1985.

A. Bayliss, B. J. Matkowsky, and A. P. Aldushin, “Dynamics of hot spots in solid fuel combustion,” Physica D, vol. 166, no. 1-2, pp. 104–130, 2002.

S. B. Margolis, H. G. Kaper, G. K. Leaf, and B. J. Matkowsky, “Bifurcation of pulsating and spinning reaction fronts in condensed two-phase combustion,” Combustion Science and Technology, vol. 43, no. 3-4, pp. 127–165, 1985.

A. Volpert, A. Volpert, S. P. Davtyan, I. N. Megrabova, and N. F. Surkov, “Threedimensional modes of unsteady solid-ame combustion,” Chaos, vol. 13, no. 1, pp. 80–86, 2003.

A. Bayliss, B. J. Matkowsky, and M. Minkoff, “Period doubling gained, period doubling lost,” SIAM Journal on Applied Mathematics, vol. 49, no. 4, pp. 1047–1063, 1989.

A. Bayliss and B. J. Matkowsky, “Fronts, relaxation oscillations, and period doubling in solid fuel combustion,” Journal of Computational Physics, vol. 71, no. 1, pp. 147–168, 1987.

A. Bayliss and B. J. Matkowsky, “Two routes to chaos in condensed phase combustion,” SIAM Journal on Applied Mathematics, vol. 50, no. 2, pp. 437–459, 1990.

S. B. Margolis, “Chaotic combustion of solids and high-density uids near points of strong resonance,” Proceedings Mathematical and Physical Sciences, vol. 433, no. 1887, pp. 131–150, 1991.

I. Brailovsky and G. I. Sivashinsky, “Chaotic dynamics in solid fuel combustion,” Physica D, vol. 65, no. 1-2, pp. 191–198, 1993.

J. Yu, L. K. Gross, and C. M. Danforth, “Complex dynamic behavior during transition in a solid combustion model,” Complexity, vol. 14, no. 6, pp. 9–14, 2009.

M. Frankel, V. Roytburd, and G. Sivashinsky, “A sequence of period doublings and chaotic pulsations in a free boundary problem modeling thermal instabilities,” SIAM Journal on Applied Mathematics, vol. 54, no. 4, pp. 1101–1112, 1994.

S. B. Margolis, “New routes to quasiperiodic combustion of solids and high-density fluids near resonant Hopf bifurcation points,” SIAM Journal on Applied Mathematics, vol. 51, no. 3, pp. 693–726, 1991.

S. B. Margolis, “An asymptotic theory of condensed two-phase flame propagation,” SIAM Journal on Applied Mathematics, vol. 43, no. 2, pp. 351–369, 1983.

B. V. Novozhilov, “Nonlinear oscillations of combustion velocity of powder,” Journal of Applied Mechanics and Technical Physics, vol. 7, no. 5, pp. 19–25, 1966.

S. B. Margolis and R. C. Armstrong, “Two asymptotic models for solid propellant combustion,” Combustion Science and Technology, vol. 47, no. 1-2, pp. 1–18, 1986.

M. Garbey, H. G. Kaper, G. K. Leaf, and B. J. Matkowsky, “Quasi-periodic waves and the transfer of stability in condensed-phase surface combustion,” SIAM Journal on Applied Mathematics, vol. 52, no. 2, pp. 384–395, 1992.

L. K. Gross and J. Yu, “Weakly nonlinear and numerical analyses of dynamics in a solid combustion model,” SIAM Journal on Applied Mathematics, vol. 65, no. 5, pp. 1708–1725, 2005.

M. R. Booty, S. B. Margolis, and B. J. Matkowsky, “Interaction of pulsating and spinning waves in condensed phase combustion,” SIAM Journal on Applied Mathematics, vol. 46, no. 5, pp. 801–843, 1986.

M. S. Reardon, Weakly Nonlinear Analysis of a Solid Propellant Combustion Model [Ph.D. thesis], University of Vermont, Burlington, Vt, USA, 2012.