Stability of diffusion flames under shear flow: Taylor dispersion and the formation of flame streets

Diffusion flame streets, observed in non-premixed micro-combustion devices, align parallel to a shear flow. They are observed to occur in mixtures with high Lewis number ($Le$) fuels, provided that the flow Reynolds number, or the Peclet number $Pe$, exceeds a critical value. The underlying mechanisms behind these observations have not yet been fully understood. In the present paper, we identify the coupling between diffusive-thermal instabilities and Taylor dispersion as a mechanism which is able to explain the experimental observations above. The explanation is largely based on the fact that Taylor dispersion enhances all diffusion processes in the flow direction, leading effectively to anisotropic diffusion with an effective (flow-dependent) Lewis number in the flow direction which is proportional to $1/Le$ for $Pe\gg 1$. Validation of the identified mechanism is demonstrated within a simple model by investigating the stability of a planar diffusion flame established parallel to a plane Poiseuille flow in a narrow channel. A linear stability analysis, leading to an eigenvalue problem solved numerically, shows that cellular (or finite wavelength) instabilities emerge for high Lewis number fuels when the Peclet number exceeds a critical value. Furthermore, for Peclet numbers below this critical value, longwave instabilities with or without time oscillations are obtained. Stability regime diagrams are presented for illustrative cases in a $Le$-$Pe$ plane where various instability domains are identified. Finally, the linear analysis is supported and complemented by time dependent numerical simulations, describing the evolution of unstable diffusion flames. The simulations demonstrate the existence of stable cellular structures and show that the longwave instabilities are conducive to flame extinction.


Introduction
Formation of diffusion flame streets in non-premixed micro-combustion channels, first reported in [1][2][3], has been shown to exhibit peculiar features that are not observed in large scale burners.First of all, they are found to occur readily when using heavier fuels but not lighter fuels such as hydrogen.In the latter case, a continuous diffusion flame is observed [1,3] unless the inlet mixture is sufficiently diluted with helium which is known to increase the effective Lewis number.In summary, flame streets are found to occur in nonpremixed combustion when the effective Lewis number is sufficiently large [3].Moreover, they are observed only if the flow Reynolds number (or, the flow rate) exceeds a critical value [1][2][3][4][5]; sufficiently above this critical value the streets are found to be stable.When the Reynolds number is lower than the critical Reynolds number, typically the quenched state is observed.However, depending on the conditions prevailing at the channel exit, sometimes repetitive extinction-ignition events starting from the exit region and involving edge flame propagation are observed.Computational studies investigating these flame streets [6][7][8][9] have qualitatively reproduced the experimental observations.The critical Reynolds number can depend on a number of relevant controlling parameters such as the mixture strength, heat loss mechanism, type of the channel cross-section, etc; for example, the critical Reynolds number is found to decrease with increasing wall temperature [5].Regardless, the important point to note is that the existence of a critical Reynolds number provides strong evidence for a flow-induced effect causing the formation of a flame street.However, the underlying physical mechanisms behind the experimental observations have not been yet fully identified.In particular, despite numerous efforts, the following two issues still lack clear explanations: 1. Why planar diffusion flames aligned with the shear flow do not exist for heavier fuels when effective Lewis number of the inlet mixtures is large?2. What justifies the existence of a critical Reynolds number (or, equivalently a critical Peclet number) above which stable diffusion flame streets are encountered?The suggestion in [2,3] that the observed patterns are the result of an instability appears not to have been pursued yet.In the present contribution, we shall show that the two aforementioned issues may be explained as being caused by well-known diffusive-thermal flame instabilities, provided account is made of Taylor's dispersion mechanism [10].
Although Taylor dispersion played a fundamental role in a wide range of practical applications involving unidirectional flows, its importance in combustion systems has not been recognized until recently.Its influence on the structure and stability of premixed flames are investigated in [11][12][13][14][15] and the structure of Burke-Schumann diffusion flame is studied in [16,17] .An important result revealed by the analyses of [12,16] is the presence of a flow-dependent effective Lewis number in the flow direction, say the x-direction, given by Here Le is the fuel Lewis number, Pe the Peclet number and γ a constant which is determined by the velocity profile.The significance of this formula is clearly elucidated by evaluating it for the following three values of Pe: Thus, the flow-dependent effective Lewis number, Le x , varies from being purely determined by molecular diffusion and equal to Le in the absence of flow to being equal to 1/Le when Pe ≫ 1.In other words, the effect of a shear flow is such that a weakly diffusing scalar appears effectively as a strongly diffusing scalar and a strongly diffusing scalar appears effectively as a weakly diffusing one if Pe > 1/ √ γLe.This peculiar feature of Taylor dispersion is expected to produce interesting effects on flame characteristics.For instance, when the diffusion flame lies perpendicular to the fluid stream [16], subadiabatic flame temperatures are predicted for low Lewis number fuels and superadiabatic flame temperatures for large Lewis number fuels if Pe > 1/ √ γLe.Although formula (1) is not strictly valid when thermal expansion effects are taken into account, the two limiting forms in (2) for Pe = 0 and Pe ≫ 1 are still correct as shown in [17].Finally, we can note from formula (1) that Le x = 1 for all values of Pe if Le = 1.It is important to emphasize that the enhancement of diffusion is in the longitudinal (or, flow direction) only and not in the transverse direction.This leads to anisotropic diffusion characterised by Le x in the x-direction and by Le in the transverse direction.It follows that the flame alignment with respect to the flow is an essential aspect to consider when attempting to explain the observations in nonpremixed micro-combustion devices.

