<figure><img src="image_1.png" /></figure>

HAL
open science

# Three-dimensional diffusive-thermal instability of flames propagating in a plane Poiseuille flow

Aiden Kelly, Prabakaran Rajamanickam, Joel Daou, Julien R Landel

► To cite this version:

Aiden Kelly, Prabakaran Rajamanickam, Joel Daou, Julien R Landel. Three-dimensional diffusive-thermal instability of flames propagating in a plane Poiseuille flow. Proceedings of the Combustion Institute, 2024, 40 (1-4), pp.105258. 10.1016/j.proci.2024.105258. hal-04651782

HAL Id: hal-04651782

https://hal.science/hal-04651782v1

Submitted on 17 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Three-dimensional diffusive-thermal instability of flames propagating in a plane Poiseuille flow

Aiden Kelly<sup>a</sup>, Prabakaran Rajamanickam<sup>a</sup>, Joel Daou<sup>a,*</sup>, Julien R. Landel<sup>a,b</sup>

## ARTICLE INFO

Keywords:

Diffusive-thermal instability

Poiseuille flow

Flame-flow interaction

Asymmetric flame

## ABSTRACT

The three-dimensional diffusive-thermal stability of a two-dimensional flame propagating in a Poiseuille flow is examined. The study explores the effect of three non-dimensional parameters, namely the Lewis number $Le$, the Damköhler number $Da$, and the flow Peclet number $Pe$. Wide ranges of the Lewis number and the flow amplitude are covered, as well as conditions corresponding to small-scale narrow ($Da \ll 1$) to large-scale wide ($Da \gg 1$) channels. The instability experienced by the flame appears as a combination of the traditional diffusive-thermal instability of planar flames and the recently identified instability corresponding to a transition from symmetric to asymmetric flame. The instability regions are identified in the $Le-Pe$ plane for selected values of $Da$ by computing the eigenvalues of a linear stability problem. These are complemented by two- and three-dimensional time-dependent simulations describing the full evolution of unstable flames into the non-linear regime. In narrow channels, flames are found to be always symmetric about the mid-plane of the channel. Additionally, in these situations, shear flow-induced Taylor dispersion enhances the cellular instability in $Le < 1$ mixtures and suppresses the oscillatory instability in $Le > 1$ mixtures. In large-scale channels, however, both the cellular and the oscillatory instabilities are expected to persist. Here, the flame has a stronger propensity to become asymmetric when the mean flow opposes its propagation and when $Le < 1$; if the mean flow facilitates the flame propagation, then the flame is likely to remain symmetric about the channel mid-plane. For $Le > 1$, both symmetric and asymmetric flames are encountered and are accompanied by temporal oscillations.

## 1. Introduction

Flame propagation in a Poiseuille flow is a classical problem relevant to many practical combustion devices with large and small scales. It provides a suitable framework to investigate fundamental flame-flow interaction phenomena such as flame blow-offs, flashbacks [1] and Taylor dispersion effects on flames. The two-dimensional channel configuration offers an attractive platform to investigate such flame phenomena both from the theoretical and experimental viewpoints, see e.g. [2–4]. This platform is also ideal to investigate flame instabilities, including the Darrieus–Landau and diffusive-thermal instabilities as done in [5–7]. This study will complement the findings of such investigations, by focusing specifically on the effect of the flow scale (channel width) and amplitude on the diffusive-thermal instability.

Fig. 1 shows a schematic representation of the configuration adopted in our study. Shown are symmetric flames whose propagation is either opposed or aided by the Poiseuille flow. The flow field is prescribed

and its components are given by

$$v_x^* = \frac{3U}{2} \left(1 - \frac{y^*}{h^2}\right), \quad v_y^* = v_z^* = 0. \qquad (1)$$

Here $x^*$ is the streamwise coordinate, $y^*$ the wall-normal coordinate, $z^*$ the spanwise coordinate, $U$ the mean flow speed, and $h$ the channel half-width. Denoting by $D_T$ the thermal diffusivity and by $D_F$ the fuel diffusion coefficient, we can define three relevant non-dimensional parameters, namely the Lewis number $Le$, the flow Peclet number $Pe$ and the Damköhler number $Da$ by

$$Le = \frac{D_T}{D_F}, \quad Pe = \frac{Uh}{D_T}, \quad Da = \frac{h^2/D_T}{\delta_L^2/D_T}, \qquad (2)$$

where $\delta_L$ is the laminar flame thickness. Large and small values of $Pe$ here indicate strong and weak flow intensities, respectively. Also, large values of $Da$ describe wide channels (or thin flames) while small values pertain to narrow channels (or thick flames).

