<figure><img src="image_1.png" /><figcaption>Fig. 1. Flame propagation against a unidirectional flow in a channel with width 2H.</figcaption></figure>

out a numerical study addressing the influence of Taylor dispersion on flame stability which appears to be the first in the literature. A particular effort will be devoted to reconciling apparent discrepancies between the theoretical findings, obtained in the asymptotic limit $\beta \to \infty$, and the numerical findings, obtained for finite values of the Zeldovich number $\beta$. This will highlight the importance of assessing the influence of $\beta$ on flame stability, as can be surmised from previous studies addressing this issue in the absence of Taylor dispersion [15–19].

In this investigation, heat losses will be modelled for simplicity as being volumetric, although they are mainly due in practice to non-adiabatic walls. Approximating such losses to walls by a volumetric sink term in the energy equation is justifiable, from the theoretical point of view, for narrow channels or thick flames [20], provided they are sufficiently weak, so that the walls are essentially near adiabatic. A detailed analytical and numerical study on the propagation of non-adiabatic flames in channels, addressing in particular this specific point, can be found in Daou and Matalon [21]. Although the weak heat loss assumption is admittedly difficult to satisfy from the experimental point of view as the walls thermal control poses non-trivial challenges, it is adopted herein as it is useful to advance the theory.

The paper is structured as follows. The problem formulation is given in the next section, where a two-dimensional model incorporating Taylor dispersion is presented which results in enhanced diffusion in the longitudinal direction and hence anisotropic diffusion. This is followed by a linear stability analysis of the planar flame propagating in the longitudinal direction, carried out analytically in the asymptotic limit $\beta \to \infty$ and culminating in the derivation of a dispersion relation. The implications of the dispersion relation are then investigated to identify the combined effect of Taylor dispersion and heat loss on the flame cellular and oscillatory instabilities. The paper closes with numerical simulations addressing both types of instabilities, including a comparison between theory and numerics with particular emphasis on clarifying the role played by the Zeldovich number.

## 2. Formulation

### 2.1. Thermo-diffusive model with Taylor dispersion

In order to have a simple model accounting for Taylor dispersion effects on flames, we consider the case of a channel Poiseuille flow, although other shear flows such as a Couette flow could also be used. The configuration of study is represented in Fig. 1 depicting a flame propagating in a channel of width $2H$ against a Poiseuille flow of amplitude $\hat{A}$. In a Cartesian frame of reference $(\hat{x}, \hat{y}, \hat{z})$ attached to the walls, the problem can be represented by the following depth-averaged 2D model

$$ \frac{\partial \bar{T}}{\partial \hat{t}} + \frac{2\hat{A}}{3} \frac{\partial \bar{T}}{\partial \hat{x}} = D_T (1 + \gamma \text{Pe}^2) \frac{\partial^2 \bar{T}}{\partial \hat{x}^2} + D_T \frac{\partial^2 \bar{T}}{\partial \hat{y}^2} + \frac{q}{c_p} B \bar{Y}_F e^{-E/R\bar{T}} - K(\bar{T} - T_u) \quad (1a) $$

$$ \frac{\partial \bar{Y}_F}{\partial \hat{t}} + \frac{2\hat{A}}{3} \frac{\partial \bar{Y}_F}{\partial \hat{x}} = D_F (1 + \gamma \text{Pe}^2 \text{Le}^2) \frac{\partial^2 \bar{Y}_F}{\partial \hat{x}^2} + D_F \frac{\partial^2 \bar{Y}_F}{\partial \hat{y}^2} - B \bar{Y}_F e^{-E/R\bar{T}} \quad (1b) $$

$$ \bar{T} = T_u \quad \bar{Y}_F = Y_{Fu} \quad \text{as} \quad \hat{x} \to -\infty \quad (1c) $$

$$ \frac{\partial \bar{T}}{\partial \hat{x}} = 0 \quad \bar{Y}_F = 0 \quad \text{as} \quad \hat{x} \to +\infty \quad (1d) $$

The model has been derived in Daou [6], except that it is slightly modified here to incorporate the effect of weak heat losses. These are accounted for in a simple and common way through the addition of the term $-K(\bar{T} - T_u)$ in the temperature equation where $K$ is a heat loss coefficient.

In this model, a chemical reaction following an Arrhenius law with pre-exponential factor $B$ and activation energy $E$, and heat release $q$ is adopted. $R$ is the universal gas constant, $Y_F$ the fuel mass fraction (assumed deficient), $T$ the temperature and the subscript $u$ indicates values in the unburnt mixture (as $\hat{x} \to -\infty$). The parameters Pe and Le are the Peclet and Lewis numbers given by

$$ \text{Pe} = \frac{\hat{A}H}{D_T} \quad \text{and} \quad \text{Le} = \frac{D_T}{D_F} $$

respectively, where $D_T$ is the thermal diffusivity and $D_F$ the fuel diffusion coefficient. The bar indicates quantities depth-averaged across the channel width such that

$$ \bar{T} = \bar{T}(\hat{x}, \hat{y}, \hat{t}) = \frac{1}{2H} \int_{-H}^{H} T d\hat{z} $$

Note that the term $2\hat{A}/3$ is the average of the Poiseuille flow $\hat{u} = \hat{A}(1 - \hat{z}^2/H^2)$ and that the diffusion coefficients $D_T$ and $D_F$ have enhanced effective values in the longitudinal direction $\hat{x}$. This enhancement is in agreement with Taylor-Aris dispersion formula [2] and involves the parameter $\gamma$ which is a numerical coefficient determined by the velocity profile; $\gamma = 8/945$ for the channel Poiseuille flow considered here. This result and a rigorous derivation of the model used herein can be found in appendix 1 of [6]. The derivation is obtained in the double limit $\epsilon \to 0$ with $\epsilon \text{Pe} \to 0$ where $\epsilon = H/\delta_L$ is the channel width scaled with the laminar flame thickness $\delta_L$. This double limit is more general than, and comprises the distinguished limit $\epsilon \to 0$ with $\text{Pe} = O(1)$ used in Pearce and Daou [4], Daou et al. [5].

Since the enhancement of diffusion is in the longitudinal $\hat{x}$-direction only, and not in the $\hat{y}$-direction, a notable feature of the problem is the presence of anisotropic diffusion.

Our aim is to investigate the stability of the planar flame solutions, independent of $\hat{y}$, of problem (1a)-(1d). That is, we need to revisit the thermo-diffusive instabilities of flames accounting for the consequences of Taylor dispersion and the corresponding anisotropic diffusion, as well as heat losses.
## 2.2. Non-dimensional equations in a frame attached to the unperturbed planar flame front

In a frame of reference attached to the unperturbed (planar) flame front, the problem (1a)-(1d) takes the non-dimensional form

$$ \frac{\partial \theta}{\partial t} + U \frac{\partial \theta}{\partial x} = (1 + p^2) \frac{\partial^2 \theta}{\partial x^2} + \frac{\partial^2 \theta}{\partial y^2} + \omega - \frac{\kappa \theta}{\beta} \quad (2a) $$

$$ \frac{\partial y_F}{\partial t} + U \frac{\partial y_F}{\partial x} = \frac{1}{\mathrm{Le}} (1 + p^2 \mathrm{Le}^2) \frac{\partial^2 y_F}{\partial x^2} + \frac{1}{\mathrm{Le}} \frac{\partial^2 y_F}{\partial y^2} - \omega \quad (2b) $$

$$ \theta = 0 \quad y_F = 1 \quad \text{as} \quad x \to -\infty \quad (2c) $$

