HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS

JOHN HORNIBROOK^{∗}, SANJEEVA BALASURIYA^{†}, AND STÉPHANE LAFORTUNE^{‡}

Abstract. The wavespeed and stability of wavefronts associated with a one-dimensional combustion model with Arrhenius kinetics and no heat loss are analyzed. The focus is on the singular limit of very large Lewis number, in which fuel diffusivity is small in comparison to that of heat. Many of the established results for the infinite Lewis number are recovered, and an empirical wavespeed formula of excellent accuracy is determined. An Evans function technique is used to verify that the linear operator arising from the linearization about this wavefront solution does not possess any eigenvalues of positive real part, thereby supporting well-established numerical evidence on the stability of the infinite Lewis number front. In the very large (but not infinite) Lewis number instance, a similarly detailed assessment of the wavespeed is obtained. However, the Evans function method shows that such wavefronts are inherently unstable.

Key words. Combustion waves, high Lewis number, wavespeed, stability, Evans function, Melnikov’s method

AMS subject classifications. 80A25, 35K57, 35B35, 34E10, 34C37

1. Introduction. In this article, we study the wavespeed and stability of a combustion wavefront along a one-dimensional medium. This is a fundamental idealized problem towards understanding how flame fronts propagate, and therefore has received a considerable amount of attention. There are several (non-dimensional) parameters which apparently are important: the Lewis number Le, the exothermicity parameter β, and the heat loss parameter ℓ. The first of these, the Lewis number, measures the relative importance of fuel diffusivity in comparison to that of heat. The exothermicity β is the ratio of the activation energy to the heat of reaction. The structure of the governing equations is such that the infinite Lewis number instance is considerably easier to deal with than allowing for fuel diffusivity. Many studies of this “solid” regime appear in the literature [4,6,22,27,28], and also the “gaseous” regime Le ≈ 1 has been frequently studied because of a symmetry in the equations [6,15,28].

∗School of Mathematics & Statistics, University of Sydney, NSW 2006, Australia

(johnh@maths.usyd.edu.au).

†[Corresponding author] School of Mathematics & Statistics, University of Sydney, NSW 2006,

Australia (sanjeeva@maths.usyd.edu.au).

‡Dept. of Mathematics, College of Charleston, 66 George Street, Charleston, SC 29424, USA

(lafortunes@cofc.edu)

2 HORNIBROOK, BALASURIYA AND LAFORTUNE

Usually, the heat loss is neglected in these “adiabatic” studies. In several of these articles [15, 22, 27, 28] the condition β ≫ 1 is essential to the wavespeed and stability analysis. The case β ≪ 1 has also been studied [7], in which a perturbative method is used to model the temperature. The bifurcation structure with respect to the heat loss parameter ℓ is addressed in [26], which obtains a stability diagram with respect to ℓ and the wavespeed.

We note that the limit of small fuel diffusivity (large, but not infinite, Lewis number) has not received much attention, perhaps because of the singularity of this limit in the governing equations. Yet this limit may be argued to be particularly appropriate for solid fuels. It is this limit which we study in this article, without restricting β. We do a detailed analysis of the wavespeed of combustion waves which can be supported. We also study the stability of such wavefronts using an Evans function technique.

The model we use is for a premixed fuel in one dimension, with no heat loss and with an Arrhenius law for the reaction rate. These combustion dynamics can be represented in non-dimensional form by [4, 7, 15, 22, 26–28]

