# Study on symmetric/asymmetric hydrogen flame shapes in the thickness of a Hele-Shaw burner

Ziyin Chen, Yves Ballossier, Song Zhao, Bruno Denet, Christophe Almarcha,

Pierre Boivin

► To cite this version:

Ziyin Chen, Yves Ballossier, Song Zhao, Bruno Denet, Christophe Almarcha, et al. Study on symmetric/asymmetric hydrogen flame shapes in the thickness of a Hele-Shaw burner. Combustion and Flame, 2025, 277, pp.114208. 10.1016/j.combustflame.2025.114208. hal-05086359

HAL Id: hal-05086359

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

Submitted on 27 May 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Study on symmetric/asymmetric hydrogen flame shapes in the
thickness of a Hele-Shaw burner

Ziyin Chen<sup>a,*</sup>, Yves Ballossier<sup>c</sup>, Song Zhao<sup>a</sup>, Bruno Denet<sup>b</sup>, Christophe Almarcha<sup>b</sup>, Pierre Boivin<sup>a</sup>

<sup>a</sup>Aix Marseille Univ, CNRS, Centrale Med, M2P2, Marseille, France

<sup>b</sup>Aix Marseille Univ, CNRS, Centrale Med, IRPHE, Marseille, France

<sup>c</sup>Institut de Combustion, Aérothermique, Réactivité et Environnement, CNRS-INSIS UPR 3021,
Orléans, France

**Abstract**

Premixed flame front shape, which drastically depends on boundary conditions, is
closely related to its propagation speed. In this work, we focus on the symmetry of
steady premixed hydrogen-air flames propagating in a narrow channel, like a Hele-
Shaw burner. A wide range of equivalence ratios (0.35 – 2.0) and channel widths
(1.8 mm – 4.8 mm) are analyzed by performing detailed simulations validated by
experiments.

A multiplicity of steady flame shape is found for channel widths above a cer-
tain critical value, that is related to flame cutoff wavelength. Our numerical re-
sults successfully reproduce steady flame fronts observed in experiments. Notably,
transitions from symmetric to asymmetric with equivalence ratio are reproduced.
Additionally, an increase in channel width reduces the region of symmetric solu-
tions. Furthermore, The study explores the effects of the Darrieus-Landau instabil-
ity and thermodiffusive effects on flame shapes, presenting a stability diagram for
symmetric/asymmetric flame configurations. Throughout the study, an increase
in flame area is associated with the asymmetry level of the flame front, show-
ing a trend that first increases and then decreases with the equivalence ratio. The
global consumption rate relative to laminar flame speed decreases monotonously
with increasing equivalence ratio. It is determined by the flame area increment and
stabilizing (destabilizing) effects on convex flame fronts at Lewis number greater
(smaller) than 1. This effect is quantified and proved independent of flame sym-
metry and channel width. For very lean mixtures, the differential species diffusion
significantly strengthens the consumption rate. A prediction model is established

Email address: ziyin.chen@univ-amu.fr (Ziyin Chen)

<footnote>*Corresponding author.</footnote>
to determine the flame front length given a certain equivalence ratio and channel width.

*Keywords:* Hydrogen combustion; symmetric/asymmetric flames; flame instabilities; narrow channel

# 1. Novelty and Significance Statement

This study is, to the best of our knowledge, the first attempt to recover stable premixed hydrogen-air flame shapes observed in experiments in the thickness of a Hele-Shaw burner by performing detailed simulations. The evolution of stable flame shapes has been qualitatively and quantitatively analyzed for a wide range of equivalence ratios (0.35–2.0) and channel widths (1.8 mm–4.8 mm). The impacts of hydrodynamic instability and thermodiffusive effects are investigated, and the correlation between flame shape and consumption rate is analyzed. A prediction model is established to determine the flame front length. It is significant because the prediction of flame shape and velocity is crucial to hydrogen combustion safety and enables to mimic 3D problem by 2D simulations.

# 2. Author Contributions

* Z. Chen performed research, analyzed data and wrote the paper

* Y. Ballossier designed experimental setup and corrected the paper

* S. Zhao analyzed data and proof read

* B. Denet supervised the project and corrected the paper

* C. Almarcha supervised experiments and corrected the paper

* P. Boivin supervised the project and corrected the paper

# 3. Introduction

The reduction of the carbon footprint is driving the exploration of hydrogen technology as a sustainable energy carrier. Due to the extensive flammability, significantly higher temperature and propagation speed, the safety of hydrogen combustion in terms of storage and leakage requires special consideration. In particular,
the ability to accurately simulate flame regimes in confined environments is criti-
cal, e.g. for micro-combustion devices [1, 2], where combustion instabilities sig-
nificantly influence flame dynamics. These instabilities include: 1) the Rayleigh-
Taylor (RT) instability, which occurs between two fluids with different densities
when the lighter one resides below the heavier one [3, 4]; 2) Saffman-Taylor (ST)
instability, which occurs between two fluids with different viscosity [5]; 3) ther-
moacoustic (TA) instability induced by the interactions between flame dynamics
and acoustics [6, 7]; 4) the Darrieus-Landau (DL) [8–12] instability, which arises
from the density variation across the flame front; and 5) the thermodiffusive (TD)
instability, which arises from the competition between heat and species diffusion
across the flame front, and thus could stabilize or destabilize flames depending
on the Lewis number and the curvature of the flame front [13]. Nevertheless,
premixed flames propagating from top to bottom of a narrow gap between two
vertically-placed parallel plates are unaffected by RT instability, and ST instability
does not impact flame dynamics in the thickness between the two plates. Addition-
ally, studies have demonstrated that the thermoacoustic instability can be mitigated
or suppressed in experiences [12, 14–16]. In this study, we focus on flame stability
without thermoacoustic disturbances.

The impact of Darrieus-Landau (DL) and thermodiffusive (TD) instabilities on planar flame stability has been extensively analyzed theoretically. A linear regime has been identified with an exponential growth (decrease) of the initial harmonic perturbation imposed on the planar flame. Considering an effective Lewis number close to unity, Matalon et al. [11] proposed a dispersion relation to quantify the relative growth rate of the initial perturbation of a given wavelength; Meanwhile, the dispersion model developed by Sivashinsky et al. [8] allows for large Lewis number variation, but assumes a small density contrast across the flame front. Further numerical approaches [17, 18] have shown that these theoretical models are not capable of accurately predicting premixed hydrogen/air flames, which present a significantly lower-than-unity Lewis number (for lean mixture) and a high density jump.

In addition to theoretical and numerical investigations on the dispersion rela-
tion of flame intrinsic instabilities [8–11, 17, 18], both numerical and experimental
studies of flame propagation and stability have been carried out. In particular,
flame morphology in narrow channels can be significantly affected by flame stability [19, 20], potentially leading to flashback [21, 22]. Numerical simulations have further demonstrated the existence of multiple steady solutions for flames propagating in narrow channels: symmetric and asymmetric [23].

According to [24], symmetric flame fronts are unstable and transform into asymmetric flames in the case of lean hydrogen-air mixtures. Moreover, asymmetric flames may only exist for a certain critical channel width when one-step chemistry is assumed [25, 26]. Research in [27] also highlights the dependence of symmetry stability on equivalence ratios and Lewis numbers. Recently, experimental investigations in [28] have observed the existence of symmetric and asymmetric flame shapes of premixed hydrogen-air laminar flames within the thickness of a Hele-Shaw burner. The quasi-2D assumption, in which a stationary symmetric flame shape within the thickness is usually assumed and the momentum equation within the thickness is reduced to a relation similar to Darcy's law [29], is found to be inapplicable for asymmetric flame shapes. Hence, symmetric/asymmetric flame front shapes are a key element of premixed hydrogen-air flame propagation in confined environments. Predicting the flame front shape at different equivalence ratios and channel widths is therefore essential for the safety of hydrogen combustion. However, further studies considering detailed hydrogen-air chemistry are still needed to understand flame front shapes under these conditions. Additionally, the independence of final steady solutions of initial ignition conditions calls for further studies as only simulations with one-step chemistry have been performed [25, 30].

Furthermore, the theoretical analysis in [31] has established correlations between the maximum velocity of stationary flames and their expansion ratios. However, the specificity of lean hydrogen-air mixtures requires more focused numerical investigations. In the presence of TD instability, a higher propagation speed is expected due to a larger flame area with multiscale corrugations [32], but also to a local increase in the flame speed [17, 33]. It is worth mentioning that these studies employ periodic boundary conditions instead of confined walls, and direct numerical simulations (DNS) coupled with detailed chemistry are still needed to quantify the correlation between flame morphology and propagation speed for both lean and rich laminar hydrogen-air flames in confined environments.

The current study manages to retrieve experimentally symmetric/asymmetric
flame fronts in a Hele-Shaw burner over a wide range of equivalence ratios, from lean (0.35) to rich (2.0), and several channel widths (from 1.8 mm to 4.8 mm) by conducting a series of two-dimensional simulations with detailed chemistry and confined walls. The article is organized as follows: after the description of the cutoff wavelength, numerical configuration, and general equations in Section 4, Section 5 investigates different cutoff wavelengths/widths for different boundary conditions, as well as the dependence of final steady solutions on initial ignition conditions. Then Section 6 is entirely dedicated to qualitative and quantitative analysis of flame front shapes, focusing on their symmetry/asymmetry properties and stability. The correlation between flame shapes and flame dynamics is studied in Section 7, followed by Section 8, which establishes a model to predict flame front configurations. Conclusions are presented in the last section.

