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

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 𝑢_𝑔=(1−(𝜌_𝑠/𝜌_𝑔))(𝜕Φ_𝑠/𝜕𝑡)̃ 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: 𝑈𝑥=𝜆_𝑠𝑈𝑥,𝑦=𝜆_𝑠𝑈𝑦,𝑧=𝜆_𝑠𝑧−Φ_𝑆𝑈,𝑡=²𝜆_𝑠̃𝑇𝑡,Θ=𝑇−_𝑢𝑇_𝑎−𝑇_𝑢𝑈,Φ=𝜆_𝑠Φ_𝑆.(1) 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 Δ=(1−(𝑇_𝑢/𝑇_𝑎))(𝐸_𝑠/𝑅_𝑔𝑇_𝑎)≫1, 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 𝑧=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𝜕Θ−𝜕𝑡𝜕Φ𝜕𝑡𝜕Θ=∇𝜕𝑧²Θ,𝑧<0,(2a)𝜌𝑐𝜕Θ𝜕𝑡−𝑐𝜕Φ𝜕𝑡𝜕Θ∇𝜕𝑧=𝜌𝑐𝜆²Θ,𝑧>0,(2b) the jump conditions 𝜕Φ𝜕𝑡=−1+𝜕Φ𝜕𝑥²+𝜕Φ𝜕𝑦^21/2exp𝐹Θ|_𝑧=0−𝐹(1),(2c)𝜌𝑐𝜆𝜕Θ𝜕𝑧_𝑧=0⁺−𝜕Θ𝜕𝑧_𝑧=0⁻=−1+𝜕Φ𝜕𝑥²+𝜕Φ𝜕𝑦²⁻¹(𝑐−1)Θ|_𝑧=0−𝑐𝜕Φ,𝜕𝑡(2d) and the 𝑧-boundary and continuity conditions lim_{𝑧→−∞}Θ=0,lim_𝑧→∞Θ=1,lim_𝑧→0+Θ=lim_𝑧→0−Θ.(2e) 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,𝐿𝑤}=0,for𝑤=𝑥,𝑦,𝜂=Θ,Φ.(2f)The operator ∇² appearing in (2a) and (2b) is the Laplacian in the moving coordinate system and is given by ∇²=𝜕²𝜕𝑥²+𝜕²𝜕𝑦²+1+𝜕Φ𝜕𝑥²+𝜕Φ𝜕𝑦²𝜕²𝜕𝑧²−2𝜕Φ𝜕𝜕𝑥𝜕𝑥𝜕²𝑧−2𝜕Φ𝜕𝜕𝑦²−𝜕𝜕𝑦𝜕𝑧₂Φ𝜕𝑥²+𝜕²Φ𝜕𝑦²𝜕.𝜕𝑧(3) The function 𝐹(Θ|_𝑧=0) appearing in (2c) is determined implicitly from Υexp𝜈Δ2Θ|_𝑧=0−1×𝐺ΔΘ|_𝑧=0−1−𝐹^1/2=exp(𝐹),(4) where 𝐺 is defined as 𝐺(𝑢)=^𝑢₀𝑒^−𝜈𝑤𝑤^𝑛𝑑𝑤(5) 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” 𝜆₀=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].

1.3. Linear Stability of the Basic Solution

