1. Introduction
Premixed flames are unstable to perturbations owing to gas expansion effects in what is called the Darrieus–Landau or hydrodynamic instability (DL) (Darrieus Reference Darrieus1938; Landau Reference Landau1945). There have been numerous analytical, numerical and experimental studies of premixed flames. The basic configuration corresponds to a nominally planar flame propagating against a mixture of reactants (or stabilized by an opposing flow matching the flame speed). Apart from the DL instability, flames can be unstable owing to thermo-diffusive effects (different diffusivity of heat and mass) (Sivashinsky Reference Sivashinsky1976; Clavin & Williams Reference Clavin and Williams1982; Matalon & Matkowsky Reference Matalon and Matkowsky1982; Pelcé & Clavin Reference Pelcé and Clavin1982), flow interaction (Buckmaster Reference Buckmaster1977; Matalon Reference Matalon1983; Searby & Clavin Reference Searby and Clavin1986), and heat loss to walls and gravity (Joulin & Sivashinsky Reference Joulin and Sivashinsky1994). The most amplified wavenumber depends, among other things, on the stabilizing role of curvature and strain (through diffusion) that can be parameterized in linear stability analysis by one or two length scales called the Markstein lengths (Markstein Reference Markstein1964; Clavin & Williams Reference Clavin and Williams1982; Clavin & Graña-Otero Reference Clavin and Graña-Otero2011), that are of the order of the preheat-zone flame thickness. Furthermore, after the linear phase of instability one encounters a nonlinear phase where the amplitude of the flame stabilizes (Sivashinsky Reference Sivashinsky1977, Reference Sivashinsky1990; Creta & Matalon Reference Creta and Matalon2011a ). Further reviews on laminar and turbulent premixed flames can be found in Sivashinsky (Reference Sivashinsky1983), Matalon (Reference Matalon2007), Driscoll (Reference Driscoll2008) and Lipatnikov & Chomiak (Reference Lipatnikov and Chomiak2010).
Analytical studies can start with the so-called hydrodynamic limit where the stability properties of the flame can be explained quite accurately by ignoring the internal structure of the flame and considering it as an infinitesimally thin jump in density from the high value upstream of the flame to the low value downstream of the flame front. This works well for Lewis numbers close to unity, where the Lewis number is defined as the ratio of the thermal to the mass diffusivity of the species. The effect of non-uniformities in the incoming flow has been considered in numerous studies, e.g. Clavin & Williams (Reference Clavin and Williams1982), Ashurst & Sivashinsky (Reference Ashurst and Sivashinsky1991), Berestycki & Sivashinsky (Reference Berestycki and Sivashinsky1991), Brailovsky & Sivashinsky (Reference Brailovsky and Sivashinsky1993, Reference Brailovsky and Sivashinsky1995), Aldredge (Reference Aldredge1996), Kagan, Sivashinsky & Makhviladze (Reference Kagan, Sivashinsky and Makhviladze1998), Kortsarts et al. (Reference Kortsarts, Kliakhandeler, Shtilman and Sivashinsky1998), Audoly, Berestycki & Pomeau (Reference Audoly, Berestycki and Pomeau2000), Kagan & Sivashinsky (Reference Kagan and Sivashinsky2000). These studies considered essentially the effects induced by changes in the streamwise (flame normal) aligned velocity (some in the small perturbation limit, others nonlinearly). Here, we consider a new case where the flame tangential velocity is not uniform. This can occur when a flame propagates into a region of free shear flow. The purpose of this communication is the analytical description of the simplest situation: near unity Lewis numbers; no other flow but that induced by shear; and completely linear behaviour (early time). With sufficient extensions, this configuration may be useful in understanding more complex problems in turbulent premixed flames. Recent direct numerical simulations (DNS) by Hawkes et al. (Reference Hawkes, Chatakonda, Kolla, Kerstein and Chen2012), where a set of twin premixed flames propagate into a region of turbulence created by planar shear layers, are an example of such a flow.
Sketch of flame configuration.
2. Flame model and mathematical formulation
The hydrodynamic linear stability analysis assumes that a planar flame propagates at a constant speed in the direction normal to itself, which is counterbalanced by an opposite flow of premixed reactants. The problem configuration is illustrated in figure 1. This generalization has the feature that there is a transverse pressure gradient orthogonal to the density gradient and it has a background of constant (initial) vorticity. The flame is considered as an infinitesimally thin sheet where fast reaction occurs. Effects of buoyancy are neglected. Given these assumptions, the only influence of combustion in this problem is thermal expansion of the flow and the consequential discontinuities of thermal properties. As is well known, thermal expansion deflects the downstream streamlines from the parallel upstream streamlines incoming towards the flame. The downstream streamlines are also parallel far from the flame. The present configuration is a generalization of the classical Darrieus–Landau problem. The present extension incorporates a transverse shear flow parallel to the flame plane which can be described analytically by
$v=Sx$
, where
$S$
is a constant. This flow is still divergence free. The flame front is described by the implicit relation
$G=x-F(y,t)=0$
and the flame normal and local flame velocity in the case of small perturbation are given by
$$\begin{eqnarray}\displaystyle \boldsymbol{n}=\frac{\boldsymbol{{\rm\nabla}}G}{|\boldsymbol{{\rm\nabla}}G|}=\frac{1}{|\boldsymbol{{\rm\nabla}}G|}\left(1,-\frac{\partial F}{\partial y}\right),\quad s_{f}=\frac{1}{|\boldsymbol{{\rm\nabla}}G|}\frac{\partial F}{\partial t}, & & \displaystyle\end{eqnarray}$$
where
$\boldsymbol{n}$
is the flame normal pointing towards the burned gases and
$s_{f}$
denotes the flame surface velocity. The flame front separates the neighbourhood of the flame into unburned (left) and burned (right) regions. Fluid properties of interest, including flow velocity
$u$
and
$v$
, pressure
$p$
, and density
${\it\rho}$
, are denoted with subscripts
$(\cdot )_{-}$
and
$(\cdot )_{+}$
in unburned and burned gas regions, respectively. The thermal expansion ratio
$r$
is defined as the ratio of burned gas density to unburned gas density, and it is therefore less than unity. The local flame speed,
$s_{L}$
, is assumed given by
$$\begin{eqnarray}\displaystyle s_{L}=s_{L}^{o}\left(1-\mathscr{L}{\bf\kappa}+\mathscr{L}_{s}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\right), & & \displaystyle\end{eqnarray}$$
where
$s_{L}^{o}$
denotes the burning velocity of the laminar planar flame,
${\it\kappa}$
denotes the curvature of the flame front, and
$\mathscr{L}$
and
$\mathscr{L}_{s}$
denote the Markstein lengths due to curvature and strain, respectively, see Clavin & Graña-Otero (Reference Clavin and Graña-Otero2011). The Markstein lengths have a characteristic magnitude of the order of the laminar flame thickness.
In order to reduce the problem to its simplest form, quantities are non-dimensionalized using the laminar planar flame thickness
$l_{F}$
and burning velocity
$s_{L}^{o}$
, giving the following non-dimensional quantities:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle u^{\ast }=\frac{u}{s_{L}^{o}},\quad v^{\ast }=\frac{v}{s_{L}^{o}},\quad p^{\ast }=\frac{p}{{\it\rho}_{u}s_{L}^{o2}},\quad {\it\rho}^{\ast }=\frac{{\it\rho}}{{\it\rho}_{u}},\quad \mathscr{M}=\frac{\mathscr{L}}{l_{F}},\\ \displaystyle x^{\ast }=\frac{x}{l_{F}},\quad y^{\ast }=\frac{y}{l_{F}},\quad t^{\ast }=\frac{t}{l_{F}/s_{L}^{o}},\quad S^{\ast }=\frac{Sl_{F}}{s_{L}^{o}},\quad \mathscr{M}_{s}=\frac{\mathscr{L}_{s}}{l_{F}}.\end{array}\!\!\right\} & & \displaystyle\end{eqnarray}$$
To simplify further the notation, we will drop the ‘*’ superscript since it will not induce any confusion later. The results will depend on the Markstein numbers
$\mathscr{M}$
and
$\mathscr{M}_{s}$
. Under this normalization, the flow velocity ahead of the flame is
$U_{-}=1$
and the density is
${\it\rho}_{-}=1$
, while behind the flame we have
$U_{+}$
(to be determined below from the jump conditions) and
${\it\rho}_{+}=r$
.
The flow is governed by Euler’s equations (in this inviscid approximation) at both sides of the flame sheet, such that
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial x}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{{\it\rho}}\frac{\partial p}{\partial y}. & \displaystyle\end{eqnarray}$$
The pressure and velocity at either side of the flame are related by the corresponding jump conditions (Williams Reference Williams1985), which are derived from conservation of mass and momentum across the flame, giving
$$\begin{eqnarray}\displaystyle & \displaystyle (r-1)\frac{\partial F}{\partial t}=r\biggl(u_{+}-\frac{\partial F}{\partial y}v_{+}\biggr)-\biggl(u_{-}-\frac{\partial F}{\partial y}v_{-}\biggr), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (ru_{+}-u_{-})\frac{\partial F}{\partial t}=ru_{+}\biggl(u_{+}-\frac{\partial F}{\partial y}v_{+}\biggr)-u_{-}\biggl(u_{-}-\frac{\partial F}{\partial y}v_{-}\biggr)+p_{+}-p_{-}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (rv_{+}-v_{-})\frac{\partial F}{\partial t}=rv_{+}\biggl(u_{+}-\frac{\partial F}{\partial y}v_{+}\biggr)-v_{-}\biggl(u_{-}-\frac{\partial F}{\partial y}v_{-}\biggr)-\frac{\partial F}{\partial y}(p_{+}-p_{-}). & \displaystyle\end{eqnarray}$$
$v$
terms which cannot be neglected from the outset. Using (2.7) for the flat flame reference conditions we obtain
$U_{+}=1/r$
. Then, one can obtain from the momentum conservation equations that
$p_{\pm }=-Sy+p_{0\pm }$
. Note that it is not known a priori that the shear at both sides of the flame,
$S_{\pm }$
, is the same, but one can verify that the jump conditions are solvable only if
$S_{-}=S_{+}=S$
. The quantitative relation between constants of integration
$p_{-}$
and
$p_{+}$
remains unknown. Given the equalities above, pressure terms can be used to substitute horizontal velocity terms into the jump conditions, (2.8), providing 
$$\begin{eqnarray}\displaystyle \frac{1}{r}-1+p_{0+}-p_{0-}=0. & & \displaystyle\end{eqnarray}$$
If one sets the additive constant
$p_{0-}=0$
, then
$p_{0+}$
is equal to
$(r-1)/r$
.
3. Perturbation analysis
Following the notation in Peters (Reference Peters2010) for the classical DL analysis, we introduce a perturbative expansion on a small parameter
${\it\epsilon}$
denoting the amplitude of a deformation in the original flat flame surface. The fields are now written as
$$\begin{eqnarray}\displaystyle u_{-}=1+{\it\epsilon}\hat{u} _{-}+O({\it\epsilon}^{2}),\quad u_{+}=\frac{1}{r}+{\it\epsilon}\hat{u} _{+}+O({\it\epsilon}^{2}), & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle v_{-}=Sx+{\it\epsilon}\hat{v}_{-}+O({\it\epsilon}^{2}),\quad v_{+}=Sx+{\it\epsilon}\hat{v}_{+}+O({\it\epsilon}^{2}), & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle p_{-}=-Sy+{\it\epsilon}\hat{p}_{-}+O({\it\epsilon}^{2}),\quad p_{+}=1-\frac{1}{r}-Sy+{\it\epsilon}\hat{p}_{+}+O({\it\epsilon}^{2}). & & \displaystyle\end{eqnarray}$$
Furthermore, we introduce a coordinate transformation to simplify the analysis, given by
$x={\it\xi}+F({\it\eta},{\it\tau})$
,
$y={\it\eta}$
, and
$t={\it\tau}$
. In this frame of reference, the flame surface
$G=0$
corresponds to
${\it\xi}=0$
. Finally, we write
$F={\it\epsilon}f+O({\it\epsilon}^{2})$
and take
$f$
to be a quantity of order one. Introducing these expansions in (2.4)–(2.6) results to first order in
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \hat{u} }{\partial {\it\xi}}+\frac{\partial \hat{v}}{\partial {\it\eta}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial \hat{u} }{\partial {\it\tau}}+\frac{1}{r}\frac{\partial \hat{u} }{\partial {\it\xi}}+S{\it\xi}\frac{\partial \hat{u} }{\partial {\it\eta}}+\frac{1}{r}\frac{\partial \hat{p}}{\partial {\it\xi}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial }{\partial {\it\tau}}\left(Sf+\hat{v}\right)+\frac{1}{r}\frac{\partial \hat{v}}{\partial {\it\xi}}+S\hat{u} +S{\it\xi}\frac{\partial \hat{v}}{\partial {\it\eta}}+\frac{1}{r}\frac{\partial \hat{p}}{\partial {\it\eta}}=0, & \displaystyle\end{eqnarray}$$
$\hat{u} ({\it\xi},{\it\eta},r,{\it\tau})$
,
$\hat{v}({\it\xi},{\it\eta},r,{\it\tau})$
and
$\hat{p}({\it\xi},{\it\eta},r,{\it\tau})$
can be used to represent both solutions to the left and right of
${\it\xi}=0$
:
$\hat{u} _{-}=\hat{u} ({\it\xi},{\it\eta},1,{\it\tau})$
for
${\it\xi}<0$
and
$\hat{u} _{+}=\hat{u} ({\it\xi},{\it\eta},r,{\it\tau})$
for
${\it\xi}>0$
; and similarly for the other fields. Introducing the perturbation expansion in the jump conditions, one gets to first order (after evaluating the jump condition at
${\it\xi}=0$
) 
$$\begin{eqnarray}\displaystyle & \displaystyle (r-1)\frac{\partial f}{\partial {\it\tau}}=r\hat{u} _{+}-\hat{u} _{-}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=2(\hat{u} _{+}-\hat{u} _{-})+\hat{p}_{+}-\hat{p}_{-}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle 0=\hat{v}_{+}-\hat{v}_{-}-\left(1-\frac{1}{r}\right)\frac{\partial f}{\partial {\it\eta}}. & \displaystyle\end{eqnarray}$$
These equations are coupled with the kinematic equation of the flame surface, which requires
$$\begin{eqnarray}\displaystyle s_{f}=\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{n}-s_{L}. & & \displaystyle\end{eqnarray}$$
This condition is imposed using a level set equation, given by
$$\begin{eqnarray}\displaystyle \frac{\partial G}{\partial t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}G=s_{L}|\boldsymbol{{\rm\nabla}}G|, & & \displaystyle\end{eqnarray}$$
where the flame velocity
$s_{L}$
incorporates the effects of curvature and strain, according to
$$\begin{eqnarray}\displaystyle s_{L}=1-\mathscr{M}{\bf\kappa}+\mathscr{M}_{s}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u\boldsymbol{\cdot }\boldsymbol{n}, & & \displaystyle\end{eqnarray}$$
where
${\bf\kappa}$
denotes the curvature of the flame,
${\bf\kappa}=-\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{n}$
. With all these elements, we obtain
$$\begin{eqnarray}\displaystyle \frac{\partial G}{\partial t}+u\frac{\partial G}{\partial x}+v\frac{\partial G}{\partial y}=(1-\mathscr{M}{\bf\kappa}+\mathscr{M}_{s}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}u\boldsymbol{\cdot }\boldsymbol{n})\sqrt{\left(\frac{\partial G}{\partial x}\right)^{2}+\left(\frac{\partial G}{\partial y}\right)^{2}}. & & \displaystyle\end{eqnarray}$$
Introducing the perturbation expansion, we obtain to first order in
${\it\epsilon}$
, the following boundary condition:
$$\begin{eqnarray}\displaystyle \hat{u} _{-}|_{{\it\xi}=0}=\frac{\partial f}{\partial {\it\tau}}-\mathscr{M}\frac{\partial ^{2}f}{\partial {\it\eta}^{2}}+\mathscr{M}_{s}\left(\frac{\partial \hat{u} _{-}}{\partial {\it\xi}}\bigg|_{{\it\xi}=0}-S\frac{\partial f}{\partial {\it\eta}}\right). & & \displaystyle\end{eqnarray}$$
4. Solution and dispersion relation
We now seek a general solution of the above equations and boundary condition of the form
$$\begin{eqnarray}\displaystyle {\it\psi}=\tilde{{\it\psi}}({\it\xi})\exp ({\it\sigma}{\it\tau}-\text{i}k{\it\eta}), & & \displaystyle\end{eqnarray}$$
where
${\it\psi}$
denotes any of the unburned or burned side fields (
$\hat{u} _{\pm }$
,
$\hat{v}_{\pm }$
,
$\hat{p}_{\pm }$
) and
$$\begin{eqnarray}\displaystyle f=\tilde{f}\exp ({\it\sigma}{\it\tau}-\text{i}k{\it\eta}). & & \displaystyle\end{eqnarray}$$
The parameter
${\it\sigma}$
is termed here the amplification factor and
$k$
denotes the wavenumber of the perturbation. The system of seven linear differential equations is closed for
$\tilde{u} _{\pm }$
,
$\tilde{v}_{\pm }$
,
$\tilde{p}_{\pm }$
and
$\tilde{f}$
, giving
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\tilde{u} }{\text{d}{\it\xi}}-\text{i}k\tilde{v}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle ({\it\sigma}-\text{i}kS{\it\xi})\tilde{u} +\frac{1}{r}\frac{\text{d}\tilde{u} }{\text{d}{\it\xi}}+\frac{1}{r}\frac{\text{d}\tilde{p}}{\text{d}{\it\xi}}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle S{\it\sigma}\tilde{f}+\tilde{u} S+({\it\sigma}-\text{i}kS{\it\xi})\tilde{v}+\frac{1}{r}\frac{\text{d}\tilde{v}}{\text{d}{\it\xi}}-\text{i}k\frac{1}{r}\tilde{p}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{u} _{-}|_{{\it\xi}=0}-\mathscr{M}_{s}\frac{\partial \hat{u} _{-}}{\partial {\it\xi}}\bigg|_{{\it\xi}=0}=({\it\sigma}+\mathscr{M}k^{2}+\text{i}k\mathscr{M}_{s}S)\tilde{f}. & & \displaystyle\end{eqnarray}$$
The general solution of
$\tilde{u} ({\it\xi},r)$
can be obtained from this closed system of differential equations by eliminating
$\tilde{v}$
(using (4.3)) and
$\tilde{p}$
(using (4.5)) and has the general form
$$\begin{eqnarray}\displaystyle \tilde{u} ({\it\xi},r)=c_{1}\text{e}^{-k{\it\xi}}+c_{2}\text{e}^{k{\it\xi}}+c_{3}\left\{\text{e}^{-k{\it\xi}}\text{erfi}[g_{-}({\it\xi})]-\text{e}^{k{\it\xi}+2\text{i}({\it\sigma}/S)}\text{erfi}[g_{+}({\it\xi})]\right\}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle g_{\pm }({\it\xi})=\left(1+\text{i}\right)\frac{k({\it\xi}rS\pm \text{i})+\text{i}r{\it\sigma}}{2\sqrt{krS}}, & & \displaystyle\end{eqnarray}$$
and
$\text{erfi}(z)=-\text{i}\,\text{erf}(\text{i}z)$
denotes the imaginary error function (an entire function). Note that
$r$
is retained as an independent parameter in the function
$\tilde{u}$
in order to represent both solutions
$\tilde{u} _{-}$
(when setting
$r=1$
) and
$\tilde{u} _{+}$
. The constants
$c_{1\pm }$
,
$c_{2\pm }$
and
$c_{3\pm }$
at both sides of the flame must be determined by requiring that perturbations vanish as
$|{\it\xi}|\rightarrow \infty$
. From the properties of
$\text{erfi}(z)$
for large
$z$
, we have
$$\begin{eqnarray}\displaystyle \text{erfi}(z)\sim \text{i}\;\text{sgn}(\text{Re}(z))+\frac{\text{e}^{z^{2}}}{\sqrt{{\rm\pi}}z}\quad \text{for }z\rightarrow \infty . & & \displaystyle\end{eqnarray}$$
This implies that
$c_{1-}=c_{3-}=0$
and
$c_{2+}=\text{i}\text{e}^{2\text{i}{\it\sigma}/S}c_{3+}$
. The complete solution is cumbersome to write but it can be evaluated at
${\it\xi}=0$
easily and combined with (3.7)–(3.9) to form a linear system of equations for the free parameters
$c_{2-}$
,
$c_{1+}$
and
$c_{3+}$
. In order for this system to have a non-trivial solution, its determinant must be zero, resulting in the dispersion relation
$$\begin{eqnarray}\displaystyle & & \displaystyle -\sqrt{ks}\left[r^{2}(k(\mathscr{M}_{c}-\mathscr{M}_{s}{\it\phi})+2{\it\phi}+1)+r(k(\mathscr{M}_{c}+\mathscr{M}_{s}({\it\phi}+2))+\text{i}\mathscr{M}_{s}s-1)+\text{i}\mathscr{M}_{s}s\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\,\left(k^{2}r(\mathscr{M}_{c}(r{\it\phi}+1)+\mathscr{M}_{s}({\it\phi}+1))+k\left(r^{2}{\it\phi}({\it\phi}+1)-\text{i}\mathscr{M}_{s}s(r{\it\phi}-1)\right)+2\text{i}rs{\it\phi}\right)\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{\sqrt{{\rm\pi}}(r-1)}{(1-\text{i})}\text{e}^{(\text{i}k(r{\it\phi}+1)^{2}/2s)}\left[1-\text{erf}\left(\frac{\left(1+\text{i}\right)\sqrt{k}(r{\it\phi}+1)}{2\sqrt{s}}\right)\right]=0,\end{eqnarray}$$
where
${\it\phi}={\it\sigma}/k$
denotes the phase speed,
$s=rS$
denotes a reduced shear parameter, and
$\mathscr{M}_{c}=\mathscr{M}-\mathscr{M}_{s}$
. As is the case for DL, (4.10) has two complex roots (owing to the quadratic power in
${\it\phi}$
) and we are only interested in the root with the largest real part, which determines the stability properties of the flame. The solution of (4.10) is obtained by a classical Newton method. First, one can easily verify that for
$k=0$
we have
${\it\phi}=0$
. But, this is not always the best starting point of the Newton method if we need to obtain the two roots of the solution, to ensure we are selecting the most unstable mode. To select an appropriate initial point in the iterative method, the most reliable approach is to observe that for small
$k$
, the exponential and the error function in (4.10) are both very nearly one and zero, respectively. Then, one can solve the quadratic equation for
${\it\phi}$
(valid for small
$k$
to order
$\sqrt{k}$
) and obtain expressions for the two approximate roots. The first root behaves as
${\it\phi}\sim \sqrt{k}$
and the second as
${\it\phi}\sim 1/k$
. From these initial conditions it is simple to trace the complete curves describing the two solutions of
${\it\phi}$
as a function of
$k$
. It turns out that it is always the first root that has the larger real part.
The current dispersion relationship can be compared with the DL relationship without shear but including curvature and strain, see Creta & Matalon (Reference Creta and Matalon2011b ) (they considered the case of a single Markstein number), which has solution
$$\begin{eqnarray}\displaystyle {\it\phi}\! & = & \displaystyle \!\{\!k(\mathscr{M}_{c}r+\mathscr{M}_{s})+r- ([k(\mathscr{M}_{c}r+\mathscr{M}_{s})+r]^{2}-r[k\mathscr{M}_{s}(1-r)+r+1]\nonumber\\ \displaystyle & & \displaystyle \quad \times \,[k(2\mathscr{M}_{c}-\mathscr{M}_{s}(r-3))+r-1])^{1/2}\!\}/\{r[k\mathscr{M}_{s}(r-1)-(r+1)]\}.\end{eqnarray}$$
Several distinguished limits will be discussed next to elucidate the effects of the different parameters in the solution. Then, we will consider a few prototypical situations of combustion of hydrogen, methane and propane with air, with Markstein numbers suggested by Clavin & Graña-Otero (Reference Clavin and Graña-Otero2011).
4.1. Limit of vanishing Markstein numbers
First, we consider the limit of zero Markstein numbers which is analogous to the original DL instability result for a planar flame without shear. The dispersion relationship in this case is obtained by setting
$\mathscr{M}_{s}=\mathscr{M}_{c}=0$
in (4.10), giving
$$\begin{eqnarray}\displaystyle & & \displaystyle -\sqrt{ks}\left[r^{2}(2{\it\phi}+1)-r\right]+\left[k(r^{2}{\it\phi}({\it\phi}+1))+2\text{i}rs{\it\phi}\right]\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{\sqrt{{\rm\pi}}(r-1)}{(1-\text{i})}\text{e}^{(\text{i}k(r{\it\phi}+1)^{2}/2s)}\left[1-\text{erf}\left(\frac{\left(1+\text{i}\right)\sqrt{k}(r{\it\phi}+1)}{2\sqrt{s}}\right)\right]=0.\end{eqnarray}$$
Figure 2 shows the behaviour of
$\text{Re}({\it\sigma})$
as a function of
$s$
, as well as the original DL curve (
$s=0$
limit), for
$r=1/7$
. As can be seen, the main role of shear is to stabilize the flame at low wavenumbers (large wavelengths). The cross-over wavenumber,
$k_{0}$
, at which the flame becomes unstable depends on
$s$
and
$r$
and it can be obtained in principle by setting
${\it\phi}=\text{i}{\it\phi}_{i}$
in (4.12) and solving for
$k$
. Figure 3 shows the behaviour of
$k_{0}$
for different values of
$s$
and
$r$
. It is observed that the dependence on
$s$
is apparently linear while
$r$
changes the slope of these lines.
Classical Darrieus–Landau stability curves with shear (zero Markstein numbers) for
$r=1/7$
.
Dependence of
$k_{0}$
on
$s$
for different
$r$
.
4.2. Limit of vanishing density jump
Another interesting limit is obtained by setting
$r=1-{\it\delta}$
where
$0<{\it\delta}\ll 1$
, analogously to the thermodiffusive approximation. The main item that must be considered carefully is the behaviour of
$\mathscr{M}_{s}$
and
$\mathscr{M}_{c}$
in this limit. It can be easily confirmed from Clavin & Graña-Otero (Reference Clavin and Graña-Otero2011) that both Markstein numbers remain finite as
${\it\delta}\rightarrow 0$
. Therefore, to leading order in
${\it\delta}$
, (4.10) becomes
$$\begin{eqnarray}\displaystyle & & \displaystyle -\sqrt{ks}\left[(k(\mathscr{M}_{c}-\mathscr{M}_{s}{\it\phi})+2{\it\phi}+1)+(k(\mathscr{M}_{c}+\mathscr{M}_{s}({\it\phi}+2))+\text{i}\mathscr{M}_{s}s-1)+\text{i}\mathscr{M}_{s}s\right]\nonumber\\ \displaystyle & & \displaystyle \quad +\, (\!k^{2}(\mathscr{M}_{c}({\it\phi}+1)+\mathscr{M}_{s}({\it\phi}+1))+k\left({\it\phi}({\it\phi}+1)\right.\nonumber\\ \displaystyle & & \displaystyle \quad -\,\left.\text{i}\mathscr{M}_{s}s({\it\phi}-1)\right)+2\text{i}s{\it\phi}\!)=0,\end{eqnarray}$$
which can be solved explicitly for
${\it\phi}$
, giving
$$\begin{eqnarray}\displaystyle {\it\phi}=-(k(\mathscr{M}_{s}+\mathscr{M}_{c})+\text{i}s\mathscr{M}_{s})=-(k\mathscr{M}+\text{i}s\mathscr{M}_{s}). & & \displaystyle\end{eqnarray}$$
This result is consistent with the same limit for the DL analysis without shear, indicating that the flame is stable for positive
$\mathscr{M}$
and unstable for negative
$\mathscr{M}$
, and that it is density variation that makes the flame unstable for positive
$\mathscr{M}$
. Shear affects
${\it\sigma}$
only through the imaginary term (phase speed) that plays no role in stability of the flame.
4.3. Limit of strong density jump
An interesting limit arises when the density jump across the flame is very large, as might be found in thermonuclear deflagrations. In this case,
${\it\rho}_{+}$
can be much smaller than
${\it\rho}_{-}$
, resulting in the limiting case
$0<r\ll 1$
. Note that it is possible to take the limit
$r\rightarrow 0$
in the dispersion relation (which gives a finite result) but it is not completely physical for the flow (it is a singular limit) because mass conservation will result in an infinite velocity in the burned side, see (3.1b
).
The dispersion relation for this limit can be obtained by manipulating (4.10), which after some algebraic manipulations, gives
$$\begin{eqnarray}\displaystyle {\it\sigma}=\frac{1-2k\mathscr{M}_{c}-3k\mathscr{M}_{s}-2\text{i}\mathscr{M}_{s}S}{2\mathscr{M}_{s}}. & & \displaystyle\end{eqnarray}$$
The key technical detail in deriving the above equation is to observe that for large
$z$
we have
$$\begin{eqnarray}\displaystyle \text{erf}(z)\approx 1-\frac{\text{e}^{-z^{2}}}{\sqrt{{\rm\pi}}z}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle z=\frac{(1+\text{i})\sqrt{k}(r{\it\phi}+1)}{2\sqrt{rS}}, & & \displaystyle\end{eqnarray}$$
with
$r$
approaching zero. Then the terms in the last line of (4.10) become
$$\begin{eqnarray}\displaystyle \frac{\sqrt{rS}(r-1)}{\sqrt{k}(1+r{\it\phi})}, & & \displaystyle\end{eqnarray}$$
which transforms the dispersion relation into a quadratic equation that can be solved explicitly for
${\it\phi}$
. This leads, after expansion in powers of
$r$
, to (4.15). One can also verify easily that by setting
$S=0$
in (4.15) the amplification factor agrees with the corresponding limit of (4.11) when
$r\rightarrow 0$
. The main observation from this (now simple) equation (4.15) is that shear does not affect the most unstable mode since it contributes a purely imaginary term and the amplification factor is identical to that of the DL instability without shear.
Real (a,c) and imaginary (b,d) values of
${\it\sigma}$
with the larger real part as a function of reduced wavenumber for
$\mathscr{M}=1$
,
$r=1/7$
and various
$s$
.
4.4. Zero strain case
This case corresponds to the original correction of Markstein to the DL analysis, where only curvature is accounted for in the problem,
$\mathscr{M}_{s}=0$
and
$\mathscr{M}_{c}=\mathscr{M}$
. The interest of this situation is that the dispersion relation can be rescaled by
$\mathscr{M}$
, by choosing
$k\rightarrow \mathscr{M}k$
and
$s\rightarrow rS\mathscr{M}$
, and the result is therefore generic and equivalent to
$\mathscr{M}=1$
. We focus mostly on the case
$r=1/7$
to discuss the qualitative properties of the solution. The results are similar for other sufficiently small values of
$r$
; the dependence of the solution on
$r$
is investigated in more general terms in figure 6 below. Figure 4 shows real and imaginary parts of the root of (4.10) with the larger real part for various reduced shear values and
$r=1/7$
, as well as the DL amplification factor without shear, (4.11).
As is well known, the first term in (4.11) dominates at small wavenumbers while the second (stabilizing) term dominates at high wavenumbers. When comparing the effect of shear, one can verify that
${\it\sigma}$
approaches the DL result as
$s\rightarrow 0$
(as shown in figure 4
a), which confirms that the analysis was carried out correctly. First, concentrating on the real part of
${\it\sigma}$
and as the value of
$s$
increases, the curves move to the right and immediately introduce a gap of stability (negative
$\text{Re}({\it\sigma})$
) at the smallest wavenumber, as shown in figure 2. At the same time, the peaks of the curves decrease and eventually they become completely negative at
$s\approx 0.17$
. Beyond that value of
$s$
, the flame is stable at all wavenumbers. As
$s$
continues to increase, the curve bends up and first becomes tangent to the
${\it\sigma}=0$
axis at
$s\approx 2.4$
. Beyond this value of
$s$
, we recover a bell-shaped curve with a positive peak at a location dependent on the value of
$s$
and a further decay towards negative values for growing wavenumber (a consequence of the stabilizing role of curvature), see figure 4(c). Second, figure 4(b) shows that the magnitude of the imaginary part of
${\it\sigma}$
is negative (being zero for DL in the unstable range of wavenumbers) for low
$s$
while it becomes positive for large
$s$
. The consequence of
${\it\sigma}$
not being purely real is that the solution will have a travelling wave character, see (4.1).
To determine the speed of the travelling wave character of the solution, we investigate the phase speed
$c=\text{Im}({\it\phi})$
in figure 5. It can be seen that the phase speed is negative for small
$s$
and positive for large
$s$
. The latter is consistent with the direction of the imposed shear, but the former may appear counterintuitive, being in the opposite direction. In fact, we can see in the distinguished limit
$r\rightarrow 1$
, (4.14), that the expected result is to have a negative phase speed (assuming
$S>0$
) for
$\mathscr{M}_{s}>0$
. Here, one may note the fact that the results in figure 5 refer to
$\mathscr{M}_{s}=0$
and therefore (4.14) predicts zero phase speed in this case. To address this difficulty, (5.1) below shows the result for a Kelvin–Helmholtz interface between two fluids of different densities using the same reference conditions as in the present flame. This would be the limit of high shear, where the self-propagation feature of the premixed flame is negligible in comparison with the effect of shear. Since
$r<1$
, the phase speed is indeed negative for the Kelvin–Helmholtz mode and this is irrespective of the values of Markstein numbers. It is, therefore, the positive phase speed behaviour of the flame in the presence of large shear that is unexpected, and it must be the result of the nonlinearities in the dispersion relation that couple with the finite heat release. In other words, the phase speed is determined by the sign of
$r-1$
, as in (5.1), but in a manner controlled by
$s$
.
Phase speed
$c$
as a function of wavenumber for
$\mathscr{M}=1$
,
$r=1/7$
and various
$s$
. The phase speed is only plotted in the regions where
$\text{Re}({\it\sigma})>0$
.
Most unstable wavenumber as a function of
$s$
for different
$r$
and
$\mathscr{M}=1$
.
Figure 6 shows the wavenumber of the most unstable mode,
$k^{\ast }$
(at the peak of
$\text{Re}({\it\sigma})$
), as a function of the reduced shear parameter
$s$
for different values of
$r$
relevant to typical combustion:
$1/7\leqslant r\leqslant 1/3$
(other lower values of
$r$
give similar results). Figure 6 shows that in all cases one observes that the ‘island’ of stability, whose precise boundaries depend on the value of
$r$
, is present. This island appears to exist for all practical values of
$r$
and therefore should be relevant in applications when premixed flames propagate in regions of shear since the flame may experience different stabilizing–destabilizing conditions depending on the local value of
$s$
.
Vorticity for
$r=1/7$
at the most unstable wavenumber for
$s=0$
(a) (DL with curvature effect only),
$s=0.1$
(b) and
$s=5$
(c).
To conclude, figure 7 shows isocontours of vorticity for DL, with curvature stabilization only, and
$s=0.1$
and 5 in the post-flame region (downstream); vorticity is zero upstream for all values of
$s$
. The vorticity can be calculated explicitly and results in the following formula (valid for
$\mathscr{M}=1$
):
$$\begin{eqnarray}\displaystyle \tilde{{\it\omega}} & = & \displaystyle \left(\frac{\partial \tilde{v}}{\partial {\it\xi}}+\text{i}k\tilde{u} \right)\text{e}^{-\text{i}k{\it\eta}}\nonumber\\ \displaystyle & = & \displaystyle -\frac{2(1-\text{i})\sqrt{sk}}{\sqrt{{\rm\pi}}}\text{exp}\left[-k\left(({\it\xi}+\text{i}{\it\eta})+\text{i}((1-r{\it\phi}+\text{i}s{\it\xi})^{2}/2s)\right)\right],\end{eqnarray}$$
which excludes the temporal modulation since the complete vorticity is travelling with phase speed
$c$
. Figure 7 shows the real part of (4.19). The main change induced by shear is the tilting of the isocontours in the direction of the imposed shear (upwards in these figures). The picture for low
$s$
is similar to that at high
$s$
but there is a larger density of ripples (smaller wavelengths) in the flame-normal direction while the transverse wavelength depends on
$s$
, as can be seen in figure 6. Furthermore, the complete vorticity solution will be traveling at the phase speed
$c$
(downwards or upwards depending on
$s$
) and being amplified exponentially in time (since we are only concerned with the unstable mode). Note that the vorticity decays to zero at large
${\it\xi}$
but at a rate that is slower than that of the DL solution.
Amplification factors for a prototypical hydrogen–air flame (a,b), lean methane–air flame (c,d) and lean propane–air flame (e,f).
Amplification factors for a flame with
$\mathscr{M}_{s}=2$
,
$\mathscr{M}_{c}=-5$
and
$r=1/6$
.
4.5. Behaviour in prototypical combustion conditions
We consider now the dispersion relation solution as a function of
$s$
for several Markstein numbers discussed in Clavin & Graña-Otero (Reference Clavin and Graña-Otero2011) as representative of combustion of low-Lewis-number mixtures corresponding to hydrogen–air and close to unity-Lewis-number mixtures characteristic of hydrocarbon combustion.
Figure 8 shows dispersion relation results for a prototypical rich hydrogen–air flame with
$\mathscr{M}_{s}=3$
and
$\mathscr{M}_{c}=-2.5$
(a,b), a prototypical lean methane flame with
$\mathscr{M}_{s}=3.5$
and
$\mathscr{M}_{c}=1.5$
(c,d) and a lean propane flame with
$\mathscr{M}_{s}=0.5$
and
$\mathscr{M}_{c}=1.5$
(e,f), all at
$r=1/6$
. In all cases, we reproduce qualitatively the result in the previous section. Starting at
$s=0$
and with increasing
$s$
, the flame becomes less unstable and eventually it stabilizes completely at some small value of
$s$
. But, if
$s$
continues to increase, the flame becomes unstable again within a range of wavenumbers, that grows with increasing
$s$
, as seen in all the figures.
4.6. Negative Markstein numbers
Finally, we investigate the effect of shear on flames with negative Markstein numbers. Traditionally, this regime is not often discussed because the flame is absolutely unstable (at all wave modes). Also note that for
$s=0$
and negative strain, Markstein number from equation (4.11) is singular at the finite value of
$k=(r+1)/((r-1)\mathscr{M}_{s})$
. This is presumably an intrinsic limitation of the description based on a strictly zero flame thickness when
$s=0$
.
In the present configuration, the stabilization introduced by shear implies that negative Markstein number flames can become conditionally unstable. This is clearly observed in figure 9 where the amplification factor for a flame with
$\mathscr{M}_{s}=2$
and
$\mathscr{M}_{c}=-5$
is shown. This flame has negative Markstein number for the curvature term and positive Markstein number for the strain term in the flame speed relationship. We observe a stable region for small
$s$
at low wavenumbers, as before, but a conditionally stable situation at higher
$s$
. This latter regime is analogous to the behaviour observed for positive Markstein numbers before because the stabilization effect of strain is still present. Note that the corresponding amplification factor for
$s=0$
is unconditionally unstable for this flame. Figure 10 shows the result for a flame with both Markstein numbers being negative,
$\mathscr{M}_{s}=-2$
and
$\mathscr{M}_{c}=-5$
, where it is seen that the flame is conditionally unstable. The amplification factor for small shear is very similar to that of the
$s=0$
case (small shear does not stabilize the flame). We also do not find a singularity of the amplification factor for negative Markstein numbers when
$s$
is larger than approximately
$0.32$
(for this flame) and the flame seems to be conditionally unstable always.
5. Discussion: the second region of instability
As discussed in the preceding sections, all flames are conditionally unstable for sufficiently high shear, even with negative Markstein numbers. The reasons this happens is unclear. One might conjecture that in this limit the flame propagation speed is negligible and the interface is therefore Kelvin–Helmholtz unstable with a density jump given by heat release. This situation would correspond to that of a sheared interface of two immiscible fluids of different densities. Now, a straightforward stability analysis, see Drazin & Reid (Reference Drazin and Reid1981), gives an amplification factor of
$$\begin{eqnarray}\displaystyle {\it\sigma}=\frac{2\text{i}(r-1)S}{(r+1)}, & & \displaystyle\end{eqnarray}$$
which corresponds to a neutrally stable interface. Therefore, it is not a Kelvin–Helmholtz instability that dominates at high shear, but a truly intrinsic interaction of the flame dynamics with shear.
Amplification factors for a flame with
$\mathscr{M}_{s}=-2$
,
$\mathscr{M}_{c}=-5$
and
$r=1/6$
.
6. Conclusions
A linear stability analysis of the behaviour of a laminar premixed flame subjected to a constant lateral shear has been carried out analytically. The interest of this problem resides in applications where a premixed flame can propagate into regions of free shear flow, which can alter the flame behaviour. The analysis extends the classical Darrieus–Landau instability analysis to include the stabilization effect produced by flame curvature and strain. The dispersion relation is obtained and its roots are discussed in detail. The stability analysis indicates that the flame show three regimes as a function of a reduced shear parameter. At low shear, the behaviour is analogous to Darrieus–Landau analysis, albeit with some correction, the flame becoming progressively more stable as shear increases. Further increase in shear results in unconditionally stable behaviour. Finally, as shear continues to increase, the flame becomes unstable again with higher amplification rate and at lower wavelengths than the low shear case. All solutions are stable at sufficiently high wavenumbers due to curvature and strain stabilization. This island of stability is observed at all heat release levels considered.
Acknowledgements
This work was supported in part by the Department of Energy, National Nuclear Security Administration, under award no. DE-NA0002382 and the California Institute of Technology. We also thank M. Matalon for the comments on the first draft of the manuscript.