# 4. Theoretical background and general formulation

## 4.1. *Definition of cutoff wavelength*

The cutoff wavelength $\lambda_c$ represents an efficient parameter for evaluating the flame instability: It is a critical value determined by both the hydrodynamic instability and thermodiffusive effects. Below this value, the harmonic disturbances with small amplitude have a negative growth rate and then disappear. Otherwise, the amplitude increases exponentially in the linear phase and transitions into the nonlinear phase, where the flame shape no longer varies harmonically. A smaller cutoff wavelength implies a larger range of unstable modes.

The dispersion relation describes the related growth rate as a function of the perturbation wave number during the linear phase, and can determine the cutoff wavelength. In [11] a dispersion relation is derived assuming the mixture Lewis number close to 1, which is not the case for lean H₂/air mixtures; another dispersion relation allowing for small Lewis number is derived in [8], with the assumption of a small expansion ratio, which also does not hold for H₂/air mixtures. A comparison between numerical dispersion relation and these theoretical models has been done in [18]: as the two dispersion relations are based on assumptions that do not correspond to H₂/air mixtures, the predicted cutoff wavelengths differ from the numerical ones. Therefore, in the present study, we compare our results
not with the theoretical cutoff wavelength, but with the numerical results in [18].

## 4.2. Configuration for cutoff wavelength/width determination

The simulations for cutoff wavelength/width (cutoff wavelength $\lambda_c$ for periodic boundary and cutoff width $L_{y,c}$ for slip/non-slip wall boundary) determination are conducted in a two-dimensional domain of size $L_x \times L_y$ with two open ends in the x direction, and the flame propagates from left to right, as shown in Fig. 1. $L_x$ is set to be ten times $L_y$ to avoid boundary influence. Three different boundary condition in the y direction have been tested to investigate the confinement effect on flame intrinsic instability: periodic, adiabatic slip wall, and adiabatic non-slip wall. As discussed in [34] where one-step chemistry is used, the presence of slip adiabatic walls would halve the cutoff width compared with periodic boundary conditions. Additionally, this cutoff value decrease further with non-slip walls. The cutoff wavelengths/widths at different boundary conditions are noted hereafter as: $\lambda_{c,P}$, $L_{y,c,S}$, and $L_{y,c,NS}$, with subscript P, S and NS standing for periodic, slip wall and non-slip wall, respectively. Note that in the propane-air case, we have measured $\lambda_{c,P}$ in Hele-Shaw cell experiments and used it to compare to a non linear model equation, the Sivashinsky equation [35, 36].

<figure><img src="image_2.png" /><figcaption>Figure 1: Snapshot of the initial configuration for instability analysis: a sinusoidal perturbed initial planar flame propagating from left to right, with a wavelength $\lambda = L_y$.</figcaption></figure>

The cutoff wavelength/width is measured as in [37]: a series of simulations with initially planar flames perturbed by a sinusoidal perturbation are performed to see if intrinsic instabilities would be triggered by this perturbation of wavelength $\lambda$. The perturbation amplitude used is $A = 0.4\delta_T$, where the thermal flame thickness $\delta_T = (T_b - T_\infty)/\max(\nabla T)$ is calculated from the planar flame. Unless $\lambda$ is
greater than $\lambda_c$ or $L_{y,c}$, this initial perturbation will decay into a planar flame with periodic and slip boundary conditions, while the flame would exhibit a slightly curved and symmetric shape with adiabatic non-slip walls. The obtained cutoff wavelength/width has an uncertainty related to the distance between two simulated wave numbers on either side of the cutoff wavelength/width.

The evolutions to the related final steady flame fronts for three different types of boundary conditions, at $\phi = 0.8$ and for different channel widths $h$, are illustrated in Fig. 2. In all cases, the initial condition is set as the perturbed sinusoidal planar flame, with the wavelength equal to channel width and the perturbation amplitude around the thermal flame thickness. The scales of flame fronts have been adjusted for comparison purposes.

It should be noted that the definition of cutoff wavelengths/widths depends on the boundary conditions: (i) with periodic boundaries, $\lambda_{c,P}$ is the wavelength below which the disturbance is damped and a planar flame is recovered as the final solution (Fig. 2a); (ii) with slip adiabatic walls, $L_{y,c,S}$ is the channel width below which the flame becomes flat again (Fig. 2c), otherwise the flame becomes asymmetric (half a wavelength in the domain) (Fig. 2d); (iii) with non-slip walls, the $L_{y,c,NS}$ is the limiting channel width below which the flame remains symmetric, while not strictly flat as in the cases with periodic or slip walls. Instead, a slightly curved front (concave toward fresh gases because of the zero velocity on the wall) would be established (Fig. 2e). With channel width above $L_{y,c,NS}$, an asymmetric shape is obtained (Fig. 2f).

A grid size of $\Delta x = 3 \times 10^{-5}$ m is shown sufficient for $L_y > 1$ mm, otherwise the resolution is reduced to $\Delta x = 10^{-5}$ m for better accuracy in determining flame front.

### 4.3. Configuration for symmetry analysis

The simulations for the flame front symmetry analysis are conducted within the thickness between the two plates of the Hele-Shaw cell. The detailed experimental setup is described in a previous study [28] and summarized hereafter: a Hele-Shaw cell, as shown in Fig. 3, consisting of two glass plates of size 30 mm × 150 mm separated by a variable distance $h$ (ranging from 1.5 mm to 4.2 mm), is placed vertically with an open end at the top and a closed end at the bottom. Premixed
<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of flame front from the same perturbed initial conditions (the very left curve in each case) to final steady solutions (the very right curve) at different channel widths, with $\phi = 0.8$.</figcaption></figure>

hydrogen-air flames propagate from the top to the bottom. Experimental analysis is performed far from the closed end to limit its influence [38], hence the TA instability is not considered in the present study. Flame dynamics in the thickness of the Hele-Shaw cell is simulated in a two-dimensional domain of size $L \times h$, where $L = 24$ mm, with both ends open, as illustrated in Fig. 4. The Hele-Shaw glass plates are modeled as two non-slip adiabatic walls.

The comparison between full-length channel with closed end on the right hand side ($L = 150$ mm) and reduced domain with free boundary condition on the right hand side ($L = 24$ mm) is presented in Fig. 5, both with an equivalence ratio of $\phi = 0.8$ and a channel width of $h = 2.7$ mm. For the reduced length case (Fig. 5b), the contour represents the final steady solution, while for the full length case (Fig. 5a), the temporal evolution of flame front is presented with flame propagating from left to right. The scale of the flame front contours is consistent across both cases, while the axis scale is adjusted to accommodate the full-length domain for the full picture. It is noteworthy that for the full length case, the initial burnt side is extracted from the reduced case to obtain a steady shape (prior to the onset of thermoacoustic instability). The flame front maintains a steady shape before thermoacoustic instabilities, which occur only after the flame has propagated
<figure><img src="image_4.png" /><figcaption>Figure 3: Sketch of a Hele-Shaw cell [28].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Sketch of asymmetric flame propagating from left to right in the thickness of a Hele-Shaw cell with a gap size *h*.</figcaption></figure>
halfway through the channel. Figure 5c compares the steady flame front from reduced length (red dashed curve) and one instantaneous flame front from full length before thermoacoustic instability (black solid curve). The comparison reveals no significant difference between the two cases, demonstrating that the combination of the attached reference frame and boundary conditions is capable of recovering the steady flame shape in the full-length channel before thermoacoustic instability.

<figure><img src="image_6.png" /><figcaption>Figure 5: Flame front contours evolution for full length channel ($L = 150$ mm) with right closed end (a); Final steady flame front contour for reduced length ($L = 24$ mm) with right open end (b); Comparison between steady flame front contour for reduced length (red dashed curve) and one instant for full length channel before thermoacoustic instability (black solid curve).</figcaption></figure>

To evaluate the symmetry of the flames, a dimensionless asymmetry factor $S$ is defined as

$$S = \frac{1}{\frac{h}{2} \times r_f} \int_{0}^{r_f} \int_{0}^{\frac{h}{2}} \|c(x, y) - c(x, \frac{h}{2} - y)\| \, dx dy, \quad (1)$$

where $c$ is the progress variable defined as $Y_{H_2O}/Y_{H_2O,b}$, with $Y_{H_2O}$ the mass fraction of $H_2O$ and subscript $b$ standing for the burnt gases state, and $r_f$ is the distance between the tip and the tail of the flame front, as shown in Fig. 4. $S$ varies from 0 (when the flame front is perfectly symmetric) to 1. It is worth mentioning that this parameter is specifically defined for one-single cellular flame structure.

To investigate the dependency of final steady solutions on initial ignition conditions, a planar flame perpendicular to walls, a tilted planar flame, a sinusoidal perturbed planar flame, and a gaussian hot spot of burnt gases superimposed close to the domain exit, have been set as the initial condition for some cases, as illus-
trated in Fig. 6. This dependency will be further discussed in Section 5.3.

<figure><img src="image_7.png" /><figcaption>Figure 6: Four different initial conditions</figcaption></figure>

A resolution of $\Delta x = 3 \times 10^{-5}$ m is chosen to adequately resolve the internal structure of the flame, unless specified otherwise.

### 4.4. Numerical treatment and governing equations

To investigate the flame evolution over a sufficiently long time, for both instability and symmetry analysis, the simulation is carried out in the reference frame attached to the flame moving at the flame consumption rate $U_{ref} = S_c$, and $S_c$ is defined as

