Diffusive-thermal instabilities of a planar premixed flame aligned with a shear flow

Abstract

The stability of a thick planar premixed flame, propagating steadily in a direction transverse to that of unidirectional shear flow, is studied. A linear stability analysis is carried out in the asymptotic limit of infinitely large activation energy, yielding a dispersion relation. The relation characterises the coupling between Taylor dispersion (or shear-enhanced diffusion) and the flame thermo-diffusive instabilities, in terms of two main parameters, namely, the reactant Lewis number $Le$ and the flow Peclet number $Pe$. The implications of the dispersion relation are discussed and various flame instabilities are identified and classified in the $Le$-$Pe$ plane. An important original finding is the demonstration that for values of the Peclet number exceeding a critical value, the classical cellular instability, commonly found for $Le<1$, exists now for $Le>1$ but is absent when $Le<1$. In fact, the cellular instability identified for $Le>1$ is shown to occur either through a finite-wavelength stationary bifurcation (also known as type-I${}_s$) or through a longwave stationary bifurcation (also known as type-II${}_s$). The latter type-II${}_s$ bifurcation leads in the weakly nonlinear regime to a Kuramoto-Sivashinsky equation, which is determined. As for the oscillatory instability, usually encountered in the absence of Taylor dispersion in $Le>1$ mixtures, it is found to be absent if the Peclet number is large enough. The stability findings, which follow from the dispersion relation derived analytically, are complemented and examined numerically for a finite value of the Zeldovich number. The numerical study involves both computations of the eigenvalues of a linear stability boundary-value problem and numerical simulations of the time-dependent governing partial differential equations. The computations are found to be in good qualitative agreement with the analytical predictions.

Keywords:

• Taylor dispersion
• diffusive-thermal instability
• shear flow
• transverse propagation
• anisotropic diffusion

1. Introduction

Thick premixed flames propagating in shear flows, such as flames propagating in a narrow-channel Poiseuille flow, are subject to enhanced diffusion in the longitudinal flow direction due to Taylor's dispersion mechanism, as reported in recent investigations [1–4]. Meanwhile, diffusion transport in a direction transverse to the flow direction is only associated with molecular diffusion. Therefore, diffusion appears as effectively anisotropic which in turn significantly affect flame propagation and stability. The thermo-diffusive stability of premixed flames in the presence of Taylor dispersion has been recently studied, analytically and numerically, in the case of flames propagating in the longitudinal (or flow) directions. This has been done assuming adiabatic conditions in [3] and accounting for heat losses in [5]. The equivalent stability problem for flames aligned with the shear flow, propagating in a direction transverse to the flow has not been addressed yet. This problem is important to investigate on account of the anisotropy of diffusion aforementioned, which is expected to markedly alter the previous findings. Another motivation to consider the problem of transverse premixed-flame propagation is due to the interesting findings of a recent investigation on the effect of Taylor dispersion on nonpremixed flames aligned with the direction of a shear flow [6]. This investigation demonstrates that Taylor dispersion can lead unexpectedly to cellular instability of a planar diffusion flame in mixtures with Lewis numbers above unity. This result has been argued in [6] to present a plausible mechanism to explain the formation of so-called diffusion flame streets observed experimentally in non-premixed microcombustors [7–9]. The question arises naturally therefore whether a similar cellular instability can occur for large Lewis numbers in the premixed case; namely, for flames propagating in the transverse direction, noting that such cellular instability for large Lewis numbers does not occur for flames propagating in the longitudinal direction [3, 5]. One of the main objectives of this paper is to provide a clear answer to this question; as well as to describe the instabilities encountered and identify the conditions for their occurrence in terms of the parameters, mainly the Lewis number characterising the mixture and the Peclet number characterising the flow. In order to focus the analysis on these objectives, mainly concerned with the coupling between Taylor dispersion and the thermo-diffusive flame instabilities [[10],[pp. 477–479][11]], we shall adopt the thermo-diffusive approximation of constant density and constant transport properties and neglect further the effect of heat losses.

The paper is structured as follows. The problem formulation is given in §2 within a Hele-Shaw or two-dimensional channel configuration[12–14], involving a unidirectional shear flow. In this configuration, the depth-averaged governing equations are written down and account for anisotropic diffusion, with enhanced diffusion in the longitudinal flow direction. A reformulation of the problem suitable for the limit of infinite Zeldovich number, $\text{β}\to \mathrm{\infty }$, is provided in §3. This reformulation is used in the linear stability analysis of §4, which is carried out analytically and culminates in the derivation of a dispersion relation. The implications of the dispersion relation are discussed in detail in §5, where various flame instabilities are identified and classified in the parameters space. The stability findings, which follow from the dispersion relation derived analytically in the limit $\text{β}\to \mathrm{\infty }$, are then complemented and examined numerically in §6 for a finite value of the Zeldovich number, $\text{β}=10$. The numerical study in §6 involves both the computation of the eigenvalues of a linear stability boundary-value problem and the numerical simulations of the time-dependent governing partial differential equations for illustrative cases. Concluding remarks are provided in §7.

2. Problem formulation

Consider a narrow channel as depicted in Figure 1 in which a planar premixed flame is propagating in the negative $\hat{y}$-direction with laminar flame speed $S_L$, in the presence of a shear flow. The shear flow is assumed, for simplicity, to have zero mean such as in the case of a Couette flow. Alternatively, in the case of Poiseuille flow for example, the frame of reference is chosen to be moving with the mean flow speed in the $\hat{x}$-direction. The frame is furthermore assumed to move with speed $S_L$ in the negative $\hat{y}$-direction, so that the flame may be considered steady in the unperturbed state. In other words, the flow field adopted is of the form (1) ${\hat{v}}_x=U\hat{u}\left(\hat{z}\right),{\hat{v}}_y=S_L,{\hat{v}}_z=0,$(1) where U is the flow (maximum) amplitude and $\hat{u}\left(\hat{z}\right)$ is the scaled zero-mean shear-flow profile, such that e.g. $\hat{u}=\hat{z}/H$ for a Couette flow and $\hat{u}=1/3-{\hat{z}}^2/H^2$ for a plane Poiseuille flow.

Figure 1. Planar premixed flame in a channel of width $2H$, aligned with the direction ($\hat{x}$) of a shear flow and propagating in the negative transverse direction ($\hat{y}$).