<footnote><sup>a</sup>Department of Mathematics, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote><sup>b</sup>Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, Ecole Centrale de Lyon, INSA Lyon, LMFA, UMR5509, 69622 Villeurbanne, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: joel.daou@manchester.ac.uk (J. Daou).</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 1. Schematic illustration of a premixed flames whose propagation is opposed by a two-dimensional Poiseuille flow (top) or aided by the flow (bottom). The flame propagation speed is $S_T - U$ with respect to the walls, where $S_T$ is the effective burning speed, i.e. the flame speed with respect to the mean flow.</figcaption></figure>

An investigation of the diffusive-thermal instability of flames such as the ones in Fig. 1 was undertaken by Kurdyumov [8] for various values of Da, Le and Pe. The investigation revealed interesting results, in particular, it emphasised the propensity of the flame to become asymmetric (with respect to $y^* = 0$) following the loss of its initially symmetric shape for certain values of the parameters. It can be noted, however, that the stability analysis in [8] is two-dimensional as it considers perturbations independent of the $z^*$-coordinate. In this paper, we shall discard this restriction by considering a three-dimensional analysis. This is important because, as we shall confirm, the dominant unstable mode turns out in most cases to correspond to perturbations depending on the spanwise coordinate $z^*$.

A stability analysis considering perturbations depending on the spanwise coordinate $z^*$ has been carried out in [9], but it was limited to thick flames in narrow channels, that is to cases with $Da \ll 1$. The main objective of the current paper is to remove the restrictions associated with the perturbations being $z^*$-independent or $Da \ll 1$. This will allow us to provide a more complete characterisation of the flame diffusive-thermal instabilities, applicable for three-dimensional flames in narrow and wide channels. To this end, we shall adopt a simple constant-density adiabatic model as done in [8,9] in order to focus on the diffusive-thermal instability. The reader may refer to the literature for information about the influence of other factors pertinent to flame propagation in a Poiseuille flow which are not taken into account herein. Such factors include heat loss [10,11], equivalence ratio [12,13], complex chemistry [14] and an axisymmetric geometry [15].

## 2. Governing equations

We shall ignore the effects of thermal expansion and heat loss in order to focus on the diffusive-thermal instability, as mentioned above and as done in [8,9]. Under these conditions, the velocity field may be prescribed and is given by (1). As shown in Fig. 1, the flow direction is from left to right when $U > 0$ and from right to left when $U < 0$. A two-dimensional symmetric premixed flame, independent of $z^*$, is depicted which propagates in the negative $x^*$-direction with constant speed $S_T - U$ with respect to the walls. Here, $S_T$ is the flame effective burning speed, that is the flame speed with respect to the mean flow. To render the flame steady, a reference frame moving with the flame is used by introducing the coordinate transformation $(x^*, y^*, z^*) \mapsto (x^* + S_T t^* - Ut^*, y^*, z^*)$. In this frame, the velocity field (1) becomes

$$v_x^* = \frac{U}{2} \left( 1 - \frac{3y^*z^*}{h^2} \right) + S_T, \quad v_y^* = v_z^* = 0.$$

The analysis assumes fuel-lean conditions for the reactive mixture and adopts a one-step reaction model. The mass of fuel burnt per unit volume and unit time can then be written as $B\rho Y_F e^{-E/RT}$, where B is the pre-exponential factor, $\rho$ is the constant density, $Y_F$ is the fuel mass fraction, $E/R$ is the activation temperature and T is the temperature. We also define the heat-release parameter $\alpha$, the Zeldovich number $\beta$ and the planar laminar burning speed $S_L$ (for $\beta \gg 1$) by the usual expressions

$$\alpha = \frac{T_{ad} - T_u}{T_u}, \quad \beta = \frac{E(T_{ad} - T_u)}{RT_{ad}^2}, \quad S_L = (2Le\beta^{-2}BD_T e^{-E/RT_{ad}})^{1/2}.$$

Here, $T_{ad} = T_u + qY_{F,u}/c_p$ is the adiabatic flame temperature, q the amount of heat released per unit mass of fuel burnt, $Y_{F,u}$ the fuel mass fraction in the unburnt mixture, and $c_p$ the constant-pressure specific heat. The laminar flame thickness $\delta_L$ is defined by $\delta_L = D_T/S_L$ and the non-dimensional effective burning speed S by $S = S_T/S_L$.

We adopt the non-dimensionalised variables

$$t = \frac{D_T t^*}{\delta_L^2}, \quad (x, z) = \frac{1}{\delta_L}(x^*, z^*), \quad y = \frac{y^*}{h}, \quad y_F = \frac{Y_F}{Y_{F,u}}, \quad \theta = \frac{T - T_u}{T_{ad} - T_u}.$$