$$S_c = -\frac{1}{\rho_u Y_{D,u} h} \int_0^L \int_0^h \dot{\omega}_D dx dy \quad (2)$$

where $D$ stands for the deficient fuel in fresh mixtures and $\dot{\omega}_D$ is the fuel consumption rate. $\rho_u$ and $Y_{D,u}$ are the density and the deficient fuel mass fraction of the unburnt mixture.

Initially, the simulation is carried out in the laboratory reference frame, but once the flames reach halfway, the reference frame of the flame is adopted so the flame remains close to the center for the remainder of the simulation. The transient time before reaching a stable flame front shape varies with the equivalence ratio and channel width, but the simulation is carried out until convergence.

The mass, momentum, energy and species equations are solved under conservative form using the Lattice-Boltzmann solver ProLB [37, 39–41]. The descrip-
tion of Lattice-Boltzmann methods (LBM) are provided in Appendix A. For a detailed derivation of the numerical method, the reader is referred to [40]. The mass and momentum equations are then in following forms:

$$ 
\begin{aligned}
\frac{\partial \rho}{\partial t} + \nabla \cdot [\rho(\mathbf{u} - U_{ref})] &= 0 \\
\frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \cdot [\rho(\mathbf{u}(\mathbf{u} - U_{ref}))] &= -\nabla p + \nabla \cdot \Gamma
\end{aligned}
\qquad (3) $$

where $\rho$ is the density, $\mathbf{u}$ is the velocity in laboratory reference, $p$ is the pressure and $\Gamma$ is the viscous tensor. The conservation of species is considered as:

$$ \frac{\partial \rho Y_k}{\partial t} + \nabla \cdot [\rho(\mathbf{u} - U_{ref})Y_k] = \nabla \cdot (-\rho V_k Y_k) + \dot{\omega}_k \qquad (4) $$

where $Y_k$ is the mass fraction of species $k$, $V_k$ stands for its diffusion velocity, and $\dot{\omega}_k$ its reaction rate. The energy conservation is expressed in the form of total energy ($e_t$), as following:

$$ 
\begin{aligned}
\frac{\partial \rho e_t}{\partial t} + \nabla \cdot [\rho(\mathbf{u} - U_{ref})e_t] &= -\nabla q + \nabla \cdot [(\Gamma - pI)(\mathbf{u} - U_{ref})] \\
e_t &= \sum_{k=1}^{N} Y_k e_{t,k}, \quad e_{t,k} = \int_{T_0}^{T} C_{v,k}(T) dT - \frac{RT_0}{W_k} + \Delta h_{f,k}^{o} + \frac{1}{2}(\mathbf{u} - U_{ref})(\mathbf{u} - U_{ref})
\end{aligned}
\qquad (5) $$

with $I$ representing the identity matrix. $C_{v,k}$ is the species heat capacity at constant volume, the perfect gas constant $R = 8.314 \text{ J/molK}$, the reference temperature $T_0 = 300 \text{ K}$, $W_k$ is species atomic weight, and $\Delta h_{f,k}^{o}$ is the species mass formation enthalpy. For the heat energy flux $q$, we have:

$$ q = -\lambda \Delta T + \rho \sum_{k=1}^{N} h_k Y_k V_k \qquad (6) $$

where $\lambda$ is the thermal conductivity and $\lambda = \frac{\mu C_p}{Pr}$ with $\mu$ the dynamic viscosity following the temperature power-law $\mu = \mu_0(\frac{T}{T_0})^{0.7}$. Pr is the constant Prandtl number. Ideal gas equation of state $\frac{\rho RT}{W} = p$ is applied for thermodynamic closure.

The flow is considered as an ideal gas and a 12-step skeletal hydrogen-air mechanism is considered [42, 43]. Transport properties are calculated assuming
constant Prandtl number ($Pr = 0.75$) and species-specific Schmidt numbers for fresh gases ($Sc_k$, H$_2$: 0.21, H: 0.14, O$_2$: 0.8, OH: 0.53, O: 0.53, H$_2$O: 0.6, HO$_2$: 0.8, H$_2$O$_2$: 0.82, N$_2$: 1.0). These parameters are those recommended in Cerfacs' database [44] and validated for this mechanism. Species-specific Lewis number $Le_k$ is the ratio of $Pr$ to $Sc_k$. There are multiple definitions of effective Lewis number for blended fuels as depicted in [45]. Here we use the formula in [46], which is commonly applied by the community. The effective Lewis number is computed as a weighted average of the fuel and oxidizer Lewis numbers:

$$ Le_{eff} = \frac{Le_E + ALe_D}{1 + A}, \quad A = \begin{cases} 1 + \beta(\phi^{-1} - 1), & \phi \le 1 \\ 1 + \beta(\phi - 1), & \phi > 1 \end{cases} \quad (7) $$

where $Le_E$ and $Le_D$ are the excess and deficient reactants Lewis numbers in fresh gases, respectively. $\beta = \frac{E_a(T_b - T_\infty)}{RT_b^2}$ is the Zeldovich number, where $R = 8.314 \, Jmol^{-1}K^{-1}$ is the perfect gas constant, $T_b$ and $T_\infty$ are the temperature of burnt gases and fresh gases, respectively, and $E_a$ is the activation energy and can be estimated through the following formula [16, 18]

$$ \frac{E_a}{R} = -2 \frac{d\ln(\rho_u S_L)}{d(1/T_b)}. \quad (8) $$

The effective Lewis numbers in the present study are less than 1 if $\phi < 1$, otherwise $Le_{eff} > 1$, as shown in Fig. 7.

<figure><img src="image_8.png" /><figcaption>Figure 7: Effective Lewis number</figcaption></figure>
The domain is initially filled with premixed mixtures at equivalence ratios $\phi$ varying from 0.35 to 2.0, at a temperature $T_{\infty} = 300$ K and pressure $p = 1$ bar. A constant pressure outlet at 1 bar is applied on the burnt side, formulated as

$$x \rightarrow -\infty : \frac{\partial \mathbf{u}}{\partial x} = \frac{\partial T}{\partial x} = \frac{\partial Y_k}{\partial x} = 0 \qquad (9)$$

A non-reflecting boundary condition is applied on the fresh side, see Table 9.4: B2 in [47].

# 5. Preliminary study

## 5.1. Cutoff wavelengths/widths at different boundary conditions

We start by investigating the impact of different boundary conditions on cutoff wavelength/width: periodic, slip walls, and non-slip walls. These three different boundary conditions are applied to $\phi$ from 0.35 to 2.0 to obtain the related $\lambda_c/L_{y,c}$, as shown in Fig. 8. It is worth mentioning that for the slip and non-slip wall cases, the error bars are used to represent the uncertainty around 0.05 mm. However, for the case with periodic boundary condition the uncertainty is within 0.015 mm, for which the error bars are not presented.

Our results show a good agreement with the study in [18] for periodic cases. It can be observed that despite significant variations in equivalence ratios, $L_{y,c,S}$ consistently remains approximately half the value of $\lambda_{c,P}$, which aligns with the findings of [34]. This is because a domain with slip adiabatic walls is equivalent to half of a periodic domain. The $L_{y,c,NS}$ decreases further if non-slip walls are applied. For three cases, the variation of $\lambda_c/L_{y,c}$ follows the same variation in function of $\phi$, which is: from $\phi = 0.35$ to $\phi = 0.8$, $\lambda_c/L_{y,c}$ keeps decreasing, then increasing if $\phi > 0.8$. The variation of $\lambda_{c,P}$ adimensionlized by the diffusive flame thickness $\delta_D$ with equivalence ratios is shown in Fig. 9 with solid black circles. $\delta_D$ is defined as the ratio of thermal diffusivity $D_{th}$ to planar flame speed $S_L$ ($\delta_D = D_{th}/S_L$). The variation of $\delta_D$ with $\phi$ is represented by red squares. $\delta_D$ decreases first with increasing $\phi$ to $\phi = 1.4$, then increases slightly with $\phi > 1.4$. Together with the evolution of $\lambda_{c,P}$, a monotonic decrease of $\lambda_{c,P}/\delta_D$ is observed as $\phi$ decreases, which corresponds to the increasing thermodiffusive instability with
<figure><img src="image_9.png" /><figcaption>Figure 8: Cutoff wavelength/width at different boundary conditions: black circles for periodic boundary conditions, red stars for reference [18] with periodic boundaries, green circles for slip walls, blue circles for non-slip walls.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Variation of $\lambda_{c,P}/\delta_D$ (in black) and $\delta_D$ (in red) with equivalence ratio $\phi$.</figcaption></figure>
decreasing effective Lewis number.

## 5.2. Evolution towards one-hump shaped flame

For the Sivashinsky model equation, solved with zero slope at the wall in [48] which would correspond more to slip boundary conditions, different thresholds are observed when the width of the domain is increased: first only solutions with half a wavelength when the width is larger than $L_{y,c,S}$, then also solutions with a full wavelength, then asymmetric solutions with a cusp on each wall [48]. These different thresholds will be obtained here in direct numerical simulations with non-slip boundary conditions.

With non-slip walls, the multiplicity of final steady solutions has been observed as the channel width increases (typically exceeding $\lambda_{c,P}$), and shows a dependency on the form of the initial condition. Apart from the tulip (see Fig. 2e) and half-hump shapes (see Fig. 2d), which are considered symmetric/asymmetric solution in previous researches [25, 34], other distinct shapes have been identified: one-hump shaped toward fresh gases [49], which can also be symmetric (one wavelength in the domain) and asymmetric depending on equivalence ratios and channel widths. Unlike the tulip and half-hump flames, the leading point of one-hump shaped flame is not located on the wall, but between the walls.