To model the chemistry, a single-step irreversible Arrhenius reaction is adopted, with pre-exponential factor B, activation energy E and heat release q per unit mass of fuel consumed. For sufficiently fuel-lean conditions, the fuel burning rate per unit volume can be written as $\rho B{Y}_F{e}^{-E/RT}$, where ρ is the density assumed to be constant, R is the universal gas constant, $Y_F$ is the fuel mass fraction and T is the temperature. The adiabatic flame temperature $T_{ad}$, the Zeldovich number β and the heat release parameter α are defined by $T_{ad}=T_u+\frac{q{Y}_{F,u}}{c_p},\alpha =\frac{T_{ad}-T_u}{T_{ad}},\text{β}=\frac{E(T_{ad}-T_u)}{R{T}_{ad}^2}$where $T_u$ and $Y_{F,u}$ are the temperature and fuel mass fraction in the unburnt gas and $c_p$ is the specific heat at constant pressure.

For non-dimensionalization, we shall use the laminar flame thickness ${\delta }_L=D_T/S_L^0$ as unit length and ${\delta }_L/S_L^0$ as unit time, where $D_T$ is the thermal diffusivity (assumed constant). Here $S_L^0$ refers to the laminar flame speed (for $\text{β}\to \mathrm{\infty }$) which is given by $S_L^0=[2{\text{β}}^{-2}LeB{D}_T\exp (-E/R{T}_{ad}){]}^{1/2}$, where $Le$ is the (fuel) Lewis number. The dependent variables $Y_F$ and T are normalised by introducing (2) $y_F=\frac{Y_F}{Y_{F,u}},\theta =\frac{T-T_u}{T_{ad}-T_u}.$(2) We consider the limit $H\ll {\delta }_L$ which leads upon averaging in the $\hat{z}$-direction to a problem which is effectively two-dimensional in the first approximation, as shown in [1–4]. The non-dimensional depth-averaged equations read (3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) (4) $\begin{array}{rl}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+S\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}& =(1+p^2)\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{y}^2}+\omega \end{array}$(4) where $p=\sqrt{\gamma }Pe$ is a parameter proportional to the Peclet number $Pe$, $S=\frac{S_L}{S_L^0},\text{{Pe}}=\frac{UH}{D_T},\omega =\frac{{\text{β}}^2}{2\text{{Le}}}{y}_F\exp \left[\frac{\text{β}(\theta -1)}{1+\alpha (\theta -1)}\right].$The constant γ is a numerical coefficient, which is determined by the shear-flow profile and is given by $\gamma ={\int }_0^1\mathrm{d}{z}^{\prime }{\left[{\int }_0^{z^{\prime }}\mathrm{d}z\hat{u}\right(z\left)\right]}^2,$where $z=\hat{z}/H$, so that $\gamma =8/945$ for the Poiseuille flow and $\gamma =1/20$ for the Couette flow, introduced above. We note that the parameter p quantifies the enhancement of diffusion by Taylor dispersion and that this enhancement is in the longitudinal x-direction, but not in the transverse y-direction. The boundary conditions in the y-direction are given by (5) $\begin{array}{rl}& y_F=1,\theta =0\mathrm{a}\mathrm{s}y\to -\mathrm{\infty },\end{array}$(5) (6) $\begin{array}{rl}& y_F=0,\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(6) The main focus of this investigation is the determination of steady, one-dimensional solutions of the problem (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) )–(6(6) $\begin{array}{rl}& y_F=0,\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(6) ) independent of x and t and their stability. The stability analysis is addressed analytically in the asymptotic limit of infinitely large Zeldovich number in §3–§5 and the predictions are then examined numerically for finite value of β in §6.

3. Formulation in the near-equidiffusional flame (NEF) limit

The stability analysis will be carried out asymptotically in the limit $\text{β}\to \mathrm{\infty }$ using the so-called near-equidiffusional flame (NEF) approximation based on the assumption that the Lewis number deviates little from unity [15, p. 33]. Within this approximation, the reduced Lewis number $l\equiv \text{β}(Le-1)$ is $O\left(1\right)$ as $\text{β}\to \mathrm{\infty }$ and Equations (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) ) and (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+S\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}& =(1+p^2)\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{y}^2}+\omega \end{array}$(4) ) can be in terms of the leading-order temperature ${\theta }^0$ and $h\sim \text{β}(y_F+{\theta }^0-1)$ as (7) $\begin{array}{rl}& {\theta }_t^0+S{\theta }_y^0=(1+p^2){\theta }_{xx}^0+{\theta }_{yy}^0,\end{array}$(7) (8) $\begin{array}{rl}& h_t+S{h}_y=(1+p^2)h_{xx}+h_{yy}+l\left[\right(1-p^2){\theta }_{xx}^0+{\theta }_{yy}^0]\end{array}$(8) which are applicable outside an infinitely thin reaction sheet, given by $y=f(x,t)$ say. The equations are subject to the boundary conditions (9) $\begin{array}{rl}& {\theta }^0=0,h=0\mathrm{a}\mathrm{s}y\to -\mathrm{\infty }\end{array}$(9) (10) $\begin{array}{rl}& {\theta }^0=1,h\mathrm{i}\mathrm{s}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }\end{array}$(10) and the jump conditions (11a) $\begin{array}{rl}& \text{〚}{\theta }^0\text{〛}=0,\text{〚}h\text{〛}=0\end{array}$(11a) (11b) $\begin{array}{rl}& \text{〚}{h}_y\text{〛}+\frac{1+f_x^2(1-p^2)}{1+f_x^2(1+p^2)}l\text{〚}{\theta }_y^0\text{〛}=0\end{array}$(11b) (11c) $\begin{array}{rl}& \sqrt{1+f_x^2(1+p^2)}\text{〚}{\theta }_y^0\text{〛}=-e^{h/2}\end{array}$(11c) applicable at $y=f(x,t)$. Here we have used the notation $\text{〚}\phi \text{〛}=\phi {|}_{y=f^{+}}-\phi {|}_{y=f^{-}}$. It is worth pointing out that jump conditions (11a(11a) $\begin{array}{rl}& \text{〚}{\theta }^0\text{〛}=0,\text{〚}h\text{〛}=0\end{array}$(11a) ), which account for the presence of Taylor dispersion, are derived from an analysis of the structure of the reaction zone. Since the derivation is similar to that reported in [2], its details are not included here.

4. Linear stability analysis for $\text{β}\to \mathrm{\infty }$