and introduce the linear operator

$$\mathcal{L}_{Le} \equiv \left[ \frac{Pe(1 - 3y^2)}{2\sqrt{Da}} + S \right] \frac{\partial}{\partial x} - \frac{1}{Le} \nabla^2 \quad (3)$$

where $\nabla^2 = \partial^2/\partial x^2 + Da^{-1}\partial^2/\partial y^2 + \partial^2/\partial z^2$ and the parameters Pe, Le and Da are as given in (2). With these notations, the governing equations take the form

$$\frac{\partial y_F}{\partial t} + \mathcal{L}_{Le} y_F = -\omega, \quad \frac{\partial \theta}{\partial t} + \mathcal{L}_1 \theta = \omega \quad (4)$$

where $\mathcal{L}_1 = \mathcal{L}_{Le=1}$ and

$$\omega(y_F, \theta) = \frac{\beta^2 y_F}{2Le} \exp \left[ \frac{-\beta(1-\theta)}{1-\alpha(1-\theta)} \right]. \quad (5)$$

The corresponding boundary conditions are prescribed below.

## 3. Steadily propagating symmetric flames

The steady two-dimensional solutions whose stability will be investigated correspond to solutions of (4) which are independent of $z$ and $t$ and symmetric with respect to the $y=0$ plane. These are denoted by $y_F = Y(x, y)$ and $\theta = \Theta(x, y)$, and satisfy the equations

$$\mathcal{L}_{Le} Y = -\omega, \quad \mathcal{L}_1 \Theta = \omega, \quad (6)$$

subject to the boundary conditions

$$x \to -\infty : \quad Y - 1 = \Theta = 0, \quad (7)$$

$$x \to +\infty : \quad Y = \Theta - 1 = 0, \quad (8)$$

$$y = \pm 1 : \quad \frac{\partial Y}{\partial y} = \frac{\partial \Theta}{\partial y} = 0. \quad (9)$$

Since the solutions sought are symmetric, the problem can be solved on the domain corresponding to the upper-half infinite strip defined by the intervals $x \in (-\infty, \infty)$ and $y \in [0, 1]$ using the boundary conditions

$$\frac{\partial Y}{\partial y} = \frac{\partial \Theta}{\partial y} = 0 \quad \text{at} \quad y = 0, 1. \quad (10)$$

Selecting $\beta = 10$ and $\alpha = 0.85$ as representative values for $\beta$ and $\alpha$, the effective burning speed S and the functions Y and $\Theta$ are determined numerically for given values of Le, Pe and Da. These computations are performed using COMSOL Multiphysics as described in [11,16]. They are based on the finite-element method and use a non-uniform grid with typically 300,000 triangular elements along with local refinement around the reaction zones. A summary of the results is provided in Figs. 2 and 3.

Fig. 2(a) shows the dependence of S on Da for selected values of Pe and Le, whereas Fig. 2(b) displays its dependence on Pe for Le = 1 and selected values of Da. As can be seen, the scaled burning
<figure><img src="image_6.png" /><figcaption>Fig. 2. The burning speed S computed for selected values of Da, Le and Pe. All computations are performed with β = 10 and α = 0.85.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Reaction-rate ω-fields for selected values of Pe, Le and Da with β = 10 and α = 0.85. The x and y length scales are indicated in each subfigure. The unburnt gas is in all cases to the left of the flame.</figcaption></figure>

speed S exhibits a non-monotonic variation with respect to Da for given values of Pe and Le. This variation is associated with the fact that for small values of Da, the diffusion transport and therefore S is influenced by Taylor dispersion, whereas at large values of Da, S is strongly influenced by curvature effects, as explained in [16].

The trends observed in Fig. 2(b) are discussed in [2]. Particularly, the burning speed grows quadratically for small values of Pe and grows linearly for |Pe| ≫ 1 although with different slopes for positive and negative values of Pe, as predicted in [2].

A sample of the computed results are presented in Fig. 3 where reaction-rate fields are shown for selected values of Le, Pe and Da. The top subfigure, corresponding to Le = 1, illustrates the transition as Da is increased from Taylor dispersion-controlled thick flames to curved thin flames. The middle and bottom subfigures, pertaining to Le ≠ 1 and Da = 3, confirm the well-known fact that preferential diffusion and curvature effects determine local burning rates in thin or moderately thick flames. For example, the local reaction rate is seen to increase when Le < 1 in regions where the flame is convex towards the unburnt gas (the left side of the flame in the figure). At low values of Da, such curvature-related effects are in fact negligible although the flame is still affected by Le through Taylor dispersion [16].

## 4. Linear stability problem

The linear stability analysis for the steady flames covered in the previous section will now be pursued. To do that, we shall assume