An immediate advantage over the dimensional model is that (2a)–(2f) admits a basic solutionΘ=𝑒𝑥,𝑦,𝑧,𝑡^𝑧1Φfor𝑧≶0,(6a)𝑥,𝑦,𝑡=−𝑡(6b)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 Θ𝑥,𝑦,𝑧,𝑡=𝜃₀(𝑧)+𝜖𝜃₁,Φ𝑥,𝑦,𝑥,𝑡𝑥,𝑦,𝑡=𝜙₀(𝑡)+𝜖𝜙₁𝑥,𝑦,𝑡,(7) 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 𝐹(Θ|_𝑧=0) in powers of 𝜖 is discussed in Appendix A where the functions 𝐹_𝑖 at 𝑂(𝜖^𝑖) are expressed in terms of {𝜃_𝑗} where 𝑗≤𝑖. In particular, we note that 𝐹₀=𝐹(1)=−ln𝜆₀ and 𝐹₁=𝐶×𝜃₁ 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 𝑂(𝜖):ℒ₋𝜃₁,𝜙₁=𝜕𝜃₁+𝜕𝑡𝜕𝜃₁𝜕𝑧−∇²𝜃₁+𝑒^𝑧𝜕²𝜙₁𝜕𝑥²+𝜕²𝜙₁𝜕𝑦²−𝜕𝜙₁ℒ𝜕𝑡=0,𝑧<0,(8a)₊𝜃₁,𝜙₁=𝜌𝑐𝜕𝜃₁𝜕𝑡+𝑐𝜕𝜃₁𝜕𝑧−𝜌𝑐𝜆∇²𝜃₁ℳ𝜃=0,𝑧>0,(8b)₁,𝜙₁=𝜕𝜙₁𝜕𝑡+𝐶×𝜃_1𝑧=0𝒩𝜃=0,(8c)₁,𝜙₁=𝜌𝑐𝜆𝜕𝜃₁𝜕𝑧_𝑧=0⁺−𝜕𝜃₁𝜕𝑧_𝑧=0⁻−(𝑐−1)𝜃_1𝑧=0−𝜕𝜙₁𝜃𝜕𝑡=0,(8d)_{1𝑧=0+𝑧=0−}=0,lim_𝑧→±∞𝜃₁=0,(8e)𝜕𝜂𝜕𝑤_{𝑤=0,𝐿𝑤}=0,for𝑤=𝑥,𝑦,𝜂=𝜃,𝜙,(8f)where the differential operator ∇² denotes the traditional 3D Laplacian: (𝜕²/𝜕𝑥²)+(𝜕²/𝜕𝑦²)+(𝜕²/𝜕𝑧²).