We examine herein the stability of the steady planar flame solution, denoted by an overbar, which satisfies Equations (7(7) $\begin{array}{rl}& {\theta }_t^0+S{\theta }_y^0=(1+p^2){\theta }_{xx}^0+{\theta }_{yy}^0,\end{array}$(7) )–(11a(11a) $\begin{array}{rl}& \text{〚}{\theta }^0\text{〛}=0,\text{〚}h\text{〛}=0\end{array}$(11a) ) with $\mathrm{\partial }/\mathrm{\partial }t=0$ and $\mathrm{\partial }/\mathrm{\partial }x=0$ and is given by (12) $S=1,\overline{f}=0,\overline{\theta }=\{\begin{array}{ll}{e}^y& \mathrm{f}\mathrm{o}\mathrm{r}y<0\\ 1& \mathrm{f}\mathrm{o}\mathrm{r}y>0\end{array},\overline{h}=\{\begin{array}{ll}-ly{e}^y& \mathrm{f}\mathrm{o}\mathrm{r}y<0\\ 0& \mathrm{f}\mathrm{o}\mathrm{r}y>0\end{array}.$(12) To this solution, we add infinitesimal normal-mode disturbances such that (13) $\begin{array}{r}\left[\begin{array}{c}f\\ {\theta }^0\\ h\end{array}\right]=\left[\begin{array}{c}0\\ \overline{\theta }\left(y\right)\\ \overline{h}\left(y\right)\end{array}\right]+e^{\sigma t+ikx}\left[\begin{array}{c}\tilde{f}\\ \tilde{\theta }\left(y\right)\\ \tilde{h}\left(y\right)\end{array}\right]\end{array}$(13) where k denotes the real wavenumber and σ a constant (not necessarily real) characterising the growth rate/frequency of the perturbation.

The methodology used to obtain the dispersion relation $\phi (\sigma ,k,l,\lambda )=0$ is classical, see e.g. [10]. We begin by substituting the perturbed solutions (13(13) $\begin{array}{r}\left[\begin{array}{c}f\\ {\theta }^0\\ h\end{array}\right]=\left[\begin{array}{c}0\\ \overline{\theta }\left(y\right)\\ \overline{h}\left(y\right)\end{array}\right]+e^{\sigma t+ikx}\left[\begin{array}{c}\tilde{f}\\ \tilde{\theta }\left(y\right)\\ \tilde{h}\left(y\right)\end{array}\right]\end{array}$(13) ) into Equations (7(7) $\begin{array}{rl}& {\theta }_t^0+S{\theta }_y^0=(1+p^2){\theta }_{xx}^0+{\theta }_{yy}^0,\end{array}$(7) )–(11a(11a) $\begin{array}{rl}& \text{〚}{\theta }^0\text{〛}=0,\text{〚}h\text{〛}=0\end{array}$(11a) ). This leads to an eigenvalue boundary problem for the functions $\tilde{\theta }\left(y\right)$ and $\tilde{h}\left(y\right)$ which is given by (14) $\begin{array}{rl}& {\tilde{\theta }}_{yy}-{\tilde{\theta }}_y-[\sigma +k^2(1+p^2\left)\right]\tilde{\theta }=0,\end{array}$(14) (15) $\begin{array}{rl}& {\tilde{h}}_{yy}-{\tilde{h}}_y-[\sigma +k^2(1+p^2\left)\right]\tilde{h}=-l[{\tilde{\theta }}_{yy}+k^2(p^2-1\left)\tilde{\theta }\right]\end{array}$(15) applicable for $y\ne 0$. These equations are subject to the boundary conditions (16) $\tilde{\theta }=0,\tilde{h}=0\mathrm{a}\mathrm{s}y\to \pm \mathrm{\infty },$(16) and the linearised jump conditions (17) $\text{〚}\tilde{\theta }\text{〛}=\tilde{f},\text{〚}\tilde{h}\text{〛}=-l\tilde{f},\text{〚}{\tilde{h}}_y\text{〛}+l\text{〚}{\tilde{\theta }}_y\text{〛}=-l\tilde{f},\text{〚}{\tilde{\theta }}_y\text{〛}=\tilde{f}-\frac{1}{2}\tilde{h}\left(0^{+}\right)$(17) applicable at y = 0, i.e. $\text{〚}\phi \text{〛}=\phi {|}_{y=0^{+}}-\phi {|}_{y=0^{-}}$.

The solution for $\mathrm{R}\mathrm{e}\{\sigma +k^2(1+p^2\left)\right\}>0$ is given by (18) $\tilde{\theta }=-\tilde{f}\{\begin{array}{l}{e}^{(1+\mathrm{\Gamma })y/2}\\ 0\end{array},\tilde{h}=\tilde{f}\{\begin{array}{ll}[1-\mathrm{\Gamma }+l(1+\chi y\left)\right]e^{(1+\mathrm{\Gamma })y/2}& \mathrm{f}\mathrm{o}\mathrm{r}y<0\\ (1-\mathrm{\Gamma })e^{(1-\mathrm{\Gamma })y/2}& \mathrm{f}\mathrm{o}\mathrm{r}y<0\end{array}$(18) with $\chi =(1+\mathrm{\Gamma }{)}^2/4\mathrm{\Gamma }+{\kappa }^2(2\lambda -1)/\mathrm{\Gamma }$, satisfying the solvability condition or the dispersion relation (19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) Here (20) $\mathrm{\Gamma }=\sqrt{1+4\sigma +4{\kappa }^2},\kappa =k\sqrt{1+p^2},\lambda =\frac{p^2}{1+p^2}\in [0,1].$(20) Equation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) can be solved for $\sigma =\sigma (\kappa ;l,\lambda )$, yielding in general three roots ${\sigma }_1$, ${\sigma }_2$ and ${\sigma }_3$, or less. Clearly, the stability of a mode with wavenumber κ is dictated by the root σ whose real part is equal to $max\left(\mathrm{R}\mathrm{e}\right\{{\sigma }_1\},\mathrm{R}\mathrm{e}\{{\sigma }_2\},\mathrm{R}\mathrm{e}\{{\sigma }_3\left\}\right)$. We shall denote this eigenvalue by (21) ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}={\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}(\kappa ;l,\lambda ).$(21) Furthermore, $\mathrm{R}\mathrm{e}\left\{{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right\}$ attains a maximum value at $\kappa ={\kappa }_m$ say, with corresponding complex growth rate $\sigma ={\sigma }_m$. These values characterise the most unstable mode and are functions of l and λ, (22) ${\sigma }_m={\sigma }_m(l,\lambda )\mathrm{a}\mathrm{n}\mathrm{d}{\kappa }_m={\kappa }_m(l,\lambda ).$(22)

5. Implications of the dispersion relation

In this section, we will examine the implications of the dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) on flame stability using the notations (20(20) $\mathrm{\Gamma }=\sqrt{1+4\sigma +4{\kappa }^2},\kappa =k\sqrt{1+p^2},\lambda =\frac{p^2}{1+p^2}\in [0,1].$(20) )–(22(22) ${\sigma }_m={\sigma }_m(l,\lambda )\mathrm{a}\mathrm{n}\mathrm{d}{\kappa }_m={\kappa }_m(l,\lambda ).$(22) ).

5.1. Illustrative cases and terminology

The function ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\kappa \right)$ assumes various forms depending on the value of l and λ, each one resembling one of the dispersion curves depicted in Figure 2. In this figure, a first category of dispersion curves can be identified for which the maximum growth rate ${\sigma }_m$ is real and occurs for a non-zero value of the wavenumber, ${\kappa }_m\ne 0$; this is the case of subfigures 2(a) to 2(d). In this category, the instability may be termed as a cellular instability since the unstable planar flame solutions are expected to evolve into cellular structures, at least near the instability onset. A second category of dispersion curves can be identified for which the instability may be characterised as being an oscillatory instability, since it corresponds to ${\sigma }_m$ having a non-zero imaginary part. This is the case of subfigures 2(e) to 2(h). The oscillatory instability appears here as a longwave instability when ${\kappa }_m=0$ as in subfigures 2(e) and 2(f) or as a finite-wavelength instability when ${\kappa }_m\ne 0$ as in subfigures 2(g) and 2(h). Note that in the case of subfigure 2(e), a finite wavelength cellular instability can occur instead of the oscillatory longwave instability, if the domain is not large enough. There are also cases such as in subfigure 2(i) where ${\kappa }_m=0$ and ${\sigma }_m$ is real, while ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\left(\kappa \right)$ has a non-zero imaginary part except close to $\kappa =0$. Solutions for such cases exhibit a non-oscillatory longwave instability near onset, at least in sufficiently large domains; in smaller domains, oscillations corresponding to non-zero wavenumbers are however expected near onset.

Figure 2. Growth rate $\mathrm{R}\mathrm{e}\left({\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)$ versus wavenumber κ for selected values of l and λ. Solid lines indicate that ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is real ($\mathrm{I}\mathrm{m}\left({\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)=0$) and dashed lines indicate that ${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}$ has a non-zero imaginary part ($\mathrm{I}\mathrm{m}\left({\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right)\ne 0$). The horizontal and vertical scales are not indicated and are chosen individually for each subfigure for the sake of clarity.

5.2. Stability regime diagram

A convenient way of summarising the various stability results is to delimit instability regions in an l-λ plane, which is equivalent to $Le$-$Pe$ plane, as done in Figure 3. The colour scale in this figure represents the wavenumber ${\kappa }_m$ of the most unstable mode (with complex growth rate ${\sigma }_m$). The solid lines in this figure separate the stable (white) region from the unstable ones, whereas the dashed lines¹ for l>0 separate the region of cellular instability (above these lines) from that of non-cellular instability. Note that the region of non-cellular instability below the dashed lines is subdivided into three subregions by the dotted lines. To the right of the vertical dotted line (l>16) below the dashed lines, we have a non-oscillatory longwave instability and to its left $(l<16)$ we have an oscillatory instability.

Figure 3. Stability regime diagram of a premixed flame propagating transversely to a shear flow in the l-λ plane, which is equivalent to the $Le$-$Pe$ plane. The figure is based on the dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ), obtained in the limit $\text{β}\to \mathrm{\infty }$. The colour scale represents the wavenumber ${\kappa }_m$ of the most unstable mode (with complex growth rate ${\sigma }_m$). The types of the bifurcation curves, represented by solid lines, are discussed in §5.3.

As is well known [11, pp. 477–479], when $\lambda =0$ (or $Pe=0$), the cellular instability emerges for l<−2 and the oscillatory instability for $l>32/3$, as can be seen in the figure. The most striking observation in the presence of Taylor dispersion ($\lambda \ne 0$) is that the cellular instability region comprises positive values of l in addition to negative values. Specifically, cellular instability can now occur, according to our theoretical analysis, for any value of l such that $|l|>2$. It is worth noting in particular that for l>2, the cellular instability can be achieved if the Peclet number is large enough. Such cellular instability appears to be more accessible in real reactive mixtures than the oscillatory instability which requires large values of l, namely $l>32/3$ (in the absence of heat loss). Another noteworthy observation pertinent to subunity Lewis number cases is that the cellular instability is hampered by Taylor dispersion when l<0 and suppressed completely for sufficiently strong shear flow, $\lambda >1/2$.

5.3. Bifurcation curves

To characterise the bifurcations from stable to unstable cases in Figure 3 occurring at the marginal condition $\mathrm{R}\mathrm{e}\left\{\sigma \right\}=0$, we shall adopt the terminology used in [[16], pp. 75–81, [17]]. According to this terminology, the types of bifurcations encountered, as denoted in Figure 3, are type-I${}_s$, type-II${}_s$, type-I${}_o$ and type-III${}_o$ bifurcations. These are illustrated schematically in Figure 4 where the real growth rate $\mathrm{R}\mathrm{e}\left\{\sigma \right(\kappa \left)\right\}$ is plotted for the eigenvalue σ which crosses during the bifurcation from the left-half to the right-half complex plane. In the terminology adopted, the subscript s refers to a stationary or non-oscillatory bifurcation, also known as a zero-eigenvalue bifurcation and corresponds to $\sigma =0$. Similarly, the subscript o refers to an oscillatory bifurcation and corresponds to $\mathrm{R}\mathrm{e}\left(\sigma \right)=0$ with $\mathrm{I}\mathrm{m}\left(\sigma \right)\ne 0$.

Figure 4. A schematic illustration of the four types of bifurcations identified in Figure 3. Plotted is $\mathrm{R}\mathrm{e}\left\{\sigma \right(\kappa \left)\right\}$ as a function of wavenumber κ for the eigenvalue σ which crosses during the bifurcation from the left-half to the right-half complex plane. Solid lines indicate that σ is real ($\mathrm{I}\mathrm{m}\left(\sigma \right)=0$) and dashed lines indicate that σ has a non-zero imaginary part ($\mathrm{I}\mathrm{m}\left(\sigma \right)\ne 0$).

The bifurcation curves in Figure 3 are determined as follows. The bifurcation curve labelled type-II${}_s$ in Figure 3 is obtained from the condition $d^2\sigma /d{\kappa }^2{|}_{\kappa =0}$, which is clear from Figure 4(a), leading to the explicit relation (23) $\lambda =\frac{l+2}{2l}\mathrm{f}\mathrm{o}\mathrm{r}l\le -2\mathrm{a}\mathrm{n}\mathrm{d}2<l\le 6.$(23) Note that this expression is applicable in the domains $l\le -2$ and $2<l\le 6$; for l>6, the boundary between stability and instability regions corresponds to different types of bifurcation curves, as indicated in Figure 3. Further details pertaining to this type-II${}_s$ bifurcation are given in §5.4 mainly dedicated to the derivation of a Kuramoto-Sivashinsky equation.

We turn now to the type-I${}_s$ bifurcation curve identified in Figure 3. The curve is obtained from the conditions $\sigma =0$ and $d\sigma /d\kappa =0$ as seen in Figure 4(b), which determine at the bifurcation the wavenumber² $\kappa ={\kappa }_m\ne 0$ and provide the relation (24) $\lambda =\frac{2}{l}(\sqrt{2l+4}-2)\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}6\le l\le 4(1+\sqrt{3}).$(24) We note that the type-I${}_s$ and type-II${}_s$ boundary curves introduced meet at a tricritical point (25) $(l,\lambda )=(6,2/3).$(25) This terminology is borrowed from phase transition theory [18, p. 493] where a tricritical point refers to a point where first-order and second-order phase transition curves meet; such transition curves are mathematically analogous to our type-I${}_s$ and type-II${}_s$ bifurcation curves.³

Furthermore, the type-III${}_o$ bifurcation curve is determined from the requirements $\mathrm{R}\mathrm{e}\left\{\sigma \right\}=0$, $d\mathrm{R}\mathrm{e}\left\{\sigma \right\}/d\kappa =0$ at $\kappa ={\kappa }_m=0$; see Figure 4(c). The curve is found to be given by the equation (26) $l=4(1+\sqrt{3})\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}\lambda \in [0.3902,0.5646],$(26) which corresponds to a vertical line segment in Figure 3. The fact that this curve is vertical is straightforward consequence of the dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) being independent of λ when $\kappa =0$.