$$ \begin{bmatrix} y_F \\ \theta \end{bmatrix} = \begin{bmatrix} Y(x, y) \\ \Theta(x, y) \end{bmatrix} + \begin{bmatrix} \hat{Y}(x, y) \\ \hat{\Theta}(x, y) \end{bmatrix} e^{ikz + \sigma t} \quad (11) $$

with |ŷ| ≪ |Y| and |θ̂| ≪ |Θ|. Here, k denotes the spanwise wavenumber which is a real number and σ is the growth rate of the perturbation,

which is in general a complex-valued eigenvalue. The equations governing Ŷ(x, y) and ˆΘ(x, y), obtained by substituting (11) into (4) and linearising, read

$$ (\mathcal{L}_{Le} + f + k^2/Le + \sigma) \hat{Y} = -g\hat{\Theta}, \quad (12) $$

$$ (\mathcal{L}_1 - g + k^2 + \sigma) \hat{\Theta} = f \hat{Y}, \quad (13) $$

where the functions f(Θ) and g(Θ, Y) are given by

$$ f = \frac{g}{\beta Y} [1 + \alpha(\Theta - 1)]^2 = \frac{\beta^2}{2Le} \exp \left[ \frac{\beta(\Theta - 1)}{1 + \alpha(\Theta - 1)} \right]. $$

The appropriate boundary conditions are given by

$$ x \to \pm\infty : \hat{Y} = \hat{\Theta} = 0, \quad (14) $$

$$ y = \pm 1 : \frac{\partial \hat{Y}}{\partial y} = \frac{\partial \hat{\Theta}}{\partial y} = 0. \quad (15) $$

The solutions to the eigenvalue problem (12)–(15) fall into two disjoint groups: even solutions (with respect to y) for which Ŷ(x, y) = Ŷ̂(x, -y) and ˆΘ(x, y) = ˆΘ̂(x, -y) and odd solutions for which Ŷ̂(x, y) = -ŷ(x, -y) and ˆΘ̂(x, y) = - ˆΘ̂(x, -y). The problem can thus be solved in the upper-half infinite strip, as done for determining the steady solutions, by imposing the conditions

$$ y = 0 : \frac{\partial \hat{Y}}{\partial y} = \frac{\partial \hat{\Theta}}{\partial y} = 0 \text{ for even solutions,} \quad (16) $$

$$ y = 0 : \hat{Y} = \hat{\Theta} = 0 \text{ for odd solutions.} \quad (17) $$

## 5. Indices for eigenfunctions and terminology

A spatial eigenfunction and its eigenvalue σ can be indexed using three numbers m, n and k as

$$ \hat{Y}_{mn}(x, y)e^{ikz}, \quad \sigma_{mn} \equiv \sigma_{mn}(k). \quad (18) $$
<figure><img src="image_8.png" /><figcaption>Fig. 4. Stable to unstable bifurcation of particular mth dispersion curve $\sigma_{mn}(k)$ as physical parameters are varied, for fixed values of $n$. Four types of bifurcations are represented, as described in the text. Solid lines pertain to eigenvalues having zero imaginary part, while dashed lines to eigenvalues having a non-zero imaginary part.</figcaption></figure>

The spanwise wavenumber $k$ is continuous (since unbounded domain in the $z$-direction is assumed) and lies in the interval $k \in [0, \infty)$. On the other hand, the indices $m$ and $n$, which characterise eigenmodes in the $x$ and $y$ directions, respectively, are discrete and can be ordered such that $m = 1, 2, 3, ...$ and $n = 0, 1, 2, 3, ...$. We note that $n$ represents the number of zero-crossings of the eigenfunctions in the $y$-direction. It is also important to emphasise that the $n$-modes are arranged such that even solutions correspond to even integer values of $n$ and odd solutions to odd integer values of $n$. For example, when $Pe = 0$, it can be checked that

$$ \hat{Y}_{mn}(x, y) = \tilde{Y}_m(x) \begin{cases} \cos(n\pi y/2), & \text{even } n \\ \sin(n\pi y/2), & \text{odd } n \end{cases} \\ \sigma_{mn}(k) = \tilde{\sigma}_m(\sqrt{k^2 + n^2\pi^2/4Da}) $$

where $\tilde{\sigma}_m(k)$ and $\tilde{Y}_m(x)$ are eigenvalues and corresponding eigenfunctions of the classical linear stability problem of the planar flame. The discreteness of the index $n$ is due to confinement in the $y$-direction, while the discreteness of $m$, as in the planar flame stability problem, results from requiring the perturbations to decay as $x \to \pm\infty$. Furthermore, since the even eigenfunctions do not affect the symmetry of the base solution about the $y=0$ plane, we shall refer to them as *symmetric* eigenmodes. Similarly, we shall refer to the odd eigenfunctions as *asymmetric* eigenmodes.

