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

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.


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 reactiondiffusion 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 autooscillatory 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 2 ISRN Applied Mathematics model. Soon after, experimental results confirming the existence of auto-oscillatory combustion were produced by Merzhanov et al. in [2]. The increased popularity of selfpropagating 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][5][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][8][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][13][14] regimes. The latter is often accompanied by period doubling cascades [15][16][17][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][20][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.

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 negativẽ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̃1,̃2, 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 thẽ3 direction can then be related to the regression rate of the sublimation surface via mass continuity yielding̃= (1 − (̃/̃))(Φ /̃) which is then used to removẽ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 closed system of equations can then be derived using thermodynamic, hydrodynamic, and pyrolysis laws which relate the unknown reagent concentration functioñ (̃,̃,̃,̃), the sublimation surface location, and the temperature distributioñ(̃,̃,̃,̃) and is given by equations (2.23)-(2.31) in [20]. The temperature of the reagent for̃≪ Φ is given bỹ, while the temperature of the products approaches the adiabatic flame temperaturẽ= (̃/̃)(̃+ ) for̃≫Φ . The parameters and represent the gas and solid phase heat capacities, respectively whilẽ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Φ /ĩ n 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.

The Asymptotic
Model. Margolis and Armstrong [20] then consider the nondimensional model in the case of a large modified nondimensional activation energy Δ = (1 − (̃/̃))(̃/̃̃) ≫ 1, wherẽis the solidphase activation energy and̃is the gas constant. In this well-known limit, the combustion zone collapses to an asymptotically thin sheet located at = 0. A stretched coordinate is introduced and the nondimensional model is divided into an outer nonreacting region for | | ≫ 0 and an inner combustion region for ≈ 0. The functions , Θ, and Φ as well as the nondimensional model equations are then expanded in powers of 1/Δ 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 = 0. 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 lim 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 conditions =0, The operator∇ 2 appearing in (2a) and (2b) is the Laplacian in the moving coordinate system and is given bỹ The function (Θ| =0 ) 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-tosolid ratios of density, thermal diffusivity, heat capacity, and activation energy, respectively. Note that the dimensional propagation velocitỹis dependent on the "burn rate eigenvalue" 0 = exp(− (1)) 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].

Linear Stability of the Basic Solution. An immediate advantage over the dimensional model is that (2a)-(2f) admits a basic solution
which describes a laminar combustion wave propagating into the reagent at a constant velocity Φ/ = −1. 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 ≪ 1 whose physical significance will be determined later. Let our perturbed solution be given by where 0 and 0 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 (Θ| =0 ) in powers of is discussed in Appendix A where the functions at ( ) are expressed in terms of { } where ≤ . In particular, we note that 0 = (1) = − ln 0 and 1 = × 1 where the constant is given by (A.2). The (1) equations are satisfied by the basic solution, and we obtain the following set of equations at ( ): where the differential operator ∇ 2 denotes the traditional 3D Laplacian: Separable solutions to (8a)-(8c) and (8e)-(8f) are given by where , ∈ R, ∈ Z, = / and Γ = 1 − ( / ). The arbitrary constants ∈ Z and ∈ R 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 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 < 0 corresponds to exponential growth of the linear solutions, while for > 0, the linear solutions decay.
The requirement = 0 along with the two dispersion relation equations Re(D) = Im(D) = 0 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 > 0 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 = min modes are the most linearly unstable.

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 = Ω ( min ) + 2 , ≪ 1. The model geometry is chosen so that min is an allowable wave number of the linear solution while all other modes remain stable. Upon expanding the dispersion relation about ( min , Ω( min )) one obtains = Ω + 2 2 + ( 3 ), 2 ∈ Z. 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.

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 T = 2 . We consider = ( , ,T) and replace the differential operator / in the asymptotic model equations (2a)-(2f) with ( / ) + ( / ) + 2 ( / T). We also expand Θ( , , , , T) and Φ( , , , T) asymptotically in powers of : where 0 and 0 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 (1) are satisfied while the equations at ( ) for ≥ 1 are given by The linear operators L ± , M, and N 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 ̸ = min decay exponentially in fast time, we will neglect these modes in our weakly nonlinear solution form and include only the = min mode. We also replace the "true" (computed with = Ω ( min ) + 2 ) with Ω (computed with = Ω ( min )) 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 Note that the constants = /Ω and Γ = 1 − 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 ).