Figure 10 illustrates the evolution towards a steady one-single hump shaped flame at $\phi = 0.43$ and $h = 4.2$ mm. During the early stage after ignition, there is a transient phase where multiple humps may appear, as shown in Fig. 10a. These humps propagate from left to right, with light and hard red lines representing the isothermal line at $1/3T_a$ and $4/5T_a$ respectively, where $T_a$ is the adiabatic flame temperature. The two contours give an estimate of the flame thickness. However, this double-hump structure is transient, despite both humps exceeding the corresponding cutoff wavelength (here $\lambda_{c,P} = 1.275$ mm). This is because the intrinsic hydrodynamic instability induces the tangential component of the velocity at the flame front and the induced velocity field ahead of the flame front will "push" the relatively smaller humps to the walls as the flame propagates forward [50], as shown in Fig. 10b and Fig. 10c. Time is adimensionlized by the characteristic time $\tau_c$, defined as $\tau_c = \delta_T/S_L$ with $S_L$ the planar flame speed. As the wall conditions are imposed, these humps disappear, leading to a steady configuration with a single
asymmetric hump, as shown in Fig. 10d. It should be noted that for the channel widths considered in the present studies, only the single hump steady-state solution exists.

<figure><img src="image_11.png" /><figcaption>Figure 10: Streamlines evolution during the transient phase after initialization, with red outlines representing the isotherm of $T = 1/3T_a$ (dark red) and $T = 4/5T_a$ (light red), where $T_a$ is the adiabatic flame temperature. H₂-air mixture at $\phi = 0.43$.</figcaption></figure>

### 5.3. Initial condition dependency

The existence of one-hump flame fronts depends on the channel width. Below a certain width $h_c$, there exists only tulip shape or half-hump. Figure 11 shows the related $h_c$ for different equivalence ratios. Likewise, the error bars stand for an uncertainty of 0.05 mm.

If the channel width is greater than these critical values, for instance, at $\phi = 0.8$ and $h = 2.7$ mm, the final solutions would be a half-hump shape if a planar flame is set as the initial condition (Fig. 11a); a one-hump front toward fresh gases would be the final steady solution if other types of initial conditions are applied: a tilted planar flame (Fig. 11b), a sinusoidally perturbed planar flame (Fig. 11c), and a gaussian hot spot (Fig. 11d). Figures 11e and 11f ($h = 0.3$ mm and $h = 0.4$ mm respectively) show cases with channel widths smaller than the critical value, as
<figure><img src="image_12.png" /><figcaption>Figure 11: Critical width $h_c$ for one-hump flame front shape with non-slip walls. Sub-figures (a)-(f) are illustrations of possible steady solutions at $\phi = 0.8$: (a)-(d) at $h > h_c$, and (e)-(f) at $h < h_c$.</figcaption></figure>

discussed in Section 5.1.

Overall, the results reveal the existence of four different steady flame front solutions: 1) symmetric tulip shape (Fig. 11e); 2) asymmetric half-hump shape (Fig. 11a and f); 3) asymmetric one-hump shape (Fig. 11b-d); 4) symmetric one-hump shape (not shown in the figure). Solution 1) exists only with $h$ below the cutoff width $L_{y,c,NS}$; solution 2) exists for $h$ above $L_{y,c,NS}$. The critical value of $h_c$ determines whether the steady one-hump flame front solutions, e.g. 3) and 4), exist.

For the following studies, the channel width is always greater than the critical width for each equivalence ratio. Furthermore, to ensure the one-hump shape, we use only the center-shifted gaussian hot spot as the initial condition, not only because it is commonly used in numerical simulations but also because it is similar to the initial condition from electrodes in experiments. We set four cases (C1 to C4) to confirm the independence of the steady solution from the position of hot spots, with $\phi = 0.8$, $h = 2.7$ mm and radius $r_c = 0.5$ mm:
C1: $y_c = 0.0$ mm

C2: $y_c = 0.55$ mm

C3: $y_c = 1.15$ mm

C4: $y_c = 1.35$ mm

where $y_c$ is the hot spot center position in the y direction. The variation in the hot spot position yields negligible differences in the final flame fronts. The relative differences compared to case C3 are defined as $\epsilon = (f_{Ci}(y) - f_{C3}(y))/\max(f_{C3})$ with $f_{Ci}(y)$ standing for flame front configuration of case Ci and i varying from 1 to 4, and all $\epsilon$ are lower than 6.5%. Moreover, a combination of two hot spots (case C1C1) is initially set, with one center located at $y = 0$ mm and another at $y = 2.7$ mm. The final solution also converges to the same shape. Hereafter, a center-shifted by 0.2 mm gaussian hot spot is used as the initial condition to ensure a faster transition to asymmetric flames when required and the top and bottom walls are all set as non-slip adiabatic walls.

# 6. Results and discussion

## 6.1. Asymmetry factor

Simulations are first carried out for equivalence ratios from 0.35 to 2.0 and channel widths of 2.7 mm, 3.2 mm and 4.2 mm. The steady-state flame front solutions are all one-hump shaped given channel widths greater than the critical width for the one-hump shape, as shown in Fig. 12, which qualitatively compares with experimentally observed flame fronts (numerical results are superimposed as red lines). It should be noted that we are comparing the leading tip observed in the experiments, and the impact of the transverse curvature [51] is not considered in our simulations. According to the integrated information collected by the schlieren diagnostic, which shows different cell tips visible in the front view corresponding to different lateral positions, the shape of the front in the thickness is approximately the same in most cases [28]. Numerical flame fronts are detected by the maximum gradient of density $\rho$ (corresponding to the darkest area in the experimental Schlieren images). All flames are one-hump shaped towards fresh gases as
illustrated in Fig. 4. An excellent qualitative agreement is obtained between exper-
imental and numerical fronts for the considered range of heights and equivalence
ratios. Discrepancies are observed nonetheless for a few configurations, which
could relate to neglecting heat losses at the wall and/or gravity effects (stabilizing
for downward propagating flames) in the simulation. Nevertheless, the agreement
between 2D numerical simulations and experiments in Hele-Shaw cells is good.
In large gaps, asymmetric flames are observed over a wide range of equivalence
ratios, while this range narrows when decreasing the gap size.

<figure><img src="image_13.png" /><figcaption>Figure 12: Flame fronts between walls with red outlines representing numerical results. The post-processed experimental flame front is highlighted in black lines.</figcaption></figure>

For a quantitative analysis, the variation of the asymmetry factor *S* with *φ* for different *h* is visualized in Fig. 13a, where solid circles represent numerical results while crosses are from experiments. Despite experimental uncertainties due to the front detection algorithm, the quantitative results confirm the trend observed in Fig. 12 for the evolution of *S* as a function of *φ* and *h*. Figure 13b shows the corresponding flame front shapes extracted from the simulations (also detected as the maximum of density gradient), with both coordinates normalized by the gap size *h*. Interestingly, all symmetric flame shapes (*S* = 0) on the rich side are identical, regardless of *φ* and *h*.

As the equivalence ratio increases from 0.35, the asymmetry factor increases and then decreases, with a maximum around $\phi = 0.6$. In addition, $S$ increases with the channel width for a given equivalence ratio. Therefore, the range of equivalence ratio with asymmetric shape narrows when the channel width decreases. The turn-
<figure><img src="image_14.png" /><figcaption>Figure 13: Evolution of (a) asymmetry factor *S* and (b) flame fronts with equivalence ratio at different channel widths. In (a) the circles are from the simulation and the crosses are from the experiments (uncertainty ±0.1); in both (a) and (b) the red, black and blue colors represent at *h* = 2.7 mm, 3.2 mm and 4.2 mm respectively.</figcaption></figure>
ing point from asymmetry to symmetry for $h = 2.7$ mm is between $\phi = 1.0$ and
1.1, while for $h = 3.2$ mm and $h = 4.2$ mm the flames become symmetric with
$\phi > 1.2$ and $\phi > 1.4$ respectively.

The decreasing $S$ with increasing $\phi$ when $\phi > 0.6$ is related to a higher effec-
tive Lewis number, as observed in previous studies [21]. On the other hand, the
decrease of $S$ for $\phi \le 0.6$ implies the possibility of having symmetrical flames
for very lean hydrogen mixtures despite the low Lewis number. For example, with
$\phi = 0.35$ and channel widths at 2.7 mm and 3.2 mm, the flame shapes become
symmetric. This unexpected behavior is discussed in Section 6.2.

6.2. Flame symmetry prediction