With the notation introduced above, our main task is to obtain the eigenvalues $\sigma_{mn}(k)$ and corresponding eigenfunctions. These are determined by solving numerically the two-dimensional eigenvalue problem (12)–(17) using the eigenvalue solver in COMSOL Multiphysics. In particular, for a given value of $n$, we can plot different dispersion curves $\sigma_{mn}(k)$ each corresponding to different values of $m$. These curves undergo changes when relevant parameters such as $Le$ are varied. Instability occurs when a portion of any of such dispersion curves corresponds to positive growth rates such that $Re\{\sigma_{mn}\} > 0$. Four categories of bifurcations from stable to unstable states are obtained in our problem when plotting $\sigma_{mn}(k)$ as illustrated in Fig. 4. Such bifurcations, which are encountered as a parameter is varied, are referred to as type-II$_s$, type-I$_o$, type-III$_s$ and type-III$_o$ following the terminology used in [17]. The subscript $s$ denotes stationary or non-oscillatory bifurcation whereas $o$ denotes oscillatory bifurcation. As can be seen from the schematic illustration shown in Fig. 4, type-II$_s$ bifurcation corresponds to a dispersion curve which always passes through the origin in the $k-\sigma$ plane whose concavity changes there as a control parameter is varied. On the other hand, type-I$_o$ bifurcation corresponds to a maximum of the dispersion curve crossing the horizontal axis at a wavenumber $k \neq 0$. The other two type-III bifurcations are similar to type-I bifurcation, except that the crossing occurs at $k = 0$ (planar mode). Such bifurcations were encountered in recent studies [18,19] that were dedicated to the limiting case $Da \ll 1$ with a focus on the role played by the direction of flame propagation.

The classical diffusive-thermal cellular instability in premixed flames is represented by type-II$_s$ bifurcation shown in Fig. 4(a), whereas the classical diffusive-thermal oscillatory instability is represented by type-I$_o$ bifurcation shown in Fig. 4(b); see e.g. [20, pp. 477–479]. These two instabilities arise for usual premixed flames in $Le < 1$ and $Le > 1$ mixtures, respectively. In addition to these two familiar instabilities, we

also come across the types shown in Fig. 4(c) and (d) where the most unstable mode occurs at $k = 0$ and therefore has no dependence on the $z$-coordinate. Our computations indicate that these type-III bifurcations are linked to eigenmodes with odd values of $n$ which lead to asymmetric flames.

## 6. Stability regime diagrams

A convenient way to summarise the findings of the stability analysis is to determine the boundaries of the regions of instability in the $Le-Pe$ plane for fixed values of $Da$. This is carried out in Fig. 5 for $Da = 0.01, 3, 30$ and $100$. The white regions here correspond to stable, symmetric flames, and the shaded regions to unstable flames. The solid lines represent the bifurcation curves and are labelled according to the type of bifurcations as defined in the previous section and in Fig. 4. Cellular instabilities are encountered for $Le < 1$ and oscillatory ones for $Le > 1$. It is also found that asymmetric flames can appear both in $Le < 1$ and $Le > 1$ mixtures.

Fig. 5(a), pertaining to $Da = 0.01$, exhibit the instability regions which occur for thick flames. The cellular instability is found to occur for $Le < 1$ with its boundary in the $Le-Pe$ plane being independent of $Pe$, although the dispersion curves themselves, similar to those in Fig. 4(a), do depend on $Pe$. Drastic changes however happen for $Le > 1$ mixtures since a slight increase in $Pe$ is found to completely suppress the oscillatory instability. These results caused by Taylor dispersion are consistent with the findings reported in [9] for $Da \ll 1$. Furthermore, within the range of $Pe$ and $Le$ shown in Fig. 5(a), asymmetric flames are not encountered and the instability domains appear invariant under the transformation $Pe \mapsto -Pe$.

Fig. 5(b) exhibit the instability regions for $Da = 3$. Here we can observe that the boundary (left black curve) of the cellular instability region depends on $Pe$, unlike in the case of Fig. 5(a). More importantly, we now also have a type-III bifurcation associated with the emergence of asymmetric flames for $Le < 1$ and $Pe > 0$ (blue curve). Therefore, flame propagating against a Poiseuille flow ($Pe > 0$) for subunity Lewis numbers can exhibit a cellular pattern in the spanwise $z$-direction and an asymmetric pattern in the $y$-direction which is normal to the wall. However, if the flame propagation is aided by the flow ($Pe < 0$), then only cellular symmetric flames can arise. The direction of flame propagation relative to the flow direction is therefore an important factor in determining the emergence of asymmetric flames.