Finally, the type-I${}_o$ bifurcation curve satisfies the conditions $\mathrm{R}\mathrm{e}\left\{\sigma \right\}=0$, $d\mathrm{R}\mathrm{e}\left\{\sigma \right\}/d\kappa =0$ at $\kappa ={\kappa }_m\ne 0$; see Figure 4(d). This curve is computed numerically and is found to extend from the point $(\lambda ,l)=(0,\frac{32}{3})$ to the point $(\lambda ,l)=(0.3902,4(1+\sqrt{3}\left)\right)$.

5.4. Further mathematical implications and weakly nonlinear analysis

The regime diagram of Figure 3 suggests a rich variety of mathematical behaviours in the vicinity of the bifurcation curves and in particular near the tricritical point. Indeed, in the vicinity of the bifurcation curves, weakly nonlinear analyses can be carried out accounting for the coupling between Taylor dispersion and the thermo-diffusive instabilities. For example, in the vicinity of type-II${}_s$ bifurcation curve, a Kuramoto-Sivashinsky (KS) equation can be derived and, in the neighbourhood of type-I${}_s$ bifurcation curve, a Swift-Hohenberg (SH) equation can be derived [17]. More interestingly, near the tricritical point, sixth-order (in space) partial differential equations can be obtained, describing the flame evolution in the weakly nonlinear regime. These mathematical aspects are discussed elsewhere [19]. Here we shall only address the flame dynamics near the type-II${}_s$ bifurcation curve, which is classically described by a KS equation. The KS equation in our case takes the form (27) $f_t+\frac{p^2-1}{2}(l-l_c)f_{xx}+\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}{f}_{xxxx}+\frac{p^2+1}{2}{f}_x^2=0$(27) where (28) $l_c=\frac{2(p^2+1)}{p^2-1}=\frac{2}{2\lambda -1}.$(28) Note that the nonlinear term in (27(27) $f_t+\frac{p^2-1}{2}(l-l_c)f_{xx}+\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}{f}_{xxxx}+\frac{p^2+1}{2}{f}_x^2=0$(27) ) may be derived using a semi-heuristic kinematic argument as explained in [20], and as done in [5] in the similar problem of flames propagating in a longitudinal direction. The details of the derivation will not be repeated here as they are similar to those in [5].

