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

Authors: Michael Reardon; Jun Yu
Journal: ISRN Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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 asintrinsic 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 andL* 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 negativez~ direction with z~-velocity given by ∂Φ~S/∂t~, where Φ~S(x~,y~,t~) 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-z~ direction, denoted u~g, which is assumed to vary with x~1, x~2, and t~. 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 x~3 direction can then be related to the regression rate of the sublimation surface via mass continuity yielding u~g=(1-(ρ~s/ρ~g))(∂Φ~s/∂t~) which is then used to remove u~g from the model equations. Here, ρ~g and ρ~s 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.Comparison of dimensional and nondimensional asymptotic model geometries.
(a)
The dimensional model geometry
(b)
The nondimensional asymptotic model geometryA closed system of equations can then be derived using thermodynamic, hydrodynamic, and pyrolysis laws which relate the unknown reagent concentration functionY~(x~,y~,z~,t~), the sublimation surface location, and the temperature distribution T~(x~,y~,z~,t~) and is given by equations (2.23)–(2.31) in [20]. The temperature of the reagent for z~≪Φ~s is given by T~u, while the temperature of the products approaches the adiabatic flame temperature T~a=(c~g/c~s)(T~u+β~) for z~≫Φ~s. The parameters cg and cs 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 Φ~S and T~ are then replaced by their nondimensional counterparts:
(1)x=U~λ~sx~,y=U~λ~sy~,z=U~λ~s(z~-Φ~S),t=U~2λ~st~,Θ=T~-T~uT~a-T~u,Φ=U~λsΦ~S.
In the above, U~ is equal to the yet undetermined value ∂Φ~s/∂t~ in the case of steady planar combustion. The nondimensional parameters ρ, c, λ, 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 Δ=(1-(T~u/T~a))(E~s/R~gT~a)≫1, where E~s is the solid-phase activation energy and R~g is the gas constant. In this well-known limit, the combustion zone collapses to an asymptotically thin sheet located at z=0. A stretched coordinate is introduced and the nondimensional model is divided into an outer nonreacting region for |z|≫0 and an inner combustion region for z≈0. The functions Y, Θ, and ΦS 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 Y 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 z=0. We refer to the resulting system of equations as the “asymptotic model” which is given by the piecewise PDE(2a)∂Θ∂t-∂Φ∂t∂Θ∂z=∇~2Θ,z<0,(2b)ρc∂Θ∂t-c∂Φ∂t∂Θ∂z=ρcλ∇~2Θ,z>0,
the jump conditions
(2c)∂Φ∂t=-[1+(∂Φ∂x)2+(∂Φ∂y)2]1/2exp(F(Θ|z=0)-F(1)),(2d)ρcλ∂Θ∂z|z=0+-∂Θ∂z|z=0-=-[1+(∂Φ∂x)2+(∂Φ∂y)2]-1((c-1)Θ|z=0-c)∂Φ∂t,
and the z-boundary and continuity conditions
(2e)limz→-∞Θ=0,limz→∞Θ=1,limz→0+Θ=limz→0-Θ.
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(2f)∂η∂w|w=0,Lw=0,forw=x,y,η=Θ,Φ.The operator ∇~2 appearing in (2a) and (2b) is the Laplacian in the moving coordinate system and is given by
(3)∇~2=∂2∂x2+∂2∂y2+[1+(∂Φ∂x)2+(∂Φ∂y)2]∂2∂z2-2∂Φ∂x∂2∂x∂z-2∂Φ∂y∂2∂y∂z-(∂2Φ∂x2+∂2Φ∂y2)∂∂z.
The function F(Θ|z=0) appearing in (2c) is determined implicitly from
(4)Υexp(νΔ2(Θ|z=0-1))×[G(Δ(Θ|z=0-1)-F)]1/2=exp(F),
where G is defined as
(5)G(u)=∫0u‍e-νwwndw
and Υ is a constant. The nondimensional parameters ρ, c, λ, 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 U~ is dependent on the “burn rate eigenvalue” λ0=exp(-F(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].


### 1.3. Linear Stability of the Basic Solution

An immediate advantage over the dimensional model is that (2a)–(2f) admits a basic solution(6a)Θ(x,y,z,t)={ez1forz≶0,(6b)Φ(x,y,t)=-twhich describes a laminar combustion wave propagating into the reagent at a constant velocity ∂Φ/∂t=-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
(7)Θ(x,y,z,t)=θ0(z)+ϵθ1(x,y,x,t),Φ(x,y,t)=ϕ0(t)+ϵϕ1(x,y,t),
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 F(Θ|z=0) in powers of ϵ is discussed in Appendix A where the functions Fi at O(ϵi) are expressed in terms of {θj} where j≤i. In particular, we note that F0=F(1)=-lnλ0 and F1=C×θ1 where the constant C is given by (A.2). The O(1) equations are satisfied by the basic solution, and we obtain the following set of equations at O(ϵ):(8a)ℒ-(θ1,ϕ1)=∂θ1∂t+∂θ1∂z-∇2θ1+ez(∂2ϕ1∂x2+∂2ϕ1∂y2-∂ϕ1∂t)=0,z<0,(8b)ℒ+(θ1,ϕ1)=ρc∂θ1∂t+c∂θ1∂z-ρcλ∇2θ1=0,z>0,(8c)ℳ(θ1,ϕ1)=∂ϕ1∂t+C×θ1|z=0=0,(8d)𝒩(θ1,ϕ1)=ρcλ∂θ1∂z|z=0+-∂θ1∂z|z=0--(c-1)θ1|z=0-∂ϕ1∂t=0,(8e)θ1|z=0-z=0+=0,limz→±∞θ1=0,(8f)∂η∂w|w=0,Lw=0,forw=x,y,η=θ,ϕ,where the differential operator ∇2 denotes the traditional 3D Laplacian: (∂2/∂x2)+(∂2/∂y2)+(∂2/∂z2).Separable solutions to (8a)–(8c) and (8e)–(8f) are given by(9a)θ1(x,y,z,t)={A(γez+Γerz)eiωtcos(mx)cos(ny)+CCAeRzeiωtcos(mx)cos(ny)+CCforz≶0,(9b)ϕ1(x,y,t)=iCAωeiωtcos(mx)cos(ny)+12B+CC,where m,n∈ℝ, ω∈ℤ, γ=iC/ω and Γ=1-(iC/ω). The arbitrary constants A∈ℤ and B∈ℝ determine the (complex) amplitude of the linear mode and a shift in the position of combustion front, respectively. The complex exponents r and R are given by(10)r=12(1+1+4(iω+κ2)),R=12ρλ(1-1+4ρ2λ2(iωλ+κ2)),
with κ2=m2+n2. 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
(11)𝒟(ω,κ;C)=ρcλR-iCω-r(1-iCω)+C+1-c=0.
For ω=ωR+iωI, we find that ωI<0 corresponds to exponential growth of the linear solutions, while for ωI>0, the linear solutions decay.The requirementωI=0 along with the two dispersion relation equations Re(𝒟)=Im(𝒟)=0 can be used to obtain the values of C and ωR 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 C is denoted by CΩ(κ). 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 C<CΩ(κ). 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 CΩ(κ) 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.Figure 2
A typical neutral stability curve.


## 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 negativez~ direction with z~-velocity given by ∂Φ~S/∂t~, where Φ~S(x~,y~,t~) 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-z~ direction, denoted u~g, which is assumed to vary with x~1, x~2, and t~. 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 x~3 direction can then be related to the regression rate of the sublimation surface via mass continuity yielding u~g=(1-(ρ~s/ρ~g))(∂Φ~s/∂t~) which is then used to remove u~g from the model equations. Here, ρ~g and ρ~s 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.Comparison of dimensional and nondimensional asymptotic model geometries.
(a)
The dimensional model geometry
(b)
The nondimensional asymptotic model geometryA closed system of equations can then be derived using thermodynamic, hydrodynamic, and pyrolysis laws which relate the unknown reagent concentration functionY~(x~,y~,z~,t~), the sublimation surface location, and the temperature distribution T~(x~,y~,z~,t~) and is given by equations (2.23)–(2.31) in [20]. The temperature of the reagent for z~≪Φ~s is given by T~u, while the temperature of the products approaches the adiabatic flame temperature T~a=(c~g/c~s)(T~u+β~) for z~≫Φ~s. The parameters cg and cs 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 Φ~S and T~ are then replaced by their nondimensional counterparts:
(1)x=U~λ~sx~,y=U~λ~sy~,z=U~λ~s(z~-Φ~S),t=U~2λ~st~,Θ=T~-T~uT~a-T~u,Φ=U~λsΦ~S.
In the above, U~ is equal to the yet undetermined value ∂Φ~s/∂t~ in the case of steady planar combustion. The nondimensional parameters ρ, c, λ, 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 Δ=(1-(T~u/T~a))(E~s/R~gT~a)≫1, where E~s is the solid-phase activation energy and R~g is the gas constant. In this well-known limit, the combustion zone collapses to an asymptotically thin sheet located at z=0. A stretched coordinate is introduced and the nondimensional model is divided into an outer nonreacting region for |z|≫0 and an inner combustion region for z≈0. The functions Y, Θ, and ΦS 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 Y 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 z=0. We refer to the resulting system of equations as the “asymptotic model” which is given by the piecewise PDE(2a)∂Θ∂t-∂Φ∂t∂Θ∂z=∇~2Θ,z<0,(2b)ρc∂Θ∂t-c∂Φ∂t∂Θ∂z=ρcλ∇~2Θ,z>0,
the jump conditions
(2c)∂Φ∂t=-[1+(∂Φ∂x)2+(∂Φ∂y)2]1/2exp(F(Θ|z=0)-F(1)),(2d)ρcλ∂Θ∂z|z=0+-∂Θ∂z|z=0-=-[1+(∂Φ∂x)2+(∂Φ∂y)2]-1((c-1)Θ|z=0-c)∂Φ∂t,
and the z-boundary and continuity conditions
(2e)limz→-∞Θ=0,limz→∞Θ=1,limz→0+Θ=limz→0-Θ.
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(2f)∂η∂w|w=0,Lw=0,forw=x,y,η=Θ,Φ.The operator ∇~2 appearing in (2a) and (2b) is the Laplacian in the moving coordinate system and is given by
(3)∇~2=∂2∂x2+∂2∂y2+[1+(∂Φ∂x)2+(∂Φ∂y)2]∂2∂z2-2∂Φ∂x∂2∂x∂z-2∂Φ∂y∂2∂y∂z-(∂2Φ∂x2+∂2Φ∂y2)∂∂z.
The function F(Θ|z=0) appearing in (2c) is determined implicitly from
(4)Υexp(νΔ2(Θ|z=0-1))×[G(Δ(Θ|z=0-1)-F)]1/2=exp(F),
where G is defined as
(5)G(u)=∫0u‍e-νwwndw
and Υ is a constant. The nondimensional parameters ρ, c, λ, 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 U~ is dependent on the “burn rate eigenvalue” λ0=exp(-F(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].


## 1.3. Linear Stability of the Basic Solution

An immediate advantage over the dimensional model is that (2a)–(2f) admits a basic solution(6a)Θ(x,y,z,t)={ez1forz≶0,(6b)Φ(x,y,t)=-twhich describes a laminar combustion wave propagating into the reagent at a constant velocity ∂Φ/∂t=-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
(7)Θ(x,y,z,t)=θ0(z)+ϵθ1(x,y,x,t),Φ(x,y,t)=ϕ0(t)+ϵϕ1(x,y,t),
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 F(Θ|z=0) in powers of ϵ is discussed in Appendix A where the functions Fi at O(ϵi) are expressed in terms of {θj} where j≤i. In particular, we note that F0=F(1)=-lnλ0 and F1=C×θ1 where the constant C is given by (A.2). The O(1) equations are satisfied by the basic solution, and we obtain the following set of equations at O(ϵ):(8a)ℒ-(θ1,ϕ1)=∂θ1∂t+∂θ1∂z-∇2θ1+ez(∂2ϕ1∂x2+∂2ϕ1∂y2-∂ϕ1∂t)=0,z<0,(8b)ℒ+(θ1,ϕ1)=ρc∂θ1∂t+c∂θ1∂z-ρcλ∇2θ1=0,z>0,(8c)ℳ(θ1,ϕ1)=∂ϕ1∂t+C×θ1|z=0=0,(8d)𝒩(θ1,ϕ1)=ρcλ∂θ1∂z|z=0+-∂θ1∂z|z=0--(c-1)θ1|z=0-∂ϕ1∂t=0,(8e)θ1|z=0-z=0+=0,limz→±∞θ1=0,(8f)∂η∂w|w=0,Lw=0,forw=x,y,η=θ,ϕ,where the differential operator ∇2 denotes the traditional 3D Laplacian: (∂2/∂x2)+(∂2/∂y2)+(∂2/∂z2).Separable solutions to (8a)–(8c) and (8e)–(8f) are given by(9a)θ1(x,y,z,t)={A(γez+Γerz)eiωtcos(mx)cos(ny)+CCAeRzeiωtcos(mx)cos(ny)+CCforz≶0,(9b)ϕ1(x,y,t)=iCAωeiωtcos(mx)cos(ny)+12B+CC,where m,n∈ℝ, ω∈ℤ, γ=iC/ω and Γ=1-(iC/ω). The arbitrary constants A∈ℤ and B∈ℝ determine the (complex) amplitude of the linear mode and a shift in the position of combustion front, respectively. The complex exponents r and R are given by(10)r=12(1+1+4(iω+κ2)),R=12ρλ(1-1+4ρ2λ2(iωλ+κ2)),
with κ2=m2+n2. 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
(11)𝒟(ω,κ;C)=ρcλR-iCω-r(1-iCω)+C+1-c=0.
For ω=ωR+iωI, we find that ωI<0 corresponds to exponential growth of the linear solutions, while for ωI>0, the linear solutions decay.The requirementωI=0 along with the two dispersion relation equations Re(𝒟)=Im(𝒟)=0 can be used to obtain the values of C and ωR 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 C is denoted by CΩ(κ). 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 C<CΩ(κ). 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 CΩ(κ) 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.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. LetC be perturbed by a small amount above the neutral stability boundary, that is, let C=CΩ(κ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+O(ϵ3), ω2∈ℤ. 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 thez and x dimensions, which is appropriate in the case of small Ly as compared to Lx. Next, we introduce the time scales τ=t, T=ϵt, and 𝒯=ϵ2t. We consider t=t(τ,T,𝒯) and replace the differential operator ∂/∂t in the asymptotic model equations (2a)–(2f) with (∂/∂τ)+ϵ(∂/∂T)+ϵ2(∂/∂𝒯). We also expand Θ(x,z,τ,T,𝒯) and Φ(x,τ,T,𝒯) asymptotically in powers of ϵ:
(12)Θ~∑i=0∞‍ϵiθi(x,z,τ,T,𝒯),Φ~∑i=0∞‍ϵiϕi(x,τ,T,𝒯),
where θ0 and ϕ0 are given by our basic solution with t replaced by τ. The model equations (2a)–(2f) are then expanded and separated at orders of ϵ. The equations at O(1) are satisfied while the equations at O(ϵj) for j≥1 are given by(13a)ℒ-(θj,ϕj)=l~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),z<0,(13b)ℒ+(θj,ϕj)=L~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),z>0,(13c)ℳ(θj,ϕj)=M~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),(13d)𝒩(θj,ϕj)=N~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),(13e)∂θj∂x|x=0,Lx=∂ϕ∂x|x=0,Lx=0,limz→±∞θj=0=θj|z=0-z=0+.The linear operators ℒ±, ℳ, and 𝒩 in (13a) are defined similarly to those in (8a)–(8f) with the y-derivatives removed and t replaced by τ.


### 2.2. TheO(ϵ) Problem

The equations atO(ϵ) 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 y-dependence of the linear solutions, and therefore κ=m. As we have chosen the parameters to ensure all modes with m≠κmin decay exponentially in fast time, we will neglect these modes in our weakly nonlinear solution form and include only the m=κmin mode. We also replace the “true” ω (computed with C=CΩ(κmin)+ϵ2) with Ω (computed with C=CΩ(κmin)) which removes the slow exponential growth from the linear solution. Instead, we allow the arbitrary amplitude constant A to become a function of the slow time scales. The constant B 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 O(ϵ) by(14a)Θ(x,z,τ,T,𝒯)={ez+ϵ(A(T,𝒯)(γez+Γerz)eiΩτcos(mx)+CC)1+ϵ(A(T,𝒯)eRzeiΩτcos(mx)+CC),
for z≶0 and
(14b)Φ(x,τ,T,𝒯)=-τ+ϵ(γeiΩτcos(mx)+CC+B(T,𝒯)).Note that the constants γ=iC/Ω and Γ=1-γ are defined using Ω instead of ω as in (9a) and (9b). Also, the roots r and R 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 O(ϵ2).


### 2.3. Solvability Conditions

In subsequent sections, we employ theFredholm alternative theorem (FAT) to obtain solvability conditions on the inhomogeneous multiscale equations at O(ϵj), j=2,3. These solvability conditions result in differential equations which describe the evolution of A(T,𝒯) and B(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 ℒ(u)=f 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 ker(ℒ†) can be found in [24]. The inner product in the range space of ℒ (which contains the domain of ℒ†) is defined by
(15)〈u,v〉=limt→∞1tLx∫0t‍∫-∞∞‍∫0Lx‍uv-dxdzdτ,
where we have assumed Ly→0. Separable solutions to ℒ†(θ†)=0 are of the form(16)θ†(x,y,z,τ)={er†zeR†z}eiω†τcos(mx)forz≶0,
where
(17)r†=12(-1+1+4(iω†+m2)),R†=12ρλ(-1-1+4ρ2λ2(iω†λ+m2)),
and ω†=ω-. There also exists a planar, time-independent solution with m=ω†=0 given by
(18)θ0†(z)={1forz<0e-(1/ρλ)zforz>0.
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. TheO(ϵ2) Problem

Upon inserting our solution forms into theO(ϵ2) equations, we obtain the differential equations and jump conditions for θ2 and ϕ2 as follows: (19a)ℒ-(θ2,ϕ2)=l~2=∑j,k‍l~2,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,(19b)ℒ+(θ2,ϕ2)=L~2=∑j,k‍L~2,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,(19c)ℳ(θ2,ϕ2)=M~2=∑j,k‍M~2,j,k(T,𝒯)ejiΩτcos(kmx)+CC,(19d)𝒩(θ2,ϕ2)=N~2=∑j,k‍N~2,j,k(T,𝒯)ejiΩτcos(kmx)+CC, with j,k∈{0,2} or {j,k}={1,1}.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 θ2→θ~2 such that ℳ(θ~2,ϕ2)=𝒩(θ~2,ϕ2)=0. We define our new function θ~2 by
(20)θ~2=θ2-f(z)M~2-g(z)N~2,
where f(z) and g(z) are two real valued twice-differentiable functions in L2(ℝ) satisfying the following boundary and jump conditions:
(21)limz→±∞f(z)=0,f|z=0-z=0+=0,f(0)=1C,ρcλf′(0+)-f′(0-)=c-1C,limz→±∞g(z)=0,g|z=0-z=0+=0,g(0)=0,ρcλg′(0+)-g′(0-)=1.
After inserting θ~2 into (19a)–(19d) we arrive at a new set of O(ϵ2) equations:(22a)ℒ(θ~2,ϕ2)={ℒ-(θ~2,ϕ2)=l~2-ℒ-(fM~2,0)-ℒ-(gN~2,0),z<0ℒ+(θ~2,ϕ2)=L~2-ℒ+(fM~2,0)-ℒ+(gN~2,0),z<0(22b)ℳ(θ~2,ϕ2)=N(θ~2,ϕ2)=0.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
(23)∂A∂T=0.
Requiring the RHS of (22a) and (22b) to be orthogonal to (16) leads to the solvability condition
(24)∂B∂T=β2|A(𝒯)|2,
and therefore
(25)B(T,𝒯)=β2|A(𝒯)|2T+B~(𝒯).
The constant β2 is defined in Appendix B.


### 2.5. TheO(ϵ3) Problem

After obtaining theO(ϵ2) solutions for θ2 and ϕ2 as described in Appendix C we can insert these into the problem at O(ϵ3), where we find a system of inhomogeneous equations of the form:(26a)ℒ-(θ3,ϕ3)=l~3=∑j,k‍l~3,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,z<0,(26b)ℒ+(θ3,ϕ3)=L~3=∑j,k‍L~3,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,z>0,(26c)ℳ(θ3,ϕ3)=M~3=∑j,k‍M~3,j,k(T,𝒯)ejiΩτcos(kmx)+CC,(26d)𝒩(θ3,ϕ3)=N~3=∑j,k‍N~3,j,k(T,𝒯)ejiΩτcos(kmx)+CC,where j,k∈{1,3} or {j,k}={0,0}. Proceeding as in the O(ϵ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 A(𝒯)(27)dAd𝒯=iω2A+β3A|A|2.
The constant β3 is defined in Appendix D. Applying the FAT with the time-independent adjoint solution yields the solvability condition(28)B2∂B∂T+B~′(𝒯)=(A2+A-2)β2|A(𝒯)|2.
As ∂B/∂T=β2|A(𝒯)|2, we obtain(29)B~′(𝒯)=(A2+A-2-B2)β2|A(𝒯)|2.
We note that the constants A2 and B2 that multiply the homogeneous solution to the O(ϵ2) problem must be determined at higher order which is not pursued here.


## 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 thez and x dimensions, which is appropriate in the case of small Ly as compared to Lx. Next, we introduce the time scales τ=t, T=ϵt, and 𝒯=ϵ2t. We consider t=t(τ,T,𝒯) and replace the differential operator ∂/∂t in the asymptotic model equations (2a)–(2f) with (∂/∂τ)+ϵ(∂/∂T)+ϵ2(∂/∂𝒯). We also expand Θ(x,z,τ,T,𝒯) and Φ(x,τ,T,𝒯) asymptotically in powers of ϵ:
(12)Θ~∑i=0∞‍ϵiθi(x,z,τ,T,𝒯),Φ~∑i=0∞‍ϵiϕi(x,τ,T,𝒯),
where θ0 and ϕ0 are given by our basic solution with t replaced by τ. The model equations (2a)–(2f) are then expanded and separated at orders of ϵ. The equations at O(1) are satisfied while the equations at O(ϵj) for j≥1 are given by(13a)ℒ-(θj,ϕj)=l~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),z<0,(13b)ℒ+(θj,ϕj)=L~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),z>0,(13c)ℳ(θj,ϕj)=M~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),(13d)𝒩(θj,ϕj)=N~j(θj-1,θj-2,…,θ0,ϕj-1,ϕj-2,…,ϕ0),(13e)∂θj∂x|x=0,Lx=∂ϕ∂x|x=0,Lx=0,limz→±∞θj=0=θj|z=0-z=0+.The linear operators ℒ±, ℳ, and 𝒩 in (13a) are defined similarly to those in (8a)–(8f) with the y-derivatives removed and t replaced by τ.


## 2.2. TheO(ϵ) Problem

The equations atO(ϵ) 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 y-dependence of the linear solutions, and therefore κ=m. As we have chosen the parameters to ensure all modes with m≠κmin decay exponentially in fast time, we will neglect these modes in our weakly nonlinear solution form and include only the m=κmin mode. We also replace the “true” ω (computed with C=CΩ(κmin)+ϵ2) with Ω (computed with C=CΩ(κmin)) which removes the slow exponential growth from the linear solution. Instead, we allow the arbitrary amplitude constant A to become a function of the slow time scales. The constant B 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 O(ϵ) by(14a)Θ(x,z,τ,T,𝒯)={ez+ϵ(A(T,𝒯)(γez+Γerz)eiΩτcos(mx)+CC)1+ϵ(A(T,𝒯)eRzeiΩτcos(mx)+CC),
for z≶0 and
(14b)Φ(x,τ,T,𝒯)=-τ+ϵ(γeiΩτcos(mx)+CC+B(T,𝒯)).Note that the constants γ=iC/Ω and Γ=1-γ are defined using Ω instead of ω as in (9a) and (9b). Also, the roots r and R 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 O(ϵ2).


## 2.3. Solvability Conditions

In subsequent sections, we employ theFredholm alternative theorem (FAT) to obtain solvability conditions on the inhomogeneous multiscale equations at O(ϵj), j=2,3. These solvability conditions result in differential equations which describe the evolution of A(T,𝒯) and B(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 ℒ(u)=f 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 ker(ℒ†) can be found in [24]. The inner product in the range space of ℒ (which contains the domain of ℒ†) is defined by
(15)〈u,v〉=limt→∞1tLx∫0t‍∫-∞∞‍∫0Lx‍uv-dxdzdτ,
where we have assumed Ly→0. Separable solutions to ℒ†(θ†)=0 are of the form(16)θ†(x,y,z,τ)={er†zeR†z}eiω†τcos(mx)forz≶0,
where
(17)r†=12(-1+1+4(iω†+m2)),R†=12ρλ(-1-1+4ρ2λ2(iω†λ+m2)),
and ω†=ω-. There also exists a planar, time-independent solution with m=ω†=0 given by
(18)θ0†(z)={1forz<0e-(1/ρλ)zforz>0.
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. TheO(ϵ2) Problem

Upon inserting our solution forms into theO(ϵ2) equations, we obtain the differential equations and jump conditions for θ2 and ϕ2 as follows: (19a)ℒ-(θ2,ϕ2)=l~2=∑j,k‍l~2,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,(19b)ℒ+(θ2,ϕ2)=L~2=∑j,k‍L~2,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,(19c)ℳ(θ2,ϕ2)=M~2=∑j,k‍M~2,j,k(T,𝒯)ejiΩτcos(kmx)+CC,(19d)𝒩(θ2,ϕ2)=N~2=∑j,k‍N~2,j,k(T,𝒯)ejiΩτcos(kmx)+CC, with j,k∈{0,2} or {j,k}={1,1}.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 θ2→θ~2 such that ℳ(θ~2,ϕ2)=𝒩(θ~2,ϕ2)=0. We define our new function θ~2 by
(20)θ~2=θ2-f(z)M~2-g(z)N~2,
where f(z) and g(z) are two real valued twice-differentiable functions in L2(ℝ) satisfying the following boundary and jump conditions:
(21)limz→±∞f(z)=0,f|z=0-z=0+=0,f(0)=1C,ρcλf′(0+)-f′(0-)=c-1C,limz→±∞g(z)=0,g|z=0-z=0+=0,g(0)=0,ρcλg′(0+)-g′(0-)=1.
After inserting θ~2 into (19a)–(19d) we arrive at a new set of O(ϵ2) equations:(22a)ℒ(θ~2,ϕ2)={ℒ-(θ~2,ϕ2)=l~2-ℒ-(fM~2,0)-ℒ-(gN~2,0),z<0ℒ+(θ~2,ϕ2)=L~2-ℒ+(fM~2,0)-ℒ+(gN~2,0),z<0(22b)ℳ(θ~2,ϕ2)=N(θ~2,ϕ2)=0.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
(23)∂A∂T=0.
Requiring the RHS of (22a) and (22b) to be orthogonal to (16) leads to the solvability condition
(24)∂B∂T=β2|A(𝒯)|2,
and therefore
(25)B(T,𝒯)=β2|A(𝒯)|2T+B~(𝒯).
The constant β2 is defined in Appendix B.


## 2.5. TheO(ϵ3) Problem

After obtaining theO(ϵ2) solutions for θ2 and ϕ2 as described in Appendix C we can insert these into the problem at O(ϵ3), where we find a system of inhomogeneous equations of the form:(26a)ℒ-(θ3,ϕ3)=l~3=∑j,k‍l~3,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,z<0,(26b)ℒ+(θ3,ϕ3)=L~3=∑j,k‍L~3,j,k(z,T,𝒯)ejiΩτcos(kmx)+CC,z>0,(26c)ℳ(θ3,ϕ3)=M~3=∑j,k‍M~3,j,k(T,𝒯)ejiΩτcos(kmx)+CC,(26d)𝒩(θ3,ϕ3)=N~3=∑j,k‍N~3,j,k(T,𝒯)ejiΩτcos(kmx)+CC,where j,k∈{1,3} or {j,k}={0,0}. Proceeding as in the O(ϵ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 A(𝒯)(27)dAd𝒯=iω2A+β3A|A|2.
The constant β3 is defined in Appendix D. Applying the FAT with the time-independent adjoint solution yields the solvability condition(28)B2∂B∂T+B~′(𝒯)=(A2+A-2)β2|A(𝒯)|2.
As ∂B/∂T=β2|A(𝒯)|2, we obtain(29)B~′(𝒯)=(A2+A-2-B2)β2|A(𝒯)|2.
We note that the constants A2 and B2 that multiply the homogeneous solution to the O(ϵ2) 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 Υ=1, n=1, and ν=2.


### 3.1. Limit Cycles forA(𝒯)

From the form of (27), it is evident that there is a fixed point at A=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 |A| 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
(30)A(𝒯)⟶ϱ(𝒯)eiψ(𝒯),
where ϱ,ψ:[0,∞)→ℝ, and we let ω2=ω2R+iω2I and β3=β3R+iβ3I. This form is then inserted into (27). After equating real and imaginary parts, the complex differential equation for A(𝒯) can be transformed into a pair of partially decoupled real differential equations for ϱ(𝒯) and ψ(𝒯):(31a)dϱd𝒯=-ω2Iϱ+β3Rϱ3,(31b)dψd𝒯=ω2R+β3Iϱ2.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 A=0. There is also a possibility of a second fixed point at ϱ^=ω2I/β3R. As we are given ω2I<0, this fixed point exists whenever β3R<0 in which case there is then a corresponding limit cycle with radius ϱ^ which is a global attractor for A(𝒯).Note that (31a) is real and separable for ϱ∉{0,ϱ^} for which an analytical solution is given by
(32)ϱ(𝒯)=(c1e2ω2I𝒯+β3Rω2I)-1/2.
The constant c1 is real and positive and is defined as
(33)c1=1ϱ02-1ϱ^2.
This function can then be substituted into (31b) which can be subsequently integrated to obtain the following:
(34)ψ(𝒯)=(ω2R+β3Iω2Iβ3R)𝒯-β3I2β3Rln(c1e2ω2I𝒯+β3Rω2I)+c2,
where c2 is given in terms of the initial phase ψ0=ψ(0) by
(35)c2=ψ0+β3I2β3Rln(c1+β3Rω2I)=ψ0-β3Iβ3Rln(ϱ0).
From the form of these equations, it is clear that as ϱ→ϱ^, ψ→(ω2R+β3Iω2I/β3R)𝒯+Constant=ω^𝒯+Constant. The limit cycle frequency is then shifted by O(ϵ2) from the neutrally stable frequency: Ω→Ω+ϵ2ω^.


### 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±z-direction. Comparing the values of r and R, one finds the dynamics are stretched in the positive z-direction due to the additional velocity of the gaseous exhaust moving away from the front (i.e., the decay is faster in the -z direction as |Re(R)|<|Re(r)|) and oscillations about this value caused by the imaginary parts of R and r die away. The decay rate has important implications for the numerical solution as it is necessary to greatly extend the numerical z-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 Figure3. 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Θ(x,z,100).Figure 4
Φ
(
x
,
t
) and Θ(x,0,t): numerical interface location and temperature for 90≤t≤100. 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ρ, c, 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 ϵ=C-CΩ 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) amplitudeLimit cycle frequency and radius for varyingρ, λ and c with ϵ=0.05.
(a)
λ
=
c
=
1
(b)
ρ
=
c
=
1
(c)
ρ
=
λ
=
1Figure 6
Numerical and analytical values ofΘ(0,0,t) (smaller amplitude oscillations) and |(∂Φ/∂t)(0,t)| (larger amplitude oscillations) along with their frequency power spectra for ϵ=0.05. In this parameter region, the numerical and analytical solutions are indistinguishable even for large times.Figure 7
Numerical and analytical values ofΘ(0,0,t) and |(∂Φ/∂t)(0,t)| along with their frequency power spectra for ϵ=0.1. 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 O(ϵ2).To explore the dynamics in a region where the analysis predicts a subcritical Hopf bifurcation, we choose to vary the nondimensional heat capacity ratioc. From Figure 5, we expect that as we decrease c from c=1 for fixed ρ=λ=1, 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 c=1,0.9, 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 c 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.Θ
(
0,0
,
t
) and |(∂Φ/∂t)(0,t)| for decreasing c with ρ=λ=1 and ϵ=0.05. For c=0.85, the analytical solution has become unbounded and is not plotted.
(a)
c
=
0.95
(b)
c
=
0.85


## 3.1. Limit Cycles forA(𝒯)

From the form of (27), it is evident that there is a fixed point at A=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 |A| 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
(30)A(𝒯)⟶ϱ(𝒯)eiψ(𝒯),
where ϱ,ψ:[0,∞)→ℝ, and we let ω2=ω2R+iω2I and β3=β3R+iβ3I. This form is then inserted into (27). After equating real and imaginary parts, the complex differential equation for A(𝒯) can be transformed into a pair of partially decoupled real differential equations for ϱ(𝒯) and ψ(𝒯):(31a)dϱd𝒯=-ω2Iϱ+β3Rϱ3,(31b)dψd𝒯=ω2R+β3Iϱ2.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 A=0. There is also a possibility of a second fixed point at ϱ^=ω2I/β3R. As we are given ω2I<0, this fixed point exists whenever β3R<0 in which case there is then a corresponding limit cycle with radius ϱ^ which is a global attractor for A(𝒯).Note that (31a) is real and separable for ϱ∉{0,ϱ^} for which an analytical solution is given by
(32)ϱ(𝒯)=(c1e2ω2I𝒯+β3Rω2I)-1/2.
The constant c1 is real and positive and is defined as
(33)c1=1ϱ02-1ϱ^2.
This function can then be substituted into (31b) which can be subsequently integrated to obtain the following:
(34)ψ(𝒯)=(ω2R+β3Iω2Iβ3R)𝒯-β3I2β3Rln(c1e2ω2I𝒯+β3Rω2I)+c2,
where c2 is given in terms of the initial phase ψ0=ψ(0) by
(35)c2=ψ0+β3I2β3Rln(c1+β3Rω2I)=ψ0-β3Iβ3Rln(ϱ0).
From the form of these equations, it is clear that as ϱ→ϱ^, ψ→(ω2R+β3Iω2I/β3R)𝒯+Constant=ω^𝒯+Constant. The limit cycle frequency is then shifted by O(ϵ2) from the neutrally stable frequency: Ω→Ω+ϵ2ω^.


## 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±z-direction. Comparing the values of r and R, one finds the dynamics are stretched in the positive z-direction due to the additional velocity of the gaseous exhaust moving away from the front (i.e., the decay is faster in the -z direction as |Re(R)|<|Re(r)|) and oscillations about this value caused by the imaginary parts of R and r die away. The decay rate has important implications for the numerical solution as it is necessary to greatly extend the numerical z-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 Figure3. 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Θ(x,z,100).Figure 4
Φ
(
x
,
t
) and Θ(x,0,t): numerical interface location and temperature for 90≤t≤100. 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ρ, c, 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 ϵ=C-CΩ 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) amplitudeLimit cycle frequency and radius for varyingρ, λ and c with ϵ=0.05.
(a)
λ
=
c
=
1
(b)
ρ
=
c
=
1
(c)
ρ
=
λ
=
1Figure 6
Numerical and analytical values ofΘ(0,0,t) (smaller amplitude oscillations) and |(∂Φ/∂t)(0,t)| (larger amplitude oscillations) along with their frequency power spectra for ϵ=0.05. In this parameter region, the numerical and analytical solutions are indistinguishable even for large times.Figure 7
Numerical and analytical values ofΘ(0,0,t) and |(∂Φ/∂t)(0,t)| along with their frequency power spectra for ϵ=0.1. 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 O(ϵ2).To explore the dynamics in a region where the analysis predicts a subcritical Hopf bifurcation, we choose to vary the nondimensional heat capacity ratioc. From Figure 5, we expect that as we decrease c from c=1 for fixed ρ=λ=1, 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 c=1,0.9, 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 c 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.Θ
(
0,0
,
t
) and |(∂Φ/∂t)(0,t)| for decreasing c with ρ=λ=1 and ϵ=0.05. For c=0.85, the analytical solution has become unbounded and is not plotted.
(a)
c
=
0.95
(b)
c
=
0.85


## 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ϵ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.



---

*Source: 580561-2013-05-16.xml*