Fig. 5(c) pertains to the moderately large value $Da = 30$, and shows that the instability regions for $Le < 1$ are qualitatively similar to those in Fig. 5(b). Here the asymmetric flames are more readily accessible as they are obtained even for smaller values of $Pe$. A noteworthy observation is that asymmetric flames can now also occur in $Le > 1$ mixtures. Furthermore, in addition to the $n = 0$ and $n = 1$ eigenmodes, we also observe $n = 2$ eigenmodes (red curve) when $Le > 1$ and $Pe > 0$.

Fig. 5(d), pertaining to $Da = 100$ show that the instability region possesses qualitatively similar characteristics to that in Fig. 5(c) in the $Le < 1$ domain. The instability region in the $Le > 1$ domain is however greatly enlarged and relatively more complex. In particular, the instability region moves to the left, towards smaller values of $Le$, most notably for $Pe < 0$. This makes this instability more easily accessible.

At this stage, it is worth comparing our results with those reported by Kurdyumov [8], which do not account for perturbations in the $z$-direction, and involve in fact only $k = 0$ modes and conditions corresponding to $Pe \ge 0$. It is clear from Fig. 4 that his analysis correctly predicts the type-III bifurcations identified herein since they correspond to situations where the most unstable mode occurs indeed at $k = 0$. In other words, the findings of [8] give correct predictions about the appearance of asymmetric flames, which correspond to the blue curves in Fig. 5. These findings do not provide however the boundaries characterising the onset of the cellular and oscillatory instabilities (black and red curves in Fig. 5), which are determined by dominant unstable modes with $k \neq 0$.
<figure><img src="image_9.png" /><figcaption>Fig. 5. Stability regime diagrams in the Le-Pe plane for four selected values of Da. The white regions correspond to stable, symmetric flames, whereas the shaded regions to unstable flames. The bifurcation curves are represented by solid lines, and they are labelled by the bifurcations types ($I_v$, $\Pi_s$, $\Pi_{ss}$ and $\Pi_{lv}$) described in Fig. 4. The label *n* refers to the index of the eigenfunction in the *y* direction defined in Section 5.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Reaction rate $\omega$-fields at selected times, computed for $\beta = 10$, $\alpha = 0.85$, $Da = 30$ and selected values of Pe and Le. The *x* and *y* length scale are indicated at the left bottom of each subfigure.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 7. Time evolution of the flame, represented by the iso-surface $\theta = 0.9$, for $Da = 3$, $Le = 0.7$, $Pe = -2$ (top subfigure) and $Pe = 2$ (bottom subfigure).</figcaption></figure>

## 7. Two-dimensional time-dependent numerical simulations

To complement the comprehensive findings of the linear stability analysis summarised in Fig. 5, time-dependent simulations have been performed to determine the evolution of unstable symmetric flames. Two-dimensional simulations are discussed in this section and three-dimensional ones in the following section. The two-dimensional simulations are carried out in the *x*-*y* domain $[-50, 50] \times [-1, 1]$, with

initial conditions corresponding to steady symmetric flames described in Section 3.

Computed reaction rate $\omega$-fields are shown in Fig. 6 for $Da = 30$ and selected values of Pe, Le and time $t$. The top subfigure, pertaining to $Le = 0.85$ and $Pe = 3$, illustrates the transition from symmetric to asymmetric flames. For this case, the transition results ultimately in a steadily-propagating stable asymmetric flame, whose effective burning speed is found to be higher than that of the symmetric flame [8,21]. The middle subfigure, pertaining $Le = 4$ and $Pe = 1$, also exhibits the aforementioned transition, although it leads ultimately to a quasi-periodic, rather than a steady, solution. The bottom subfigure, pertaining to $Le = 4.8$ and $Pe = 4$, shows that the flame evolution involves persistent quasi-periodic oscillations, with the solution remaining however symmetric, a behaviour which is consistent with the predictions of the linear stability analysis in Fig. 5(c) (in the region with the label “oscillatory + symmetric”).

## 8. Three-dimensional time-dependent simulations

The two-dimensional simulations discussed in the previous section are particularly significant in cases where the most unstable mode occurs at $k = 0$; such cases belong to the regions labelled “asymmetric” in Fig. 5. For situations where the most unstable mode does not occur at $k = 0$, three-dimensional simulations are important. Computations for three sample cases are presented in Figs. 7 and 8 which use
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Reaction rate $\omega$-fields at a selected value of time, $t = 400$, computed for $\beta = 10$, $\alpha = 0.85$, $Da = 30$, $Pe = 2.5$ and $Le = 5$. The three subfigures on the left correspond to sections of the domain by the planes $z = -10$, $z = 0$ and $z = 10$ whereas the right subfigure corresponds to a section by the plane $y = 0$.</figcaption></figure>