As for the linear part of (27(27) $f_t+\frac{p^2-1}{2}(l-l_c)f_{xx}+\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}{f}_{xxxx}+\frac{p^2+1}{2}{f}_x^2=0$(27) ), this can be simply obtained from the dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ). To this end, we note that Equation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) has always a real root $\sigma \left(k\right)$ such that $\sigma \left(0\right)=0$ and that the Taylor expansion of $\sigma \left(k\right)$ for small values of k is given by (29) $\sigma ={\alpha }_2{k}^2-{\alpha }_4{k}^4+\cdots$(29) where ${\alpha }_2$ and ${\alpha }_4$ are given by ${\alpha }_2=\frac{p^2-1}{2}(l-l_c)\mathrm{a}\mathrm{n}\mathrm{d}{\alpha }_4=\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}.$These are the coefficients of $f_{xx}$ and $f_{xxxx}$ in (27(27) $f_t+\frac{p^2-1}{2}(l-l_c)f_{xx}+\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}{f}_{xxxx}+\frac{p^2+1}{2}{f}_x^2=0$(27) ), since the linear part of (27(27) $f_t+\frac{p^2-1}{2}(l-l_c)f_{xx}+\frac{2(p^2+1{)}^2(p^2-2)}{p^2-1}{f}_{xxxx}+\frac{p^2+1}{2}{f}_x^2=0$(27) ) is equivalent to (29(29) $\sigma ={\alpha }_2{k}^2-{\alpha }_4{k}^4+\cdots$(29) ) in the case of normal modes $f(x,t)\propto e^{\sigma t+ikx}$.