\{ \begin{matrix}\frac{∂u}{∂t} & = & \frac{∂^{2}u}{∂x^{2}} + y e^{−1/u} \\ \frac{∂y}{∂t} & = & \frac{1}{Le} \frac{∂^{2}y}{∂x^{2}} − β y e^{−1/u} \end{matrix} . (1.1)

Here, u(x, t) is the temperature, and y(x, t) the fuel concentration, at a point x at time t. The parameters β and Le are as described earlier. We are neglecting heat loss (had we included it, an additional term −ℓ (u − u_{a}) would be necessary on the right-hand side of the u equation in (1.1)). This one-dimensional model is also applicable to combustion in cylinders [24], with u and y being cross-sectionally averaged quantities in this case. See also [5, 6, 14, 21] for closely related governing equations. The non-dimensionalization leading to (1.1) ensures that the cold boundary problem is circumvented (see [27] for a discussion). Since the Lewis number will be assumed large, set ϵ = 1/Le, and therefore 0 ≤ ϵ ≪ 1. This small ϵ limit clearly constitutes a singular perturbation in (1.1).

This article analyzes (1.1) as follows. In Section 2, we determine the wavespeed as a function of β and ϵ. We separate our analysis into two cases: (i) Le = ∞ (Section 2.1), and (ii) 1 ≪ Le < ∞ (Sections 2.2 and 2.3). The first of these is wellstudied, yet we are able to empirically determine a simple formula for the wavespeed

HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS 3

which has excellent accuracy. The case 1 ≪ Le < ∞ needs some asymptotics, based on which in Section 2.2 we obtain a method for computing the wavespeed. We show in Section 2.3 the possibility of obtaining a theoretical estimate for this wavespeed, using a dynamical systems technique called Melnikov’s method. Again, we are able to determine an empirical formula for the wavespeed, which is accurate for all β values (and not restricted to the “usual” large β limit).

The stability analysis of the wavefronts is given in Section 3. Having described the Evans function [1, 11, 18, 30] approach for stability in Section 3.1, we first investigate the Le = ∞ wavefront in Section 3.2. We find that there are no eigenvalues in the right half-plane, supporting the numerical evidence on the stability of this wavefront. In contrast, we show in Section 3.3 that the moment that ϵ is non-zero, eigenvalues pop into the right half-plane. We show that the situation worsens as ϵ increases, generating additional unstable eigenvalues. Thus, wavefronts associated with 1 ≪ Le < ∞ are unstable, and our numerics indicate that this is so for any value of β. A stable propagating front is thus apparently not supportable in (1.1) for large Lewis numbers.

2. Wavespeed analysis. We seek wavefronts which travel in time, and hence set u(x, t) = u(ξ) and y(x, t) = y(ξ), where ξ = x − c t and c is the traveling wave speed. Under this ansatz, (1.1) reduces to

\{ \begin{matrix}u′′ + c u′ + y e^{−1/u} & = & 0 \\ ϵ y′′ + c y′ − β y e^{−1/u} & = & 0 \end{matrix} . (2.1)

2.1. Wavefront for Le = ∞. Set ϵ = 0 in (2.1). Upon defining the new variable v = u′, the dynamics can be represented by a three-dimensional first-order system

\{ \begin{matrix}u′ & = & v \\ v′ & = & −c v − y e^{−1/u} \\ y′ & = & \frac{β}{c} y e^{−1/u} \end{matrix} . (2.2)

The system (2.2) possesses a conserved quantity

H_{c}(u, v, y) = β v + β c u + c y , (2.3)

since it is verifiable that dH_{c}/dξ = 0 along trajectories of (2.2). Thus, motion is confined to planes defined by H_{c}(u, v, y) = constant. Now, for a wavefront, we require that (u, v, y) → (0, 0, 1) as ξ → ∞; this corresponds to the region in which fuel

Fig. 2.1. Projection onto the (u, v)-plane of trajectories of (2.2) lying on different planes H_{c} = c. Here, β = 1, and the three curves correspond to c = 0.5 (dotted), 0.5707 (solid) and 0.7 (dashed).

is not yet burnt (and remains at its maximum non-dimensional mass of one) and the temperature (and its variation) is still zero. This point lies on H_{c}(u, v, y) = c, giving a well-known conservation relation [28]. At the other limit ξ → −∞, the fuel is completely burnt, and has reached a steady temperature, and so (u, v, y) → (u_{B}, 0, 0), where the temperature u_{B} is to be determined. Utilizing H_{c}(u_{B}, 0, 0) = c, we find that u_{B} = 1/β is necessary; see also [7,15,28] for alternative ways to obtain this value.

Thus, we seek a heteroclinic solution of (2.2), which progresses between the fixed points (1/β, 0, 0) and (0, 0, 1), and is confined to the plane β v + β c u + c y = c; i.e., the fuel concentration obeys

y = 1 − β u − \frac{β}{c} v (2.4)

at all values of ξ. Considering (2.2) under this restriction, we obtain

\{ \begin{matrix}u′ & = & v \\ v′ & = & − c v − (1 − β u − \frac{β}{c} v) e^{−1/u} \end{matrix} . (2.5)

Fig. 2.2. Variation of the wavespeed c with β: open circles (numerical results); unbroken curve (empirical curve (2.6)); dotted curve (exp(−0.5β), as obtained in [6, 22, 27]).

This is effectively a projection of the flow on the particular invariant plane H_{c}(u, v, y) = c onto the (u, v)-plane. Any value of c for which a heteroclinic connection exists between (1/β, 0) and (0, 0) is a permitted speed for the wavefront.

The unstable eigen-direction of the point (1/β, 0) is (−c, −βe^{−β}), and we determine c numerically by shooting along this direction, and attempting to match up with a trajectory approaching the origin. In Figure 2.1 we show several numerically computed trajectories of this form, for different values of c, where we have chosen β = 1. Note that this is not a standard (u, v)-phase space for (2.5), since each curve corresponds to a different value of the parameter c. Rather, it is a projection onto the (u, v)-plane of specialized trajectories from the invariant planes H_{c}(u, v, y) = c of the three-dimensional system (2.2). The one trajectory which makes the required connection lies in the invariant plane corresponding to c = 0.5707. The determination of this c value was obtained by making incremental adjustments of c until an appropriate connection is obtained.

We use this technique to numerically compute the wavespeeds for various values

Fig. 2.3. Temperature front at β = 1 (solid, left scale) and β = 3 (dashed, right scale)

of the fuel parameter β, and illustrate this dependence in Figure 2.2. The wavespeed decays with β. For fuels with larger β (poor fuels), the energy resulting from the reaction is insufficient to quickly activate combustion in nearby material, and combustion fronts propagate slowly. The data fits the exponential curve

c(β) = 0.926 e^{−0.486β} (2.6)

with correlation ρ > 0.9999. Equation (2.6) therefore provides an empirically determined formula of excellent accuracy, for the speed of a wavefront in perfectly solid adiabatic one-dimensional media. This result is close (and consistent with) a variety of sources: exp(−0.5β) is quoted in [27] for the small β limit; this same value is given as an upper bound in [6], and also implied in eq. (10) in [22] using a large β limit within a discontinuous front approximation. See Figure 2.2 for a comparison with our results.

The structure of the temperature front is illustrated in Figure 2.3 for β = 1 (solid curve, left scale) and β = 3 (dashed curve, right scale), demonstrating that smaller β fronts are narrower.

HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS 7

2.2. Wavefront for 1 ≪ Le < ∞. When the Lewis number is not infinite, but large, ϵ is small, and weak fuel diffusion needs to be permitted. This is a singular limit in (1.1) and (2.1), and as a consequence has been hardly examined in the literature. By defining v = u′ as before, but now also z = y′, the governing equations (2.1) can be represented as a four-dimensional system

\{ \begin{matrix}u′ & = & v \\ v′ & = & −c v − y e^{−1/u} \\ y′ & = & z \\ z′ & = & \frac{1}{ϵ} (−c z + β y e^{−1/u}) \end{matrix} . (2.7)

The quantity

G_{c}^{ϵ}(u, v, y, z) = β v + β c u + c y + ϵ z

can be verified to be a conserved quantity of (2.7). Hence, flow is confined to the invariant three-dimensional surfaces G_{c}^{ϵ} = constant. Now, we seek a wavefront solution which goes from (u, v, y, z) = (u_{B}, 0, 0, 0) to a value (0, 0, 1, 0), and we find that G_{c}^{ϵ}(u, v, y, z) = c, and u_{B} = 1/β as before. The three-dimensional invariant surface on which both points lie is

z = \frac{1}{ϵ} (c − β v − β c u − c y) .

The dynamics of (2.7) on this surface can be projected to the three variables (u, v, y), such that

\{ \begin{matrix}u′ & = & v \\ v′ & = & −c v − y e^{−1/u} \\ y′ & = & \frac{1}{ϵ} (c − β v − β c u − c y) \end{matrix} . (2.8)

We seek a heteroclinic connection from (0, 0, 1) to (1/β, 0, 0) in the small ϵ limit. However, there are two “time”-scales in this singularly perturbed system, and hence we adopt the standard trick of defining a new independent variable η = ξ/ϵ to elucidate motion in the fast “time” η. With a dot denoting the rate of change with respect to

Fig. 2.4. The hyperbolic invariant manifold (a) S_{0} for (2.10), and (b) S_{ϵ} for (2.9)

η, (2.8) becomes

\{ \begin{matrix}u̇ & = & ϵ v \\ v̇ & = & ϵ (−c v − y e^{−1/u}) \\ ẏ & = & c − β v − β c u − c y \end{matrix} . (2.9)

In the ϵ = 0 limit, the system collapses to

\{ \begin{matrix}u̇ & = & 0 \\ v̇ & = & 0 \\ ẏ & = & c − β v − β c u − c y \end{matrix} , (2.10)

in which it is clear that the plane S_{0} defined by c − β v − β c u − c y = 0 consists entirely of fixed points. This is the same plane as defined through H_{c}(u, v, y) = c for equation (2.2), on which the interesting behavior occurred for perfectly solid fuels. Each fixed point has a one-dimensional stable manifold (in the y-direction), and a two-dimensional centre manifold, which is S_{0}. Thus the plane S_{0} is invariant and normally hyperbolic with respect to (2.10); there is exponential contraction towards it as illustrated in Figure 2.4(a).

Upon switching on ϵ and considering the dynamics (2.9), S_{0} perturbs to an invariant curved entity S_{ϵ}, which is order ϵ away from S_{0}. This is because of the structural stability of normally hyperbolic sets [12], which also implies that normal hyperbolicity is preserved for small ϵ. Therefore, there is exponential decay of trajectories towards S_{ϵ} on “time”-scales of order η, as shown in Figure 2.4(b). Motion on S_{ϵ} occurs at a

HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS 9

slower rate (on “time”-scales of order ξ), and hence it is termed a ‘slow manifold’. The heteroclinic connection we seek lies on S_{ϵ}, from (u, v, y) = (1/β, 0, 0) to (0, 0, 1). Since S_{ϵ} is invariant, two-dimensional and not parallel to the y-axis, it therefore makes sense to project the motion to the (u, v)-plane in order to describe behavior. To elucidate this motion, we need to once again return to the original “time”-scale ξ – the slow time associated with motion on the slow manifold.

Return to the relationship G_{c}^{ϵ} (u(ξ), v(ξ), y(ξ), z(ξ)) = c, which upon differentiation yields

β v′ + β c u′ + c y′ + ϵ z′ = 0 ,

and since u′ = v and y′ = z,

z = −\frac{β}{c} v′ − β v − \frac{ϵ}{c}z′ .

Substituting back into G_{c}^{ϵ}(u, v, y, z) = c, we obtain

β v + β c u + c y + ϵ (−\frac{β}{c} v′ − β v + O(ϵ)) = c ,

and thus

y = 1 − \frac{β}{c} v − β u + ϵ \frac{β}{c^{2}} v′ + ϵ \frac{β}{c} v + O(ϵ^{2}) .

Substitution into the v′ equation in (2.7) or (2.8) gives

v′ (1 + ϵ \frac{β}{c^{2}} e^{−1/u}) = −c v − (1 − \frac{β}{c} v − β u + ϵ \frac{β}{c} v + O(ϵ^{2})) e^{−1/u} .

Therefore

v′ = (1 − ϵ \frac{β}{c^{2}} e^{−1/u}) [−c v − (1 − \frac{β}{c} v − β u + ϵ \frac{β}{c} v) e^{−1/u}] + O(ϵ^{2})
= −c v − (1 − \frac{β}{c} v − β u) e^{−1/u} + ϵ \frac{β}{c^{2}} (1 − \frac{β}{c} v − β u) e^{−2/u} + O(ϵ^{2}) .

Retaining only O(ϵ) terms, we obtain the (u, v)-projected approximate equations on the slow manifold

\{ \begin{matrix}u′ & = & v \\ v′ & = & −c v − (1 − \frac{β}{c} v − β u) e^{−1/u} + ϵ \frac{β}{c^{2}} (1 − \frac{β}{c} v − β u) e^{−2/u} \end{matrix} . (2.11)

Choosing different c values would change the slow manifold S_{ϵ}. So when we seek heteroclinic solutions to (2.11) from (1/β, 0) to (0, 0) by varying c, we are effectively adopting the same attitude as in the Section 2.1. The dynamics are on a

Fig. 2.5. Wavespeed variation with ϵ for different values of β: β = 0.85 (triangles), β = 1 (circles) and β = 1.16 (stars). The solid line is the theoretical approximation for β = 1 obtained from the methods of Section 2.3.

two-dimensional surface, which varies with c, from which we project the trajectories to the (u, v)-plane. This surface in this case is not the plane y = 1 − βu − βv/c (or H_{c}(u, v, y) = c), but is a nonlinear surface O(ϵ) close to this. However, for each chosen small value of ϵ, we can use exactly the same strategy as in Section 2.1 to determine the wave speeds c at which an appropriate heteroclinic connection exists.

The wavespeed speed was calculated for various ϵ and β values (see Figure 2.5). We notice that c decreases as we increase ϵ, that is, when we decrease the Lewis number. Now, in dimensional form Le = κ/(ρ c_{p} D), where ρ, κ, c_{p} and D are the respectively the density, thermal conductivity, specific heat capacity and molecular diffusivity of the fuel [5,7,24,26,28]. Increasing ϵ is equivalent to increasing the relative importance of D, ρ and c_{p} in relation to κ. Reducing κ obviously decreases the ability of heat to move, and hence the combustion speed. Higher densities result in increased fuel mass in each location, which means more heat is needed in a given area to ignite all of the fuel before the wave moves on. Fuels with increased c_{p} require more heat to increase the temperature by the a specified amount. Finally, increasing D increases

HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS 11

the transport of burnt fuel into the unburnt region and vice-versa, interfering with front propagation.

We computed the changes to the wavefront profile (akin to Figure 2.3) when ϵ is changed (not shown). We verified the obvious physical conclusion that the fuel concentration front becomes less steep when ϵ is increased from zero.

2.3. Perturbative formula for wavespeed. Here, we derive and numerically study a formula for the wavespeed correction in going from Le = ∞ to finite Lewis number. Let

c (β, ϵ) = c_{0} (β) + ϵ c_{1} (β) + O (ϵ^{2}) , (2.12)

where c_{0} is the wavespeed associated with the infinite Lewis number (ϵ = 0) combustion wavefront. In the spirit of perturbation analysis, we obtain a formula for the correction c_{1}(β) purely in terms of the unperturbed (ϵ = 0) wave, using a nontraditional application of “Melnikov’s method” [23] from dynamical systems theory.

Melnikov’s method is applied most commonly to area-preserving two-dimensional systems under time-periodic perturbations [3, 16, 29]. (Here, “time” is used loosely to mean the independent variable, in our case ξ.) Our system (2.11) turns out to be not area-preserving, and has a perturbation which is independent of the temporal variable. Under these conditions, we describe the method applied to the system

z′ = f (z) + ϵ g (z) . (2.13)

When ϵ = 0, suppose this system possesses a heteroclinic connection between the two saddle fixed points a and b as shown in Figure 2.6(a). A heteroclinic connection of this sort occurs when a branch of the one-dimensional unstable manifold of a coincides with a branch of the stable manifold of b. This heteroclinic trajectory can be represented as a solution z = ẑ(ξ) to (2.13) with ϵ = 0.

Now, for small ϵ > 0 in (2.13), the fixed points a and b perturb by O(ϵ), and retain their stable and unstable manifolds [12]. However, these need no longer coincide. Figure 2.6(b) shows how they can split apart, with the dashed curve showing the original manifold. Let d(ξ, ϵ) be a distance measure between these manifolds, measured along a perpendicular to the unperturbed heteroclinic drawn at ẑ(−ξ). The variable ξ can thus be used to identify the position along the heteroclinic curve (cf. “heteroclinic coordinates” of Section 4.5 in [29]). Since d(ξ, 0) = 0 for all ξ, this distance is Taylor

Fig. 2.6. Manifold structure for the Melnikov approach: (a) ϵ = 0, (b) ϵ ≠ 0

expandable in ϵ in the form

d(ξ, ϵ) = ϵ \frac{M(ξ)}{|f (ẑ(−ξ))|} + O(ϵ^{2}) ,

where the scaling factor |f (ẑ(−ξ))| in the denominator represents the unperturbed trajectory’s speed at the location ξ. The quantity M(ξ) is the “Melnikov function”, for which an expression turns out to be

(2.14) M (ξ) = ∫_{−∞}^{∞} exp [− ∫_{−ξ}^{μ} ∇ · f(ẑ(s)) ds] f(ẑ(μ)) ∧ g(ẑ(μ)) dμ ,

where the wedge product between two vectors is defined by (a_{1}, a_{2})^{T} ∧ (b_{1}, b_{2})^{T} = a_{1}b_{2} − a_{2}b_{1}. Obtaining the version (2.14) requires two adjustments to the standard Melnikov approaches [3, 16, 29]: incorporating the non area-preserving nature of the unperturbed flow of (2.13), and representing the distance in terms of heteroclinic coordinates. We need to ensure the persistence of a heteroclinic trajectory in (2.11) for ϵ > 0, and thus require d(ξ, ϵ) = 0 for all ξ. For this to happen for all small ϵ, we therefore need to set M(ξ) ≡ 0.

To apply this technique to our system, we begin by writing (2.11) in the form (2.13). Using the expansion (2.12), and utilizing binomial expansions for 1/(c_{0} + ϵc_{1}), we get

\{ \begin{matrix}u′ & = & v \\ v′ & = & −c_{0} v − e^{−1/u}Υ_{uv} + ϵ (−c_{1}v − \frac{βc_{1}e^{−1/u}}{c_{0}^{2}} v + \frac{βe^{−2/u}}{c_{0}^{2}} Υ_{uv}) \end{matrix} , (2.15)

where higher-order terms in ϵ have been discarded, and

Υ_{uv} = 1 − β u − \frac{β}{c_{0}} v .

HIGH LEWIS NUMBER COMBUSTION WAVEFRONTS 13

By appropriately identifying f and g from (2.15) through comparison with (2.13), we see that

(f ∧ g) (u, v) = v (−c_{1}v − \frac{βc_{1}e^{−1/u}v}{c_{0}^{2}} + \frac{βe^{−2/u}}{c_{0}^{2}} Υ_{uv})

and ∇ · f = −c_{0} + βe^{−1/u}/c_{0}. Substituting into the Melnikov formula (2.14), and setting it equal to zero, we obtain

∫_{−∞}^{∞} exp [∫_{−ξ}^{μ} (c_{0} − \frac{β}{c_{0}} e^{−1/u(s)}) ds] v (−c_{1}v − \frac{βvc_{1}e^{−1/u}}{c_{0}^{2}} + \frac{βe^{−2/u}}{c_{0}^{2}} Υ_{uv}) dμ = 0 ,

where each of u(μ) and v(μ) is evaluated along the ϵ = 0 combustion wave. Notice, however, that for this infinite Lewis number combustion wave, (2.4) tells us that the fuel concentration y(μ) = Υ_{uv}(μ) for all μ. Therefore

c_{1}(β) = \frac{β ∫_{−∞}^{∞} exp [∫_{−ξ}^{μ} (c_{0} − \frac{β}{c_{0}} e^{−1/u(s)}) ds] v y e^{−2/u} dμ}{∫_{−∞}^{∞} exp [∫_{−ξ}^{μ} (c_{0} − \frac{β}{c_{0}} e^{−1/u(s)}) ds] v^{2}(c_{0}^{2} + β e^{−1/u}) dμ} , (2.16)

where u(μ), v(μ) and y(μ) in the integrands are obtained from the ϵ = 0 combustion wave discussed in Section 2.1. The apparent dependence of c_{1} on the wave coordinate ξ is spurious: if I is an anti-derivative of the inner integrals in (2.16), a multiplicative term exp [−I(−ξ)] emerges in both the numerator and denominator, which therefore cancels. Hence, any convenient value for ξ can be chosen in (2.16), for example 0.

We note that v < 0 for the infinite Lewis number wavefront, as is clear from the phase-portrait Figure 2.1. Alternatively, u is smaller at the front of the wave, where fuel is yet to be burnt, and is therefore a decreasing function of μ, leading to v = u′ < 0. Based on this, (2.16) immediately displays that c_{1} < 0, proving the property that the wavespeed decreases when fuel diffusivity is included. This is in agreement with the numerical observations in Section 2.2.

Equation (2.16) provides an explicit perturbative formula on how the wavespeed varies through the inclusion of the finiteness of the Lewis number, expressed entirely in terms of the infinite Lewis number combustion wave. This result is used to compute the solid line in Figure 2.5, which is the theoretical wavespeed 0.5707 − 0.1552ϵ obtained by using (2.16) and (2.12) when β = 1. When ϵ is small, it forms an excellent approximation to the numerically obtained wavespeed, as described in Section 2.2. Indeed, Figure 2.5 show that the theoretical line is tangential to the curve formed by the closed circles.

Fig. 2.7. The perturbing wavespeed as a function of β

The perturbation wavespeed c_{1} as a function of β appears in Figure 2.7. There is a β value (around 2) at which the absolute influence of the finiteness of the Lewis number is greatest. Nevertheless, since c_{0} is itself a function of β, it may make sense to investigate the relative influence c_{1}/c_{0} of the perturbative term. Such is presented in the numerically computed Figure 2.8. The graph is virtually linear and has zero intercept. We therefore arrive at the empirical approximation

c(β, ϵ) = c_{0}(β) [1 − 0.267 ϵ β] = c_{0}(β) [1 − 0.267 \frac{β}{Le}] , (2.17)

for large Lewis numbers, with excellent validity across all β, and with c_{0}(β) also known through (2.6). The simplicity of this expression is remarkable.

3. Stability analysis. Here, we analyze the stability of the wavefronts we have discovered in Section 2. To do so, we perform numerical calculation of the Evans function. In cases when the dimensions of the stable and unstable spaces are greater than one, such a calculation is usually performed in an exterior algebra space (see for example [2, 8, 9, 15, 26]). For our purpose, we use the more direct method introduced in [20] in which the Evans function is defined as a determinant.

3.1. Evans function. In general, the linear stability of a localized traveling wave solution to a system of PDEs is obtained by studying the eigenvalue problem

Fig. 2.8. Relative size of perturbative wavespeed as a function of β

Lw = λ w, (3.1)

where the matrix differential operator L arises from the linearization of the PDEs. The traveling solution is said to be linearly stable if the spectrum of L lies in the closed left half-plane.

The system (3.1) can be turned into a linear dynamical system of the form

X′ = A(ξ, λ) X (3.2)

where A is a n × n square matrix depending on ξ = x − c t and the spectral parameter λ. Under some technical conditions (see Appendix A for details), it can be shown that the asymptotic behavior of the solutions to (3.2) is determined by the matrices

A^{±∞}(λ) = lim_{ξ→±∞} A(ξ, λ)

in the following sense. If μ_{+} (resp. μ_{−}) is an eigenvalue of A^{+∞} (resp. A^{−∞}) with eigenvector v_{+} (resp. v_{−}), then there exists a solution w_{+} (resp. w_{−}) to (3.2) with the property that

(3.3)  lim_{ξ→∞} w_{+}e^{−μ_{+} ξ} = v_{+} (resp. lim_{ξ→−∞} w_{−}e^{−μ_{−} ξ} = v_{−}) .

16 HORNIBROOK, BALASURIYA AND LAFORTUNE

To study the linear stability, one should consider both the essential and point spectrum of L. The essential spectrum of L consists of the values of λ for which A^{∞} or A^{−∞} has purely imaginary eigenvalues [17]. The point spectrum can be studied by means of the Evans function. Let Ω denote a domain of the complex λ plane with no intersection with the essential spectrum and let n_{s} and n_{u} denote, respectively, the number of eigenvalues of A^{∞} with negative real part and the number of eigenvalues of A^{−∞} with positive real part in Ω. We assume that n_{s} + n_{u} = n. Let w_{i}^{+}(λ, ξ), i = 1, 2, ..., n_{s} (resp. w_{i}^{−}(λ, ξ), i = 1, 2, ..., n_{u}) be linearly independent solutions to (3.2) converging to zero as ξ → ∞ (resp. ξ → −∞) which are analytic of λ in Ω. Clearly, a particular value of λ belongs to the point spectrum of L if (3.2) admits a solution that is converging to zero for both ξ → ±∞, that is if the space of solutions generated by the w_{i}^{+} intersects with the one generated by the w_{i}^{−}. To detect such values of λ in Ω, we define the Evans function as

E(λ) = det (w_{1}^{+}, w_{2}^{+}, ...w_{n_{s}}^{+}, w_{1}^{−}, w_{2}^{−}, ...w_{n_{u}}^{−}),

in which the w_{i}^{±} are evaluated at ξ = 0. This function is analytic in Ω, is real for real values of λ and the locations of the zeros of E(λ) correspond to eigenvalues of L. There is a great deal of arbitrariness in the definition of the Evans function due to the non-uniqueness of the w_{i}^{±}. However, the locations of the zeros of the Evans function are not affected by a redefinition.

3.2. Le = ∞ wavefront. In the case Le = ∞, we denote the wavefront solution by (u, y) = (u_{0}(ξ), y_{0}(ξ)) with ξ = x − c t and consider a perturbation of the form

u = u_{0}(ξ) + U(ξ) e^{λ t}, y = y_{0}(ξ) + Y (ξ) e^{λ t}. (3.4)

At first order, U and Y satisfy an eigenvalue problem of the form (3.1) which can be turned into the linear dynamical system

( \begin{matrix}U \\ V \\ Y \end{matrix} )′ = ( \begin{matrix}0 & 1 & 0 \\ λ − \frac{y_{0}}{u_{0}^{2}} e^{−1/u_{0}} & −c & −e^{−1/u_{0}} \\ \frac{β y_{0}}{c u_{0}^{2}} e^{−1/u_{0}} & 0 & \frac{λ}{c} + \frac{β}{c} e^{−1/u_{0}} \end{matrix} ) ( \begin{matrix}U \\ V \\ Y \end{matrix} ) . (3.5)

Using the fact that (u_{0}, y_{0}) → (0, 1) as ξ → ∞ and (u_{0}, y_{0}) → (1/β, 0) as ξ → −∞, one finds the limit matrices

A^{∞} = ( \begin{matrix}0 & 1 & 0 \\ λ & −c & 0 \\ 0 & 0 & \frac{λ}{c} \end{matrix} ) , A^{−∞} = ( \begin{matrix}0 & 1 & 0 \\ λ & −c & −e^{−β} \\ 0 & 0 & \frac{λ}{c} + \frac{β}{c} e^{−β} \end{matrix} ) .

It is straightforward to verify that A^{+∞} has a purely imaginary eigenvalue if λ is purely imaginary, and that all eigenvalues of A^{±∞} have a nonzero real part if ℜ(λ) > 0. The essential spectrum thus does not intersect the open right half complex plane and does include the imaginary axis. The domain Ω can therefore be taken in the open right half-plane. Furthermore, it can be verified that, for ℜ(λ) > 0, A^{−∞} has two eigenvalues with positive real parts (n_{u} = 2), μ_{i}^{−}, i = 1, 2, and A^{∞} has one with negative real part (n_{s} = 1), μ_{1}^{+}. These eigenvalues, together with their respective eigenvectors, are given by

μ_{1}^{−} = \frac{λ + β e^{−β}}{c} , v_{1}^{−} = ( \begin{matrix}−e^{−β} \\ −e^{−β}μ_{1}^{−} \\ μ_{1}^{−}^{2} + β e^{−β} \end{matrix} ) ;

μ_{2}^{−} = −\frac{c}{2} + \frac{1}{2}\sqrt{c^{2} + 4 λ}, v_{2}^{−} = ( \begin{matrix}1 \\ μ_{2}^{−} \\ 0 \end{matrix} ) ;

μ_{1}^{+} = −\frac{c}{2} − \frac{1}{2}\sqrt{c^{2} + 4 λ}, v_{1}^{+} = ( \begin{matrix}1 \\ μ_{1}^{+} \\ 0 \end{matrix} ) .

For any given λ, the solution w_{1}^{+} of (3.5) which converges to zero as ξ → ∞ is obtained numerically, by choosing the initial condition

w_{1}^{+}(L) = v_{1}^{+} (3.6)

for L sufficiently large. Similarly, the two solutions w_{1}^{−} and w_{2}^{−} which converge as ξ → −∞ are obtained numerically with the initial conditions

w_{1}^{−}(−L) = v_{1}^{−}, w_{2}^{−}(−L) = v_{2}^{−}. (3.7)

The Evans function is then defined as

E(λ) = det (w_{1}^{+}(0), w_{1}^{−}(0), w_{2}^{−}(0)). (3.8)

Fig. 3.1. Values of the Evans function E(λ) in the case Le = ∞ as λ varies on the closed right half of the circle centered at (0.001, 0) with radius 0.35 (β = 1, ϵ = 0, and c = 0.5707).

Once the solutions w_{1}^{+} and w_{i}^{−}, i = 1, 2 are obtained numerically, the Evans function can be calculated. In order to check the accuracy of our numerical results, we use the procedure described in [20] which takes advantage of the fact that the evolution of the Wronskian associated to a fundamental set of solutions of a linear dynamical system such as (3.2) is known.

Figures 3.1 and 3.2 show how the Evans function behaves as λ varies on the closed right half of the circle centered at (0.001, 0) with radius 0.35 (with β = 1, ϵ = 0, and c = 0.5707). Once the values of the Evans function on such a closed path are known, it is possible to determine the number of zeros inside the contour by numerically evaluating the winding number. In the case displayed in Figure 3.1 (and for the other values of β we have considered), the winding number was found to be zero.

Furthermore, it is possible to perform an analytic continuation of the Evans function across the essential spectrum using the Gap Lemma [13,19]. For values of λ with negative real parts, one uses the same definition of the Evans function (3.8) with initial conditions (3.6) and (3.7). However, one should notice that, when ℜ (λ) < 0, then,

Fig. 3.2. Zoom in on Figure 3.1

since ℜ (μ_{2}^{+}) < 0, the solution w_{2}^{−} is now divergent as ξ → −∞. The Gap lemma guarantees that there exists a region that extends beyond the essential spectrum and on which the Evans function can be analytically continued.

It is thus possible to consider a semi-circle which slightly intersects the left halfside of the complex plane. For instance, we translated the same half circle used before to generate the data of Figure 3.1 to the right and moved the center to (−0.001, 0). It is known [25] that the zero of the Evans at the origin function is at least of order one due to the translational invariance of (1.1). We recover this fact in our numerical computations, by obtaining a winding number of one for this displaced semi-circle, in contrast with zero for the semi-circle in the right half-plane. Thus, the only zero inside the closed curve is indeed at the origin.

It is also possible to graph the Evans function on the real λ line. Figure 3.3 shows such a plot, in which a scaled version of the Evans function is displayed to compensate for its rapid growth with λ. This figure backs up our winding number calculations which show an order one zero at the origin, and no zeros in the right half-plane.

Fig. 3.3. The Evans function (scaled) along the positive real axis when ϵ = 0 and β = 1 (c = 0.5707 here).

3.3. Le < ∞ wavefront. In the case 0 < Le = 1/ϵ, we consider the system (1.1) and its wavefront denoted (u, y) = (u_{0}(ξ), y_{0}(ξ)). With a perturbation of the form (3.4), we find an eigenvalue problem giving rise to the following 4-dimensional linear dynamical system

( \begin{matrix}U \\ V \\ Y \\ Z \end{matrix} )′ = ( \begin{matrix}0 & 1 & 0 & 0 \\ λ − \frac{y_{0}}{u_{0}^{2}} e^{−1/u_{0}} & −c & −e^{−1/u_{0}} & 0 \\ 0 & 0 & 0 & 1 \\ \frac{β y_{0}}{ϵ u_{0}^{2}} e^{−1/u_{0}} & 0 & \frac{λ}{ϵ} + \frac{β}{ϵ} e^{−1/u_{0}} & −\frac{c}{ϵ} \end{matrix} ) ( \begin{matrix}U \\ V \\ Y \\ Z \end{matrix} ) .

Again, using the fact that (u_{0}, y_{0}) → (0, 1) as ξ → ∞ and (u_{0}, y_{0}) → (1/β, 0) as ξ → −∞, one easily computes the limit matrices and finds that n_{u} = n_{s} = 2. The corresponding eigenvectors and eigenvalues are given by

μ_{1}^{−} = \frac{−c + \sqrt{c^{2} + 4 ϵ (λ + β e^{−β})}}{2 ϵ},

with v_{1}^{−} = ( \begin{matrix}−e^{−β} \\ −e^{−β}μ_{1}^{−} \\ (λ − c μ_{1}^{−}) (1 − \frac{1}{ϵ}) − \frac{β e^{−β}}{ϵ} \\ −\frac{c}{ϵ} (λ + β e^{−β}) (1 − \frac{1}{ϵ}) + μ_{1}^{−} ((λ + \frac{c^{2}}{ϵ} ) (1 − \frac{1}{ϵ}) − \frac{β e^{−β}}{ϵ} ) \end{matrix} ) ;

μ_{2}^{−} = −\frac{c}{2} + \frac{1}{2}\sqrt{c^{2} + 4 λ}, with v_{2}^{−} = ( \begin{matrix}1 \\ μ_{2}^{−} \\ 0 \\ 0 \end{matrix} ) ;

μ_{1}^{+} = −\frac{c}{2} − \frac{1}{2}\sqrt{c^{2} + 4 λ}, with v_{1}^{+} = ( \begin{matrix}1 \\ μ_{1}^{+} \\ 0 \\ 0 \end{matrix} ) ;

μ_{2}^{+} = −\frac{c + \sqrt{c^{2} + 4 λ ϵ}}{2 ϵ}, with v_{2}^{+} = ( \begin{matrix}0 \\ 0 \\ 1 \\ μ_{2}^{+} \end{matrix} ) .

The essential spectrum can again be shown not to have an intersection with the open right half-plane.

The Evans function can again be calculated numerically, and whether it possesses any zeros with positive real part can be investigated using the ideas of Nyquist diagrams and winding numbers. We will omit these diagrams, and present the results. We find that the moment ϵ becomes nonzero, a zero of the Evans function jumps off the origin onto the positive real axis. This appears to happen for any β. In fact, as ϵ is further increased, this zero progresses along the real axis, and at a certain value of ϵ, another real positive zero appears. When ϵ is increased further, a third zero bifurcates out onto the real positive axis, and this process apparently continues. We show the behavior of these zeros in Figure 3.4, computed when β = 0.4. While the detail of this is interesting, we remark that the main conclusion from this is that the wavefront for ϵ > 0 is unstable due to the presence of at least one such zero. The instability mechanism is through a linear instability mode which grows monotonically (since the zeros are on the real axis).

We illustrate the (scaled) Evans function along the real positive axis when ϵ = 0.2

Fig. 3.4. The bifurcating structure of the real positive zeros of the Evans function (here, β = 0.4)

in Figure 3.5, with β = 1 chosen for comparison with Figure 3.3. This is an instance in which there are two real positive zeros of the Evans function, associated with two different unstable modes. The wavefront will not be able to progress while maintaining its structure. Based on our computations, we conclude that the wavefront is unstable for large Lewis numbers, even though the singular limit of infinite Lewis number produces stable wavefronts. The effect of fuel diffusion is apparently destabilizing. Irrespective of how small it is, as long as it is non-zero, our results indicate that the wavefront is unstable.

Appendix A. Analyticity of the Evans function.

In this appendix, we discuss the technical details about the analyticity of the solutions of (3.2).

Consider a linear n−dimensional system of the form (3.2). Assume that A is integrable in ξ, that A^{∞}(λ) is diagonalizable, and that

∫_{1}^{∞} |A − A^{∞}| dξ (A.1)

converges. Then, if μ_{+} is an eigenvalue of A^{+∞} with eigenvector v_{+}, there exists a

Fig. 3.5. The Evans function (scaled) along the positive real axis when ϵ = 0.2 and β = 1 (c = 0.5347 here).

solution w_{+} to (3.2) satisfying (3.3) (see [10] page 104).

Furthermore, if we consider a domain Ω of the complex plane in which (A.1) converges uniformly in λ in any compact subset of Ω, the eigenvalues and eigenvectors of A^{∞} are analytic functions of λ, and the number of eigenvalues with distinct real parts does not change, then w_{+} can be chosen to be an analytic function of λ [10,25].

REFERENCES

[1] J. Alexander, R. Gardner, and C. Jones, A topological invariant arising in the stability analysis of travelling waves, J. Reine Angew. Math., 410 (1990), pp. 167–212.

[2] L. Allen and T. J. Bridges, Numerical exterior algebra and the compound matrix method, Numer. Math., 92 (2002), pp. 197–232.

[3] D. K. Arrowsmith and C. M. Place, An Introduction to Dynamical Systems, Cambridge University Press, Cambridge, 1990.

[4] A. Bayliss and B. Matkowsky, Two routes to chaos in condensed phase combustion, SIAM J. Appl. Math., 50 (1990), pp. 437–459.

[5] , From traveling waves to chaos in combustion, SIAM J. Appl. Math., 54 (1994), pp. 147– 174.

[6] J. Billingham, Phase plane analysis of one-dimensional reaction diffusion waves with degenerate reaction terms, Dyn. Stab. Sys., 15 (2000), pp. 23–33.

[7] J. Billingham and G. N. Mercer, The effect of heat loss on the propogation of strongly exothermic combustion waves, Combust. Theory Modelling, 5 (2001), pp. 319–342.

[8] T. J. Bridges, G. Derks, and G. Gottwald, Stability and instability of solitary waves of the fifth-order KdV equation: a numerical framework, Phys. D, 172 (2002), pp. 190–216.

[9] L. Q. Brin, Numerical testing of the stability of viscous shock waves, Math. Comp., 70 (2001), pp. 1071–1088.

[10] E. A. Coddington and N. Levinson, Theory of ordinary differential equations, McGraw-Hill, New York, 1955.

[11] J. W. Evans, Nerve axon equations. IV. The stable and the unstable impulse, Indiana Univ. Math. J., 24 (1974/75), pp. 1169–1190.

[12] N. Fenichel, Persistence and smoothness of invariant manifolds for flows, Indiana Univ. Math. J., 21 (1971), pp. 193–226.

[13] R. Gardner and K. Zumbrun, The Gap Lemma and geometric criteria for instability of viscous shock profiles, Commun. Pure Appl. Math., LI (1998), pp. 797–855.

[14] B. Gray, S. Kalliadasis, A. Lazarovich, C. Macaskill, J. Merkin, and S. Scott, The suppression of exothermic branched-chain flame through endothermic reaction and radical scavenging, Proc. R. Soc. Lond. A, 458 (2002), pp. 2119–2138.

[15] V. Gubernov, G. N. Mercer, H. S. Sidhu, and R. O. Weber, Evans function stability of combustion waves, SIAM J. Appl. Math., 63 (2003), pp. 1259–1275.

[16] J. Guckenheimer and P. Holmes, Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields, Springer, New York, 1983.

[17] D. Henry, Geometric theory of semilinear parabolic equations, vol. 840 of Lecture Notes in Mathematics, Springer-Verlag, Berlin, 1981.

[18] C. K. R. T. Jones, Stability of the travelling wave solution of the FitzHugh-Nagumo system, Trans. Amer. Math. Soc., 286 (1984), pp. 431–469.

[19] T. Kapitula and B. Sandstede, Stability of bright solitary-wave solutions to perturbed nonlinear Schrödinger equations, Phys. D, 124 (1998), pp. 58–103.

[20] S. Lafortune, J. Lega, and S. Madrid-Jaramillo, Instability of local deformations of an elastic rod: numerical evaluation of the Evans function, Preprint, (2005).

[21] S. Margolis and F. Williams, Diffusion/thermal instability of solid propellant flame, SIAM J. Appl. Math., 49 (1989), pp. 1390–1420.

[22] B. Matkowsky and G. Sivashinsky, Propagation of a pulsating reaction front in solid fuel combustion, SIAM J. Appl. Math., 35 (1978), pp. 465–478.

[23] V. K. Melnikov, On the stability of the centre for time-periodic perturbations, Trans. Moscow Math. Soc., 12 (1963), pp. 1–56.

[24] G. N. Mercer and R. O. Weber, Combustion waves in two dimensions and their onedimensional approximation, Combust. Theory Modelling, 1 (1997), pp. 157–165.

[25] R. L. Pego and M. I. Weinstein, Eigenvalues and instabilities of solitary waves, Phil. Trans.: Phys. Sc. and Eng., 340 (1992), pp. 47–94.

[26] P. Simon, S. Kalliadasis, J. Merkin, and S. Scott, Evans function analysis of the stability

of non-adiabatic flames, Combust. Theory Modelling, 7 (2003), pp. 545–561.

[27] F. Varas and J. Vega, Linear stability of a plane front in solid combustion at large heat of reaction, SIAM J. Appl. Math., 62 (2002), pp. 1810–1822.

[28] R. O. Weber, G. N. Mercer, H. S. Sidhu, and B. F. Gray, Combustion waves for gases (Le = 1) and solids (Le → ∞), Proc. R. Soc. Lond. A, 453 (1997), pp. 1105–1118.

[29] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, New York, 1990.

[30] E. Yanagida, Stability of fast travelling pulse solutions of the FitzHugh-Nagumo equations, J. Math. Biol., 22 (1985), pp. 81–104.