To generalize the investigation, the range of $h$ is extended to 1.8 mm for the
lean side and 4.8 mm for the rich side, exceeding the critical width required for
a stable one-hump flame shape. Figure 14 provides a complete map showing the
flame asymmetry factor $S$ at different equivalence ratios and channel widths, with
circles for the symmetric one-hump shape, and crosses for asymmetric ones with
colors indicating values of the asymmetry factor (referring to the color map for
$S$). The critical channel widths $h_c$ are represented by black circles with error bars.
For the region above these critical $h_c$ (where one-hump solutions exist), within the
equivalence ratio range of this study, there is a minimum channel width (1.8 mm)
below which only one-hump symmetric flames exist. It is worth mentioning that
for this channel width, while the flame can be well resolved with $\Delta x = 3 \times 10^{-5}$ m,
a finer mesh size at $\Delta x = 2 \times 10^{-5}$ m is applied for a better accuracy in determin-
ing the flame front configuration and flame length calculation in post-processing.
Moreover, 4.8 mm is the maximum value above which symmetric solutions no
longer exist. Meanwhile, cases around bifurcation of symmetry and asymmetry
require a more asymmetric initial condition to trigger steady asymmetric solutions,
given the relatively strong stability of symmetric flames. These special cases are
represented by the superposition of circles and crosses. For example, at $\phi = 0.9$
and $h = 2.4$ mm, with an off-centered gaussian hot spot as the initial condition,
the final steady solution is symmetric. However, if a tilted planar flame is used for
initialization, the final steady solution is asymmetric.
<figure><img src="image_15.png" /><figcaption>Figure 14: Distribution of flame asymmetry factor in function of equivalence ratios and channel widths, with circles for symmetric one-hump shape and crosses for asymmetric ones with colors standing for values of asymmetry factor. Black circles with error bars represent the critical width for one-hump shape $h_c$.</figcaption></figure>

Figure. 14 reveals general variations of flame shape with equivalence ratio and channel width, which, in summary, are: i) as the equivalence ratio increases for $\phi > 0.6$ (decrease for $\phi < 0.6$), symmetric flames are more likely; and ii) for a given equivalence ratio, asymmetric flames could appear as the channel width increases. This trend suggests the existence of a possible criterion for determining flame symmetry. This section is devoted to the analysis of the factors controlling the flame shape stability in order to predict the flame symmetry.

* Cutoff wavelength

As shown in Fig. 8, $\lambda_{c,P}$ first decreases and then increases with increasing equivalence ratios. The turning point is around $\phi = 0.8$. On the other hand, the first increasing and then decreasing values of $S$ suggest that asymmetric flames are more likely to occur for small $\lambda_{c,P}$.

* Diffusive flame thickness

Previous research by [52] investigated the stability of the convex flame front to small-scale disturbances. A characteristic length scale is the flame thickness. These perturbations along the flame fronts move with the flame toward
the fresh gases and are simultaneously directed towards the sidewalls due to the tangential component of the velocity at the flame fronts, as shown in Fig. 10. In addition, the perturbation wavelength continues to increase because of the stretching effect on the flame front. Overall, the stability of the one-hump flame front depends on the interactions between two mechanisms: i) the stretching of the disturbances by the flame front, which increases the wavelength of the disturbances; ii) the hydrodynamic induced velocity field, which guides the disturbances to the wall. A one-hump flame front can remain stable as long as the disturbances are either directly dumped onto the flame front or directed to the wall faster than the time required for the disturbances to grow [50, 52]. The diffusive flame thickness $\delta_D$ is chosen to synthesize these mechanisms, as it incorporates both the intrinsic disturbances and the velocity.

The channel width $h$ is used to non-dimensionalize the two aforementioned parameters. Figure 15 provides a stability diagram of symmetric/asymmetric flame fronts as a function of $h/\lambda_{c,P}$, and Peclet number $Pe = h/\delta_D$.

Different symbol types represent different $h$ values, with colors corresponding to the value of asymmetric factor (referring to color map), and the increasing sizes represent increasing equivalence ratios from 0.35 to 2.0. Increasing $h/\lambda_{c,P}$ and decreasing $Pe$ both increase the value of the asymmetry factor, which corresponds to the two mechanisms depicted above: at a given channel width, a higher $h/\lambda_{c,P}$ implies a smaller cutoff wavelength allowing perturbation growth, and a lower $Pe$ signifies a smaller induced flow velocity to guide disturbances to the wall. For a given channel width $h$, on the lean side with $\phi \le 0.8$, both $\lambda_{c,P}$ and $\delta_D$ decrease monotonously with the equivalence ratio, hence both $h/\lambda_{c,P}$ and $Pe$ increase with the symbol sizes, e.g. $\phi$, and the asymmetry factor increases till the highest point for each symbol set; on the other hand, $\lambda_{c,P}$ and $\delta_D$ increase with $\phi$ on the rich side with $\phi > 1.4$, hence both $h/\lambda_{c,P}$ and $Pe$ decreases with symbol size increasing, and flame shape converts into symmetric. For $\phi$ between 0.8 and 1.4, $\lambda_{c,P}$ increases while $\delta_D$ decreases, giving the flex of symbols with increasing $\phi$.

A dashed black line is assigned to the bifurcation limit between symmetry and asymmetry, and special cases around bifurcation showing dependency on the initial
<figure><img src="image_16.png" /><figcaption>Figure 15: Stability diagram for symmetric/asymmetric flame front based on $h/\lambda_{c,P}$ and $Pe$. Symbol types stand for different channel widths, and symbol sizes indicate equivalence ratio values, with colors standing for asymmetry factor values. Black dashed line represents the bifurcation of symmetric/asymmetric shapes. Symbols superposed with blue circles are cases whose final steady flame configuration depends on initial condition's asymmetry level.</figcaption></figure>
condition asymmetry level (superposition of circles and crosses) are represented by the superposition of blue circles and symbols. In addition, the bifurcation limit line with a slope much smaller than 1 indicates that at a certain $\phi$, the increase of channel width would lead to a more asymmetric shape due to a longer pathway allowing perturbation growth.

## 7. Flame shape and consumption rate

Variations in flame shape caused by asymmetries lead to strong variations in flame consumption rate $S_c$, through both changes in flame front length and the related Lewis effect. Hence, we define the normalized flame front length as $l_f/h$, the ratio of the flame front length $l_f$ to the gap size. Finally, to quantify the related Lewis effect on flame curvature, let us define the amplification factor $I_0$ as the ratio of normalized consumption rate to flame front length

$$I_0 = \frac{S_c h}{S_L l_f}, \qquad (10)$$

following [13]. $I_0$ distinguishes whether the increase in consumption speed originates from the increase in the flame surface area or from the varying reactivity along the flame front. Near the flammability limit, the definition of $I_0$ is controversial: planar flames do not exist for $\phi < 0.26$ ($S_L = 0$). However, researches in [19, 53] have unveiled that flames can be sustained in the form of stable circular fronts, double-cell fronts, and multiplicity of the two shapes, as well as splitting fronts. These flame structures exhibit different propagation speeds, and the increment of flame surface is also hard to define as flame fronts are isolated. Nevertheless, the definition of $I_0$ in Eq. 10 is valid for the scope of the present study.

The dependence of the normalized consumption rate, flame front length, and amplification factors is presented in Fig. 16, and discussed hereafter.

The normalized consumption rate $S_c/S_L$ decreases monotonically with $\phi$, for all gap sizes. The normalized flame front length initially increases then decreases with $\phi$, with a maximum in the asymmetric region (see Fig. 13). Combining both, it is observed that the amplification factor $I_0$ is monotonically decreasing with $\phi$. More interestingly, $I_0$ appears to be independent of the gap size and the flame
<figure><img src="image_17.png" /><figcaption>Figure 16: (a) Ratio of flame front length to channel width in black circles; ratio of flame consumption rate to planar flame speed in red triangles; (dotted lines for $h = 2.7$ mm; dash-dotted lines for $h = 3.2$ mm; dashed lines for $h = 4.2$ mm; blue crystals from [31]); (b) Amplification factor $I_0$ in all cases with crosses standing for asymmetric flame shapes and circles for symmetric ones.</figcaption></figure>
symmetry factor S (Fig. 16b). This was confirmed beyond the three gap sizes represented in Figs. 12, 13, and 16a by a number of additional simulations, as indicated in the subplot of Fig. 16b. There, each simulation is represented by a circle (symmetric flame), or cross (asymmetric flame), colored by the asymmetry factor S. When projected in the $(\phi, I_0)$, all points collapse very close to a single amplification factor curve, regardless of h and S. This universality suggests that the amplification factor can well quantify the Lewis effect. This trend aligns well with the experimental results [28], although the latter exhibit scattering due to assumptions, heat losses, and variations in transverse curvature. Nevertheless, it is mentionable that this universality of amplification factor could not be generalized to all gap widths.

These trends are further discussed into three sub-regimes: i) symmetric rich flames, ii) asymmetric flames, and iii) symmetric lean flames.

## 7.1. Symmetric rich flames

For large equivalence ratios, the effective Lewis number is close to or larger than unity. As pointed out when discussing Fig. 13b, the normalized flame shape is constant at $l_f/h \approx 1.22$. The normalized consumption rate is slightly smaller and decreases with increasing $\phi$, resulting in $I_0 \lesssim 1$.

It ensues that the consumption rate is not very sensitive to the gap size, with relative differences within 5%, as seen in the subplot of Fig. 16a. Finally, it is worth noting that the normalized consumption rate agrees reasonably well with the formula proposed by [31]

$$S_c/S_L \approx 1 + \frac{(\sigma^2 - 1)^2}{32\sigma^3}, \qquad (11)$$

with $\sigma$ standing for the expansion ratio defined as the ratio of the fresh gases density to the burnt gases density, and obtained assuming the flame velocity increase to be only due to the front curvature. This agreement is depicted in Fig. 16a using blue symbols.
## 7.2. Asymmetric flames

When decreasing the equivalence ratio from the rich symmetric flame regime the amplification factor starts to increase due to the decrease in the effective Lewis number. Equivalently, the consumption rate $S_c/S_L$ increases faster than the normalized flame front length.