periodic boundary conditions in the z-direction. These computations are carried out using COMSOL Multiphysics, as for the two-dimensional simulations, and are performed on non-uniform tetrahedral grids.

The iso-surface $\theta = 0.9$, taken as a representative of the flame surface, is shown in Fig. 7 at selected values of time $t$. The top subfigure, corresponding to $Le = 0.7$, $Pe = -2$ and $Da = 3$, illustrates the formation of a cellular pattern in the spanwise z-direction, while remaining symmetric at all times about the $y = 0$ plane. In this case, the flame is found to evolve in time into a stable steady state. The bottom subfigure, corresponding to $Le = 0.7$, $Pe = 2$ and $Da = 3$, demonstrates the formation of a cellular pattern in the z-direction, but also the transition to asymmetric flames. In this case, however, the flame is not found to approach a steady state. The complete flame dynamics for these two cases can be better appreciated from the animated time evolution included in the supplementary material.

Computational results for $Le = 5$, $Pe = 2.5$ and $Da = 30$ are presented in Fig. 8. Shown are the reaction rate $\omega$-fields in the $y = 0$ plane (top view) and in the planes defined by $z = -10$, $z = 0$ and $z = 10$, at a particular value of time, $t = 400$. The flame evolves into a time-dependent structure involving travelling waves along the z-direction and time oscillations of a symmetric flame in planes parallel to the $x-y$ plane. Please refer to the supplementary material for a time animation of this case.

## 9. Conclusions

In this study, the diffusive-thermal instability of flames propagating in a Poiseuille flow is investigated within a three-dimensional framework. A wide range of variation of the physical parameters is considered which covers a variety of practical applications ranging from sub-millimetre micro-burners to relatively large channel burners and experiments involving flame propagation in wide or narrow channels and phenomena such as flame flashbacks.

The flame is found, depending on the parameters, to undergo either an oscillatory instability, or a cellular instability in the spanwise direction (the $z^*$-direction in Fig. 1) combined with a transition instability from symmetric to asymmetric shapes in planes parallel to the $x^*-y^*$ plane. The instabilities encountered are determined by the system size (or, the Damköhler number $Da$), the flow strength (or, the Peclet number $Pe$) and by preferential diffusion (or, the Lewis number $Le$). The results of a linear stability problem, based on the computation of its eigenvalues, are conveniently summarised in the form of stability regime diagrams drawn in the $Le-Pe$ plane for selected values of $Da$, see Fig. 5. Two- and three-dimensional time-dependent computations are also carried out, confirming the prediction of the linear stability analysis for sample values of the parameters, and describing the full evolution of unstable flames into the non-linear regime.

In narrow channels, it is found that the flames are always symmetric about the channel mid-plane and can be strongly influenced by Taylor dispersion; in particular, Taylor dispersion is found to promote the cellular instability and to inhibit the oscillatory instability. In a moderate size or large channels, on the other hand, the flame has stronger

propensity for $Le < 1$ to become asymmetric when it is opposed by the flow ($Pe > 0$), but not when it is aided by the flow ($Pe < 0$). In addition, the instability scenarios for $Le > 1$ are more complex and involve both symmetric and asymmetric flames; see e.g. Fig. 5(c) and (d).

Since the conclusions drawn above are based on the diffusive-thermal model under adiabatic conditions, it is desirable to extend this investigation to account for thermal-expansion and heat-loss effects in the future. For instance, in the presence of heat losses, the boundaries of the instability domains drawn in Fig. 5 are anticipated to shift towards $Le = 1$ [22], making the instability more readily accessible. In the presence of density or viscosity variations, in addition to the diffusive-thermal instability addressed herein, hydrodynamic instabilities such as the Darrieus–Landau and the Saffman–Taylor instabilities are also expected to contribute to the instability of flames [23–25].

## Novelty and significance statement

The novelty of this research is due to it being the first three-dimensional investigation of the diffusive-thermal instability of premixed flames propagating in a Poiseuille flow. The investigation includes a linear stability analysis as well as time-dependent numerical simulations for flames whose propagation are either aided or opposed by the flow. It is significant because it covers wide range of Damköhler numbers, applicable to flame propagation from large-scale wide to small-scale narrow channels. It is also significant because it shows that the flame instabilities are determined by a combination of the traditional diffusive-thermal instabilities and a transition instability from symmetric to asymmetric flames. It also clarifies the role played by Taylor dispersion in narrow channels.

## CRediT authorship contribution statement