Separable solutions to (8a)–(8c) and (8e)–(8f) are given by𝜃₁=𝐴𝑥,𝑦,𝑧,𝑡𝛾𝑒^𝑧+Γ𝑒^𝑟𝑧𝑒^𝑖𝜔𝑡cos(𝑚𝑥)cos𝑛𝑦+𝐶𝐶𝐴𝑒^𝑅𝑧𝑒^𝑖𝜔𝑡()𝜙cos𝑚𝑥cos𝑛𝑦+𝐶𝐶for𝑧≶0,(9a)₁𝑥,𝑦,𝑡=𝑖𝐶𝐴𝜔𝑒^𝑖𝜔𝑡()+1cos𝑚𝑥cos𝑛𝑦2𝐵+𝐶𝐶,(9b)where 𝑚,𝑛∈ℝ, 𝜔∈ℤ, 𝛾=𝑖𝐶/𝜔 and Γ=1−(𝑖𝐶/𝜔). 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 by1𝑟=21+1+4𝑖𝜔+𝜅²,1𝑅=2𝜌𝜆1−1+4𝜌²𝜆²𝑖𝜔𝜆+𝜅²,(10) 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 𝒟(𝜔,𝜅;𝐶)=𝜌𝑐𝜆𝑅−𝑖𝐶𝜔−𝑟1−𝑖𝐶𝜔+𝐶+1−𝑐=0.(11) 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(𝒟)=Im(𝒟)=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.

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 𝐶=𝐶_Ω(𝜅_min)+𝜖², 𝜖≪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 𝜔=Ω+𝜖²𝜔₂+𝑂(𝜖³), 𝜔₂∈ℤ. 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 𝜖: Θ∼^∞_𝑖=0𝜖^𝑖𝜃_𝑖(𝑥,𝑧,𝜏,𝑇,𝒯),Φ∼^∞_𝑖=0𝜖^𝑖𝜙_𝑖(𝑥,𝜏,𝑇,𝒯),(12) 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 𝑂(1) are satisfied while the equations at 𝑂(𝜖^𝑗) for 𝑗≥1 are given byℒ₋𝜃_𝑗,𝜙_𝑗=𝑙_𝑗̃𝜃_𝑗−1,𝜃_𝑗−2,…,𝜃₀,𝜙_𝑗−1,𝜙_𝑗−2,…,𝜙₀ℒ,𝑧<0,(13a)₊𝜃_𝑗,𝜙_𝑗=𝐿_𝑗𝜃_𝑗−1,𝜃_𝑗−2,…,𝜃₀,𝜙_𝑗−1,𝜙_𝑗−2,…,𝜙₀ℳ𝜃,𝑧>0,(13b)_𝑗,𝜙_𝑗=𝑀_𝑗𝜃_𝑗−1,𝜃_𝑗−2,…,𝜃₀,𝜙_𝑗−1,𝜙_𝑗−2,…,𝜙₀,𝒩𝜃(13c)_𝑗,𝜙_𝑗=𝑁_𝑗𝜃_𝑗−1,𝜃_𝑗−2,…,𝜃₀,𝜙_𝑗−1,𝜙_𝑗−2,…,𝜙₀,(13d)𝜕𝜃_𝑗𝜕𝑥_{𝑥=0,𝐿𝑥}=𝜕𝜙𝜕𝑥_{𝑥=0,𝐿𝑥}=0,lim_𝑧→±∞𝜃_𝑗=0=𝜃_𝑗^𝑧=0+_𝑧=0−.(13e)The 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 𝑚≠𝜅_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)+𝜖²) 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Θ()=𝑒𝑥,𝑧,𝜏,𝑇,𝒯_𝑧+𝜖𝐴(𝑇,𝒯)𝛾𝑒^𝑧+Γ𝑒^𝑟𝑧𝑒^𝑖Ω𝜏()𝐴()𝑒cos𝑚𝑥+𝐶𝐶1+𝜖𝑇,𝒯^𝑅𝑧𝑒^𝑖Ω𝜏cos(𝑚𝑥)+𝐶𝐶,(14a) for 𝑧≶0 and Φ(𝑥,𝜏,𝑇,𝒯)=−𝜏+𝜖𝛾𝑒^𝑖Ω𝜏cos(𝑚𝑥)+𝐶𝐶+𝐵(𝑇,𝒯).(14b)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.3. 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 𝐴(𝑇,𝒯) 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 ker(ℒ^†) can be found in [24]. The inner product in the range space of ℒ (which contains the domain of ℒ^†) is defined by ⟨𝑢,𝑣⟩=lim_𝑡→∞1𝑡𝐿_𝑥^𝑡₀^∞_−∞^𝐿𝑥₀𝑢𝑣d𝑥d𝑧d𝜏,(15) where we have assumed 𝐿_𝑦→0. Separable solutions to ℒ^†(𝜃^†)=0 are of the form𝜃^†=𝑒𝑥,𝑦,𝑧,𝜏^𝑟†𝑧𝑒^𝑅†𝑧𝑒^{𝑖𝜔†𝜏}cos(𝑚𝑥)for𝑧≶0,(16) where 𝑟^†=12−1+1+4𝑖𝜔^†+𝑚²,𝑅^†=12𝜌𝜆−1−1+4𝜌²𝜆²𝑖𝜔^†𝜆+𝑚²,(17) and 𝜔^†=𝜔. There also exists a planar, time-independent solution with 𝑚=𝜔^†=0 given by 𝜃^†₀(𝑧)=𝑒1for𝑧<0_{−(1/𝜌𝜆)𝑧}for𝑧>0.(18) 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: ℒ₋𝜃₂,𝜙₂=𝑙₂̃=_𝑗,𝑘𝑙̃_2,𝑗,𝑘(𝑧,𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()ℒcos𝑘𝑚𝑥+𝐶𝐶,(19a)₊𝜃₂,𝜙₂=𝐿₂=_𝑗,𝑘𝐿_2,𝑗,𝑘(𝑧,𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()ℳ𝜃cos𝑘𝑚𝑥+𝐶𝐶,(19b)₂,𝜙₂=𝑀₂=_𝑗,𝑘𝑀_2,𝑗,𝑘(𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()𝒩𝜃cos𝑘𝑚𝑥+𝐶𝐶,(19c)₂,𝜙₂=𝑁₂=_𝑗,𝑘𝑁_2,𝑗,𝑘(𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}cos(𝑘𝑚𝑥)+𝐶𝐶,(19d) with 𝑗,𝑘∈{0,2} or {𝑗,𝑘}={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 𝜃₂→𝜃₂ such that ℳ(𝜃₂,𝜙₂)=𝒩(𝜃₂,𝜙₂)=0. We define our new function 𝜃₂ by 𝜃₂=𝜃₂−𝑓(𝑧)𝑀₂−𝑔(𝑧)𝑁₂,(20) where 𝑓(𝑧) and 𝑔(𝑧) are two real valued twice-differentiable functions in 𝐿²(ℝ) satisfying the following boundary and jump conditions: lim_𝑧→±∞𝑓(𝑧)=0,𝑓|^𝑧=0+_𝑧=0⁻(0)=1=0,𝑓𝐶,𝜌𝑐𝜆𝑓^′0⁺−𝑓^′(0⁻)=𝑐−1𝐶,lim_𝑧→±∞𝑔(𝑧)=0,𝑔|^𝑧=0+_𝑧=0⁻=0,𝑔(0)=0,𝜌𝑐𝜆𝑔^′0⁺−𝑔^′(0⁻)=1.(21) After inserting 𝜃₂ into (19a)–(19d) we arrive at a new set of 𝑂(𝜖²) equations:ℒ𝜃₂,𝜙₂=ℒ₋𝜃₂,𝜙₂=𝑙₂̃−ℒ₋𝑓𝑀₂,0−ℒ₋𝑔𝑁₂,ℒ,0𝑧<0₊𝜃₂,𝜙₂=𝐿₂−ℒ₊𝑓𝑀₂,0−ℒ₊𝑔𝑁₂,ℳ𝜃,0𝑧<0(22a)₂,𝜙₂=𝑁𝜃₂,𝜙₂=0.(22b)

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 𝜕𝐴𝜕𝑇=0.(23) Requiring the RHS of (22a) and (22b) to be orthogonal to (16) leads to the solvability condition 𝜕𝐵𝜕𝑇=𝛽₂|𝐴(𝒯)|²,(24) and therefore 𝐵(𝑇,𝒯)=𝛽₂|𝐴(𝒯)|²𝑇+𝐵(𝒯).(25) 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:ℒ₋𝜃₃,𝜙₃=𝑙₃̃=_𝑗,𝑘𝑙̃_3,𝑗,𝑘(𝑧,𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()ℒcos𝑘𝑚𝑥+𝐶𝐶,𝑧<0,(26a)₊𝜃₃,𝜙₃=𝐿₃=_𝑗,𝑘𝐿_3,𝑗,𝑘(𝑧,𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()ℳ𝜃cos𝑘𝑚𝑥+𝐶𝐶,𝑧>0,(26b)₃,𝜙₃=𝑀₃=_𝑗,𝑘𝑀_3,𝑗,𝑘(𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}()𝒩𝜃cos𝑘𝑚𝑥+𝐶𝐶,(26c)₃,𝜙₃=𝑁₃=_𝑗,𝑘𝑁_3,𝑗,𝑘(𝑇,𝒯)𝑒^{𝑗𝑖Ω𝜏}cos(𝑘𝑚𝑥)+𝐶𝐶,(26d)where 𝑗,𝑘∈{1,3} or {𝑗,𝑘}={0,0}. 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 𝐴(𝒯)𝑑𝐴𝑑𝒯=𝑖𝜔₂𝐴+𝛽₃𝐴|𝐴|².(27) The constant 𝛽₃ is defined in Appendix D. Applying the FAT with the time-independent adjoint solution yields the solvability condition𝐵₂𝜕𝐵+𝐵𝜕𝑇^′(𝒯)=𝐴₂+𝐴₂𝛽₂|𝐴(𝒯)|².(28) As 𝜕𝐵/𝜕𝑇=𝛽₂|𝐴(𝒯)|², we obtain𝐵^′(𝒯)=𝐴₂+𝐴₂−𝐵₂𝛽₂|𝐴(𝒯)|².(29) 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 Υ=1, 𝑛=1, and 𝜈=2.

3.1. Limit Cycles for 𝐴(𝒯)

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(𝜔₂)<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 𝐴(𝒯)⟶𝜚(𝒯)𝑒^{𝑖𝜓(𝒯)},(30) where 𝜚,𝜓∶[0,∞)→ℝ, 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 𝜓(𝒯):𝑑𝜚𝑑𝒯=−𝜔^𝐼₂𝜚+𝛽^𝑅₃𝜚³,(31a)𝑑𝜓𝑑𝒯=𝜔^𝑅₂+𝛽^𝐼₃𝜚².(31b)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 𝜚̂=𝜔^𝐼₂/𝛽^𝑅₃. As we are given 𝜔^𝐼₂<0, this fixed point exists whenever 𝛽^𝑅₃<0 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 𝜚∉{0,𝜚̂} for which an analytical solution is given by 𝜚(𝒯)=𝑐₁𝑒^{2𝜔𝐼2𝒯}+𝛽^𝑅₃𝜔^𝐼₂^−1/2.(32) The constant 𝑐₁ is real and positive and is defined as 𝑐₁=1𝜚²₀−1𝜚̂².(33) This function can then be substituted into (31b) which can be subsequently integrated to obtain the following: 𝜓(𝒯)=𝜔^𝑅₂+𝛽^𝐼₃𝜔^𝐼₂𝛽^𝑅₃𝛽𝒯−^𝐼32𝛽^𝑅₃ln𝑐₁𝑒^{2𝜔𝐼2𝒯}+𝛽^𝑅₃𝜔^𝐼₂+𝑐₂,(34) where 𝑐₂ is given in terms of the initial phase 𝜓₀=𝜓(0) by 𝑐₂=𝜓₀+𝛽^𝐼₃2𝛽^𝑅₃ln𝑐₁+𝛽^𝑅₃𝜔^𝐼₂=𝜓₀−𝛽^𝐼₃𝛽^𝑅₃ln𝜚₀.(35) From the form of these equations, it is clear that as 𝜚→𝜚̂, 𝜓→(𝜔^𝑅₂+𝛽^𝐼₃𝜔^𝐼₂/𝛽^𝑅₃)𝒯+Constant=𝜔𝒯+Constant. 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 |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 Figure 4.

Figure 3 

Numerical solution for Θ(𝑥,𝑧,100).

Figure 4 

Φ(𝑥,𝑡) and Θ(𝑥,0,𝑡): numerical interface location and temperature for 90≤𝑡≤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 𝜌, 𝑐, 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) 𝜆=𝑐=1

(b) 𝜌=𝑐=1

(c) 𝜌=𝜆=1

Figure 5 

Limit cycle frequency and radius for varying 𝜌, 𝜆 and 𝑐 with 𝜖=0.05.

Figure 6 

Numerical and analytical values of Θ(0,0,𝑡) (smaller amplitude oscillations) and |(𝜕Φ/𝜕𝑡)(0,𝑡)| (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,𝑡) and |(𝜕Φ/𝜕𝑡)(0,𝑡)| 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 𝑂(𝜖2).

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 Figures 7, 8(a), and 8(b) which correspond to 𝑐=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 𝑐 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) 𝑐=0.95

(b) 𝑐=0.85

Figure 8 

Θ(0,0,𝑡) and |(𝜕Φ/𝜕𝑡)(0,𝑡)| for decreasing 𝑐 with 𝜌=𝜆=1 and 𝜖=0.05. For 𝑐=0.85, 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 𝜖²≪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.

Appendices

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 𝐹₁=𝐶×𝜃_1𝑧=0,(A.1) where the parameter grouping 𝐶 is defined as𝐶=Δ𝜈+Υ²ln𝜆₀^𝑛𝜆₀^2−𝜈2+Υ²ln𝜆₀^𝑛𝜆₀^2−𝜈.(A.2) Note that constant 𝐶 is proportional to the nondimensional activation energy Δ. In the special case when 𝜈=2, we have, 𝐶=Δ. At 𝑂(𝜖²), we have 𝐹₂=𝐶𝜃_2𝑧=0+𝐶₂𝜃_1𝑧=0²,(A.3) where 𝐶₂=−4𝐶²+Δ²𝜈²4+2Υ²ln𝜆₀^𝑛𝜆₀^2−𝜈+𝑛(Δ−𝐶)²ln𝜆₀^𝑛−1+𝜈Δ²−𝐶²ln𝜆₀^𝑛Υ²𝜆₀^2−𝜈4+2Υ²ln𝜆₀^𝑛𝜆₀^2−𝜈.(A.4) At 𝑂(𝜖³) we find that the function 𝐹₃ is given by 𝐹₃=𝐶𝜃₃|_𝑧=0+2𝐶₂𝜃₁𝜃₂|_𝑧=0+𝐶₃𝜃₁|_𝑧=0³,(A.5) where 𝐶₃=16ln𝜆₀^𝑛−2𝑛−𝑛²(𝐶−Δ)³−ln𝜆₀^𝑛𝐶𝐶₂1𝜈+6𝜈²𝐶³−Δ³+ln𝜆₀^𝑛−1𝑛𝐶₂()+1𝐶−Δ6𝜈𝑛2𝐶³−3𝐶²Δ+Δ³Υ²𝜆₀^2−𝜈−4𝐶𝐶₂−43𝐶³+16𝜈³Δ³×12+Υ²ln𝜆₀^𝑛𝜆₀^2−𝜈.(A.6) Note that for 𝜈=2, 𝐶₂=𝐶₃=0.

B. Terms from the 𝑂(𝜖²) Problem

The constant 𝛽₂ appearing in (25) is defined 𝛽₂=𝐶𝑐𝜈_2,0,0−1𝑟𝑢_2,0,0−𝑈_2,0,0+𝜇1/𝜌𝜆−𝑅^2,0,0+𝐶𝐶(B.1) with 𝑢_2,𝑗,𝑘=12𝑟𝛾Γ(−1)^{(𝑗−𝑘)/2}(𝑘+1)𝑚²,()𝑢+𝑖ΩB.2a_2,𝑗,𝑘=12(−1)^(2−𝑗)/2𝛾²𝑘𝑚²+𝑗2,()𝑈𝑖ΩB.2b_2,𝑗,𝑘=12𝜌𝑐𝜆𝛾𝑅𝑖Ω+()𝜌𝜆−1^{(𝑗−𝑘)/2}(𝑘+1)𝑚²,()𝜇B.2c_2,𝑗,𝑘=14−𝐶²−2𝐶₂+𝛾²(−1)^{(𝑗−𝑘)/2}𝑚²,()𝜈B.2d_2,𝑗,𝑘=12−(𝑐−1)𝑖Ω𝛾−𝛾²(−1)^{(𝑗−𝑘)/2}𝑚².(B.2e)

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: 𝜃_2,𝑃=_𝑗∑_,𝑘𝑎_𝑗,𝑘𝑒^𝑟𝑧+𝑎_𝑗,𝑘𝜒_𝑗,𝑘(𝑧)+𝑎̌_𝑗,𝑘𝑒^̌𝑟𝑧𝐴^(2+𝑗)/2𝐴^(2−𝑗)/2×𝑒^{𝑗𝑖Ω𝜏}()∑cos𝑘𝑚𝑥+𝐶𝐶_𝑗,𝑘𝛼_𝑗,𝑘𝑒^𝑅𝑧+𝛼̌_𝑗,𝑘𝑒^̌𝑅𝑧𝐴^(2+𝑗)/2𝐴^(2−𝑗)/2×𝑒^{𝑗𝑖Ω𝜏}()()𝜙cos𝑘𝑚𝑥+𝐶𝐶for𝑧≶0,C.1a_2,𝑃=_𝑗,𝑘𝑏_𝑗,𝑘𝐴^(2+𝑗)/2𝐴^(2−𝑗)/2𝑒^{𝑗𝑖Ω𝜏}cos(𝑘𝑚𝑥)+𝐶𝐶,(C.1b) where 𝑗,𝑘∈{0,2} and𝜒_𝑗,𝑘(𝑧)=𝑧𝑒^𝑧𝑒for𝑗=𝑘=0,_𝑧otherwise.(C.1) The functions 𝑒^{𝑗𝑖Ω𝜏}𝑒^{̌𝑞𝑗,𝑘𝑧}cos(𝑘𝑚𝑥) lie in the null space of ℒ_± with 𝑞̌_𝑗,𝑘=𝑟̌_𝑗,𝑘 or 𝑅̌_𝑗,𝑘 given by 𝑟̌_𝑗,𝑘=121+1+4𝑗𝑖Ω+𝑘²𝑚²,𝑅̌_𝑗,𝑘=12𝜌𝜆1−1+4𝜌𝜆²𝑗𝑖Ω𝜆+𝑘²𝑚².(C.2) As 𝑅̌_0,0=0, we must choose 𝛼̌_0,0=0. Furthermore, 𝑟̌_0,0=1 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 𝑧=0. Upon inspection, we find that 𝑏_0,0 does not enter into any of these equations, though as the functional form 𝑏_0,0𝑒^{(0)𝑖Ω𝜏} is identical to that of the homogeneous solution, we may choose 𝑏_0,0=0. Equating terms in the ℒ_± equations yields 9 solutions and 3 equations: 𝑎_𝑗,𝑘=𝑢_2,𝑗,𝑘𝑠_𝑗,𝑘,𝛼_𝑗,𝑘=𝑈_2,𝑗,𝑘𝑆_𝑗,𝑘{},𝑎for𝑗,𝑘∈0,2_0,0=𝑢_2,0,0𝑠̂_0,0=−𝑢_2,0,0,(C.3) where the {2,𝑗,𝑘} constants without tildes are then defined similarly to their counterparts in (B.2a)–(B.2e) as: 𝑢_2,𝑗,𝑘=𝑢_2,𝑗,𝑘,𝑢_2,𝑗,𝑘=𝑢_2,𝑗,𝑘+𝛿_𝑗,𝑘,0𝛽₂2,𝜇_2,𝑗,𝑘=𝜇_2,𝑗,𝑘−𝛿_𝑗,𝑘,0𝛽₂2,𝑈_2,𝑗,𝑘=𝑈_2,𝑗,𝑘,𝜈_2,𝑗,𝑘=𝜈_2,𝑗,𝑘+𝛿_𝑗,𝑘,0𝛽₂2,𝑠_𝑗,𝑘=𝑗𝑖Ω+𝑟−𝑟²+𝑘²𝑚²,𝑠̂_𝑗,𝑘=−𝛿_𝑗,𝑘,0+𝑗𝑖Ω+𝑘²𝑚²,𝑆_𝑗,𝑘=𝜌𝑐𝑗𝑖Ω+𝑐𝑅−𝜌𝑐𝜆𝑅²+𝜌𝑐𝜆𝑘²𝑚².(C.4) We also obtain 3 equations: 𝑎_𝑗,𝑘−𝑏_𝑗,𝑘=𝑢_2,𝑗,𝑘𝑠̂_𝑗,𝑘,for𝑗,𝑘={0,2},{2,0},or{2,2}.(C.5) As 𝛼̌_0,0=0 is known, we find 𝑎̌_0,0=𝛼_0,0−𝑎_0,0 for 𝑗=𝑘=0, otherwise, we obtain three equations from the continuity condition: 𝛼̌_𝑗,𝑘−𝑎̌_𝑗,𝑘−𝑎_𝑗,𝑘=𝑎_𝑗,𝑘−𝛼_𝑗,𝑘.(C.6) The ℳ and 𝒩 equations yield six more equations: 𝑗𝑖Ω𝑏_𝑗,𝑘+𝐶𝛼̌_𝑗,𝑘=𝜇_2,𝑗,𝑘−𝐶𝛼_𝑗,𝑘𝑅̌(C.7)𝜌𝑐𝜆_𝑗,𝑘−(𝑐−1)𝛼̌_𝑗,𝑘−𝑎_𝑗,𝑘−𝑟̌_𝑗,𝑘𝑎̌_𝑗,𝑘−𝑗𝑖Ω𝑏_𝑗,𝑘=𝜈_2,𝑗,𝑘+(𝑐−1)−𝜌𝑐𝜆𝑅𝛼_𝑗,𝑘+𝑟𝑎_𝑗,𝑘,(C.8) for {𝑗,𝑘}≠{0,0}. 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 𝛽₃𝑢=−^3,1,1+𝑈_3,1,1+𝜌𝑐𝜆𝑅^†−𝑟^†/𝐶𝜇_3,1,1−𝜈_3,1,1𝜌𝑐/𝑅+𝑅^†−Γ/𝑟+𝑟^†−𝛾𝜌𝑐𝜆𝑅^†−𝑟^†/𝐶+1(D.1) where 𝑢_3,1,1(𝑧)=𝑟𝛽₂−12𝑟²𝛾²𝑚²Γ𝑒^𝑟𝑧+𝑖Ω𝑟2𝑏_2,0+𝑏_2,2+14𝛾²𝑟²𝑚²Γ𝑒^𝑟𝑧+2𝛾𝑚²𝑏_0,2+𝑏_0,2−𝑏_2,2−38𝛾²−𝑖Ω𝛾2𝑏_2,0+𝑏_2,2+𝛾𝛽₂𝑒^𝑧+_𝑗,𝑘𝛾2^𝑘/2(−1)^{(𝑗−𝑘)/2}(𝑘+1)𝑚²×+𝑖Ω𝑟𝑎_𝑗,𝑘𝑒^𝑟𝑧+𝑟̌_𝑗,𝑘𝑎̌_𝑗,𝑘𝑒^{̌𝑟𝑗,𝑘𝑧}+𝑎_𝑗,𝑘1+𝛿_𝑗,𝑘,0𝑧𝑒^𝑧+_𝑗,𝑘𝛿_𝑗,0𝛾2^𝑘/2(−1)^{(𝑗−𝑘)/2}(𝑘+1)𝑚²×+𝑖Ω𝑎_𝑗,𝑘𝑟𝑒^𝑟𝑧+𝑎̌_𝑗,𝑘𝑟̌_𝑗,𝑘𝑒^{̌𝑟𝑗,𝑘𝑧}+𝑎_𝑗,𝑘1+𝛿_𝑗,𝑘,0𝑧𝑒^𝑧,()𝑈D.2a_3,1,1(𝑧)=𝑅𝑐𝛽₂−12𝜌𝑐𝜆𝛾²𝑅²𝑚²𝑒^𝑅𝑧+𝑖𝑐Ω𝑅2𝑏_2,0+𝑏_2,2+14𝜌𝑐𝜆𝛾²𝑅²𝑚²𝑒^𝑅𝑧+_𝑗,𝑘𝛾2^𝑘/2𝑅𝛼_𝑗,𝑘𝑒^𝑅𝑧+𝛿_𝑗,0𝛼_𝑗,𝑘𝑅𝑒^𝑅𝑧×𝜌𝑐𝜆(−1)^{(𝑗−𝑘)/2}(𝑘+1)𝑚²+𝑖𝑐Ω+_𝑗,𝑘𝛾2^𝑘/2𝑅̌_𝑗,𝑘𝛼̌_𝑗,𝑘𝑒^{̌𝑅𝑗,𝑘𝑧}+𝛿_𝑗,0𝛼̌_𝑗,𝑘𝑅̌_𝑗,𝑘𝑒^{𝑅𝑗,𝑘𝑧}×𝜌𝑐𝜆(−1)^{(𝑗−𝑘)/2}(𝑘+1)𝑚²,()𝜇+𝑖𝑐ΩD.2b_3,1,1=18𝛾²𝐶+𝛾𝑏_2,2−𝑏_0,2−𝑏_0,2𝑚²−38𝐶³−94𝐶𝐶₂−94𝐶₃−𝐶²+2𝐶_2𝑗,𝑘2^−𝑘/2𝛼_𝑗,𝑘+𝛼̌_𝑗,𝑘+𝛿_𝑗,0𝛼_𝑗,𝑘+𝛼̌_𝑗,𝑘,()𝜈D.2c_3,1,1=(1𝑐−1)4𝛾²𝑚²−𝛽₂−𝑖Ω2𝑏_2,0+𝑏_2,2+2𝛾𝑚²𝑏_0,2+𝑏_0,2−𝑏_2,2+18𝑖Ω𝛾²−𝑖Ω𝛾(𝑐−1)_𝑗,𝑘2^−𝑘/2𝛼_𝑗,𝑘+𝛼̌_𝑗,𝑘+𝛿_𝑗,0𝛼_𝑗,𝑘+𝛼̌_𝑗,𝑘.(D.2d) In these summations, 𝑗,𝑘∈{0,2}.