It is worth noting that the onset of the (type-II${}_s$) instability corresponds to the condition ${\alpha }_2=0$ which is equivalent to the condition $d^2\sigma /d{k}^2{|}_{k=0}=0$ used earlier in §5.3. This is so, provided the coefficient ${\alpha }_4$ of the fourth-derivative term is positive. Indeed, if ${\alpha }_4$ were negative, then the critical wavenumber $k_c$ at the onset of instability would be necessarily non-zero, which contradicts the requirement $k_c=0$ of a type-II${}_s$ bifurcation. The condition ${\alpha }_4>0$ requires $p>\sqrt{2}$ or p<1, that is equivalently, $\lambda >2/3$ or $\lambda <1/2$. Therefore the critical condition given by (28(28) $l_c=\frac{2(p^2+1)}{p^2-1}=\frac{2}{2\lambda -1}.$(28) ) is only applicable in this range of λ (or p), which determines the black type-II${}_s$ bifurcation curve in Figure 3. For $\lambda <2/3$ and l>0, cellular flames still exist but the transition from the stable to the unstable region occurs through a type-I${}_s$ bifurcation, represented by the red curve in Figure 3.

6. Computational results for finite Zeldovich number β

The theoretical results discussed so far are all based on the dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) obtained in the asymptotic limit $\text{β}\to \mathrm{\infty }$. In this section, we shall carry out computations with a finite value of β in order to assess the applicability of theoretical predictions, at least qualitatively, and also to examine the nonlinear evolution of unstable flames. To this end, we consider the problem consisting of governing Equations (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) ) and (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+S\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}& =(1+p^2)\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{y}^2}+\omega \end{array}$(4) ) and boundary conditions (5(5) $\begin{array}{rl}& y_F=1,\theta =0\mathrm{a}\mathrm{s}y\to -\mathrm{\infty },\end{array}$(5) )–(6(6) $\begin{array}{rl}& y_F=0,\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(6) ), adopting the numerical values $\text{β}=10$ and $\alpha =0.85$, while varying the parameters $Le$ and $p^2$ (or equivalently $\lambda =p^2/(1+p^2)$).

6.1. Preliminary considerations

Before presenting the finite-β numerical results in the next two subsections, we make a few preliminary remarks which will facilitate the discussion and interpretation of the results. We begin by noting upon examining Equations (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) ) and (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+S\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}& =(1+p^2)\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{y}^2}+\omega \end{array}$(4) ) that the presence of a flow-dependent effective Lewis number in the x-direction, as noted in [2, 4], namely, (30) ${\text{{Le}}}_x=\frac{\text{{Le}}(1+\gamma {\text{{Pe}}}^2)}{1+\gamma {\text{{Pe}}}^2{\text{{Le}}}^2}=\frac{\text{{Le}}(1+p^2)}{1+p^2{\text{{Le}}}^2}=\frac{\text{{Le}}}{1+\lambda ({\text{{Le}}}^2-1)},$(30) where use has been made of the definitions $p^2=\gamma {Pe}^2$ and $\lambda =p^2/(1+p^2)$.

An important implication of formula (30(30) ${\text{{Le}}}_x=\frac{\text{{Le}}(1+\gamma {\text{{Pe}}}^2)}{1+\gamma {\text{{Pe}}}^2{\text{{Le}}}^2}=\frac{\text{{Le}}(1+p^2)}{1+p^2{\text{{Le}}}^2}=\frac{\text{{Le}}}{1+\lambda ({\text{{Le}}}^2-1)},$(30) ) is that ${Le}_x\to Le$ as $Pe\to 0$ and $Le\to 1/Le$ as $Pe\to \mathrm{\infty }$. In other words, the effect of a shear flow is such that a weakly diffusing reactant ($Le>1$) appears effectively as strongly diffusing (in the x-direction, ${Le}_x<1$) and a strongly diffusing reactant ($Le<1$) appears effectively as weakly diffusing (${Le}_x>1$) provided the Peclet number is large enough, more precisely when $Pe$ exceeds the value $1/\sqrt{\gamma Le}$. In fact, this value is determined from the fact that ${Le}_x=1$ irrespective of $Le$ when $Pe=1/\sqrt{\gamma Le}$, which follows from (30(30) ${\text{{Le}}}_x=\frac{\text{{Le}}(1+\gamma {\text{{Pe}}}^2)}{1+\gamma {\text{{Pe}}}^2{\text{{Le}}}^2}=\frac{\text{{Le}}(1+p^2)}{1+p^2{\text{{Le}}}^2}=\frac{\text{{Le}}}{1+\lambda ({\text{{Le}}}^2-1)},$(30) ).

6.2. Linear stability analysis based on eigen-boundary value problem