**Aiden Kelly:** Performed research, Analyzed data. **Prabakaran Rajamanickam:** Performed research, Analyzed data, Wrote the paper.

**Joel Daou:** Designed research, Analyzed data, Wrote the paper. **Julien**

**R. Landel:** Analyzed data, Reviewed the paper.

## 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

This work was supported by the UK EPSRC through grant EP/V004840/1.

## Appendix A. Supplementary data

Three animated time evolution of three-dimensional numerical simulations are included as video files.

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.proci.2024.105258.
References

[1] B. Lewis, G. Von Elbe, Combustion, Flames and Explosions of Gases, Academic Press, London, UK, 1987.

[2] J. Daou, M. Matalon, Flame propagation in poiseuille flow under adiabatic conditions, Combust. Flame 124 (3) (2001) 337–349.

[3] 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.

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

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

[6] F. Veiga-López, D. Martínez-Ruiz, E. Fernández-Tarrazo, M. Sánchez-Sanz, Experimental analysis of oscillatory premixed flames in a Hele-Shaw cell propagating towards a closed end, Combust. Flame 201 (2019) 1–11.

[7] E. Al Sarraf, C. Almarcha, J. Quinard, B. Radisson, B. Denet, P. Garcia-Ybarra, Darrieus–Landau instability and Markstein numbers of premixed flames in a Hele-Shaw cell, Proc. Combust. Inst. 37 (2) (2019) 1783–1789.

[8] V.N. Kurdyumov, Lewis number effect on the propagation of premixed flames in narrow adiabatic channels: Symmetric and non-symmetric flames and their linear stability analysis, Combust. Flame 158 (7) (2011) 1307–1317.

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

[10] V.N. Kurdyumov, C. Jiménez, Propagation of symmetric and non-symmetric premixed flames in narrow channels: Influence of conductive heat-losses, Combust. Flame 161 (4) (2014) 927–936.

[11] J. Daou, A. Kelly, J. Landel, Flame stability under flow-induced anisotropic diffusion and heat loss, Combust. Flame 248 (2023) 112588.

[12] M. Sánchez-Sanz, D. Fernández-Galisteo, V.N. Kurdyumov, Effect of the equivalence ratio, Damköhler number, Lewis number and heat release on the stability of laminar premixed flames in microchannels, Combust. Flame 161 (5) (2014) 1282–1293.

[13] D. Fernández-Galisteo, C. Jiménez, M. Sánchez-Sanz, V.N. Kurdyumov, Effects of stoichiometry on premixed flames propagating in narrow channels: symmetry-breaking bifurcations, Combust. Theory Model. 21 (6) (2017) 1050–1065.

[14] D. Fernández-Galisteo, C. Jiménez, M. Sánchez-Sanz, V.N. Kurdyumov, The differential diffusion effect of the intermediate species on the stability of premixed flames propagating in microchannels, Combust. Theory Model. 18 (4–5) (2014) 582–605.

[15] V.N. Kurdyumov, C. Jiménez, Structure and stability of premixed flames propagating in narrow channels of circular cross-section: Non-axisymmetric, pulsating and rotating flames, Combust. Flame 167 (2016) 149–163.

[16] P. Rajamanickam, J. Daou, A thick reaction zone model for premixed flames in two-dimensional channels, Combust. Theory Model. 27 (4) (2023) 487–507.

[17] M.C. Cross, P.C. Hohenberg, Pattern formation outside of equilibrium, Rev. Modern Phys. 65 (3) (1993) 851.

[18] J. Daou, P. Rajamanickam, Diffusive-thermal instabilities of a planar premixed flame aligned with a shear flow, Combust. Theory Model. 28 (1) (2024) 20–35.

[19] J. Daou, P. Rajamanickam, Premixed flame stability under shear-enhanced diffusion: Effect of the flow direction, Phys. Rev. Fluids 8 (12) (2023) 123202.

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

[21] D. Rodríguez-Gutiérrez, R. Gómez-Miguel, E. Fernández-Tarrazo, M. Sánchez-Sanz, Characterization of symmetric to non-symmetric flamefront transition in slender microchannels, Proc. Combust. Inst. 39 (2) (2023) 1813–1821.

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

[23] G. Joulin, G.I. 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.

[24] Y. Han, M. Modestov, D.M. Valiev, Effect of momentum and heat losses on the hydrodynamic instability of a premixed equidiffusive flame in a Hele–Shaw cell, Phys. Fluids 33 (10) (2021).

[25] D. Fernández-Galisteo, A. Dejoan, J. Melguizo-Gavilanes, V.N. Kurdyumov, A three-dimensional study of the influence of momentum loss on hydrodynamically unstable premixed flames, Proc. Combust. Inst. 39 (2) (2023) 1545–1554.