$$ \theta_x = 0 \quad y_F = 0 \quad \text{as} \quad x \to +\infty \quad (2d) $$

in terms of

$$ y_F = \frac{\bar{Y}_F}{Y_{Fu}} \qquad \theta = \frac{\bar{T} - T_u}{T_{ad} - T_u} $$

Here, $U$ is the non-dimensional propagation speed of the unperturbed front with respect to the gas. Furthermore, $T_{ad} = T_u + qY_{Fu}/c_p$ is the adiabatic flame temperature ($c_p$ being the mixture's heat capacity) and $\beta = E(T_{ad} - T_u)/RT_{ad}^2$ is the Zeldovich number. For non-dimensionalisation, we have chosen $S_L$ as unit speed, $\delta_L$ as unit length, and $\delta_L/S_L$ as unit time where

$$ S_L = \sqrt{\frac{2}{\beta^2} \mathrm{Le} D_T \mathrm{Be}^{-E/R T_{ad}}} \quad \text{and} \quad \delta_L = \frac{D_T}{S_L} $$

are the laminar flame speed (for $\beta \gg 1$) and the planar flame thickness in the absence of Taylor dispersion and heat losses ($p=0$ and $\kappa=0$). Throughout the paper, $p^2 = \gamma \mathrm{Pe}^2$ represents the Taylor dispersion coefficient and $\kappa = \beta KD_T/S_L^2$ the non-dimensional heat loss coefficient. Finally, $\omega$ is the reaction rate given by

$$ \omega = \frac{\beta^2}{2\mathrm{Le}} y_F \exp \left( \frac{\beta (\theta - 1)}{1 + \alpha_h (\theta - 1)} \right) $$

where $\alpha_h = (T_{ad} - T_u)/T_{ad}$.

The analysis is most consistently carried out analytically in the limit $\beta \to \infty$ using the so-called near-equidiffusional flame (NEF) approximation based on the assumption that the Lewis number deviates little from unity [22, p. 33]. Within this approximation, the reduced Lewis number $l = \beta(\mathrm{Le} - 1)$ is $O(1)$ and Eqs. (2a) and (2b) can be written in terms of the leading order temperature $\theta^0 \sim \theta$ and $h \sim \beta(\theta + y_F - 1)$. The reformulated problem is given by the equations

$$ \theta_t^0 + U \theta_x^0 = (1 + p^2) \theta_{xx}^0 + \theta_{yy}^0 \quad (3) $$

$$ h_t + U h_x = (1 + p^2) h_{xx} + h_{yy} + l \left[ (1 - p^2) \theta_{xx}^0 + \theta_{yy}^0 \right] - \kappa \theta^0 \quad (4) $$

which are applicable outside an infinitely thin reaction sheet, given by $x = f(y, t)$ say, subject to the boundary conditions

$$ \theta^0 = 0 \quad h = 0 \quad \text{as} \quad x \to -\infty \quad (5) $$

$$ \theta^0 = 1 \quad h \text{ has no exp. growth as } x \to +\infty \quad (6) $$

and the jump conditions

$$ \llbracket \theta^0 \rrbracket = 0 \qquad \llbracket h \rrbracket = 0 \quad (7a) $$

$$ \llbracket h_x \rrbracket + \frac{1 + f_y^2 - p^2}{1 + f_y^2 + p^2} l \llbracket \theta_x^0 \rrbracket = 0 \quad (7b) $$

$$ (1 + f_y^2 + p^2)^{\frac{1}{2}} \llbracket \theta_x^0 \rrbracket = - \exp \left( \frac{h}{2} \right) \quad (7c) $$

at $x = f(y, t)$. Here we have used the notation $\llbracket \phi \rrbracket = \phi(x = f^+) - \phi(x = f^-)$.

The jump conditions (7), accounting for anisotropic diffusion, are fully derived in Daou [6] following a methodology commonly used in premixed combustion (see e.g., Buckmaster and Ludford [22, p. 39]). They reduce to the well known jump conditions of premixed flames [23, p. 527] when $p^2 = 0$, that is when Taylor dispersion is absent.

# 3. Linear stability analysis

## 3.1. The planar flame solution

The planar flame solution (denoted by a tilde) whose stability is being investigated is governed by Eqs. (3)-(7) with $\partial/\partial t = 0$ and $\partial/\partial y = 0$. The solution is given by

$$ \tilde{f} = 0 \quad (8a) $$

$$ \tilde{\theta} = \begin{cases} e^{\frac{U_x}{1+p^2}} & (x < 0) \\ 1 & (x > 0) \end{cases} \quad (8b) $$

$$ \tilde{h} = \begin{cases} \left[ \frac{-2\kappa(1+p^2)}{U^2} + \left( \frac{\kappa}{U} - \frac{l(1-p^2)U}{(1+p^2)^2} \right) \right] x & (x < 0) \\ \left[ \frac{-2\kappa(1+p^2)}{U^2} - \frac{\kappa x}{U} \right] & (x > 0) \end{cases} \quad (8c) $$

where the propagation speed $U$ is determined by

$$ \frac{U^2}{1 + p^2} \ln \left( \frac{U^2}{1 + p^2} \right) = -2\kappa \quad (8d) $$

The last equation can be written as

$$ U_*^2 \ln U_* = -\kappa \quad (9) $$

where

$$ U_* \equiv \frac{U}{(1 + p^2)^{\frac{1}{2}}} \quad (10) $$

represents the planar flame propagation speed $U = U(\kappa, p)$ scaled by its adiabatic value $(1 + p^2)^{1/2}$. A plot of $U_*$ versus $\kappa$ is shown in Fig. 2 where a classical inverse-C shaped curve is obtained exhibiting an extinction point labelled C with coordinates

$$ \kappa^{\text{ext}} \equiv \frac{1}{2e} \quad U_*^{\text{ext}} \equiv \frac{U^{\text{ext}}}{\sqrt{1 + p^2}} = e^{-\frac{1}{2}} \quad (11) $$

## 3.2. Normal modes analysis

A normal-mode stability analysis can now be applied to the basic solution (8) by considering perturbations of the form

$$ \begin{bmatrix} f \\ \theta^0 \\ h \end{bmatrix} = \begin{bmatrix} \tilde{f} \\ \tilde{\theta}(x) \\ \tilde{h}(x) \end{bmatrix} + \delta e^{st+iky} \begin{bmatrix} 1 \\ \hat{\theta}(x) \\ \hat{h}(x) \end{bmatrix} \quad (12) $$

where $\delta$ is a small number representing the amplitude of the perturbation, and $k$ and $s$ are a real and a complex numbers representing its wavelength and its growth rate, respectively. The main aim is to derive a dispersion relation expressing the dependence of $s$ on $k$ and the parameters $l, \kappa$ and $p^2$. Instability will correspond of course to situations for which $Re(s) > 0$ occurs.
<figure><img src="image_2.png" /><figcaption>Fig. 2. Scaled propagation speed $U_*$ versus $\kappa$. The stability of the planar solutions corresponding to the points labelled A, B, C and D will be examined in Fig. 3 based on the dispersion relation (13). All points on the lower branch have been confirmed to be unstable, as expected.</figcaption></figure>

Note that the speed $U$ appearing in the governing Eqs. (3)-(7) and the basic solution (8) is a known quantity, determined by Eq. (8d) in terms of $\kappa$ and $p$, and corresponds to the propagation speed of the unperturbed planar front. The methodology used to derive the dispersion relation is classical, see e.g., [24,25], and the derivation is given in Appendix A. Briefly, when perturbed solutions of the form (12) are substituted into Eqs. (3)-(7), an eigen-boundary value problem for the functions $\hat{\theta}(x)$ and $\hat{h}(x)$ is obtained in the linear approximation corresponding to small values of $\delta$. The dispersion relation is found to be given by

$$2\Gamma^2(\Gamma - 1) + 2(\Gamma^2 - 1)\ln U_* + \frac{l}{1 + p^2}[(\Gamma - 1 - 2s_*)(1 - p^2) + 4p^2k_*^2] = 0 \quad (13)$$

where

$$\Gamma = (1 + 4(s_* + k_*^2))^\frac{1}{2}$$

Here the starred variables $s_*$ and $k_*$ are rescaled versions of $s$ and $k$ defined by

$$s_* \equiv \frac{s}{U_*^2}, \quad k_* \equiv \frac{k}{U_*} \quad (14)$$

where $U_*$ is the rescaled propagation speed, introduced in (10), which is a bivalued function of $\kappa$ given by (9) and represented in Fig. 2.

Once a value of $U_*$ has been selected, corresponding to a point on either the upper or lower branch of Fig. 2 and hence to a specific value of $\kappa$, the stability of the corresponding solution can be investigated for arbitrary values of $p$ by turning to Eq. (13). In summary, the dispersion relation (13) encapsulates all the information characterising the stability of the flame to small perturbations for arbitrary values of $l$, $\kappa$ and $p$. Before exploring its implications, we note that (13) reduces as it should to the dispersion relation derived by Sivashinsky [24] when $p=0$ and $\kappa=0$, to that derived by Joulin and Clavin [25] when $p=0$, and to that derived by Daou [6] when $\kappa=0$.

# 4. Implications of the theoretical results

## 4.1. Stability regions

We note that the dispersion relation (13) incorporates the effect of Taylor dispersion and the resulting anisotropy of diffusion on the stability of non-adiabatic flames.

The overall picture provided by the dispersion relation is illustrated in Fig. 3 where bifurcation curves are plotted in the $l-k_*$ plane for selected values of the parameters $p$ and $\kappa$ (or rather $U_*$ which is a bivalued function of $\kappa$). The first row in Fig. 3 corresponds to point A in Fig. 2, the second row to point B, the third row to the extinction point C, and the fourth row to point D on the lower branch of Fig. 2.

When $p=0$ (first column) and $\kappa=0$ ($U_*=1$, first row) we recover the classical stability diagram pertaining to an adiabatic flame. Here we have a stationary bifurcation for $l<-2$ (leading to a cellular instability) and a Hopf bifurcation for $l>32/3$ (oscillatory instability).

The effect of $\kappa$ on the thermo-diffusive instability as studied by Joulin and Clavin [25] corresponds to the first column (where the dispersion coefficient $p=0$). It is seen, in agreement with the findings of [25], that a decrease in the value of $U_*$ destabilises the flame, promoting both the cellular and oscillatory instabilities. In particular, the last row confirms that point D on the lower branch is unstable, and so is in fact any other point on this lower branch (and this is found to hold for any value of $p$).

The first row summarises the effect of $p$ in the adiabatic case; as $p$ is increased the Hopf bifurcation curve is displaced to the right and disappears for $p \ge 1$. This indicates that enhanced diffusion strongly impedes the oscillatory instability. Also, the size of the stability domain in the left half-plane (delimited from the left by the stationary bifurcation curve) is reduced by an increase in $p$ which indicates that Taylor dispersion may somewhat promote the cellular instability in that it destabilises some wave numbers which are otherwise stable. This suggests that flames with finite or restricted extent in the $y$-direction (so that $k$ must be larger than a minimum non-zero value) may be stable when $p=0$ but unstable for non-zero values of $p$. This is confirmed by numerical simulations below; see Fig. 7 where such possibility is illustrated in two cases corresponding to $\kappa=0$ and $\kappa \ne 0$. In general, Fig. 3 demonstrates the ability of Taylor dispersion and heat loss to combine so as to significantly affect the stability of the flame. For example, while the oscillatory instability is promoted by an increase in $\kappa$ (decrease in $U_*$), it is hampered by an increase in $p$. This is clearly illustrated in Fig. 4, where the stability regions are determined in the $\kappa-l$ plane for $p=0$ (left) and $p=0.5$ (right).

## 4.2. The longwave cellular instability

We note that the dispersion relation (13) has always a real root $s_*(k_*)$ such that $s(0)=0$ and that the Taylor expansion of $s_*(k_*)$ for small values of $k_*^2$ is given by

$$s_* = -\frac{l+2+4\ln U_*}{2+4\ln U_*} k_*^2 - \frac{l^2(6-l+(6+l)p^2)}{(1+p^2)(2+4\ln U_*)^3} k_*^4 + O(k_*^6)$$

provided $U_* > U_*^{\text{ext}} \equiv e^{-\frac{1}{2}}$ (see (11)). It is this root that is at the origin of the cellular instability, a longwave instability whose onset corresponds to the coefficient of $k_*^2$ in the expansion becoming zero, which corresponds to $l=l_c$ where

$$l_c = -2 - 4 \ln U_* \quad (15)$$

The critical value $l_c$ given by (15) coincides with that derived in Joulin and Clavin [25] when $p=0$, but it is found here to be applicable irrespective of the value of $p$, provided the dependence on $p$ is incorporated in $U_*$ as expressed in (10).
<figure><img src="image_3.png" /><figcaption>Fig. 3. Bifurcation curves in the $l-k$ plane for selected values of parameters $p$ and $U_*$, The parameter $p$ increases through the values 0, 0.5 and 1 from left to right. $U_*$ decreases from top to bottom, taking the values corresponding to points A, B, C and D in Fig. 2; that is $U_* = 1$ (corresponding to $\kappa = 0$), $U_* = 0.7$ ($\kappa \approx 0.174$), $U_* = e^{-1/2}$ ($\kappa = \kappa^{\text{ext}} \approx 0.184$) and $U_* = 0.5$ ($\kappa \approx 0.173$).</figcaption></figure>

In the weakly unstable regime near the onset of instability, $l \lesssim l_c$, the expansion above can be written as

$$s_* = \frac{l-l_c}{l_c} k_*^2 - \left( \frac{1-p^2}{1+p^2} - \frac{6}{l_c} \right) k_*^4 + \dots \quad (16)$$

For weakly unstable flames, $l \lesssim l_c$, the unstable modes are those with wavenumbers $k = U_*k_*$ such that

$$0 < k_*^2 < -\frac{(1+p^2)(l-l_c)}{(6-l_c) + (6+l_c)p^2}$$

and the most amplified mode is characterized by a wavenumber $k_0 = U_*k_{0*}$ where $k_{0*}^2$ is the midpoint of the interval above, that is

$$\frac{k_0^2}{U_*^2} = -\frac{1}{2} \frac{(1+p^2)(l-l_c)}{(6-l_c) + (6+l_c)p^2}$$

Therefore, using (15), $k_0$ and the corresponding growth rate $s_0$ are given by

$$k_0^2 = -\frac{U_*^2 (1+p^2)(l+2+4\ln U_*)}{8 (2+\ln U_* + (1-\ln U_*)p^2)} \quad (17)$$

$$s_0 = \frac{U_*^2}{32} \frac{(1+p^2)(l+2+4\ln U_*)^2}{(1+2\ln U_*)(2+\ln U_*+(1-\ln U_*)p^2)} \quad (18)$$
<figure><img src="image_4.png" /><figcaption>Fig. 4. Bifurcation curves and stability regions in the $\kappa - l$ plane for $p = 0$ (left) and $p = 0.5$ (right). The left subfigure is in agreement with Fig. 2 of Joulin and Clavin [25].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 5. Growth rate versus the wavenumber $k$ for selected values of $p$ in the adiabatic case $\kappa = 0$ ($U_* = 1$) with the reduced Lewis number $l = -2.2$ (left subfigure) and $l = -3$ (right subfigure). The solid curves are based on the asymptotic formula (16) and the dashed curves are obtained numerically from the dispersion relation (13).</figcaption></figure>

These relations imply that an increase in $p^2$ widens the range of unstable wave numbers and raises the maximum growth rate and has therefore a destabilising effect. These conclusions are confirmed in Fig. 5, where the growth rate is plotted versus $k$ for selected values of $p$ in the illustrative near-onset case $l = -2.2$ (left subfigure); the solid curves are based on the asymptotic formula (16) and the dashed curves are obtained numerically from the dispersion relation (13). The right subfigure of Fig. 5 is similar except that it pertains to the case $l = -3$ which is further from onset. In this case, it is seen that the solid curves deviate significantly from the dashed curves from the quantitative point of view, although the deviations are perhaps insignificant from the qualitative point of view. Still, when comparing with the finite-$\beta$ numerical simulations to be presented below, we shall use growth rates based on the dispersion relation (13) rather than formula (16), and on the dashed, rather than solid, curves.

We note that (17) can be used to provide a length scale $\lambda_0 = 2\pi/k_0$ which may plausibly characterize the size of cells appearing in the weakly unstable regime $l \lesssim l_c$. This length scale is determined here in terms of $p$ and $\kappa$ (or $U_*$) and is given by

$$ \lambda_0 = \frac{4\pi \sqrt{2}}{U_*} \sqrt{\frac{2 + \ln U_* + (1 - \ln U_*)p^2}{(1 + p^2)|l + 2 + 4 \ln U_*|}} \quad (19) $$

The formula indicates that when $p$ increases the cells size predicted by $\lambda_0$ should decrease, with the ratio $\lambda_0(p, U_*)/\lambda_0(p = 0, U_*)$ tending to $\sqrt{(1 - \ln U_*)/(2 + \ln U_*)}$ as $p \to \infty$; hence more cells should appear per unit transverse length if $p$ is larger (except very close to extinction conditions<fnref n="1" />). Such qualitative predictions will be roughly confirmed in the numerical simulations below. They should be treated however with caution, from the quantitative point of view at least, due to several limitations. An obvious limitation is that the flame extent in the transverse direction is necessarily restricted in computations as well as in experiments, so that wavelengths close to $\lambda_0$ may not be allowed in moderately large domains. Another more important limitation regarding the theoretical results in general is that there are based on the assumption that the Zeldovich number $\beta$ is infinite (used e.g., in the jump conditions (7)), while the numerical computations are carried out with finite values of $\beta$. We shall examine in more detail below this limitation related to the Zeldovich number, whose importance was noted in a number of studies on flame stability, e.g., Denet and Haldenwang [15], Lasseigne et al. [16], Sharpe [17,18], Brailovsky et al. [19]. Finally, it is worth restating that an-

<footnote n="1">¹ Near extinction $\lambda_0$ becomes independent of $p$ as implied by formula (19) in the limit $U_* \to e^{-1/2}$.</footnote>
other obvious limitation of our study is the predominance of the Landau-Darrieus instability for real flames, an aspect which is not considered herein.

### 4.3. The Kuramoto-Sivashinsky equation

Eq. (16) implies the following linear evolution equation for the flame front $f(y, t) \propto e^{st+iky}$:

$$f_t = \left(1 - \frac{l}{l_c}\right) f_{yy} - \left(\frac{1 - p^2}{1 + p^2} - \frac{6}{l_c}\right) U_*^{-2} f_{yyyy}$$

where use has been made of (14). We note that the coefficient of the fourth derivative $f_{yyyy}$ remains negative for all values of $l_c$ in its range [-2, 0) and p in [0, $\infty$), so that the last term is always stabilising as in Sivashinsky's original equation [24].

This linear equation need to be supplemented in the weakly unstable regime by a non-linear term which is required in order to saturate the instability. Such term may be determined by a semi-heuristic kinematic argument as explained in Sivashinsky [10]. The resulting non-linear evolution equation is then given by

$$f_t = \left(1 - \frac{l}{l_c}\right) f_{yy} - \left(\frac{1 - p^2}{1 + p^2} - \frac{6}{l_c}\right) U_*^{-2} f_{yyyy} - \frac{U_*}{2(1 + p^2)^{\frac{1}{2}}} f_y^2 \quad (20)$$

The last term in this equation may be plausibly justified as follows. For the unperturbed planar flame the propagation speed (with respect to the gas) in the negative x direction given by (8d) reads

$$U = U_* (1 + p^2)^{\frac{1}{2}}$$

Consider now that the flame is perturbed such that its local slope is $f_y$. Assuming that the structure of such flame is unchanged from that of a planar tilted flame with slope $f_y$, its local propagation speed in the negative x-direction with respect to the gas can be determined using the governing equations and jump conditions (7) to be

$$U' = U_* (1 + f_y^2 + p^2)^{\frac{1}{2}}$$

From the definition of f it then follows that

$$f_t = U - U' = U_* (1 + p^2)^{\frac{1}{2}} - U_* (1 + p^2 + f_y^2)^{\frac{1}{2}}$$

A Taylor expansion for $|f_y| \ll 1$ then yields

$$f_t = - \frac{U_*}{2(1 + p^2)^{\frac{1}{2}}} f_y^2 + \dots$$

which suggests the non-linear term of Eq. (20).

It is worth noting when $p=0$ and $\kappa=0$ (so that $U_* = 1$ and $l_c = -2$) that Eq. (20) reduces to the Kuramoto-Sivashinsky equation

$$f_t - \left(1 + \frac{1}{2} l\right) f_{yy} + 4 f_{yyyy} + \frac{1}{2} f_y^2 = 0$$

which is essentially Eq. (7) of Sivashinsky [10].

## 5. Numerical simulations and discussion

In this section, we present numerical results based on the time dependent simulations of problem (2) for finite values of $\kappa$, in contrast to the theoretical results presented above based on the limit $\kappa \to \infty$. The simulations are carried out using the finite-element package Comsol Multiphysics. This has been extensively tested in combustion applications, including in our publications [4,5], where more detailed descriptions of the numerical procedure can be found. The domain is covered by a grid of approximately 200,000 triangular elements, with local refinement around the reaction zone. Solutions have been tested to be independent of

the mesh and the size of the domain. In the simulations, periodic boundary conditions in the transverse y direction are imposed, as it is common in studies on the thermo-diffusive flame instabilities or the Kuramoto-Sivashinsky equation [15,26], and the transverse domain size, $L_y$ say, is explicitly specified.

An important point to note is that although U in Eq. (2) refers to the propagation speed of the unperturbed planar front which is time independent, it is considered in the simulations as a time dependent eigenvalue. This is determined at each time step so that the flame remains anchored to the origin of the computational domain, by imposing the constraint $\theta = 0.5$ at the origin. Thus U shall characterize below the instantaneous flame propagation speed at $y=0$ in the negative x-direction. Another useful quantity which will be referred to below, is the effective (or average) flame propagation speed in the negative x-direction, $U_T$ say. This is defined at any given time t as the total burning rate per unit transverse flame area, which in our 2D simulations is given by

$$U_T = \frac{1}{L_y} \iint \omega dxdy \quad (21)$$

In all simulations, the stability of the planar flame is examined by solving an initial value problem, where the initial conditions correspond to the planar flame profiles (obtained numerically as stationary solutions of a time independent problem) to which random perturbations of small amplitude are added.<fnref n="2" />

In presenting the numerical results below, the main focus will be on the cellular instability to which the next three sections are dedicated. This is followed by a fourth section addressing the oscillatory instability.

### 5.1. The cellular instability in small and moderately large domains

We begin with Fig. 6 which illustrates the cellular instability for $L_e = 0.7$ and $\kappa = 10$ for a moderately large transverse domain size, $L_y = 40$, and selected values of p and $\kappa$ specified in the captions. The top row pertains to the adiabatic case $\kappa = 0$ and the bottom row to a non-adiabatic case. Shown are reaction rate contours, which represent stationary cells to which the flame front has evolved starting from an initial profile corresponding to the stationary planar flame solution to which random perturbations of small amplitude are added as described above. We note that more cells are obtained as p is increased, which is in qualitative agreement with our comments regarding the expected cells size following Eq. (19).

We turn now to Fig. 7, which is similar to the previous figure with a smaller transverse size, $L_y = 10$. The figure shows that the flame front is stable when $p=0$, but that it is destabilised by Taylor dispersion, evolving to stationary cells when $p \neq 0$. Thus, the simulations demonstrate that Taylor dispersion is able to destabilise an otherwise thermo-diffusively stable planar front. It turns out, as we shall discuss below, that the finiteness of the Zeldovich number $\kappa$ plays an essential part in the flame destabilisation in the cases of this figure, for which the flame should be stable according to the asymptotic analysis.

### 5.2. The cellular instability in large domains

We now consider the case of large transverse sizes, $L_y > 100$ say, which is known to typically lead to a chaotic front dynam-

<footnote n="2">² For illustration, the initial condition adopted for the mass fraction field is $y_F = \bar{y}_F(x) + \epsilon \text{rand}(x, y)$ where $\bar{y}_F(x)$ corresponds to the stationary planar solution, $\text{rand}(x, y)$ a uniformly distributed random function with mean zero and unit amplitude, and $\epsilon = 10^{-4}$. Numerically, the perturbation term added to $\bar{y}_F$ is replaced by zero when $\bar{y}_F$ is not in the interval $(\epsilon, 1-\epsilon)$ to avoid nonphysical values of $y_F$ lying outside the interval $(0,1)$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Reaction rate fields for $Le = 0.7$, $\beta = 10$ and selected values of $p$ after a longtime evolution of a planar flame front; $p = 0$ (left), $p = 1$ (middle) and $p = 3$ (right). The first row corresponds to $\kappa = 0$ and the second row to $\kappa = 0.11$. The domain transverse size is $Ly = 40$. The cells appearing are stationary. We note that as $p$ is increased, the size of the cells decrease.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Reaction rate fields for $Le = 0.7$, $\beta = 10$ and selected values of $p$ after a longtime evolution of a planar flame front; $p = 0$ (left), $p = 1$ (middle) and $p = 3$ (right). The first row corresponds to $\kappa = 0$ and the second row to $\kappa = 0.11$. The domain transverse size is $Ly = 10$. The cells appearing are stationary. We note that the flame is stable for $p = 0$, but it is destabilised by an increase in $p$, both for $\kappa = 0$ and $\kappa = 0.11$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Temperature fields (left) and reaction rate fields (right) for $Le = 0.7$, $\beta = 10$ and specific values of $p$ at selected times; $p = 0$ (top row), $p = 1$ (middle) and $p = 3$ (bottom). The domain transverse size is $Ly = 120$ and $\kappa = 0$.</figcaption></figure>

ics, at least in the absence of Taylor dispersion [10,15,26]. Two il-
lustrative sets of simulations corresponding to $L_y = 120$ and three
selected values of $p$ are shown in Fig. 8, pertaining to the adia-
batic case $\kappa = 0$, and in Fig. 9, pertaining to a non-adiabatic case.
The figures represent at several values of time $t$ snapshots of the
temperature and reaction rate fields of initially planar fronts sub-
ject to small random perturbations at $t = 0$. In each row, the sec-
ond figure from the left corresponds to a time close to the insta-
bility development time where a cellular structure emerging from
the planar front becomes discernable. The overall picture of the
front dynamics is further clarified in Fig. 10 where the propa-
gation speed $U$ is plotted versus time for the adiabatic cases of
Fig. 8 (top row) and the non-adiabatic cases of Fig. 9 (bottom
row).

Several remarks are worth making in connection with the last
three figures:

(1) The evolution of the planar flame into a front with an appar-
ently chaotic spatiotemporal dynamics is in line with the lit-
erature [10,15,26], and clearly involves the mechanism of tip-
splitting and cells merging emphasized in [15] which seems
to persist in the presence of moderate Taylor dispersion, e.g.,
when $p = 1$.

(2) In general, the effect of an increase in $p$ is seen to have a destabilising effect on the planar front, by shortening e.g., the instability development time observed in Fig. 10 (roughly the initial duration for $U$ or $U_T$ to deviate significantly from its planar value). Furthermore, the typical size of the cells which first become discernable is seen to decrease with increasing $p$. Both these remarks are in line qualitatively with the theoretical predictions of Section 4.2.

(3) An interesting new outcome revealed by the numerical sim-
ulations is that the chaotic time-dependent behaviour fol-
lowing the destabilisation of the planar front is suppressed
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Temperature fields (left) and reaction rate fields (right) for $L_y = 0.7$, $\beta = 10$ and specific values of $p$ at selected times; $p = 0$ (top row), $p = 1$ (middle) and $p = 3$ (bottom). The domain transverse size is $L_y = 120$ and $\kappa = 0.135$.</figcaption></figure>

if $p$ is large enough. This is observed e.g., for $p = 3$ where time-independent stable cellular structures are obtained for large times, as seen in Figs. 8 and 9. Furthermore, the occurrence of such stable cellular structures is found to be quite robust. Indeed, we have checked numerically that similar time-independent behaviours are obtained by changing the periodic lateral boundary conditions from periodic to Neumann conditions, the domain size (to $L_y = 150$ e.g.), the value of $p$ (to $p = 2$ or $p = 4$ e.g.), or the Lewis number (to $L_e = 0.75$ or $L_e = 0.6$ e.g.).

(4) The development of the cellular instability observed for $p = 3$ seems to be somewhat in contradiction with the theoretical predictions of Section 4.2 when analysed carefully, even from the qualitative point of view. For example, the stable stationary structures obtained for large times just mentioned are not obtained when solving the Kuramoto–Sivashinsky Eq. (20) (calculations not shown herein). Furthermore the size of the cells first emerging from the destabilised planar front is quite different from that predicted by the linear stability theory. It turns out, as discussed in the next para-

graph, that the dependence on the Zeldovich number is at the origin of such discrepancies between the numerical and asymptotic studies, as can be surmised from similar studies on flame instabilities [15–19].

### 5.3. Effect of the Zeldovich number on the cellular instability

In this section we examine the effect of the Zeldovich number $\beta$ which will allow us to reconcile apparent discrepancies between the numerical results carried out so far with $\beta = 10$ and the theoretical stability analysis carried out in the asymptotic limit $\beta = \infty$. We shall also identify flame behaviours obtained for finite realistic values of $\beta$ which can be missed in an asymptotic analysis. For simplicity, we shall restrict the discussion to the adiabatic case $\kappa = 0$.

We begin by focusing our attention on the case $p = 3$ of Fig. 8. In this case, the figure reveals that the planar front subject to random perturbations at $t = 0$ evolves into a cellular structure which is first discernable with 12 cells emerging at time $t \approx 80$, suggesting that the fastest linearly growing mode under this confinement
<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Propagation speed *U* versus time *t* for selected values of *p* (solid line). The first row pertains to *κ* = 0 and the cases of Fig. 8. The second row pertains to *κ* = 0.135 and the cases of Fig. 9. The dashed red curves in the first two figures represent the effective propagation speed *U*<sub>*T*</sub> defined in (21). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Growth rate versus the wavenumber *k* for selected values of *p*, *κ* = 0, and reduced Lewis number *l* = −3. The dashed curves are based on the dispersion relation (13) as in Fig. 5(b). The vertical lines represent the discrete wave numbers allowed for periodic solutions in a domain with transverse size *L*<sub>*y*</sub> = 120 given by *k*<sub>*n*</sub> = 2π*n*/*L*<sub>*y*</sub> with *n* = 1, 2, 3, ... For *p* = 3 the most unstable mode is seen to correspond to *k*<sub>7</sub>, indicating that 7 cells are expected to appear near the instability onset for *β* ≫ 1. Numerically, this is found to be true when *β* > 25 approximately, as seen in Fig. 12.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Temperature fields for $l = -3$, $p = 3$ and selected values of $\beta$. The domain transverse size is $L_y = 120$ and $\kappa = 0$.</figcaption></figure>

has wavelength $L_y/12$. Ultimately, the flame is seen to evolve into a stationary 10-cell structure. In principle, the emergence of the 12-cell structure should be predicted by the dispersion relation (13) of the linear stability analysis. To this end, the growth rate versus $k$ based on (13) is reported in Fig. 11 for the case $l = \beta(Le - 1) = -3$ of Fig. 8. Also plotted on the figure are vertical lines corresponding to the discrete values $k_n = 2n\pi/L_y$ ($n = 1, 2, 3...$) of the wavenumber $k$ allowed by the confinement when $L_y = 120$. For $p = 3$ the figure indicates that the most unstable mode corresponds to $k_7$ suggesting that 7 cells should first emerge from the unstable planar front according to the theory, rather than 12 cells as found numerically. To examine the influence of $\beta$ on the findings, the numerical simulations are repeated in this case for various values of $\beta$ with $l = \beta(Le - 1)$ held fixed with $l = -3$. The results are reported in Fig. 12 where snapshots of the temperature field are given at selected times and prescribed values of $\beta$. The figure reveals that the number of cells in the cellular structure first emerging from the unstable planar front decreases with increasing $\beta$; this number is equal to 15 when $\beta = 7$, 12 when $\beta = 10$, 9 when $\beta = 15$, 8 when $\beta = 20$ (not shown), and 7 when $\beta = 25$. Thus, the number of emerging cells predicted by the stability theory agrees with that found numerically provided that $\beta$ is sufficiently large, possibly unrealistically large ($\beta > 25$ approximately in this case). It is worth pointing out that the long-time behaviour when $\beta = 7$, 10, 15 and 20 correspond to a stationary cellular structure (with 10, 10, 9 and 8 cells, respectively). In contrast, when $\beta = 25$ the long-time behaviour corresponds to an unsteady cellular structure with typically 7 cells, with an apparently chaotic dynamics and this behaviour is in line with the predictions based on the Kuramoto-

Sivashinsky Eq. (20). These remarks are confirmed in Fig. 13 where the effective propagation speed $U_T$ is plotted versus $t$ for the cases of Fig. 12; in particular, the apparently chaotic dynamics is clear for $\beta = 25$. Finally, we note that the deviations between the theoretical and numerical results associated with the size of $\beta$ under discussion appear to be most pronounced for larger values of the dispersion coefficient $p$. For example, when $p = 0$, Fig. 11 implies that $k_4$ is the most unstable mode which suggests that a 4-cell-structure should first emerge from the unstable planar front; this is indeed what is observed in the numerical simulations with $\beta = 10$ in Fig. 8.

To close this section, we note that Fig. 11, where the vertical lines associated with the discrete values $k_n$ ($n = 1, 2, 3...$) represent the allowed values of the wavenumber $k$ when $L_y = 120$, can also be used for comparison with the numerical results of Fig. 6 where $L_y = 40$, or with those of Fig. 7 where $L_y = 10$. In the case $L_y = 40$, of course, the allowed values of $k$ are restricted to $k_3, k_6, k_9, k_{12}$, etc., and in the case $L_y = 10$, the allowed values of $k$ are $k_{12}, k_{24}, k_{36}$ and so on. Let us briefly comment on the $L_y = 10$ case. This case is indeed interesting, because according to Fig. 11 all allowed modes ($k_{12}, k_{24}$, etc.) have negative growth rates and are thus stable and thus the planar flame is expected to be stable for $p = 0, 1$ or $3$; according to the numerical simulations in Fig. 7 the planar flame is however unstable when $p = 1$ or $3$. This demonstrates again that the finite-$\beta$ numerical results can show qualitative differences with the theoretical predictions, and such differences are more pronounced for larger values of $p$. Luckily, it is still possible to reconcile the numerical and theoretical findings, by taking $\beta$ to be large enough in the numerical simulations. For ex-
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Effective propagation speed $U_T$ versus time for $l = -3$, $p = 3$ and selected values of $\beta$. The domain transverse size is $L_y = 120$. All curves shown settle asymptotically to stationary values, but the curve corresponding to $\beta = 25$ remains time independent for large times.</figcaption></figure>

ample, when $p = 1$ we have found (calculations not shown) that the planar flame is indeed stable when $\beta > 13$, approximately, and unstable when $\beta < 13$. Similarly, when $p = 3$ we have found that the planar flame is stable when $\beta > 43$ (an unrealistically large value), and stable otherwise. In summary, the asymptotic theory is robust, although the discrepancies aforementioned limit somewhat our confidence in its predictive ability for realistic values of $\beta$, especially for larger values of $p$.

## 5.4. The flame oscillatory instability

The numerical simulations above have focused on the flame cellular instability. For sake of completeness, we now address briefly the flame oscillatory instability. This is the subject of the last two figures, with Fig. 14 pertaining to the adiabatic case $\kappa = 0$ and $Le = 2.5$, and Fig. 15 to the non adiabatic case $\kappa = 0.11$ and $Le = 1.8$. Both figures confirm that a small amount of dispersion can stabilize an otherwise oscillatory unstable flame, in agreement with our theoretical predictions. Note also that the oscillations observed in Fig. 15 are suppressed if $\kappa$ is set to zero for the cases of this figure (calculations not shown), confirming the destabilising role of $\kappa$ suggested by the theory.

As discussed in connection with the cellular instability above, it is worth examining the effect of $\beta$ on the oscillatory instability. To this end, the calculations of Fig. 14 pertaining to $\beta = 15$ have been repeated with a value of $\beta = 10$ while maintaining $l = \beta(Le - 1)$ fixed. It is then found that the planar flame is stable in all cases. A similar outcome is reached by repeating the calculations with $Le$ being held fixed. Thus, we can conclude that a decrease in $\beta$ has always a stabilizing effect on the oscillatory instability. This conclusion is true both in the presence and in the absence of Taylor dispersion; in the latter case ($p = 0$), this observation is in line with the findings of [16].

<figure><img src="image_14.png" /><figcaption>**Fig. 14.** Propagation speed versus time for selected values of $p$; $p = 0$ (top) $p = 0.12$ (middle), $p = 0.25$ (bottom). $\beta = 15$, $Le = 2.5$, $\kappa = 0$ and $L_y = 20$. We note that when $p = 0$ or $p = 0.12$, the flame is unstable (oscillatory instability) but it is stable for $p = 0.25$. The origin of time in this and next figure is chosen after a long transient regime.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 15.** Propagation speed versus time for selected values of *p*; *p* = 0 (top) *p* = 0.1 (middle), *p* = 0.2 (bottom). *β* = 15, *Le* = 1.8, *κ* = 0.11 and *Ly* = 20. We note that when *p* = 0 and *p* = 0.1, the flame is unstable but it is stable for *p* = 0.2.</figcaption></figure>

**6. Conclusion**

In this paper, two original contributions have been made: (1) an analytical investigation of the effect of Taylor dispersion on the thermo-diffusive instabilities of non-adiabatic flames, and (2) an extensive numerical study addressing the influence of Taylor dispersion on flame stability which appears to be the first in the literature.

The analytical investigation consists of a linear stability analysis carried out in the asymptotic limit $\beta \to \infty$, where $\beta$ is the Zeldovich number. This has led to the dispersion relation (13), encapsulating the effects of the three parameters $l$ (the reduced Lewis number), $p$ (the dispersion coefficient) and $\kappa$ (the volumetric heat loss coefficient). The dispersion relation derived reduces as it should to the dispersion relation derived by Sivashinsky [24] when $p=0$ and $\kappa=0$, to that derived by Joulin and Clavin [25] when $p=0$, and to that derived by Daou [6] when $\kappa=0$. A thorough discussion of the implications of the analytical findings on the flame cellular and oscillatory instabilities has been provided. This has included the determination of bifurcation curves and domains of stability instability in the parameters space (see Figs. 3 and 4) and the derivation of the Kuramoto-Sivashinsky type Eq. (20) involving the parameters $l$, $p$ and $\kappa$ and characterising the flame dynamics in the weakly non-linear regime near the onset of the cellular instability. The theoretical results demonstrate the ability of Taylor dispersion and heat loss to significantly affect the stability of the flame. For example, while the oscillatory instability is promoted by an increase in $\kappa$, it is hampered by an increase in $p$. On the other hand, both $p$ and $\kappa$ have a destabilising effect in connection with the cellular instability, and the typical size of cells first emerging from such instability can be predicted (near the instability onset) by a simple formula given by (19). In general, this typical size is an decreasing function of $p$, suggesting that more cells should appear by unit transverse length for larger valued of $p$.

A numerical investigation has also been conducted which largely confirms and significantly extends the theoretical findings. In particular, the Zeldovich number $\beta$ is found to be at the origin of apparent discrepancies between the theoretical and numerical results since $\beta$ is assumed to be infinitely large in the asymptotic analysis while taking finite values, $\beta \approx 10$, numerically. Such discrepancies are more pronounced for larger values of $p$ and can have in fact a qualitative nature, in addition to being quantitative. For example, when the domain lateral size is relatively large as in the case of Fig. 13 corresponding to $p=3$, the long-time behaviour of an initially planar flame developing a cellular instability is found to be typically stationary for realistically large values of $\beta$, say $\beta < 25$, and chaotic for larger values or in the limit $\beta \to \infty$. Similarly, when the domain lateral size is relatively small as in the case of Fig. 7, a flame which is stable according to the asymptotic theory as $\beta \to \infty$, may turn out to be unstable for finite large values of $\beta$. Luckily, the asymptotic theory is found to be robust in the sense that its predictions are recovered numerically if $\beta$ is taken large enough, although such predictions should be questioned for realistic values of $\beta$.

Finally, here is a compact summary of the asymptotic and numerical findings: an increase in $p$ or a decrease in $\beta$ have a destabilising effect in connection with the cellular instability, and a stabilizing effect in connection with the oscillatory instability; both instabilities are promoted by an increase in $\kappa$.

We close the paper by highlighting that this investigation has adopted the simplifying assumption of constant density, which is quite common in studies on the thermo-diffusive flame instabilities [15,25,27] and Clavin and Searby [28, p. 77–81, 473–480]. This has been done deliberately for sake of analytical tractability, and also in order to switch off the Darrieus-Landau (DL) instability so as to be able to focus on the thermo-diffusive instabilities as these are directly affected by Taylor dispersion. Of course, the DL instability is prevalent for real flames, and the thermo-diffusive instabilities are found to be superimposed on them in suitable mixtures such as propane rich or hydrogen lean mixtures, as observed in photographs of real flames [10,11–28], p. 78]. It is interesting to examine how Taylor dispersion affects the coupling between the DL and the thermo-diffusive instabilities, an in particular how gas expansion may influence the findings of this paper. This is as yet an open and challenging question, which we shall address in the near future as a natural follow up of the current work, both analytically and numerically. From the analytical point of view, we anticipate that the stability analysis will be facilitated by adopting first a depth-averaged model which generalises the one used herein by incorporating the coupling between variable density and Taylor dispersion as in Pearce and Daou [4], Daou et al. [5]. As a framework for such future analytical and numerical studies, it is convenient to continue adopting the Hele-Shaw channel configuration which has been the focus of recent experiments and numerical simulations such as those reported in the publications [12,13]. In these publications, the effect of forced convection on flame stability has not been considered; this ingredient is essential however to reveal the effect of Taylor dispersion and corresponding diffusion anisotropy on the flame. In fact, forced convection may be incorporated in theoretical and numerical models by prescribing an incoming (Poiseuille) flow, which corresponds for example to maintaining an upwards flow in the vertical Hele-Shaw burner used in Al Sarraf et al. [12]. Conceptually, a Couette shear flow may also be envisaged, which would correspond to the channel walls moving in opposite directions, although this may be more challenging to the experimentalist. Furthermore, the direction of the shear flow may be taken, in theory, to be arbitrary with respect to the direction of propagation of the unperturbed flame, and this may provide a deeper understanding of the effect of shear induced anisotropic diffusion on flame propagation.
## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

JD is grateful to the EPSRC for financial support under grant [EP/V004840/1].

## Appendix A. Derivation of the dispersion relation

This appendix is dedicated to the derivation of the dispersion relation (13). We begin by substituting the perturbed solutions (12) into Eqs. (3)–(7). This leads to an eigen-boundary value problem for the functions $\hat{\theta}(x)$ and $\hat{h}(x)$ which is given by

$$ (1 + p^2)\hat{\theta}'' - U\hat{\theta}' - (s + k^2)\hat{\theta} = 0 \quad (A.1) $$

$$ (1 + p^2)\hat{h}'' - U\hat{h}' - (s + k^2)\hat{h} = -l\left\{(1 - p^2)\hat{\theta}'' - k^2\hat{\theta}\right\} + \kappa\hat{\theta} \quad (A.2) $$

for $x \neq 0$, subject to the boundary conditions

$$ \hat{\theta} = 0 \quad \hat{h} = 0 \quad \text{as} \quad x \to \pm\infty \quad (A.3) $$

and the jump conditions

$$ \|\hat{\theta}\| = \frac{U}{1 + p^2} \qquad \|\hat{h}\| = -\frac{U l (1 - p^2)}{(1 + p^2)^2} \quad (A.4) $$

$$ \|\hat{h}'\| + \frac{1 - p^2}{1 + p^2} l \|\hat{\theta}'\| = - \frac{U^2 l (1 - p^2)}{(1 + p^2)^3} \quad (A.5) $$

$$ \|\hat{\theta}'\| = \frac{U^2}{(1 + p^2)^2} - \frac{U}{1 + p^2} \left( \frac{\hat{h}(0^+)}{2} - \frac{\kappa}{2U} \right) \quad (A.6) $$

at $x=0$. It is to be noted that the jump conditions at the reaction sheet $x=f(y,t)$ have been transferred to $x=0$, using Taylor expansions for small values of $\delta$, since $f=O(\delta)$ according to (12).

The solution of the linear second order differential Eqs. (A.1) and (A.2) subject to all conditions except (A.5) is found to be

$$ \hat{\theta} = \begin{cases} \frac{U}{1+p^2} \exp(r^+ x), & r \to 0 \\ \left[ \frac{U(1-\Gamma)}{1+p^2} + \frac{\kappa}{U} + \frac{U l (1-p^2)}{(1+p^2)^2} + \chi x \right] \exp(r^+ x), & r > 0 \end{cases} $$

where

$$ \Gamma = \left(1 + \frac{4}{U^2}(1 + p^2)(s + k^2)\right)^{\frac{1}{2}}, \quad r^{\pm} = \frac{U(1 \pm \Gamma)}{2(1 + p^2)} \quad \text{and} $$

$$ \chi = \frac{l \left\{ \left[ U^2 (\Gamma + 1) + 2(1 + p^2)s \right] (1 - p^2) - 4(1 + p^2)p^2 k^2 \right\}}{2\Gamma(1 + p^2)^3} - \frac{\kappa}{\Gamma(1 + p^2)} $$

In deriving the solution above, unstable modes characterised by $Re(s) \ge 0$ are assumed; this assumption implies that $Re(\Gamma^+) \ge 1$, and hence $Re(r^+) \ge 0$, and $Re(r^- \le 0)$. Use of these expressions in the jump condition (A.5) leads, after a few algebraic manipulations, to the dispersion relation (13).

## References

[1] G.I. Taylor, Dispersion of soluble matter in solvent flowing slowly through a tube, Proc R Soc Lond. Ser A 219 (1137) (1953) 186-203.

[2] R. Aris, On the dispersion of a solute in a fluid flowing through a tube, Proc R Soc Lond. Ser A 235 (1200) (1956) 67-77.

[3] H. Brenner, D. Edwards, Macrotransport Processes, Butterworth-Heinemann, 1993.

[4] P. Pearce, J. Daou, Taylor dispersion and thermal expansion effects on flame propagation in a narrow channel, J. Fluid Mech. 754 (2014) 161.

[5] J. Daou, P. Pearce, F. Al-Malki, Taylor dispersion in premixed combustion: questions from turbulent combustion answered for laminar flames, Phys. Rev. Fluids 3 (2) (2018) 023201.

[6] J. Daou, Effect of Taylor dispersion on the thermo-diffusive instabilities of flames in a Hele-Shaw burner, Combust. Theor. Model. 25 (4) (2021) 765-783.

[7] A. Liñán, P. Rajamanickam, A.D. Weiss, A.L. Sánchez, Taylor-diffusion-controlled combustion in ducts, Combust. Theor. Model. 24 (6) (2020) 1054-1069.

[8] P. Rajamanickam, A.D. Weiss, Effects of thermal expansion on Taylordiffusion-controlled diffusion flames, Combust. Theor. Model. 26 (1) (2022) 50-66.

[9] J. Buckmaster, The structure and stability of laminar flames, Annu. Rev. Fluid Mech. 25 (1) (1993) 21-53.

[10] G.I. Sivashinsky, Instabilities, pattern formation, and turbulence in flames, Annu. Rev. Fluid Mech. 15 (1) (1983) 179-199.

[11] M. Matalon, Intrinsic flame instabilities in premixed and nonpremixed combustion, Annu. Rev. Fluid Mech. 39 (2007) 163-191.

[12] E. Al Sarraf, C. Almarcha, J. Quinard, B. Radisson, B. Denet, Quantitative analysis of flame instabilities in a Hele-Shaw burner, Flow, Turbul. Combust. 101 (3) (2018) 851-868.

[13] D. Fernández-Galisteo, V.N. Kurdyumov, P.D. Ronney, Analysis of premixed flame propagation between two closely-spaced parallel plates, Combust. Flame 190 (2018) 133-145.

[14] G. Joulin, G. Sivashinsky, Influence of momentum and heat losses on the large-scale stability of quasi-2D premixed flames, Combust. Sci. Technol. 98 (1-3) (1994) 11-23.

[15] B. Denet, P. Haldenwang, Numerical study of thermal-diffusive instability of premixed flames, Combust. Sci. Technol. 86 (1-6) (1992) 199-221.

[16] D. Lasseigne, T. Jackson, L. Jameson, Stability of freely propagating flames revisited, Combust. Theor. Model. 3 (4) (1999) 591.

[17] G. Sharpe, Linear stability of planar premixed flames: reactive Navier-Stokes equations with finite activation energy and arbitrary lewis number, Combust. Theor. Model. 7 (1) (2003) 45.

[18] G.J. Sharpe, Thermal-diffusive instability of premixed flames for a simple chain-branching chemistry model with finite activation energy, SIAM J. Appl. Math. 70 (3) (2009) 866-884.

[19] I. Brailovsky, P.V. Gordon, L. Kagan, G. Sivashinsky, Diffusive-thermal instabilities in premixed flames: stepwise ignition-temperature kinetics, Combust. Flame 162 (5) (2015) 2077-2086.

[20] J. Daou, J. Dold, M. Matalon, The thick flame asymptotic limit and Damköhler's hypothesis, Combust. Theor. Model. 6 (1) (2002) 141-153.

[21] J. Daou, M. Matalon, Influence of conductive heat-losses on the propagation of premixed flames in channels, Combust. Flame 128 (4) (2002) 321-339.

[22] J. Buckmaster, G. Ludford, Lectures on Mathematical Combustion, Society for Industrial Mathematics, Philadelphia, 1987.

[23] G. Joulin, P. Vidal, An introduction to the instability of flames, shocks and detonations, in: C. Godreche, P. Manneville (Eds.), Hydrodynamics and Non-linear Instabilities, Collection Alea-Saclay: Monographs and Texts in Statistical Physics, Cambridge University Press, Cambridge, 1998.

[24] G. Sivashinsky, Nonlinear analysis of hydrodynamic instability in laminar flames—I. Derivation of basic equations, AcAau 4 (11) (1977) 1177-1206.

[25] G. Joulin, P. Clavin, Linear stability analysis of nonadiabatic flames: diffusion-thermal model, Combust. Flame 35 (1979) 139-153.

[26] J.M. Hyman, B. Nicolaenko, The Kuramoto-Sivashinsky equation: a bridge between PDE's and dynamical systems, Physica D 18 (1-3) (1986) 113-126.

[27] G. Sivashinsky, Diffusional-thermal theory of cellular flames, Combust. Sci. Technol. 15 (3-4) (1977) 137-145.

[28] P. Clavin, G. Searby, Combustion Waves and Fronts in Flows: Flames, Shocks, Detonations, Ablation Fronts and Explosion of Stars, Cambridge University Press, 2016.