Oxidizer
In the present paper, we shall show that the issues mentioned above can be explained as being the consequence of diffusional-thermal instabilities occurring in a medium where diffusion is effectively anisotropic due to Taylor dispersion.This is demonstrated by investigating the stability of a planar diffusion flame aligned with the shear flow accounting for the presence of effectively two distinct Lewis numbers, Le x and Le, in the longitudinal and transverse directions respectively.It is important to emphasize that this problem, which addresses the stability of the diffusion flame under the influence of shear-enhanced diffusion, is both novel and fundamental and therefore deserves investigation on its own right.To this end, and in order to highlight the new instability mechanism involved, we shall purposely adopt suitable simplifying assumptions.These include ignoring complications related to entrance and exit conditions by focusing on the fully-developed flow region as well as ignoring other effects such as the presence of heat losses and density variations.The model configuration is sketched in Fig. 1 which depicts a diffusion flame aligned parallel to a two-dimensional Poiseuille flow.The instabilities of the underlying one-dimensional diffusion flames have been studied previously, in the absence of longitudinal flows, in [7,[18][19][20][21][22].
The paper is structured as follows.The problem is formulated in §2.The linear stability of the planar diffusion flame is characterized by an eigenvalue problem derived in §3.The solution of the eigenvalue problem is presented in §4, where the emerging instabilities are discussed and delineated in a Le-Pe plane.Numerical simulations complementing the linear stability findings are presented in §5 for selected cases, illustrating the nonlinear time evolution of unstable flames.

Governing equations and boundary conditions
We consider a simple constant-density model to investigate the stability of a planar flame aligned with the direction of a shear flow, as sketched in Fig. 1.Thermal expansion and heat loss effects are ignored so as to demonstrate that account of Taylor dispersion is sufficient to obtain the instabilities under consideration.The chemistry is modeled by a single step irreversible Arrhenius reaction for which the mass of fuel burnt per unit volume per unit time is given by Here B is the preexponential factor, ρ is the constant density, Y * i is the mass fraction of species i, T * is the temperature, E is the activation energy and R is the universal gas constant.
The temperatures at both sides y = ±L are assumed equal and denoted by T u .It is also assumed that the concentrations of fuel and oxidizer are maintained fixed at y = ±L, namely, at the oxidizer side y = +L.Although such conditions may be delicate to fully achieve experimentally [23,24], they are adopted herein for sake of simplicity, as it is done in many other theoretical investigations involving diffusion flames such as [7,[19][20][21][22][25][26][27][28][29][30].
The Zeldovich number β, the heat release parameter α, the stoichiometry parameter S and the Damköhler number Da are defined by where T ad = T u + qY F,F /c p (S + 1) is the adiabatic flame temperature at the stoichiometric location.Further, s and q denote the mass of oxygen consumed and the amount of heat released per unit mass of fuel burnt, c p the constant-pressure specific heat and D T the constant thermal diffusivity.The oxidizer Lewis number Le O is taken equal to unity and the fuel Lewis number Le F is simply denoted as Le.
The problem is analyzed in a frame of reference moving with the flow mean velocity U.In this frame, the velocity components are given by In terms of the non-dimensional variables the governing equations read where ǫ = h/L, Pe = Uh/D T and In the limit ǫ → 0, the dependent variables tend to be uniform in the z direction and the convective terms on the left-side of (4)-( 6) lead upon depth averaging of the equations to enhanced longitudinal diffusion as derived in [11,12,16,17].The resulting problem is two dimensional and is given by where p = √ γPe and γ = 2/105.The corresponding lateral boundary conditions are As for the boundary condition in the x-direction, we shall adopt, for simplicity, periodic conditions, which will allow us in particular to avoid complications associated with entrance and exit conditions.This is sufficient for our main purpose, which is to identify the coupling between shear-enhanced diffusion and the instability of the planar diffusion flame, the focus being on the region where the flow is fully developed.For investigations that incorporate the entrance region conditions and other effects such as heat loss, the reader is referred to [6], where a two-dimensional depth-averaged model is derived and solved numerically, and to [7][8][9], which are based on three-dimensional computations.
Because Le O is taken to be unity, a simplified formulation of the above system of equations can be obtained [31] by introducing the mixture fraction Z and the excess enthalpy variable H, which are normalized such that Z − 1 = H = 0 at y = −1 and Z = H = 0 at y = +1.A linear combination of ( 8)-( 10) yields the chemistry free equation which upon integration leads to Consequently the problem can be described in terms of two dependent variables only, which are selected here to be H and Y F .These are governed by where the scaled variable X = x/ 1 + p 2 is introduced for convenience and where ω(H, Y F ) is obtained by substituting the equations which follow from ( 13) and ( 14) into the expression of ω given by (7).It is worth noting in ( 15)-( 16) the presence of the longitudinal Lewis number Le x defined in (1).When Le = 1, we have Le x = 1 and H = 0 and one needs to solve only for Y F .

The base solution
The base solution whose stability is under investigation corresponds to the steady planar diffusion flame aligned with the flow direction.Since this solution, denoted by an overbar, depends only on the y coordinate, it is governed by ( 15)-( 16) with ∂/∂t = ∂/∂X = 0, that is by equations whereas the second equation is integrated numerically.Specifically, the nonlinear differential equation is solved numerically using COMSOL Muliphysics software, which determines Da as an unknown parameter along with the solution for values of the temperature prescribed at a given location, sat at y = 0.A sample of the numerical results is given in Fig. 2. Shown are plots of the maximum of the temperature field, θ max , and its location, y flame , versus the Damköhler number Da for selected values of Le.In the limit Da → ∞, these quantities are seen to approach their Burke-Schumann values The solid curves correspond to near-equilibrium diffusion flames and the dashed curves to flames in the so-called partial burning regime [32] which are typically unstable.For each Le, the turning point Da = Da ext corresponds to extinction.Since Da ext varies significantly with Le as evident from the figure, it is convenient for later reference to use as the relative measure of the Damköhler number.

Perturbed state
To investigate the stability of the planar diffusion flame, we consider infinitesimal normalmode perturbations superimposed on the base solution such that where σ is the complex growth rate and κ is the real wave number.Similarly, one may define the perturbations Y O (y) and θ(y), which according to ( 17)-( 18) are related to Y F (y) and H(y) by Substituting the perturbed dependent variables in ( 15)-( 16) and linearizing, we obtain the eigenvalue problem 1 Le where .
The quantities H and Y F satisfy homogeneous Dirichlet boundary conditions, namely, This eigenvalue problem possesses a discrete spectrum of eigenvalues σ n indexed by an integer n.For example, when δ → ∞, the eigenvalues are given by a simple formula, derived in Appendix A. In general, the eigenvalue with the largest real part (real growth rate), max n [Re{σ n }], determines the stability characteristics.Henceforth, σ shall refer to this eigenvalue and not to the entire spectrum.

Preliminary remarks and notations
The main object of our stability analysis is to describe the relationship between σ (the eigenvalue with the largest real part), the wavenumber κ and the parameters Le, p and δ.All eigenvalues are computed by solving the eigenvalue problem ( 20)-( 22) numerically with a tolerance of 10 −6 using the Chebfun package [33], which is based on a spectral collocation method using Chebyshev polynomials.
In doing so, computations are restricted to the strongly burning, near-equilibrium diffusion flames, the top branches in the left plot of Fig. 2. As the scaled Damköhler number δ defined in (19) tends to infinity, the corresponding Burke-Schumann flame is known to be unconditionally stable when p = 0 [19].This conclusion remains true for p = 0 as shown in Appendix A. As δ is decreased the flame becomes unstable for values of δ smaller than a critical value δ c , determined by the marginal stability condition Re{σ} = 0. Planar adiabatic near-equilibrium diffusion flames are thus unstable in the range 0 < δ < δ c .
The flame stability for p = 0 has been investigated in [19] and [20] through asymptotic analyses in the so-called slowly varying flame and near-equidiffusional flame limits [34], respectively.In particular, it has been shown in [20,30] that diffusion flames are unstable for Le > Le c , where Le c > 1 is a critical value depending on δ.The ensuing instability was found to be a longwave instability with or without pulsations in time.
In our problem allowing for non-zero values of p, we also find that no instability occurs for Le < Le c where Le c = Le c (δ, p) > 1.However, unlike in the p = 0 case, the instability need not always be a longwave instability, i.e., an instability for which the most unstable mode corresponds to κ = 0. Indeed, we also find for a certain parametric range, finite wavelength instability for which the most unstable mode corresponds to κ = κ m = 0, as we shall confirm below.
To facilitate the discussion, it is convenient to introduce the notation where the values σ 0 and σ m correspond to the condition dσ/dκ = 01 .In fact, we may regard σ 0 as the growth rate (possibly complex) characterising the development of the longwave instability.In this case, the instability will be appropriately referred to as a non-oscillatory longwave instability if σ 0 is real and as an oscillatory longwave instability otherwise.Similarly, σ m may be regarded as the growth rate (found always to be real) characterising the development of the finite wavelength (cellular) instability.

Results for cases with δ varying and fixed values of p and Le
The effect of varying δ on flame stability for the top-branches in the left subfigure of Fig. 2 is illustrated in Fig. 3 pertaining to Le = 2 and p = 0.5 (left subfigure) and p = 3 (right subfigure).Plotted is the real part of σ versus the wave number κ for selected values of δ.
The corresponding curves (dispersion curves) are drawn with solid lines when σ(κ) is real and dashed lines when σ(κ) has a non-zero imaginary part.The left subfigure corresponding to p = 0.5 is found to be qualitatively similar to the p = 0 case (not shown); here the most unstable mode always corresponds to κ = 0, indicating the presence of a longwave instability.This longwave instability is non-oscillatory when 0 < δ < δ * and oscillatory when δ * < δ < δ c , where δ * is some critical value (δ * ≈ 0.0168 when Le = 2 irrespective of the value of p; see below).
We turn now to the right subfigure which corresponds to p = 3.Here the notable feature is that the most unstable mode corresponds to a non-zero value κ m of κ.Therefore we are in the presence of a finite wavelength (cellular) instability.Furthermore, it is worth pointing out that the dispersion curve is associated with non-oscillatory modes characterized by real growth rates for 0 < δ < δ * (solid lines in the subfigure), but that unstable oscillatory modes also exist when δ * < δ < δ c (dashed lines).
Before closing this section, we note that the parameter δ * introduced above corresponds to σ 0 defined in (24) changing, as δ is increased, from being a real number to being a complex number with nonzero imaginary part.Since σ 0 is independent of p (see footnote 1), δ * , a function of Le, is also independent of p.This function is plotted in Fig. 4 and the curve thus computed characterizes the transition from oscillatory to non-oscillatory longwave instabilities.

Results for cases with p varying and fixed values of δ and Le
In this section, we investigate the effects of varying p for a fixed value of δ.Variations in p can in fact be achieved conceptually by adjusting the flow rate as done in experiments [2,3].An illustrative case is shown in Figure 5 corresponding to δ = 0.01 and Le = 1.6 (left subfigure) and Le = 2 (right subfigure).Plotted is the real part of σ versus the wave number κ for selected values of p.As in figure 3, the sections of the dispersion curves drawn with solid lines correspond to σ(κ) being real and dashed sections to σ(κ) having a non-zero imaginary part.It can be seen that σ 0 ≡ σ(κ = 0) (as defined in (24)) is independent of p, in line with the remark made in footnote 1.
Focusing first on the left subfigure, pertaining to Le = 1.6, we note that the dispersion curve for p = 0 is entirely complex (plotted as a dashed line) and as its maximum occurs at κ = 0, it indicates the presence of an oscillatory longwave instability.On the other hand, for p = 1, the dispersion curve consists of a complex branch (dashed line) which meets a real branch (solid line) at a finite value of κ where a cusp is encountered.Although for this value of p, the real branch is still in the stable region, it gradually crawls up on the complex branch as p is increased, thus narrowing the range of oscillatory modes.The maximum of the real branch (κ m , σ m ) enters the unstable region as p exceeds a certain value of p corresponding to σ m = 0.When p is increased further, a critical value p = p A is encountered for which σ m = Re{σ 0 }.For p > p A , σ m > Re{σ 0 } and therefore the most unstable mode occurs at κ = κ m = 0.This indicates the presence of a finite wavelength (cellular) instability.The transition between the oscillatory longwave instability (occurring for p < p A ) and the finite wavelength instability (occurring for p > p A ) may be viewed as a jump transition, in the sense that a cellular pattern is predicted to suddenly emerge as soon as p exceeds p A .We turn now to the right subfigure to examine how the findings above are affected by adopting a higher value of the Lewis number, Le = 2. Here, the dispersion curve for p = 0 implies a non-oscillatory longwave instability, unlike the corresponding curve in the left subfigure where an oscillatory longwave instability occurs when p = 0.As can be seen, an increase in p leads the concavity of the dispersion curve at κ = 0 to change from being negative to positive as p crosses a critical value p = p B , corresponding to the condition d 2 σ/dκ 2 | κ=0 = 0.For p > p B , the dispersion curve has a maximum at κ m = 0, implying the presence of a finite wavelength instability.For p < p B , the dispersion curve peaks at κ = 0 and no maximum occurs at a finite values of κ and the instability is a non-oscillatory longwave one.Note however that the transition at p = p B between the non-oscillatory longwave instability and the finite wavelength instability is continuous in the sense that κ m → 0 as p decreases towards p B .

Stability regime diagram
Based on our discussion above of the various types of instabilities and transitions, we can synthesize the findings in the form of a regime diagram in the Le-p plane.This is computed and shown in Fig. 6 in the illustrative case δ = 0.01.Solid lines in this diagram separate regions possessing different stability characteristics, whereas the dash-dotted curve represents the equation p = 1/ √ Le which is equivalent, according to (2), to the condition Le x = 1.We note that in the domain above this dash-dotted curve, we have Le x < 1 if Le > 1 and Le x > 1 if Le < 1.Four regions can be identified in the diagram: a stable region on the left (white region) and three unstable regions on the right (shaded regions), corresponding to the three types of instabilities encountered.This section corresponds to a bifurcation associated with a non-oscillatory (or stationary) finite wavelength instability, often referred to as a type-I s instability [35].
We now briefly comment on the boundaries between the three unstable regions aforementioned, which are seen to intersect at a cusp point (Le, p) = (1.85,0.76).First, the vertical line at Le ≈ 1.85 for p ∈ [0, 0.76] is obtained from Fig. 4 using the condition δ * = 0.01.The curve, labeled p A , corresponding to the critical condition p = p A introduced in the previous subsection, separates the regions of oscillatory longwave instability and finite wavelength instability.Similarly, the curve, labeled p B , corresponds to the critical condition p = p B and separates the regions of non-oscillatory longwave instability and finite wavelength instability.An important observation related to this figure is that the finite wavelength instability region lies entirely in the region Le x < 1 (above the dash-dotted curve).This demonstrates that a necessary condition for the occurrence of the finite wavelength (or cellular) instability is that Le x < 1 provided the fuel Lewis number Le > 1.
It is also interesting to examine how a change in δ affects the stability regime diagram of Fig. 6, by carrying out similar computations for δ = 10 −4 .The results are summarized in Fig. 7.The figure shows that the four regions identified earlier persist in this case, except that the unstable regions are extended now further to the left (Le c ≈ 1.17), thereby narrowing the stable region.Moreover, the extent of the oscillatory longwave instability region is found now to be narrower in comparison with that of Fig. 6.
To close this section, we note that our results have been computed for the stoichiometrically balanced case, S = 1.However, it is possible to provide qualitative comments relevant to other values of S, since the controlling Lewis number pertains to the deficient reactant which corresponds to the fuel when S ≪ 1 and to the oxidizer when S ≫ 1.Furthermore, since the influence of Taylor dispersion on the flame instability is stronger when the Lewis number is further from unity, this influence is more apparent when S ≪ 1 than when S ≫ 1 because Le O ≈ 1.In particular, we can infer that cellular instabilities are easier to realize in mixtures with Le F = Le > 1 and S ≪ 1 and more difficult in mixtures with Le F = Le < 1 and S ≪ 1, a fact which is consistent with experiments [3].

Time dependent numerical simulations
The stability regime diagram of §4.4,based on the computations of the eigenvalues of the linear stability problem (20)- (22), summarizes the main conclusions that one could infer from the linear analysis.To describe the eventual fate of infinitesimal disturbances, we need to solve the full nonlinear problem (15)-( 16).This will allow, in particular, to answer the important question about the existence of cellular structures and their stability.The solution of equations ( 15)-( 16) is carried out numerically and selected cases in each of the three unstable regions of diagram 6 are presented herein.The computations are performed using COMSOL Multiphysics, a package which has been extensively tested in combustion applications, see e.g.[11][12][13][14][15].The computational domain, taken to be (−15, 15) × (−1, 1), is covered by a grid consisting of approximately 400,000 triangular elements and includes local refinement in the reaction zone.Solutions have been tested to be independent of the mesh and of the time step ∆t (whose maximum value is set to be ∆t = 0.005).Dirichlet boundary conditions corresponding to (11) and (12) are applied in the y-direction.As far as the X-direction, we adopt periodic boundary conditions.These boundary conditions are suitable for our idealized simple model, given that our main concern is to study the stability of a planar diffusion flame aligned with the direction of the flow.In particular, we disregard complicating effects associated with heat loss and with special inlet/outlet conditions in finite domains which are encountered in practice.The two types of longwave instability, namely oscillatory and non-oscillatory, identified through our linear stability analysis, are always found in our computations to ultimately result in flame quenching2 .The only difference between these two types is found to lie in their transient evolution to extinction.This is illustrated in Fig. 8 for two cases corresponding to Le = 1.6 (left subfigure) and Le = 2 (right subfigure) with selected values of p. Shown in this figure is the maximum temperature θ max versus time t.We note that the solid curves correspond to two-dimensional computations, while the dashed curves, pertaining to κ = 0, are computed as solutions of the one-dimensional problem obtained from ( 15)-( 16) by dropping the X-derivatives.The two-dimensional computations are carried out with the horizontal domain size equal to 30, as mentioned above, and therefore the minimum wavenumber allowed numerically is κ = 2π/30 = 0.2094.We note that once θ max (t) drops below the steady-state temperature θ max corresponding to the static extinction condition δ = 0 (the turning points in Fig. 2), θ max cannot recover and decays monotonically to zero.

Finite wavelength instability
In this section, we present computations for parameter values which fall in the region labelled 'finite wavelength instability' in Fig. 6.The computations starting from the initial condition corresponding to a planar diffusion flame are performed in a large number of cases, but we shall in our presentation focus on four illustrative cases.Two of these cases correspond to the points (Le = 1.6, p = 2) and (Le = 1.6, p = 3) in Fig. 6; the first point is selected to be close to the boundary p = p A in Fig. 6 and the second further from this boundary.The two other cases correspond to the points (Le = 2, p = 1.5) and (Le = 2, p = 2) in Fig. 6; the first point is selected to be close to the boundary p = p B in Fig. 6 and the second further from this boundary.The time evolution of the flame in these four cases are described in figures 9-12, showing the fields of the reaction rate ω for selected values of time.Note that in these figures, the coordinate x (which measures the horizontal distance scaled by the mixing layer thickness L) is used instead of X = x/ 1 + p 2 in order to better appreciate the size of the emerging cells.corresponds to (Le = 2, p = 2), shown in Fig. 9.Here we observe that the planar diffusion flame initially splits into a series of strongly burning regions, separated by locally quenched gaps.The strongly burning regions, each of which terminates at two edge-flames, are found to readily adjust the gaps among themselves and settle into an apparently stable steady state.At the onset of instability occurring at t ≈ 3.9, we note the emergence of 8 cells or strongly burning regions.This number of cells, which is observed to persist in time, is in good agreement with that predicted by the linear stability analysis of §4.Indeed the maximum growth rate according to the linear analysis correspond to κ = κ m = 1.65, which for our horizontal domain x = 30 1 + p 2 ≈ 67 yields 7.86 ≈ 8 cells.Similar observations can be made for unstable flames when p ≫ p A and p ≫ p B as our extensive set of computations (not shown herein) confirm.In particular, it is observed that unstable diffusion flames evolve into apparently stable steady cells and the number of these cells is rather correctly predicted by the linear stability analysis.
We now comment on the other three cases under consideration where a more complex dynamics is obtained.Figure 10 illustrates the flame evolution when (Le = 1.6, p = 3).Here after the initial splitting of the planar flame, the strongly burning regions continue to shrink to form flame spots at t ≈ 15.8.Some of these spots eventually quench, leaving a wide gap of fresh unburnt reactant mixture (t ≈ 16.7).This wide gap is invaded by propagating edgeflames from opposite ends and gradually shrinks to disappear at t ≈ 21.5.The resulting nearly planar flame, replacing the wide gap, becomes itself unstable and splits into a series of cells, which are found to settle into an apparently stable stationary cellular structure.
Turning now to Fig. 11 pertaining to the case (Le = 1.6, p = 2), we observe that the development of the flame instability is roughly similar in its initial and final stages to that in the previous figure.In particular, an apparently stable steady structure is ultimately obtained, as in Fig. 10, despite peculiar dynamics in intermediate stages.
We now discuss the case (Le = 2, p = 1.5) which is represented in Fig. 12.Here the most notable new feature is that the unstable diffusion flame does not evolve into a steady state.It evolves instead into an irregular, apparently chaotic state, where several physical mechanisms are at play including edge flame propagation and merging, flame splitting, etc.The full dynamics of the flame is best appreciated by examining the animated time evolution which is included as a supplementary material.Note that the chaotic behaviour observed seems to occur when p is close to p B , corresponding to the continuous transitions discussed in §4.3.
Finally, to complement our discussion of figures 9-12, we plot the function θ max (t) for the cases of these figures in Fig. 13.The figure confirms the asymptotic evolution of the unstable diffusion flame into a steady state, except for the case (Le = 2, p = 1.5) of Fig. 12, where an irregular time dependent behaviour is obtained.

Additional comments on unstable flames when p p A
Consideration of the left subfigure of Fig. 5 shows that the planar diffusion flame is expected to disintegrate near instability onset into a cellular structure if p > p A , as confirmed in our two-dimensional computations just presented.In fact, such cellular structures once formed can persist if p is decreased to a certain extent below p A , which may have some implications on the formation of diffusion flame streets in experiments.This statement can be confirmed by numerical computations (not shown) adopting as initial conditions cellular structures or an array of hot spots, as done in [37,38].For example, we have performed computations for the case Le = 1.6 and p = 1.5, for which the most unstable mode corresponds to the oscillatory longwave mode (κ = 0), according to Fig. 5 (left), although the finite wavelength mode is also unstable in this case.It is found that if the initial condition is taken to be a cellular structure (say, the steady solution corresponding to Le = 1.6 and p = 3), then this structure evolves into another stable cellular structure.In contrast, if the initial condition is taken to be the planar diffusion flame solution, then an oscillatory longwave instability leading to flame quenching is observed.

Concluding remarks
In this paper, we have examined the stability of a planar diffusion flame aligned with the direction of a shear flow within a simple narrow channel model.Particular attention has been focused on the effect of the flow Peclet number on flame stability.A linear stability analysis, supported by time dependent numerical simulations, confirms that the diffusion flame can be unstable for Le sufficiently larger than one and that the nature of the emerging instability depends critically on the (scaled) Peclet number p. Specifically, it is shown that a longwave instability with or without time oscillations is obtained for small values of p, while a finite wavelength (cellular) instability is obtained for p above a critical value p c .The instability domains are clearly delimited for the illustrative examples of Fig. 6 and 7, where p c corresponds to the curves labelled p A (jump transition) and p B (continuous transition).
The longwave instability is found to typically lead to flame quenching.The more crucial finding is however that a critical value p c of the Peclet number is shown to exist, above which the instability is cellular; this is consistent with the experimental observations made in the original studies [2,3].Our findings provide an original explanation to the issues related to the two questions raised in the introduction, namely, the existence of a critical Peclet number and the disintegration of the planar flame when Le > 1.Our explanation is largely based on Taylor's dispersion mechanism which leads to the presence of effectively two distinct Lewis numbers, Le x and Le, in the longitudinal and transverse directions respectively, as discussed in the introduction.The necessary condition for the occurrence of the finite wavelength (or cellular) instability is found to be Le x < 1 provided the fuel Lewis number Le > 1.
To close this section, we make a few remarks regarding possible extensions of the current study.The first useful extension is to account for heat losses which were deliberately neglected herein in order to isolate the effect of Taylor dispersion as being on its own a driving mechanism of the cellular flame instability.Heat loss effects can be significant, however, in experiments such as those involving methane diffusion flames [2,3].Such diffusion flames, characterised by a (fuel) Lewis number close to unity, are expected to be stable according to our adiabatic investigation, but are found to disintegrate into flame streets in practice, which may be explained by accounting for heat losses [6].Another useful extension of this study is to account for the influences of stoichiometry by considering values of S = 1 since flammability limits can be greatly extended due to the formation of stable cellular structures in the range δ < 0 when S ≪ 1 [37,38].It would also be interesting to have experimental verification of some key findings of the current investigation such as the existence of cellular structures near the jump transition at p A identified above, both for p p A and p > p A .Finally, it is worth investigating whether the cellular instability for large Lewis numbers induced by Taylor dispersion identified herein for diffusion flames can also be found for premixed flames aligned with the direction of a shear flow.Such investigations would enrich our appreciation of premixed flame instabilities in confined systems considered in recent studies [39][40][41][42].

Figure 1 :
Figure 1: Model configuration showing a diffusion flame aligned parallel to a shear (Poiseuille) flow.The fuel and oxidizer reach the flame by transverse diffusion.

Figure 2 :
Figure 2: Maximum temperature and flame location of the unperturbed planar diffusion flame as functions of Da for selected values of Le; here and elsewhere α = 0.85, β = 10 and S = 1.

Figure 5 :
Figure 5: Dispersion curves representing the maximum growth rate, Re{σ}, versus the wave number κ for selected values of p.The left subfigure summarizes calculations for Le = 1.6 and δ = 0.01 and the right subfigure for Le = 2 and δ = 0.01.All computations are carried out with β = 10, α = 0.85 and S = 1.Solid lines represent real branches (Re{σ} = σ) and dashed lines complex branches (Im{σ} = 0).The values of critical Peclet numbers are p A ≈ 1.56 and p B ≈ 1.00.

Figure 8 :
Figure 8: Maximum temperature θ max as a function of time t, calculated for Le = 1.6 (left) and Le = 2 (right) with δ = 0.01.The dashed lines pertain to one-dimensional computations (independent of p), whereas the solid lines pertain to two-dimensional computations for selected values of p.All computations are carried out with β = 10, α = 0.85 and S = 1.

Figure 13 :
Figure 13: Maximum temperature θ max as a function of time t, calculated for Le = 1.6 (left) and Le = 2 (right) with δ = 0.01 and selected values of p.All computations are carried out with β = 10, α = 0.85 and S = 1.The cases correspond to those of figures 9-12.All curves asymptote to a steady state except when Le = 2, p = 1.5.