Solvability Conditions.
In subsequent sections, we employ the Fredholm alternative theorem (FAT) to obtain solvability conditions on the inhomogeneous multiscale equations at ( ), = 2,3. These solvability conditions result in differential equations which describe the evolution of ( , T) and ( , T) and preclude the appearance of secular terms in the solution at each order. The FAT implies that a solution to the linear operator equation L( ) = exists if and only if the RHS is orthogonal to the kernel of the adjoint operator L † . In the current context, L( , ) 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 L as well as the derivation of L † and ker(L † ) can be found in [24]. The inner product in the range space of L (which contains the domain of L † ) is defined by where we have assumed → 0. Separable solutions to L † ( † ) = 0 are of the form † ( , , , ) = { † † } † cos ( ) for ≶ 0, (16) where † = 1 2 (−1 + √ 1 + 4 ( † + 2 )) , † = 1 2 (−1 − √ 1 + 4 2 2 ( † and † = . There also exists a planar, time-independent solution with = † = 0 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.
As the domain of definition of our differential operator L contains functions which satisfy homogeneous boundary and jump conditions, we are unable to directly apply the FAT. Therefore, we seek a change of variables 2 →̃2 such that M(̃2, 2 ) = N(̃2, 2 ) = 0. We define our new functioñ2 bỹ2 where ( ) and ( ) are two real valued twice-differentiable functions in 2 (R) satisfying the following boundary and jump conditions: After inserting̃2 into (19a)-(19d) we arrive at a new set of ( 2 ) equations: We are now in a position to apply the FAT, and upon inserting the RHS of (22a) and (22b) into the inner product Requiring the RHS of (22a) and (22b) to be orthogonal to (16) leads to the solvability condition and therefore The constant 2 is defined in Appendix B.

The ( 3 ) Problem.
After obtaining the ( 2 ) solutions for 2 and 2 as described in Appendix C we can insert these into the problem at ( 3 ), where we find a system of inhomogeneous equations of the form: , , ( , , T) Ω cos ( ) + , < 0, , , ( , , T) Ω cos ( ) + , where , ∈ {1, 3} or { , } = {0, 0}. Proceeding as in the ( 2 ) 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 (T) The constant 3 is defined in Appendix D. Applying the FAT with the time-independent adjoint solution yields the solvability condition 2 +̃(T) = ( 2 + 2 ) 2 | (T)| 2 .
As / = 2 | (T)| 2 , we obtaiñ We note that the constants 2 and 2 that multiply the homogeneous solution to the ( 2 ) problem must be determined at higher order which is not pursued here.

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 Υ = 1, = 1, and ] = 2.

Limit Cycles for (T).
From the form of (27), it is evident that there is a fixed point at = 0. The linearization at this fixed point determines that it is unstable for Im( 2 ) < 0 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 , : [0, ∞) → R, and we let 2 = 2 + 2 and 3 = 3 + 3 . This form is then inserted into (27). After equating real and imaginary parts, the complex differential equation for (T) can be transformed into a pair of partially decoupled real differential equations for (T) and (T): 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 = 0 corresponding to the aforementioned fixed point at = 0. There is also a possibility of a second fixed point at̂= √ 2 / 3 . As we are given 2 < 0, this fixed point exists whenever 3 < 0 in which case there is then a corresponding limit cycle with radiuŝwhich is a global attractor for (T).
Note that (31a) is real and separable for ∉ {0,̂} for which an analytical solution is given by

ISRN Applied Mathematics
The constant 1 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 2 is given in terms of the initial phase 0 = (0) by From the form of these equations, it is clear that as →̂, → ( 2 + 3 2 / 3 )T + Constant =̂T + Constant. The limit cycle frequency is then shifted by ( 2 ) from the neutrally stable frequency: Ω → Ω + 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 | Re( )| < | Re( )|) 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

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 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 = 1 for fixed = = 1, the analytical and numerical solutions will experience increased amplitudes and frequency shifts. This is what we see in 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.

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 2 ≪ 1 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.

A. The Expansion of (4)
The jump condition (2d) involves the unknown function (Θ| =0 ) 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 Θ| =0 = 1 and obtain equivalent expressions in terms of { | =0 }, ≤ . 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 ] = 2, we have, = Δ. At ( 2 ), we have At ( 3 ) we find that the function 3 is given by .