This effect is independent of $S$, as shown in Fig. 16b. But as $Le_{eff}$ decreases, asymmetric flames solutions appear [23], and $S$ increases, leading to a significant increase in both normalized flame front length and consumption rate.

Asymmetric flames appear at a critical equivalence ratio, which increases with the gap size (subplot of Fig. 16b). This transition is characterized by a sudden increase in flame area and consumption rate. $S$ then keeps increasing with lower equivalence ratio until its maximum.

## 7.3. Lean symmetric flames

From its maximum value around $\phi = 0.6$, the symmetry factor decreases with $\phi$. However, the amplification factor $I_0$ increases, indicating an important influence of the thermodiffusive effect.

Figure 17a illustrates the non-unity Lewis effect on the flame for $\phi = 0.35$ and $h = 2.7$ mm. Contours of an atom-based equivalence ratio

$$ \phi_{mol} = \frac{Y_H/Y_O}{(Y_H/Y_O)_{st}}, $$

are represented. The formula is identical to the classical equivalence ratio definition in the fresh gases, but would remain constant throughout a 1D flame due to atom conservation. Values over 0.35 therefore indicate an excess of fuel compared to the fresh gases, and an associated faster reaction rate. This is also illustrated by isolines (in purple) of heat release rate (HRR) normalized by the maximum HRR of the corresponding 1D planar flame, clearly indicating enhanced reaction in the center region. The yellow line indicates the flame as $T = 1000$ K.

At $\phi = 0.35$, $Le_{eff} = 0.3$, fresh reactants diffuse faster than heat, hence the local reaction would be strengthened (weakened) by convex (concave) curvature. The central part of the flame front has a strong convex curvature, and the differ-
<figure><img src="image_18.png" /><figcaption>Figure 17: At $\phi = 0.35$, $h = 2.7$ mm, (a) the distribution of normalized heat release rate $HRR_{norm}$ (red contours), local equivalence ratio $\phi$ (blue color map), and isothermal line at $T = 1000$ K (yellow line); (b) the distribution of $HRR_{norm}$ (black line) and normalized mass fraction of hydrogen $Y_{H_2, norm}$ (red line) along the isothermal line $T = 1000$ K.</figcaption></figure>
ential species diffusion leads to a higher local molecular equivalence ratio up to 0.42. Therefore, the $HRR_{norm}$ increases significantly in this region compared to the concave regions near the walls where $\phi_{mol}$ is around 0.19 and flame quenching occurs.

Figure 17b shows the evolution of $HRR_{norm}$ and the fuel mass fraction normalized by its fresh gas value, along the 1000K isoline identified in Fig. 17a.

Similar to the observations in [33], due to the highly concave curvature near the walls, the depletion of hydrogen leads to a much lower heat release rate compared to the convex part. Overall, the intense reaction in the convex region leads to a much higher global consumption rate.

## 8. Flame front shape parametrization

Figure 16 reveals that the amplification factor does not depend on flame shapes or channel widths, but on the equivalence ratio. To predict the flame consumption rate, the missing puzzle piece remains the flame front length at a given equivalence ratio and channel width. Rich premixed hydrogen-air flame front shapes perfectly overlap when they are symmetric and one-hump toward fresh gases, as shown in Fig. 13. The asymmetric flames also exhibit similar profile despite the difference in amplitude. This section is dedicated to establishing a model to reconstruct flame front configurations at given channel widths and equivalence ratios. The cases presented in Fig. 13 are used for studies in Sections 8.1 and 8.2.

### 8.1. Polynomial fitting

For the sake of parametrization, the y-axis (transversal direction along the channel width) is normalized as $\tilde{y} = 2 \times (y/h - 0.5)$, so that $\tilde{y}$ ranges from -1 to 1 for all channel widths. The flame front configuration $f(y)$ (longitudinal direction along the channel length) is also normalized as $F(\tilde{y}) = (f(y) - f(\frac{h}{2})) / h$ for general comparison between cases. According to the expression of flame front configuration in [54], fourth-order polynomials $P(\tilde{y})$ are proposed to fit both sym-
metric and asymmetric flame shapes. These polynomials are of form

$$P_{j,h}(\tilde{y}) = \sum_{i=1}^{i=4} a_{i,j,h} \tilde{y}^i \quad (12)$$

with j standing for equivalence ratio varying from 0.35 to 2, and h representing channel width. a₀ = 0 for all j and h given all flame fronts pass by (0, 0) after normalization. The fitting results prove that flame fronts can be perfectly recovered by fourth order polynomials, as shown in Fig. 18.

<figure><img src="image_19.png" /><figcaption>Figure 18: Flame front shapes with fitting curves, with black hollow circles standing for numerical results and solid curves for fitting results.</figcaption></figure>

Given the perfect overlap of symmetric flame front shapes, the even and uneven parts in this polynomial are separated into:

$$P_{j,h,even}(\tilde{y}) = \sum_{i=2,4} a_{i,j,h} \tilde{y}^i, \quad P_{j,h,uneven}(\tilde{y}) = \sum_{i=1,3} a_{i,j,h} \tilde{y}^i \quad (13)$$

responsible for even and uneven parts, respectively. It is worth noting that for symmetric cases, $P_{uneven}(\tilde{y}) = 0$.
After regularization of both $P_{even}$ and $P_{uneven}$ by their respective amplitudes, the configurations of $\tilde{P}_{even}$ and $\tilde{P}_{uneven}$ are shown in Fig. 19. It should be mentioned that very lean cases with $\phi < 0.5$ are not shown in the figure, as the flames in these cases exhibit significant weaker reactions near the walls, leading to differences compared to other cases. Despite the variation in channel width and equivalence ratio, the regularized even and uneven parts show perfect overlap. The even part can be described by one fourth order polynomial

$$ \tilde{P}_{even}(\tilde{y}) = \sum_{i=2,4} \tilde{a}_i \tilde{y}^i, \qquad (14) $$

which no longer depends on channel width and equivalence ratio. The uneven part is represented by one third order polynomial

$$ \tilde{P}_{uneven}(\tilde{y}) = \sum_{i=1,3} \tilde{a}_i \tilde{y}^i, \qquad (15) $$

likewise, independent of $h$ and $\phi$.

<figure><img src="image_20.png" /><figcaption>Figure 19: Left: $\tilde{P}_{even}$ of flame fronts at $\phi \ge 0.5$; right: $\tilde{P}_{uneven}$ for asymmetric flame fronts at $\phi \ge 0.5$.</figcaption></figure>

## 8.2. Amplitudes quantification

To rebuild the flame front, the last missing key here is the amplitudes of $P_{even}$ and $P_{uneven}$. Hereafter, we use $A$ for this amplitude, $A_{even} = \|\min(P_{even})\|$ and $A_{uneven} = \|\max(P_{uneven}) - \min(P_{uneven})\|$.

* Symmetric cases
For symmetric flames, $A_{even}$ remains approximately constant regardless of channel width and equivalence ratio, as shown in Fig. 20. This is because the normalized symmetric flame fronts already show the same amplitude, as demonstrated in Figs. 13b and 18.

<figure><img src="image_21.png" /><figcaption>Figure 20: $A_{even}$ (crosses) and $A_{uneven}$ (circles) for symmetric flames at $\phi \ge 0.5$. Colors of symbols stand for different channel widths, with blue for $h = 4.2$ mm, black for $h = 3.2$ mm, and red for $h = 2.7$ mm.</figcaption></figure>

### Asymmetric cases

Given that the asymmetry property depends on $A_{uneven}$ and considering the dependency of $S$ on $h/\lambda_{c,P}$ (adimensional parameter used in Section 6.2), here we associate the amplitude $A_{uneven}$ with $h/\lambda_{c,P}$ to study the correlation. Relations

$$g_1(\phi) = A_{even} - \frac{A_{uneven}}{2}, \quad g_2(\phi) = \frac{A_{uneven}}{(h/\lambda_{c,P})^2} \quad (16)$$

are proposed to make the dependency on channel width implicit, as shown in Fig. 21, with the data represented by crosses and hollow circles, respectively. For asymmetric flames, both $g_1$ and $g_2$ can be well recovered by equations of the same form $Fit_i = m_i e^{-t_i \phi^2} + b_i$ ($i = 1, 2$), respectively, both with coefficient of determination greater than 0.997.

The exact values of these constants $\tilde{a}_{1-4}$, $m_i$, $t_i$, and $b_i$ in fitting equations $Fit_i = m_i e^{-t_i \phi^2} + b_i$ ($i = 1, 2$) are given in table 1.
<figure><img src="image_22.png" /><figcaption>Figure 21: Relations $g_1$ (crosses) and $g_2$ (circles) to relate $A_{even}$ and $A_{uneven}$ for asymmetric flames, with $g_1$ and $g_2$ defined in Equs. 16. Colors of symbols stand for different channel widths, with blue for $h = 4.2$ mm, black for $h = 3.2$ mm, and red for $h = 2.7$ mm. Fitting equations $Fit_1$ and $Fit_2$ are of form $me^{-t\phi^2} + b$, with different constants $m$, $t$ and $b$.</figcaption></figure>