We begin by examining the linear stability of the planar flame solution with the finite value $\text{β}=10$ adopted. This base solution satisfies Equations (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) ) and (4(4) $\begin{array}{rl}\frac{\mathrm{\partial }\theta }{\mathrm{\partial }t}+S\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}& =(1+p^2)\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2\theta }{\mathrm{\partial }{y}^2}+\omega \end{array}$(4) ) with $\mathrm{\partial }/\mathrm{\partial }t=\mathrm{\partial }/\mathrm{\partial }x=0$. The scaled burning speed S, which is equal to unity in the limit $\text{β}\to \mathrm{\infty }$, is now computed numerically and plotted as a function of the Lewis number $Le$ in Figure 5. If we denote the dependent variables of the base state as ${\overline{y}}_F\left(y\right)$ and $\overline{\theta }\left(y\right)$, then the stability of the base solution to small perturbations can be studied by introducing (31) $\left[\begin{array}{c}{y}_F\\ \theta \end{array}\right]=\left[\begin{array}{c}{\overline{y}}_F\\ \overline{\theta }\end{array}\right]+e^{ik\sqrt{1+p^2}x+\sigma t}\left[\begin{array}{c}{\tilde{y}}_F\\ \tilde{\theta }\end{array}\right]$(31) into the problem (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) )–(6(6) $\begin{array}{rl}& y_F=0,\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(6) ). The linear stability is then described by the eigen-boundary value problem (32) $\begin{array}{rl}& \frac{1}{\text{{Le}}}\frac{{\mathrm{d}}^2{\tilde{y}}_F}{\mathrm{d}{y}^2}-S\frac{\mathrm{d}{\tilde{y}}_F}{\mathrm{d}y}-\frac{k^2}{{\text{{Le}}}_x}{\tilde{y}}_F-\{{\tilde{y}}_F+\frac{\text{β}{\overline{y}}_F\tilde{\theta }}{[1-\alpha (1-\overline{\theta }){]}^2}\}\frac{{\text{β}}^2}{2\text{{Le}}}\exp \left[\frac{\text{β}(\overline{\theta }-1)}{1+\alpha (\overline{\theta }-1)}\right]=\sigma {\tilde{y}}_F,\end{array}$(32) (33) $\begin{array}{rl}& \frac{{\mathrm{d}}^2\tilde{\theta }}{\mathrm{d}{y}^2}-S\frac{\mathrm{d}\tilde{\theta }}{\mathrm{d}y}-k^2\tilde{\theta }+\{{\tilde{y}}_F+\frac{\text{β}{\overline{y}}_F\tilde{\theta }}{[1-\alpha (1-\overline{\theta }){]}^2}\}\frac{{\text{β}}^2}{2\text{{Le}}}\exp \left[\frac{\text{β}(\overline{\theta }-1)}{1+\alpha (\overline{\theta }-1)}\right]=\sigma \tilde{\theta }\end{array}$(33) with (34) $\begin{array}{rl}& {\tilde{y}}_F=0,\tilde{\theta }=0\mathrm{a}\mathrm{s}y\to -\mathrm{\infty }\end{array}$(34) (35) $\begin{array}{rl}& {\tilde{y}}_F=0,\frac{\mathrm{d}\tilde{\theta }}{\mathrm{d}y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(35) It is worth noting the presence in the equations of the longitudinal effective Lewis number ${Le}_x$ defined in (30(30) ${\text{{Le}}}_x=\frac{\text{{Le}}(1+\gamma {\text{{Pe}}}^2)}{1+\gamma {\text{{Pe}}}^2{\text{{Le}}}^2}=\frac{\text{{Le}}(1+p^2)}{1+p^2{\text{{Le}}}^2}=\frac{\text{{Le}}}{1+\lambda ({\text{{Le}}}^2-1)},$(30) ), in addition to the Lewis number $Le$. The problem possesses a discrete spectrum of eigenvalues σ, in which the eigenvalue with the maximum real part (growth rate) determines the stability of the base solution. The corresponding computational results are summarised in Figure 6, where the stability-instability regions are delimited in the $Le$-λ plane. The solid lines in the figure separates the stable (white) region from the unstable (grey) regions. The dash-dotted line represents the condition ${Le}_x=1$, which takes the form $\lambda =1/(Le+1)$ on using $\lambda =p^2/(1+p^2)$ in (30(30) ${\text{{Le}}}_x=\frac{\text{{Le}}(1+\gamma {\text{{Pe}}}^2)}{1+\gamma {\text{{Pe}}}^2{\text{{Le}}}^2}=\frac{\text{{Le}}(1+p^2)}{1+p^2{\text{{Le}}}^2}=\frac{\text{{Le}}}{1+\lambda ({\text{{Le}}}^2-1)},$(30) ).

Figure 5. The scaled planar-flame burning speed $S=S_L/S_L^0$ versus the Lewis number $Le$, for $\text{β}=10$ and $\alpha =0.85$.

Figure 6. Stability regime diagram of a premixed flame propagating transversely to a shear flow in the $Le$-λ plane. The figure is based on the computation of the eigenvalues of problem (32(32) $\begin{array}{rl}& \frac{1}{\text{{Le}}}\frac{{\mathrm{d}}^2{\tilde{y}}_F}{\mathrm{d}{y}^2}-S\frac{\mathrm{d}{\tilde{y}}_F}{\mathrm{d}y}-\frac{k^2}{{\text{{Le}}}_x}{\tilde{y}}_F-\{{\tilde{y}}_F+\frac{\text{β}{\overline{y}}_F\tilde{\theta }}{[1-\alpha (1-\overline{\theta }){]}^2}\}\frac{{\text{β}}^2}{2\text{{Le}}}\exp \left[\frac{\text{β}(\overline{\theta }-1)}{1+\alpha (\overline{\theta }-1)}\right]=\sigma {\tilde{y}}_F,\end{array}$(32) )–(35(35) $\begin{array}{rl}& {\tilde{y}}_F=0,\frac{\mathrm{d}\tilde{\theta }}{\mathrm{d}y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(35) ) with $\text{β}=10$ and $\alpha =0.85$. The grey region corresponds to unstable flames and the white to stable flames. Compare with Figure 3, where analogous notations are used.

Comparing Figure 6, computed for $\text{β}=10$, with Figure 3 corresponding to $\text{β}\to \mathrm{\infty }$, the following conclusions can be drawn. First, we note the good qualitative agreement between the two figures, notably, regarding the bifurcation curves separating the stability from instability regions. In particular note the presence of a tricritical point at $(Le,\lambda )\approx (1.94,0.47)$, identified in Figure 3. Second, the computations confirm the existence of the cellular instability region for $Le>1$, an original finding of the asymptotic analysis. More precisely, two cellular instability regions appear in Figure 6, whose locations are compatible with the conditions $Pe<1/\sqrt{\gamma Le}$ when $Le<1$ and $Pe>1/\sqrt{\gamma Le}$ when $Le>1$, or equivalently ${Le}_x<1$ for any $Le$, as argued in §6.1. Further the flame is found to be stable for all values of λ when $Le\in (0.84,1.19)$ in the case $\text{β}=10$, while the equivalent stability range predicted by the asymptotic analysis, $l\in (-2,2)$, provides the stability range $Le=1+l/\text{β}\in (0.8,1.2)$ for $\text{β}=10$.

6.3. Time-dependent numerical simulations

Finally, we present time-dependent numerical simulations simply to illustrate the occurrence of the cellular instability for $Le>1$ and the flame long-time evolution. The simulations are based on the numerical solution of the time-dependent problem (3(3) $\begin{array}{rl}\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }t}+S\frac{\mathrm{\partial }{y}_F}{\mathrm{\partial }y}& =\frac{1}{\text{{Le}}}\left[\right(1+p^2{\text{{Le}}}^2)\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{y}_F}{\mathrm{\partial }{y}^2}]-\omega ,\end{array}$(3) )–(6(6) $\begin{array}{rl}& y_F=0,\frac{\mathrm{\partial }\theta }{\mathrm{\partial }y}=0\mathrm{a}\mathrm{s}y\to +\mathrm{\infty }.\end{array}$(6) ), with periodic conditions in the x-direction, starting from an initial condition corresponding to a steady, planar premixed flame. The computations are carried out using COMSOL Multiphysics software as described in [4, 5]. Two cases are considered corresponding to $(Le,\lambda )=(1.94,0.6)$ and $(Le,\lambda )=(1.94,0.9)$, whose results are shown in Figures 7 and 8, respectively. Shown are reaction-rate $\left(\omega \right)$ fields at selected values of time t. These illustrate the development of the instability and the flame evolution into a cellular structure, which settles into an apparently stable state at large times.