<table><caption>Table 1: Constants of fitting equations $Fit_1$ and $Fit_2$.</caption><thead><tr><th>$\tilde{a}_1$</th><th>$\tilde{a}_2$</th><th>$\tilde{a}_3$</th><th>$\tilde{a}_4$</th><th>$m_1$</th><th>$m_2$</th><th>$t_1$</th><th>$t_2$</th><th>$b_1$</th><th>$b_1$</th></tr></thead><tbody><tr><td>0.378</td><td>-0.480</td><td>0.122</td><td>-0.520</td><td>0.206</td><td>0.047</td><td>5.895</td><td>3.147</td><td>0.211</td><td>0.011</td></tr></tbody></table>

Overall, the correlation between $A$ and $h$ is made implicit through the relations $g_1$ and $g_2$, and the correlation with $\phi$ can be depicted by fitting equations only as function of $\phi$. Once $\phi$ and $h$ are given, $A$ can also be determined.

## 8.3. Validation

With the decomposition of flame front configuration into even and uneven parts, and the quantification on their forms and amplitudes, the flame front can be re-established given a certain equivalence ratio and channel width. The process is visualized by the flow chart in Fig. 22.

To validate this model, flame front configurations for a large range of $h$ (cases shown in Section 6.2) are predicted and compared to numerical simulations. As shown in Fig. 23, the flame front length $l_f$ and the ratio to channel width $l_f/h$ show a satisfactory agreement between prediction and simulation. The relative discrepancy is less than 5%.
<figure><img src="image_23.png" /><figcaption>Figure 22: Flow chart to rebuild flame front.</figcaption></figure>

## 9. Conclusions

In this work, the variation of flame front shapes within narrow channels of a Hele-Shaw burner is qualitatively and quantitatively analyzed over a wide range of equivalence ratios and channel widths by performing detailed simulations and comparing with experiments. The cutoff wavelengths/widths for symmetric/asymmetric flames with different boundary conditions have been determined, as well as the critical widths for one-hump flame fronts, above which exists the multiplicity of steady solutions. The numerical results successfully reproduce flame fronts observed from experimental, which consistently exhibit a one-hump shape toward fresh gases throughout this work when the channel widths exceed the critical widths for the one-hump shape. As the channel width increases, the range of symmetric solutions progressively narrows. For a given channel width, flame fronts tend to become symmetric for equivalence ratios above $\phi = 0.6$. Eventually, flames with rich mixtures exhibit a self-similar shape.

The variation in flame front length increment is associated with the level of asymmetry, but the consumption rate also depends on an amplification factor, and we demonstrated that this amplification factor is solely dependent on the equiva-
<figure><img src="image_24.png" /></figure>

<table><caption>Figure 23: Comparison of flame front length $l_f$ (top left), ratio of $l_f$ to channel width $h$ (top right), both with blue color standing for predictions from parametrization model and red for results from numerical simulations. The relative error of prediction (bottom). Symbol types in three figures stand for channel widtds.</caption><tr><td>○</td><td>h = 4.8 mm</td><td>▷</td><td>h = 3.6 mm</td><td>◁</td><td>h = 3.1 mm</td><td>◇</td><td>h = 2.7 mm</td><td>◇</td><td>h = 2.3 mm</td></tr><tr><td>▽</td><td>h = 4.5 mm</td><td>○</td><td>h = 3.4 mm</td><td>○</td><td>h = 3.0 mm</td><td>◇</td><td>h = 2.6 mm</td><td>▽</td><td>h = 2.1 mm</td></tr><tr><td>△</td><td>h = 4.2 mm</td><td>□</td><td>h = 3.3 mm</td><td>○</td><td>h = 2.9 mm</td><td>★</td><td>h = 2.5 mm</td><td>△</td><td>h = 1.8 mm</td></tr><tr><td>◁</td><td>h = 3.7 mm</td><td>○</td><td>h = 3.2 mm</td><td>◇</td><td>h = 2.8 mm</td><td>◎</td><td>h = 2.4 mm</td><td></td><td></td></tr></table>
lence ratio for all gap widths investigated. Rich mixtures with $\phi > 1.5$ exhibit an amplification factor $I_0 < 1$, resulting in a global consumption rate increment smaller than the flame front length increment. Conversely, with for lean mixtures with $\phi < 0.5$, the global consumption rate continuously increases with decreasing effective Lewis number, despite the decrease in flame front length increment, due to $I_0 > 1$. For cases with $Le$ close to 1, the flame speed is entirely driven by the flame front length increment, as $I_0 = 1$. The significant increase in flame speed observed for very lean mixtures is of particular concern due to the extraordinarily enhanced differential species diffusion.

Flame front configurations are perfectly parameterized and recovered using fourth-order polynomials. A prediction model based on the parametrization has been developed and validated to determine the flame front length at a given equivalence ratio and channel width. In future work, the combination of the amplification factor and the prediction model for flame front length could be applied to correct flame propagation velocity and mimic 3D problems by 2D simulations.

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

# 11. Acknowledgments

This research was supported by IMI Institut Mécanique et Ingénierie funding and the ECOSAFE ANR project ANR-21-CE05-0028. Centre de Calcul Intensif d'Aix-Marseille and GENCI-TGCC/CINES (Grant 2023-A0152B11951) are acknowledged for granting access to their high-performance computing resources. Carmen Jiménez is gratefully acknowledged for fruitful discussion.

# References

[1] Y. Ju, K. Maruta, Microscale combustion: Technology development and fundamental research, Prog. Energy Combust. Sci. 37 (2011) 669–715.
[2] N. S. Kaisare, D. G. Vlachos, A review on microcombustion: Fundamentals, devices and applications, Prog. Energy Combust. Sci. 38 (2012) 321–359.

[3] D. H. Sharp, An overview of rayleigh-taylor instability, Physica D 12 (1984) 3–18.

[4] P. M. J. TREVELYAN, C. ALMARCHA, A. DE WIT, Buoyancy-driven instabilities of miscible two-layer stratifications in porous media and hele-shaw cells, J. Fluid Mech. 670 (2011) 38–65.

[5] H. Thomé, M. Rabaud, V. Hakim, Y. Couder, The saffman-taylor instability: From the linear to the circular geometry, Phys. Fluids A. 1 (1989) 224–240.

[6] G. Searby, Acoustic instability in premixed flames, Combust Sci Technol 81 (1992) 221–231.

[7] A. Petchenko, V. Bychkov, V. Akkerman, L.-E. Eriksson, Flame-sound interaction in tubes with nonslip walls, Combust. Flame 149 (2007) 418–434.

[8] G. I. Sivashinsky, Diffusional-thermal theory of cellular flames, Combust. Sci. Technol. 15 (1977) 137–145.

[9] B. DENET, On non-linear instabilities of cellular premixed flames, Combust. Sci. Technol. 92 (1993) 123–144.

[10] B. Denet, P. Haldenwang, A numerical study of premixed flames darrieus-landau instability, Combust. Sci. Technol. 104 (1995) 143–167.

[11] M. Matalon, C. Cui, J. K. Bechtold, Hydrodynamic theory of premixed flames: effects of stoichiometry, variable transport coefficients and arbitrary reaction orders, J. Fluid Mech. 487 (2003) 179–210.

[12] C. Almarcha, B. Denet, J. Quinard, Premixed flames propagating freely in tubes, Combust. Flame 162 (2015) 1225–1233.

[13] L. Berger, K. Kleinheinz, A. Attili, H. Pitsch, Characteristic patterns of thermodiffusively unstable premixed lean hydrogen flames, Proc. Combust. Inst. 37 (2019) 1879–1886.

[14] B. Radisson, J. Piketty-Moine, C. Almarcha, Coupling of vibro-acoustic waves with premixed flame, Phys. Rev. Fluids. 4 (2019) 121201.

[15] M. Rubio-Rubio, F. Veiga-López, D. Martínez-Ruiz, E. Fernández-Tarrazo, M. Sánchez-Sanz, Suppression of thermoacoustic instabilities by flame-structure interaction, Proc. Combust. Inst. 39 (2023) 1577–1585.

[16] J. Beeckmann, R. Hesse, S. Kruse, A. Berens, N. Peters, H. Pitsch, M. Matalon, Propagation speed and stability of spherically expanding hydrogen/air flames: Experimental study and asymptotics, Proc. Combust. Inst. 36 (2017) 1531–1538.

[17] C. E. Frouzakis, N. Fogla, A. G. Tomboulides, C. Altantzis, M. Matalon, Numerical study of unstable hydrogen/air flames: Shape and propagation speed, Proc. Combust. Inst. 35 (2015) 1087–1095.

[18] L. Berger, A. Attili, H. Pitsch, Intrinsic instabilities in premixed hydrogen flames: Parametric variation of pressure, equivalence ratio, and temperature. part 1 - dispersion relations in the linear regime, Combust. Flame 240 (2022) 111935.

[19] F. Veiga-López, M. Kuznetsov, D. Martínez-Ruiz, E. Fernández-Tarrazo, J. Grune, M. Sánchez-Sanz, Unexpected propagation of ultra-lean hydrogen flames in narrow gaps, Phys. Rev. Lett. 124 (2020) 174501.

[20] M. Kuznetsov, J. Grune, Experiments on combustion regimes for hydrogen/air mixtures in a thin layer geometry, Int. J. Hydrogen Energy 44 (2019) 8727–8742.
[21] A. Dejoan, V. N. Kurdyumov, Thermal expansion effect on the propagation of premixed flames in narrow channels of circular cross-section: Multiplicity of solutions, axisymmetry and non-axisymmetry, Proc. Combust. Inst. 37 (2019) 1927-1935.

[22] C. Jiménez, D. Fernández-Galisteo, V. N. Kurdyumov, DNS study of the propagation and flashback conditions of lean hydrogen-air flames in narrow channels: Symmetric and non-symmetric solutions, Int. J. Hydrogen Energy 40 (2015) 12541-12549.