Figure 7. Fields of reaction rate ω at selected values of time t, computed for $(Le,\lambda )=(1.94,0.6)$, $\text{β}=10$ and $\alpha =0.85$. The initial condition for the time-dependent calculations corresponds to a steady, planar premixed flame, computed numerically.

Figure 8. Fields of reaction rate ω at selected values of time t, computed for $(Le,\lambda )=(1.94,0.9)$, $\text{β}=10$ and $\alpha =0.85$. The initial condition for the time-dependent calculations corresponds to a steady, planar premixed flame, computed numerically.

7. Conclusions

In this paper, we have examined the effect of Taylor dispersion, or shear-enhanced diffusion on the stability of a premixed flame, propagating in a direction transverse to that of unidirectional shear flow. A simple thermo-diffusive model is adopted corresponding to a Hele-Shaw channel, whose walls are assumed to be adiabatic and closely spaced. In this configuration, a simple shear flow is prescribed, corresponding for example to a Couette flow. Upon depth-averaging, the problem is governed by two-dimensional transport equations in which diffusion is anisotropic with shear-enhanced diffusion coefficients in the longitudinal flow direction. A linear stability analysis is carried out in the limit $\text{β}\to \mathrm{\infty }$, where β is the Zeldovich number. A simple dispersion relation (19(19) $2{\mathrm{\Gamma }}^2(1-\mathrm{\Gamma })+l(1-\mathrm{\Gamma }+2\sigma +4\lambda {\kappa }^2)=0.$(19) ) is derived analytically involving the Lewis number $Le$ and the Peclet number $Pe$ as parameters. A stability regime diagram (Figure 3) is constructed in the $Le$-$Pe$ plane, which categorises the various instabilities and bifurcations encountered.

A remarkable finding is the demonstration that the classical cellular instability, which is usually expected to occur in $Le<1$ mixtures, exists now only for $Le>1$ mixtures when the Peclet number exceeds a critical value. The necessary conditions needed to observe cellular flames in $Le>1$ mixtures include that the geometry must be slender such as in the case of a Hele-Shaw channel, that the Peclet number $Pe$ must exceed a critical value above $1/\sqrt{\gamma Le}$ and that the flame must be aligned with the flow direction. Furthermore, the cellular instability identified for $Le>1$ is shown to occur either through a finite-wavelength stationary bifurcation (also known as type-I${}_s$) or through a longwave stationary bifurcation (also known as type-II${}_s$). In the weakly nonlinear regime, a Kuramoto-Sivashinsky equation in the vicinity of type-II${}_s$ bifurcation is derived. Moreover, it is found that the oscillatory instability, usually encountered in $Le>1$ mixtures persists under the influence of Taylor dispersion if the Peclet number is below a critical value and disappears above this value.

The stability results aforementioned, which follow from the dispersion relation obtained in the limit $\text{β}\to \mathrm{\infty }$, are complemented by numerical computations carried out for a finite value of the Zeldovich number. The computations involve the determination of the eigenvalues of a linear stability boundary-value problem and numerical simulations of the time-dependent governing partial differential equations. The numerical results are found to be in good qualitative agreement with the analytical predictions. In particular, a stability regime diagram (Figure 6) computed for $\text{β}=10$ is found to be consistent with the stability regime diagram of Figure 3, corresponding to $\text{β}\to \mathrm{\infty }$.

To close this paper, we note that the stability results obtained for the case of transverse flame propagation differ markedly from those obtained for flames propagating in the longitudinal direction [3]. For example, the appearance of cellular instability in $Le>1$ mixtures when the Peclet number is above a critical value, does not occur in the case of longitudinal propagation. In fact, in the latter case, only stable flames are encountered when $Le>1$ and $Pe\gg 1$. A natural follow-up of the current work is to extend the stability analysis to flames propagating in an arbitrary direction relative to that of the shear flow. It is also worthwhile to investigate the influence of heat losses on our predictions, following [5], since such influence is significant in practice. Interestingly, the cellular instability identified herein for $Le>1$ mixtures has been recently shown to also take place for diffusion flames aligned with the direction of the shear flow [6]. This unexpected instability has been proposed in [6] as a plausible mechanism for the formation of diffusion flame streets observed in experiments [7–9]. It is desirable to have similar experiments in the premixed case to test the presence of the cellular instability identified herein for $Le>1$ mixtures in the presence of a strong shear flow.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This research was supported by the UK Engineering and Physical Sciences Research Council (EPSRC) through [grant number EP/V004840/1].

Notes

1 The dashed curves meet at a cusp located at $(l,\lambda )=(16,\frac{1}{2})$. The crossing from the non-cellular to the cellular region as λ is increased involves a discontinuous jump in ${\kappa }_m$ for l<16 and a continuous transition for l>16. Similar transitions have been observed in diffusion flames which have been described with more detail in [6].

2 An explicit relation for ${\kappa }_m\left(l\right)$ can be obtained by substituting the expression for λ in (24(24) $\lambda =\frac{2}{l}(\sqrt{2l+4}-2)\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}6\le l\le 4(1+\sqrt{3}).$(24) ) into the equation ${\kappa }_m^2=\left[\right(\lambda l+2+\sqrt{(\lambda l+2{)}^2-6l})/12{]}^2-1/4$, which follows from the condition $d\sigma /d\kappa =0$. From this relation, we can deduce that ${\kappa }_m\sim \sqrt{l-6}$ as $l\to 6^{+}$.

3 As in phase transition theory [18, p. 496], the curves of type-I${}_s$ (24(24) $\lambda =\frac{2}{l}(\sqrt{2l+4}-2)\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{l}\mathrm{i}\mathrm{c}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{f}\mathrm{o}\mathrm{r}6\le l\le 4(1+\sqrt{3}).$(24) ) and type-II${}_s$ (23(23) $\lambda =\frac{l+2}{2l}\mathrm{f}\mathrm{o}\mathrm{r}l\le -2\mathrm{a}\mathrm{n}\mathrm{d}2<l\le 6.$(23) ) are continuous and have continuous first derivatives but discontinuous second and higher-order derivatives at the tricritical point.