[23] 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 (2011) 1307-1317.

[24] C. Jiménez, V. N. Kurdyumov, Propagation of symmetric and non-symmetric lean hydrogen-air flames in narrow channels: Influence of heat losses, Proc. Combust. Inst. 36 (2017) 1559-1567.

[25] A. Dejoan, C. Jiménez, V. N. Kurdyumov, Critical conditions for non-symmetric flame propagation in narrow channels: Influence of the flow rate, the thermal expansion, the lewis number and heat-losses, Combust. Flame 209 (2019) 430-440.

[26] C. Jiménez, D. Fernández-Galisteo, V. N. Kurdyumov, Flame-acoustics interaction for symmetric and non-symmetric flames propagating in a narrow duct from an open to a closed end, Combust. Flame 225 (2021) 499-512.

[27] 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 (2014) 1282-1293.

[28] Y. Ballossier, P. Boivin, C. Almarcha, Three dimensional shapes of hydrogen-air flames within millimetric hele shaw cells, Int. J. Hydrogen Energy 60 (2024) 333-341.

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

[30] 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 (2022).

[31] K. A. Kazakov, On-shell description of stationary flames, Phys. Fluids 17 (2005).

[32] F. Creta, P. E. Lapenna, R. Lamioni, N. Fogla, M. Matalon, Propagation of premixed flames in the presence of Darrieus-Landau and thermal diffusive instabilities, Combust. Flame 216 (2020) 256-270.

[33] L. Berger, M. Grinberg, B. Jürgens, P. E. Lapenna, F. Creta, A. Attili, H. Pitsch, Flame fingers and interactions of hydrodynamic and thermodiffusive instabilities in laminar lean hydrogen flames, Proc. Combust. Inst. 39 (2023) 1525-1534.

[34] K. Bhairapurada, B. Denet, P. Boivin, A Lattice-Boltzmann study of premixed flames thermo-acoustic instabilities, Combust. Flame 240 (2022) 112049.

[35] C. Almarcha, B. Radisson, E. Al Sarraf, E. Villermaux, B. Denet, J. Quinard, Interface dynamics, pole trajectories, and cell size statistics, Phys. Rev. E. 98 (2018) 030202.

[36] B. Radisson, B. Denet, C. Almarcha, Nonlinear dynamics of premixed flames: from deterministic stages to stochastic influence, J. Fluid Mech. 903 (2020) A17.

[37] M. Tayyab, S. Zhao, P. Boivin, Lattice-boltzmann modeling of a turbulent bluff-body stabilized flame, Phys. Fluids 33 (2021) 031701.
[38] 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.

[39] M. Tayyab, S. Zhao, Y. Feng, P. Boivin, Hybrid regularized lattice-boltzmann modelling of premixed and non-premixed combustion processes, Combust. Flame 211 (2020) 173-184.

[40] G. Farag, S. Zhao, T. Coratger, P. Boivin, G. Chiavassa, P. Sagaut, A pressure-based regularized lattice-boltzmann method for the simulation of compressible flows, Phys. Fluids 32 (2020) 066106.

[41] S. Zhao, K. Bhairapurada, M. Tayyab, R. Mercier, P. Boivin, Lattice-boltzmann modelling of the quiet and unstable preccinsta burner modes, Comput. Fluids. 260 (2023) 105898.

[42] P. Boivin, C. Jiménez, A. Sánchez, F. Williams, An explicit reduced mechanism for h2-air combustion, Proc. Combust. Inst. 33 (2011) 517-523.

[43] P. Boivin, A. Dauptain, C. Jiménez, B. Cuenot, Simulation of a supersonic hydrogen-air autoignition-stabilized flame using reduced chemistry, Combust. Flame 159 (2012) 1779-1790.

[44] http://www.cerfacs.fr/cantera/mechanisms/hydro.php#.

[45] J. Gaucherand, D. Laera, C. Schulze-Netzer, T. Poinsot, Intrinsic instabilities of hydrogen and hydrogen/ammonia premixed flames: Influence of equivalence ratio, fuel composition and pressure, Combust. Flame 256 (2023) 112986.

[46] J. Bechtold, M. Matalon, Hydrodynamic and diffusion effects on the stability of spherically expanding flames, Combust. Flame 67 (1987) 77-90.

[47] T. Poinsot, D. Veynante, Theoretical and numerical combustion, RT Edwards, Inc., 2005.

[48] B. Denet, Stationary solutions and neumann boundary conditions in the sivashinsky equation, Phys. Rev. E. 74 (2006) 036303.

[49] M. A. Liberman, M. F. Ivanov, O. E. Peil, D. M. Valiev, L.-E. Eriksson, Numerical studies of curved stationary flames in wide tubes, Combust. Theor. Model. 7 (2003) 653.

[50] M. A. Liberman, Introduction to physics and chemistry of combustion: explosion, flame, detonation, Springer Science & Business Media, 2010.

[51] 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 (2023) 1545-1554.

[52] A. B. Zel'dovich, A. G. Istratov, N. I. Kidin, V. B. Librovich, Flame propagation in tubes: Hydrodynamics and stability, Combust. Sci. Technol. 24 (1980) 1-13.

[53] R. Palomeque-Santiago, A. Domínguez-González, D. Martínez-Ruiz, M. Rubio-Rubio, E. Fernández Tarrazo, M. Sánchez-Sanz, Unveiling the bi-stable character of stealthy hydrogen-air flames, Phys. Fluids 36 (2024) 087137.

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

[55] G. Farag, T. Coratger, G. Wissocq, S. Zhao, P. Boivin, P. Sagaut, A unified hybrid lattice-Boltzmann method for compressible flows: Bridging between pressure-based and density-based methods, Phys. Fluids 33 (2021).

[56] J. Jacob, O. Malaspinas, P. Sagaut, A new hybrid recursive regularised Bhatnagar-Gross-Krook collision model for Lattice Boltzmann method-based large eddy simulation, J. Turbul. 5248 (2018) 1-26.

[57] M. Taha, S. Zhao, A. Lamorlette, J.-L. Consalvi, P. Boivin, Large eddy simulation of fire-induced flows using Lattice-Boltzmann methods, Int. J. Therm. Sci. 197 (2024) 108801.
## Appendix A. The lattice-Boltzmann method

The unified density/pressure-based Lattice-Boltzmann method (LBM) for compressible and multi-species flows [40, 55] is used for the mass and momentum transport in the current study. In practice, the time stepping of these two fields is achieved as

$$ 
\begin{aligned}
\rho(t + \Delta t, \mathbf{x}) &= \sum_{i=1}^{Q} f_i^{\text{col}}(t, \mathbf{x} - c_i \Delta t) + \rho(t, \mathbf{x})(1 - \theta(t, \mathbf{x})), \\
(\rho u_{\alpha})(t + \Delta t, \mathbf{x}) &= \sum_{i=1}^{Q} c_{i\alpha} f_i^{\text{col}}(t, \mathbf{x} - c_i \Delta t).
\end{aligned} 
\quad (A.1) $$

where $Q = 19$ corresponding to the D3Q19 lattice utilized in this study, $\theta$ is the normalized temperature defined as $\theta \equiv p/\rho/c_s^2$ with $c_s \equiv \Delta x/\Delta t/\sqrt{3}$ as the lattice sound speed. $c_i$ is the lattice speed of the $i^{th}$ population, and $f_i^{\text{col}}$ is the collide population which is evaluated locally using the hybrid recursive regularized method [56] as

$$ f_i^{\text{col}} = f_i^{\text{eq}} + \left(1 - \frac{\Delta t}{\bar{\tau}}\right) f_i^{\text{neq}} + \frac{\Delta t}{2} F_i^E, \quad (A.2) $$

where the equilibrium distribution $f_i^{\text{eq}}$ and non-equilibrium distribution $f_i^{\text{neq}}$ are calculated by

$$ 
\begin{aligned}
f_i^{\text{eq}} &= \omega_i \left[ \rho\theta + \frac{\mathcal{H}_{i,\alpha}^{(1)}}{c_s^2} \rho u_{\alpha} + \frac{\mathcal{H}_{i,\alpha\beta}^{(2)}}{2c_s^4} \rho u_{\alpha} u_{\beta} + \frac{\mathcal{H}_{i,\gamma}^{(3r)}}{6c_s^6} a_{\gamma}^{(3),\text{eq}} \right], \\
f_i^{\text{neq}} &= \omega_i \left[ \frac{\mathcal{H}_{i,\alpha\beta}^{(2)}}{2c_s^4} a_{\alpha\beta}^{(2),\text{neq}} + \frac{\mathcal{H}_{i,\gamma}^{(3r)}}{6c_s^6} a_{\gamma}^{(3r),\text{neq}} \right],
\end{aligned} 
\quad (A.3) $$

The parameter $\bar{\tau}$ in the above equation is the relaxation time in LBM, which is connected to the viscosity of the fluid

$$ \bar{\tau} = \frac{\mu}{\rho c_s^2} + \frac{\Delta t}{2}. $$

Detailed expressions of the lattice weights $\omega_i$, the Hermite polynomial $\mathcal{H}$, the third order equilibrium moments $a^{(3r),\text{eq}}$, and the non-equilibrium moments $a^{(2),\text{neq}}$ and $a^{(3r),\text{neq}}$ as well as the forcing term $F^E$ can be found in [57].