# Large Eddy Simulation of large-scale hydrogen deflagrations using the Thickened Flame Model with stretch sensitivity adaptation and thermo-diffusive instability modeling

Cédric Mehl, Sandy Poncet, Karine Truffin, Olivier Colin

► To cite this version:

Cédric Mehl, Sandy Poncet, Karine Truffin, Olivier Colin. Large Eddy Simulation of large-scale hydrogen deflagrations using the Thickened Flame Model with stretch sensitivity adaptation and thermo-diffusive instability modeling. International Journal of Hydrogen Energy, 2024, 93, pp.457-468. 10.1016/j.ijhydene.2024.10.169. hal-04877452

HAL Id: hal-04877452

https://ifp.hal.science/hal-04877452v1

Submitted on 9 Jan 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
<figure><img src="image_2.png" /></figure>

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

# Large Eddy Simulation of large-scale hydrogen deflagrations using the Thickened Flame Model with stretch sensitivity adaptation and thermo-diffusive instability modeling

Cédric Mehl*, Sandy Poncet, Karine Truffin, Olivier Colin

IFPEN, 1-4 av. du Bois Preau, Rueil-Malmaison 92852, France

Institut Carnot IFPEN Transports Energie, France

## ARTICLE INFO

**Keywords:**

Large Eddy Simulation

Hydrogen deflagrations

Thickened Flame Model

Flame stretch

Thermo-diffusive instabilities

## ABSTRACT

Due to the high costs associated to experimental testing of hydrogen ($H_2$) explosions, numerical simulations play an essential role in the design of safe industrial systems. In this work, we investigate the use of the Thickened Flame Model (TFM) to perform simulations of $H_2$ deflagrations. Specific issues arise when considering hydrogen: (i) due to differential diffusion, the flame is very sensitive to stretch; (ii) instabilities of the flame front are generated and lead to an increased flame propagation speed. Both of these phenomena are affected by flame thickening. First, the flame sensitivity to stretch is artificially increased by the thickening. The recently developed Ma-TFM model proposes a correction of diffusive fluxes to reach the correct Markstein length and is selected in this work. Secondly, thermo-diffusive and hydrodynamic flame instabilities are also altered by thickening. By computing the dispersion relationship for a hydrogen/air mixture, we show that the instability scales are proportional to the thickening factor $F$. A simple subgrid scale (SGS) model, based on the computations of laminar unstable 2-D and 3-D flames, is then considered to take into account the effect of instabilities. The Ma-TFM and SGS instability models are applied together for the first time on a semi-industrial vented lean hydrogen deflagration, leading to promising results. In particular, it is found that considering 3-D unstable flames is necessary to accurately predict the experimental pressure peaks.

## 1. Introduction

The increasing use of hydrogen ($H_2$) as an energy carrier in the industry poses a threat to the security of installations. $H_2$ is indeed more volatile and has larger flammability limits than standard hydrocarbons. As experimental testing of $H_2$ explosions is complex and costly, numerical simulations, and more specifically Computational Fluid Dynamics (CFD), play an essential role in the design of safe systems. Due to the large size of industrial installations, the whole range of turbulent flow scales can however not be solved. Two techniques have been considered in the context of safety CFD to tackle this issue, namely Unsteady Reynolds Averaged Simulations (URANS) [1–5] and Large Eddy Simulation (LES) [6–12].

In both LES and URANS, computational cell sizes are not fine enough to solve the flame front. A model is required to (i) properly describe the flame front; (ii) take into account the unresolved flame wrinkling, which has a significant impact on the turbulent flame speed. In most studies on safety CFD, a geometrical interpretation of the flame front is considered [13]. The flame is interpreted as a surface delimiting

fresh and burnt gases and is described by a filtered progress variable $\tilde{c}$. Several models have been proposed to ensure the correct propagation of the flame. The Zimont model [14] consists in writing the progress variable source term such that the flame propagates at the turbulent speed $S_t$, which depends on the laminar flame speed $S_l^0$ and the turbulent velocity fluctuation $u'$. The Zimont model has been extensively used in URANS [4,15,16]. It however suffers from two important drawbacks: (i) the model does not degenerate to the laminar flame speed when $u' = 0$; (ii) the model assumes an instantaneous equilibrium between the flame and turbulence. An extension of the Zimont model has been proposed to tackle these issues but still relies on simplifying assumptions and many model parameters [17]. An alternative is to consider Flame Surface Density (FSD) models. In this case the progress variable source term is computed from a flame surface density $\Sigma$ which is either transported or computed algebraically. FSD models are commonly used in safety URANS [1,2] and LES [9,10,12,18]. They give the correct flame speed when $u' = 0$ and, if $\Sigma$ is transported, take into account the flame response time to turbulence. Nevertheless, as both FSD and Zimont

<footnote>* Corresponding author.
E-mail address: cedric.mehl@ifpen.fr (C. Mehl).</footnote>
models only consider a progress variable to model the flame front, they do not capture complex chemical effects such as flame response to stretch (strain and curvature) or super-adiabatic temperatures in the burnt gases, which are important for hydrogen combustion.

An alternative to tackle this issue is to use direct integration to solve chemistry. In this context, the resolution issue (i) can be treated by thickening the flame front by a factor F using the Thickened Flame Model [19]. The unresolved flame wrinkling is then included through the use of an efficiency function, which can be computed algebraically [19,20] or dynamically [21,22]. The TFM model has been successfully applied for hydrocarbons on flame/obstacles configurations [6,23], highlighting its capacity to tackle safety problems. Issues however arise when considering the application of TFM to H₂ safety problems. First, although the flame is affected by stretch, the thickening leads to an incorrect response of the flame [24–26]. Corrections of TFM have thus been proposed in the literature [25,27]. These corrections were however developed and validated on simple academic cases at low to moderate values of F, and the behavior of these models at the large thickening factors encountered in safety CFD is yet unclear. Secondly, strong differential diffusion effects present with H₂ generate thermo-diffusive flame instabilities [28,29] which tend to increase the flame surface by a significant amount. As the flame structures generated by instabilities are, at least in part, smaller than typical LES cell sizes, subgrid scale (SGS) models are needed. Most approaches proposed so far to take into account the SGS contribution of TD instabilities involve an enhancement of the flame propagation speed, either through additional terms in the modeling of Sᵢ in the case of Zimont and FSD models [7,15,16,30] or through an enhancement of the efficiency factor in TFM [31,32]. Early models to estimate the flame speed increase rely on simple algebraic relations [15,16] or include the acceleration directly in the laminar flame speed estimation from spherical flame measurements [7,30]. More recently, Aniello et al. [31] used 2-D Direct Numerical Simulation (DNS) to estimate the wrinkling due to fully developed TD instabilities and exploited that information to improve the TFM efficiency function. A similar model was proposed by Hok et al. [32], which nevertheless used experimental data to inform the efficiency function model and were able to predict the early growth of TD instabilities.

The objective of the present work is to assess the capability of TFM, enhanced with models for the response to stretch and the SGS instabilities, to predict accurately large scale lean hydrogen deflagrations. The response of the thickened flame to stretch is corrected using the recently developed Ma-TFM model [25]. The impact of TD instabilities at the subgrid scale is included using the model of Aniello et al. [31]. The modeling framework is first presented in Section 2. Section 3 focuses on the assessment of the Ma-TFM model for spherical flames with a large thickening factor. Section 4 discusses the impact of thickening on TD instabilities and proposes a strategy to include their effects in LES. Finally, a semi-industrial H₂ deflagration is computed in Section 5.

## 2. Modeling and numerics

### 2.1. Thickened Flame Model (TFM)

#### 2.1.1. TFM-AMR

The Thickened Flame Model (TFM), originally proposed by Colin et al. [19], is selected in the present work to solve the propagating front. The flame is artificially thickened by a factor F and the impact of the subgrid scale (SGS) structures is modeled using an efficiency factor E [19,21]. A dynamic formulation of TFM is retained here [33]. In this formalism, thickening is restricted to the flame region, defined by S = 1, where S is a flame sensor function [34]. More specifically, we consider here the TFM-AMR extension proposed in [35], where adaptive mesh refinement (AMR) is used to resolve the flame front

(when S > 0). The transport equation for a given species mass fraction $\tilde{Y}_k$ reads:

$$ \frac{\partial \bar{\rho} \tilde{Y}_k}{\partial t} + \frac{\partial \bar{\rho} \tilde{\mu}_i \tilde{Y}_k}{\partial x_i} = \frac{\partial}{\partial x_i} \left( \left[ \mu_k^{eff} \frac{\partial \tilde{Y}_k}{\partial x_i} - \bar{\rho} \tilde{Y}_k V_i^c \right] \right) + \frac{E}{F} R_c(\gamma) \bar{\omega}_k \quad (1) $$

In Eq. (1), the left-hand side represents time variation and convection while terms on the right-hand side represent species diffusion and chemical source term. In this equation the effective diffusivity of species k, $\mu_k^{eff}$, depends on the thickening factor F, the efficiency function E, the molecular mass diffusivity $D_k$, the turbulent Schmidt number $Sc_t$, the turbulent viscosity $\mu_t$ and the sensor function S:

$$ \mu_k^{eff} = E F \frac{\bar{\rho} D_k}{\gamma} + (1-S) \frac{\mu_t}{Sc_t} \quad (2) $$

$\gamma$ and $R_c(\gamma)$ are stretch sensitivity correction parameters further detailed in Section 3.1. The correction velocity $V_i^c$ in Eq. (1) is introduced to ensure that the sum of species diffusion fluxes equals zero [13]:

$$ \bar{\rho} V_i^c = \sum_k \mu_k^{eff} \frac{\partial \tilde{Y}_k}{\partial x_i} \quad (3) $$

The thickening factor is computed as:

$$ F = 1 + (F_t - 1) S \quad (4) $$

where $F_t$ is a target thickening factor. In order to satisfy this required target thickening, the appropriate AMR refinement level n is applied [35] to reach a cell size $\Delta_x^{flame} = \Delta_x 2^{-n}$ that satisfies a prescribed resolution $n_{res}$ (here equal to 5) in the flame front:

$$ F_t \delta_t^0 (\phi, p, T_u) = n_{res} \Delta_x^{flame} \quad (5) $$

#### 2.1.2. TFM-AMR-E efficiency function

In the case of a flame front solved with AMR, the turbulence cascade inside the flame front is only partial, and the standard efficiency models, based on local values of turbulent velocities, may fail [35]. The TFM-AMR-E model proposed in [35] is here considered to tackle this effect. An effective SGS turbulent velocity $u'_{eff}$ and an effective scale $\Delta_{eff}$ are constructed by mimicking the partial turbulent cascading in the flame. This is performed by solving relaxation equations for $\Delta_{eff}$ and $u'_{eff}$ towards a local filter size $\Delta = \gamma \Delta_x$ representative of the largest resolved scales and the local turbulent velocity at scale $\hat{\Delta}$, denoted as $\hat{u}$. The relaxation time is assumed to scale with the largest turbulent scales as $\tau_t = \alpha \Delta_{eff} / \hat{u}_{eff}$. The efficiency is finally computed as [35]:

$$ E = \frac{\Xi_{\Delta_{eff}}^{\text{tot}} \left( \frac{\Delta_{eff}}{\delta_t^0}, \frac{u'_{eff}}{S_t^0} \right)}{\Xi_{\Delta_{eff}}^{\text{res}} \left( \frac{\Delta_{eff}}{\bar{\rho} \delta_t^0}, \frac{u'_{eff}}{S_t^0} \right)} \quad (6) $$

Where $\Xi_{\Delta_{eff}}^{\text{tot}}$ and $\Xi_{\Delta_{eff}}^{\text{res}}$ are terms approximating the total and resolved flame wrinkling up to scale $\Delta_{eff}$. Both terms are here computed using the flame wrinkling function $\Xi$ from Charlette et al. [21] with exponent $\beta = 0.5$. TFM-AMR-E model parameters are here set to $\gamma = 3$ and $\alpha = 1$, as suggested in [35]. Eq. (6) only tackles the SGS interactions between flame and turbulence. An extension of the model to take into account the SGS contribution of TD instabilities will be detailed in Section 4.2.

### 2.2. Chemistry and diffusive transport

$H_2$ combustion is described by the mechanism of Connaire et al. [36]. An appropriate modeling of thermal and mass diffusion coefficients is essential to recover properly the differential diffusion effects of $H_2$. An accurate approach is the mixture averaged model [13], which is however expensive in CFD because both the computation of the mass diffusion coefficients and the viscosity involve a double loop
on chemical species for each computational cell [37]. An alternative considered here is to tabulate the molecular viscosity as a function of the temperature using values in air, i.e. $\mu = \mu^{air}(T)$, and estimate the molecular mass diffusivities using a constant Schmidt number $Sc_k$ assumption:

$$D_k = \frac{\mu}{Sc_k} \qquad (7)$$

Due to its large conductivity, $H_2$ cannot be adequately described using a constant Prandtl number assumption. The conductivity is thus calculated from a mixture-averaged assumption [13]:

$$\lambda = \frac{1}{2} \left[ \sum_{k=1}^{N_S} X_k \lambda_k + \left( \sum_{k=1}^{N_S} \frac{X_k}{\lambda_k} \right)^{-1} \right] \qquad (8)$$

Where $X_k$ and $\lambda_k$ are the molar fraction and conductivity of species k, respectively. Contrarily to viscosity and diffusivities obtained from a mixture-averaged assumption, the computation of Eq. (8) only necessitates one loop over species and is therefore computationally cheaper. A validation of the current approach against the mixture-averaged model is provided in Appendix A.

## 2.3. Numerics

All the simulations presented in the following sections are performed with the CONVERGE 3.0 CFD solver [38]. It features in particular a second-order centered finite volume spatial discretization and second-order Crank-Nicholson time discretization. Chemistry integration is performed using the CVODE solver [39].

# 3. Analysis of the Ma-TFM model to correct the thickened flame response at large $F$

Flame stretch strongly affects the propagation speed of premixed hydrogen flames. In LES, the actual flame structure is however not resolved and the correct response of the flame to stretch is not guaranteed. In the specific context of TFM, the thickening leads to an alteration of the flame response, worsening as $F$ increases [24]. Several models have been proposed to improve the behavior of TFM [24–27,40]. The Ma-TFM model developed by Poncet et al. [25], which enforces the preservation of the Markstein length for the thickened flame, is here selected and briefly summarized below. It is then tested a posteriori on laminar spherical flames with $F = 100$.

## 3.1. Presentation of the Ma-TFM model

**Flame response to stretch.** The flame stretch K corresponds to the rate of change of a flame surface element subjected to strain and curvature. It affects in particular the local propagation speed of the flame. We write $S_d^b$ the flame displacement speed in the burned gases, and $\hat{S}_d^b = \sigma^{-1}S_d^b$ its normalized counterpart, where $\sigma = \rho_u/\rho_b$ is the flame expansion ratio. At low stretch, the linear relationship $\hat{S}_d^b = S_d^0 - L_d^b K$, where $L_d^b$ is the displacement Markstein length in the burned gases, is considered [41]. The corresponding Markstein number is $Ma_d^b = L_d^b/\delta_l^0$. For a flame thickened by a factor $F$, the Markstein length becomes $FL_d^b$ [24]. For hydrogen, $L_d^b < 0$ and the flame speed is thus over-estimated for positive stretch rates.

**Ma-TFM model.** Most models proposed in the literature to correct the response of thickened flames to stretch rely on a modification of mass diffusivities of species, with only the methodology to determine this correction changing [24–27,40]. The Ma-TFM model proposed in the work of Poncet et al. [25] features a correction of species diffusivities ensuring that the correct Markstein length is preserved for the thickened flame. The correction reads:

$$D_k^* = \frac{D_k}{\gamma} \Leftrightarrow Sc_k^* = \gamma Sc_k \qquad (9)$$

<figure><img src="image_4.png" /><figcaption>Fig. 1. Evolution of the correction factor $\gamma_F$ with the thickening factor $F$.</figcaption></figure>

Where $\gamma$ is a model parameter set to preserve $Ma_d^b$. This parameter depends on $F$ as well as the local thermo-chemical conditions ($p, T_0, \phi$). As the introduction of $\gamma$ modifies species diffusivities, the laminar flame speed $S_d^0$ is also modified. For this reason, a correction factor $R_\gamma(\gamma)$ is applied to all species reaction rates. For a given thickening factor $F$, the correction factor $\gamma_F$ is obtained by solving the following equation:

$$Ma_d^b(\gamma = \gamma_F) = \frac{Ma_d^b(\gamma = 1)}{F} \qquad (10)$$

where the function $Ma_d^b(\gamma)$ is obtained by computing a series of spherically expanding flames for varying $\gamma$. An interpolation of $Ma_d^b(\gamma)$ between selected $\gamma$ values is then necessary to solve Eq. (10). More details about the model can be found in [25].

## 3.2. Analysis of the model on spherical flames

Previous studies in the literature focused on the computation of stretched flames at small to moderate thickening factors and showed a good performance of the correction models [25–27]. In this section, we assess the performance of stretch correction with a large thickening factor ($F = 100$), representative of industrial security cases.

**Setup description.** Spherically expanding flames are used to determine the optimal $\gamma_F$ factors, and are also considered as a first test case for the stretch sensitivity correction model. The spherical flame equations are solved in one dimension along the radial axis. $Ma_d^b$ is obtained from a linear regression of the scaled flame displacement speed $\hat{S}_d^b$ at low stretch. We consider lean $H_2$ combustion at $\phi = 0.47$, $p = 101325$ Pa and $T_0 = 300$ K. The flame speed is then $S_d^0 = 0.4$ m/s, the thickness $\delta_l^0 = 0.49$ mm and the characteristic flame time $\tau_F = \delta_l^0/S_d^0 = 1.23$ ms.

**Determination of $\gamma_F$.** A series of spherical flames with varying $\gamma$ and $F = 1$ is first computed to determine $\gamma_F$. The evolution of $\gamma_F$ with $F$ is shown in Fig. 1. $\gamma_F$ rises steeply for $F$ up to around 20 and then slowly converges towards a constant value. The sensitivity of the flame to $\gamma_F$ is thus expected to increase with $F$.

**A posteriori spherical flame simulations.** Spherical flames are computed a posteriori using the obtained $\gamma_F$ factors. Two cases are here considered: (i) a moderate thickening factor $F = 8$; (ii) a large thickening factor $F = 100$. The robustness of the model is assessed by varying $\gamma_8$ and $\gamma_{100}$ by $\pm 1\%$ and analyzing the impact on $L_d^b$.

We first consider the case $F = 8$, for which the obtained correction factor is $\gamma_8 = 1.61$. The sensitivity study is performed by varying $\gamma_8$ by $\pm 1\%$, leading to $\gamma_8^- = 1.5939$ and $\gamma_8^+ = 1.6261$. Table 1 presents the Markstein lengths obtained for the non-thickened and thickened ($F =$
<caption>**Table 1**</caption>

Markstein lengths computed from spherical flames for the reference ($F = 1$) and thickened ($F = 8$ and $F = 100$) cases.

<table><thead><tr><th></th><th>L<sup>b</sup><sub>d</sub> (mm)</th><th>L<sup>b</sup><sub>d</sub>/L<sup>b,ref</sup><sub>d</sub></th></tr></thead><tbody><tr><td>F = 1</td><td>-0.372</td><td>1</td></tr><tr><td>F = 8, γ = 1</td><td>-3.05</td><td>8</td></tr><tr><td>F = 8, γ<sub>8</sub> = 1.61</td><td>-0.51</td><td>1.35</td></tr><tr><td>F = 8, γ<sub>8</sub><sup>-</sup> = 1.5939</td><td>-0.55</td><td>1.44</td></tr><tr><td>F = 8, γ<sub>8</sub><sup>+</sup> = 1.6261</td><td>-0.47</td><td>1.25</td></tr><tr><td>F = 100, γ = 1</td><td>-37.2</td><td>100</td></tr><tr><td>F = 100, γ<sub>100</sub> = 1.7844</td><td>-1.9</td><td>5.12</td></tr><tr><td>F = 100, γ<sub>100</sub><sup>-</sup> = 1.7666</td><td>-2.3</td><td>6.18</td></tr><tr><td>F = 100, γ<sub>100</sub><sup>+</sup> = 1.8022</td><td>-1.52</td><td>4.09</td></tr><tr><td>F = 100, L<sub>e</sub><sub>k</sub> = 1</td><td>6.98</td><td>-18.8</td></tr></tbody></table>

8) flames, as well as the ratio with the reference value ($F = 1$). The reference Markstein length in the present conditions is $L_d^{b,ref} = -0.372$ mm. As expected from theory, the non-corrected ($\gamma = 1$) simulation gives a Markstein length multiplied by 8. Although not perfectly accurate, the Ma-TFM model leads to an improvement in the prediction of $L_d^b$ compared to the uncorrected TFM. In this case, the influence of a small variation of $\gamma$ by $\pm 1\%$ is moderate, as the resulting $L_d^b$ varies by $\pm 7.8\%$. This variation is nevertheless larger than the initial $\pm 1\%$ applied on $\gamma$. As stated above, Eq. (10) requires an interpolation of $Ma_d^b(\gamma)$ to be solved, and the accuracy of the results depends on the interpolation errors. Combined with the fact that the determination of $\gamma$ becomes sensitive as $F$ increases (see Fig. 1), it can lead to non-negligible error on $L_d^b$. This error could be reduced by choosing a finer $\gamma$ grid for the computation of the spherical flames, as was done in the work of Poncet et al. [25].

The results for the $F = 100$ case are also provided in Table 1. In this case, we find a correction factor $\gamma_{100} = 1.7844$. The sensitivity study is also performed by varying $\gamma_{100}$ by $\pm 1\%$, leading to $\gamma_{100}^{-} = 1.7666$ and $\gamma_{100}^{+} = 1.8022$. An additional simulation with unity Lewis numbers is here also performed. Similarly to the $F = 8$ case, the non-corrected ($\gamma = 1$) simulation leads to a Markstein length multiplied by $F = 100$. Using the Ma-TFM model, the Markstein length is much closer to the exact value, although it mispredicts it by a factor 5.12, which is larger than the error made on the $F = 8$ case. Varying $\gamma_{100}$ by $\pm 1\%$ leads to a more substantial variation of $L_d^b$ by $\pm 20\%$, highlighting the fact that the sensitivity increases with increasing $F$. Finally, $L_d^b$ obtained with a unity Lewis assumption is positive and with a large magnitude. This is due to the fact that unity Lewis flames have a non-zero, albeit small, positive displacement Markstein length [41].

The results are graphically illustrated in Figs. 2(a) and 2(b) for $F = 8$ and $F = 100$, respectively. The displacement speed $\dot{S}_d^b$ is represented as a function of the stretch rate $K$. The $K$ scale is different in both graphs, since the actual stretch rate applied on the $F = 100$ spherical flame is much larger due to the large flame thickness. As stretch values close to 0 were not reached in the $F = 1$ simulation, the $F = 1$ line in the $F = 100$ figure corresponds to the low stretch linear extrapolation. By writing the scaled stretch rate as $K^* = K F \tau_f$, we see that both flames are actually facing similar $K^*$ values ($K^* \approx 0.1 - 0.2$). In any case, the behavior of the flame is drastically improved using the corrected TFM. Even though the error on the actual $L_d^b$ is large for $F = 100$, it is accurate on the range of curvature values to which the spherical flame is subjected. At these low $K$, the sensitivity of the model to $\gamma$ is also dampened. The performance of the Ma-TFM model is therefore depending on the resolved stretch rate applied to the flame front. This will be investigated further on an industrial application in Section 5.

## 4. TFM and thermo-diffusive instabilities

In the case of negative Markstein lengths, the flame sensitivity to stretch ultimately leads to the generation of cellular structures on the flame front [28,42–45]. The nature of these instabilities is either

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Evolution of the scaled displacement speed $\dot{S}_d^b$ as a function of stretch rate $K$ for laminar spherical $H_2$ flames. Reference ($F = 1$) cases are compared to thickened flames with and without the correction model of Section 3.1. Top: $F = 8$. Bottom: $F = 100$. Envelopes including the variation of $\gamma$ by $\pm 1\%$ are added for the Ma-TFM results (dashed lines).</figcaption></figure>

hydrodynamic or thermo-diffusive (TD), as evidenced by many authors. The cells have scales which can be much smaller than the LES cell size, and must therefore be considered in the subgrid-scale (SGS) turbulent combustion model [31,46,47]. In this section, we first investigate the impact of flame thickening on the instabilities, and then propose a simple SGS model based on the Aniello et al. model [31].

### 4.1. Artificial flame thickening and instabilities

The development of flame instabilities follows two stages [48]: (i) a linear phase, where small perturbations of the flame front grow exponentially; (ii) a non-linear phase, where the flame front is composed of cells with different sizes interacting with each other. The linear phase (i) can be investigated in spectral space by assuming that the frequencies present in the initial perturbation evolve independently. This phase can then be quantitatively studied by considering a dispersion relationship of the form $\omega(k)$, where $k$ is the wavenumber of a given perturbation and $\omega$ its growth rate. While most previous studies consider separate simulations to compute $\omega$ for each $k$ [42,48], Al Kassar et al. [49] recently proposed to determine $\omega(k)$ using one simulation initialized with a sum of sinusoids of different wavelengths.
<figure><img src="image_6.png" /><figcaption>Fig. 3. Dispersion relationships for the $\phi = 0.47$ atmospheric $H_2$ flame. The normalized growth rate $\dot{\omega} = \omega F \tau_f$ is plotted as a function of the normalized wavenumber $\hat{k} = kF \delta_l^0$. Square symbols correspond to the non-thickened ($F=1$) and circles to the thickened cases ($F=4$), respectively. Plain lines: with differential diffusion. Dashed lines: unity Lewis. Dotted lines: differential diffusion with Ma-TFM model.</figcaption></figure>

This method is selected in the present paper. The initial perturbation $y_f$ is prescribed as:

$$y_f(x, t_0) = \sum_{n=1}^{N} A_n(t_0) \sin \left( \frac{2\pi n x}{L_x} \right) \quad (11)$$

where $A_n(t_0)$ are the amplitudes of each sinusoid and $L_x$ the width of the simulation domain. The number of sinusoids $N$ is set at a value $N = L_x/(2\Delta_x)$, where $\Delta_x$ is the computational cell size. This initialization therefore includes harmonics from wavenumber $k_1 = 2\pi/L_x$ to $k_{max} = \pi/\Delta_x$. The determination of the dispersion relationship is here performed for a 2-D domain. It is obtained by performing a Discrete Fourier Transform (DFT) of an iso-contour of the temperature $T$ and extracting the exponential growth $\omega$ for each wavenumber $k$. A detailed explanation of the procedure may be found in [49].

Atmospheric $H_2$ combustion at $\phi = 0.47$ is considered. The domain width is set to $L_x = 3846^0$ and the initial sinusoids amplitudes are set to the same value $A_n(t_0) = 2.3 \times 10^{-6} \delta_l^0$ $\forall n$. These values, also proposed in the work of Al Kassar et al. [49], are found to give smooth dispersion relationships with the current numerical setup. The DFT is performed for an iso-contour $T = 1100$ KL, corresponding approximately to the location of the maximum heat release rate. $\Delta_x$ is here selected so that the flame thickness $\delta_l^0$ is resolved with 10 points. Simulations are performed for the reference flame ($F=1$) and thickened flames. To evaluate $\omega(k)$ for thickened flames, the domain size and initial perturbation amplitude are multiplied by $F$. We choose here to analyze the impact of the diffusion model by running unity Lewis and differential diffusion simulations, using the diffusivities exposed in Section 2.2, at $F=1$ and $F=4$, and the Ma-TFM model for $F=4$. An additional simulation with differential diffusion and $F=100$ is also performed to investigate the effect of a large thickening factor on the dispersion relationship.

Dispersion functions obtained for thickened and non-thickened flames are shown in Fig. 3. The normalized growth rate $\dot{\omega} = \omega F \tau_f$ is plotted as a function of the normalized wavenumber $\hat{k} = kF \delta_l^0$. The theoretical Darrieus–Landau (D–L) dispersion line is also added on the plot (dashed line). We verify, accordingly to theory [50], that the obtained dispersion functions are tangent to the D-L line as $\hat{k} = 0$. For a given diffusion model, either constant Schmidt or unity Lewis numbers, we observe that the $F=1$ and $F=4$ (without Ma-TFM correction) plots are close. This is also the case for the $F=100$ case computed with differential diffusion, showing that this observation still holds at large thickening factors. This suggests that the function $\dot{\omega}(\hat{k})$ is independent of $F$. It follows from the expressions of $\hat{k}$ and $\dot{\omega}$ that: $\omega(F) = \omega(F=1)/F$ and $k(F) = k(F=1)/F$. In other words, when thickening, the size

of the instability cells is multiplied by $F$, and the growth rate of these instabilities is divided by $F$. The unity Lewis instabilities are larger ($k$ smaller) and with a lower growth rate than the instabilities generated with differential diffusion. This is also expected from theory, as hydrodynamic instabilities are present at large scales and stabilized by diffusion effects for larger $k$ [50,51]. Thickening therefore strongly alters the nature and dynamics of flame instabilities. The conclusions drawn here are also supported by theory, as a multiplication of $\delta_l^0$ by $F$ leads to the same conclusions for both Sivashinsky and Matalon models for $\omega(k)$ [45,52].

The dispersion relationship obtained using the Ma-TFM model has a peak which is lower than the standard constant $S_{ck}$ simulations and larger than the unity Lewis cases. This is consistent with the effect of the model, which corrects the over-prediction of $\hat{S}_d^b$ by reducing the sensitivity to stretch. A decrease of the instabilities growth rate is then expected as they are driven by stretch sensitivity. The model however does not permit to recover the correct dispersion relationship, which, as argued above, would be $(\hat{k}(F=1)/F, \dot{\omega}(F=1)/F)$ in the normalized space.

The consequence of the results presented here is that for industrial security cases, where $F$ is very large (typically $> 100$), the instabilities will naturally evolve with a much reduced growth rate and much larger scales. We therefore assume in this work that the effects of instabilities only takes place at SGS scales.

## 4.2. SGS wrinkling model for instabilities

The instability cells lead to a significant flame wrinkling and therefore acceleration [28,53]. These effects may be superimposed to the effects of turbulence, which is the main source of flame wrinkling for usual hydrocarbon fuels. In the context of TFM, the unresolved flame wrinkling is modeled through the use of the efficiency factor $E$ as detailed in Section 2. We consider in this work the model of Aniello et al. [31], which we write in the context of the TFM-AMR-E model. This approach specifically assumes a decoupling of turbulence and instability in terms of both time and length scales, making it practically valid only for flames with low Karlovitz numbers. We would like to stress that this assumption of decoupling is probably not exact because to our knowledge no DNS study could show the presence of small scale laminar instabilities embedded in larger turbulent corrugations. But this assumption is at the present time one of the only ones proposed in the literature that accounts for the increase of the local consumption speed (through the stretch factor $I_0$) in unstable flames. Neglecting this effect leads to strongly under-predicted turbulent flame speeds. This hypothesis is therefore necessary in the absence of a more advanced method to address the complex interaction of physical scales. However, it is important to note that large-scale explosions often involve high turbulence levels implying large Karlovitz numbers. A proper modeling is thus needed for both low and high Karlovitz regimes.

The flame propagation speed $S_f$ is here expressed as $S_f = \Xi S_c$, where $\Xi$ is the flame wrinkling factor and $S_c = I_0 S_l^0$ is the local consumption speed, given as the product of the stretch factor $I_0$ and the laminar flame speed $S_l^0$. The TFM-AMR-E efficiency factor can be rewritten as:

$$E = \frac{I_0 \Xi_{\Delta_{eff}}^{tot}}{I_0^{res} \Xi_{\Delta_{eff}}^{res}} \quad (12)$$

where $I_0^{res}$ is the stretch factor corresponding to the resolved stretch. Since we assume that the effect of instabilities and turbulence act at separated scales, the total wrinkling can then be split as $\Xi_{\Delta_{eff}}^{tot} = \Xi_{\Delta_{eff}}^{turb} \Xi_{\Delta_{eff}}^{instab}$, where $\Xi_{\Delta_{eff}}^{turb}$ and $\Xi_{\Delta_{eff}}^{instab}$ are the contributions from turbulence and instabilities, respectively. The stretch factor $I_0$ is split in SGS and resolved contributions: $I_0 = I_0^{sgs} I_0^{res}$. The model of Aniello et al. [31] used here can be interpreted in this way: the SGS stretch factor caused by turbulence is assumed to be the same as the one caused by instabilities, which simply leads to $I_0^{sgs} = I_0^{instab}$ for both laminar
and turbulent conditions. This assumption may be considered over-simplified, but a proper modeling of this term for turbulent flows and for transitioning flows from a laminar unstable flame to a turbulent flame us still not available in the literature and thus requires further research which is out of the scope of the present paper.

The efficiency factor finally becomes:

$$E = I_{0}^{\text{instab}} \frac{\Xi_{\text{instab}}^{\text{turb}}}{\Xi_{\text{eff}}^{\text{res}}} \quad (13)$$

The term $\Xi_{\text{instab}}^{\text{turb}}/\Xi_{\text{eff}}^{\text{res}}$ is computed using Eq. (6). On the other hand, the term $E_{\text{instab}} = I_{0}^{\text{instab}}\Xi_{\text{instab}}^{\text{res}}$ will be assumed constant in the present work, and estimated using direct numerical simulations (DNS) of laminar unstable flames. This choice is motivated by the fact that ignition is not modeled here, and the initialization of simulations will be performed by adding a sphere of burned gases, whose radius is large due to the large thickening applied. It is therefore reasonable to assume that instabilities have already developed and reached a saturated flame acceleration.

### 4.3. Estimation of $E_{\text{instab}}$ using DNS of laminar unstable flames

An option to compute $E_{\text{instab}}$ is to extract this value from laminar flames [31]. This approach neglects the possible interactions between instabilities and turbulence [53], and is selected here due to its simplicity and the lack of data about instabilities/turbulence interactions, which is still a research problem under active research.

**Setup description.** A 2-D rectangular domain with height $L_y = n_y\delta_l^0$ and width $L_x = n_x\delta_l^0$ is initialized with a perturbed planar flame front at its center. The flame structure is extracted from a 1-D laminar premixed flame computed with the Cantera solver [54]. We consider a large $L_y$ value, here $n_y = 1500$, to avoid any influence from inlet and outlet boundary conditions on the results. The initial perturbation is $\psi(x) = \psi_0\sin(\frac{2\pi x}{\lambda})$, with $\psi_0 = \delta_l^0/2$ and $\lambda = 10\delta_l^0$. The flame front is refined using AMR with level $n = 3$ to reduce computational costs. A dependency of $S_f$ on the domain width $L_x$ has been reported in several studies [28,55]. Although Berger et al. [28] find an independence for very large domains, there is no theoretical bound on $n_x$ to reach independence, as unstable flame simulations are sensitive to physical models and numerics [55]. We choose here to perform simulations with $n_x = 100$ and 800 to observe the impact of the domain size on the propagation.

Recent work has shown that the 3-D nature of instability cells plays an important role on the quantitative determination of flame acceleration [56]. A 3-D simulation is therefore also performed here to assess the importance of considering three-dimensional structures to estimate $E_{\text{instab}}$. A third dimension $L_z = n_x\delta_l^0$ is thus considered. The 3-D initial perturbation is $\psi(x, z) = \psi_0\sin(\frac{2\pi x}{\lambda})\sin(\frac{2\pi z}{\lambda})$. Due to the large computational costs, only the $n_x = 100$ case is performed for the 3-D configuration.

**Results.** An instantaneous field of temperature is shown for the 2-D case with $n_x = 800$ in Fig. 4. As in the work of Berger et al. [28], we distinguish perturbations of the flame over a wide range of length scales. The smallest scales correspond to cells of unequal sizes distributed all along the flame front, typical of thermo-diffusive effects. Flame fingers, composed of several small cells are then formed and are presumably generated by hydrodynamic (Darrieus–Landau) instabilities. Fluctuations at larger scales, composed of several fingers, are also observed. A zoomed view of the flame front is also shown in Fig. 4 to highlight local variations of the heat release rate (HRR). The HRR is enhanced in regions concave towards burned gases while the flame is almost extinguished for strong curvatures in convex zones. This is a marker of strong differential diffusion effects. An iso-contour of temperature, colored by the stretch rate $K$, is provided in Fig. 5 for the 3-D case.

The same qualitative behavior, with small cells on top of larger cells, is observed.

The evolution of the wrinkling factor $\Xi$, the stretch factor $I_0$ and the normalized flame speed $S_f/S_l^0$ in time are represented in Figs. 6–8, respectively. Results are provided for the 2-D cases with $n_x = 100$ (black solid line) and $n_x = 800$ (blue line), as well as the 3-D case with $n_x = 100$ (red line). $S_f$ is here estimated by considering the integral of the $H_2$ source term $\dot{\omega}_{H_2}$ in the domain:

$$S_f = -\frac{1}{\rho_u Y_{H_2,u} A_0} \int \dot{\omega}_{H_2} dV \quad (14)$$

where $Y_{H_2,u}$ is the $H_2$ mass fraction in fresh gases and $A_0$ the projected surface of the flame. In the 2-D case, $A_0 = L_x$ and $dV = dxdy$ while in the 3-D case, $A_0 = L_xL_z$ and $dV = dxdydz$.

All cases reach a statistically stationary state after a sufficiently long time ($t \approx 100\tau_f$ for the 2-D cases and $t \approx 50\tau_f$ for the 3-D case). The 3-D case is stopped earlier than the 2-D cases as it is computationally expensive. As flame wrinkling is only generated by instabilities, we have here $E_{\text{instab}} = S_f/S_l^0$. In practice, $E_{\text{instab}}$ is computed by averaging $S_f/S_l^0$ once the steady state is reached. We see in Fig. 7 that the stretch factor $I_0$ is higher than 1, which is explained by the negative Markstein length in the present lean $H_2$ conditions. Moreover, $I_0$ is not affected by the domain width, and its value is therefore governed by the small cells. However, $I_0$ is significantly larger for the 3-D case with respect to 2-D cases. This is easily explained by a larger curvature value in 3-D, where one additional space dimension is considered. As illustrated in Fig. 6, the presence of the cells leads to a significant wrinkling of the flame ($\Xi > 1$). Once again, the wrinkling factor is larger for the 3-D case, as the ratio between a surface and its projected area is higher with 3 dimensions compared to 2 for a similar shape. Overall, this leads to a much larger propagation speed in the 3-D case, as shown in Fig. 8.

We observe in practice that the 2-D $n_x = 100$ simulation features one large unstable cell in the whole domain, which controls the value of the wrinkling. This suggests that the results at $n_x = 100$ might still depend on the width of the domain. The same phenomenon is observed for the 3-D simulation at $n_x = 100$. Due to this possible domain dependency, the consumption speed obtained for the 3D case is subject to uncertainties. However, we see that the $n_x = 800$ and $n_x = 100$ values are close for the 2D case. Recent work [56] also suggests that 3D unstable flames consumption speed converges above $n_x = 100$. We therefore admit here that the $E_{\text{instab}}$ value obtained for the 3D case is a good estimate of the expected value. Further work will nevertheless be required to obtain this value accurately at a moderate cost.

### 4.4. Extension of the model to complex cases

The model proposed in this paper is valid for a system in given conditions ($\phi, T_0, p$). The extension of the strategy to a more complex system with varying conditions, such as hydrogen explosions in a stratified mixture, would require the generation of a table $E_{\text{instab}}(\phi, T_0, p)$. As argued above, the cost of 3D DNS are prohibitive for a table generation. Further research efforts are therefore needed to extend the model presented in this paper to multi-conditions systems.

## 5. Application: industrial $H_2$ deflagration case

In this final section, the Ma-TFM model and the SGS instability model are applied together on a semi-industrial deflagration case.

### 5.1. Setup description

A $H_2$ deflagration setup experimented at INERIS is selected [57]. A photograph of the configuration is shown in Fig. 9. It consists in a 4 m³ chamber filled with premixed $H_2$/air in atmospheric conditions at an equivalence ratio $\phi = 0.47$. In the considered case, the mixture is ignited at the back of the box, and as the pressure increases, the fresh $H_2$/air mixture is vented out of the box. The experiment is
<figure><img src="image_7.png" /><figcaption>Fig. 4. Instantaneous temperature field for the 2-D simulation with $n_x = 800$. The zoom represents the heat release rate for a portion of the flame front.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Temperature iso-surface for the 3-D simulation with $n_x = 100$. The surface is colored by the stretch rate $K$ ($s^{-1}$).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 6. Evolution of the wrinkling factor $Ξ$ as a function of normalized time $t/τ_f$.</figcaption></figure>

equipped with several pressure sensors, inside and outside the box, as illustrated in Fig. 9. Several teams have simulated this setup and reported their results in the paper of Vyazmina et al. [58]. Although

<figure><img src="image_10.png" /><figcaption>Fig. 7. Evolution of the stretch factor $I_0$ as a function of normalized time $t/τ_f$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. Evolution of the normalized turbulent flame speed $S_f/S_i^0$ as a function of normalized time $t/τ_f$. Dashed lines correspond to the averaged $E_{instab}$.</figcaption></figure>

satisfying results were obtained, the simulations were based on many simplifying assumptions and ad-hoc parameters. The approach tested in the present study aims at predicting this setup with a more advanced
<figure><img src="image_12.png" /><figcaption>Fig. 9. Photograph of the INERIS setup with associated pressure probes positions.</figcaption></figure>

<table><caption>**Table 2**
Summary of the simulations performed in this study.</caption><thead><tr><td>Simulation</td><td>Diffusion</td><td>γ</td><td>E<sub>instab</sub></td></tr></thead><tbody><tr><td>A</td><td>D<sub>k</sub> = μ/γ Sc<sub>k</sub></td><td>1.7844</td><td>1</td></tr><tr><td>B</td><td>D<sub>k</sub> = μ/γ Sc<sub>k</sub></td><td>1.7844</td><td>2.83</td></tr><tr><td>C</td><td>D<sub>k</sub> = μ/γ Sc<sub>k</sub></td><td>1.7844</td><td>4.84</td></tr><tr><td>D</td><td>D<sub>k</sub> = μ/Sc<sub>k</sub></td><td>–</td><td>4.84</td></tr><tr><td>E</td><td>Le = 1</td><td>1.0</td><td>4.84</td></tr></tbody></table>

description of stretched flames and instability physics, in an attempt to be more widely applicable.

A set of simulations is performed to assess the importance of both Ma-TFM and SGS instabilities modeling. As for the spherical and unstable flames, the CONVERGE solver is used. The ignition of the hydrogen/air mixture is promoted by introducing a hemisphere of burned gases in the back of the box at the initial instant of the simulation. The TFM-AMR method is applied with a target thickening factor $\Gamma_t = 100$.

The selected simulations, numbered from A to E, are summarized in Table 2. A simulation with the Ma-TFM model and no model for SGS instabilities is first performed (simulation A). The impact of adding SGS instability modeling is then assessed by performing two simulations with $E_{instab} = 2.83$ (simulation B) and $E_{instab} = 4.84$ (simulation C), respectively. The former uses the $E_{instab}$ value obtained from the laminar 2-D simulation while the latter mimics an enhancement due to 3-D effects, as explained in Section 4.3. A simulation with $E_{instab} = 4.84$ and without the stretch correction model is also carried out to investigate the impact of the Ma-TFM model (simulation D). Finally, a computation using unity Lewis numbers $Le = 1$ is proposed (simulation E).

## 5.2. Pressure peaks predictions

The accuracy of the simulations is assessed by comparing the results to experimental pressures obtained at two locations: (i) P1, located in the box and (ii) $P(x = 2m)$, located outside the box. As ignition is not modeled in the present work, the timing of the pressure peak is not reproduced numerically. Since engineers are mostly interested by the pressure peak value, results are represented for shifted profiles, where the peak at P1 is set to be at $t = 0$ s, for both experimental and simulated pressures.

Fig. 10 shows the results obtained for the five simulations at probes P1 (top) and $P(x = 2m)$ (bottom). Before leaving the chamber, the propagation of the flame is mostly laminar and the wrinkling is thus due to instabilities. Turbulence is then generated as gases are pushed through the vent. We see that when neglecting the influence of SGS instability wrinkling (simulation A), the laminar flame propagation is very slow and unable to generate a visible pressure peak. This suggests that considering instabilities is important in this phase. By using the

<figure><img src="image_13.png" /><figcaption>Fig. 10. Comparison of the pressure evolution obtained from experiments and numerical simulations. Top: sensor P1. Bottom: sensor $P(x = 2m)$.</figcaption></figure>

SGS instability model proposed in the present work, we observe an increase of the pressure peak. A good agreement is found for $E_{instab} = 4.84$ (simulation C) regarding the magnitude of the pressure peaks for probes P1 and $P(x = 2m)$. A slight time shift is however observed for probe $P(x = 2m)$, where the computational pressure peak is delayed compared to experiments. For $E_{instab} = 2.83$ (simulation B), we see an under-prediction of the pressure peak at both probes. The 2-D instability wrinkling is thus too low, and taking into account 3-D effects is probably necessary.

The influence of the Ma-TFM model is then assessed separately by comparing simulations C and D. We see an over-prediction of the pressure peak when TFM is not corrected. This is explained by an over-prediction of $I_0^{res}$ as the Markstein length is multiplied by $\Gamma = 100$. This shows that even with the relatively low curvatures found in industrial scale systems, a high thickening factor influences significantly the propagation of the flame front. Finally, we observe that simulation E, performed with unity Lewis numbers, leads to results similar to the Ma-TFM simulation C. This is explained by the fact that using unity Lewis also decreases the flame sensitivity to stretch significantly (see Table 1). This is investigated further by focusing on stretch distribution at several instants in the next section.

## 5.3. Stretch distribution and displacement speed

An in-depth comparison of the TFM diffusion models is performed by analyzing the distribution of flame stretch and estimating the flame displacement speed for each model. Figs. 11 and 12 represent the distribution of the resolved stretch K for two instants during the deflagration, namely an instant when the flame front is inside the box and mostly laminar (Fig. 11) and a second instant when the flame has left the box and has interacted with turbulence (Fig. 12). The data is extracted from simulation C, which leads to the best match with experimental data as shown in Section 5.2. The stretch is estimated on
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Probability Density Function (PDF) of the stretch *K* in the initial laminar phase of the deflagration. Estimation of the flame propagation speed as a function of stretch using the linear relation $S_l^u - L_b K$ is added to the plot. Top: iso-surface of the progress variable $c^* = 0.8$.</figcaption></figure>

a progress variable *c* iso-surface, for $c^* = 0.8$. The progress variable is here defined as $c = (Y_{H_2O} - Y_{H_2O}^u)/(Y_{H_2O}^b - Y_{H_2O}^u)$, where $Y_{H_2O}^u$ and $Y_{H_2O}^b$ are the $H_2O$ mass fractions in the fresh and burnt gases, respectively. Similar *K* distributions were obtained for $c^* = 0.5$ (not shown here), showing that the conclusions drawn here are not dependant on the choice of $c^*$. A representation of the iso-surface $c^* = 0.8$ is provided on the top in Figs. 11 and 12. For the first instant, in Fig. 11, the stretch distribution is narrow around a slightly positive mean value $\bar{K} = 16.9 \text{ s}^{-1}$ and a standard deviation $\sigma = 43.1 \text{ s}^{-1}$. On the contrary, the distribution for the second instant, in Fig. 12, is more spread out. The mean is here $\bar{K} = 111.2 \text{ s}^{-1}$ and the standard deviation $\sigma = 177.5 \text{ s}^{-1}$. This may be explained by the influence of turbulence, which generates a wider range of strain and curvature on the flame front.

The evolution of the flame displacement speed as a function of the stretch using the linear Markstein relation $S_l^0 - L_b K$ is shown on the top of the *K* distribution in Figs. 11 and 12. Although the linear relationship $\dot{S}_d^b = S_l^0 - L_d^b K$ is not valid for high stretch rates, it gives a useful indication about the magnitude of the error between the reference and thickened flames. The results are represented for four estimations of the Markstein length $L_b$, depending on the diffusion model: (i) *F* = 1 with constant Schmidt numbers (reference value); (ii) *F* = 100 and constant Schmidt without Ma-TFM; (iii) *F* = 100 and constant Schmidt with Ma-TFM; (iv) *F* = 100 and unity Lewis. The corresponding displacement Markstein length values can be found in Table 1. In Fig. 11, where the stretch distribution is narrow, we see that the propagation speed for the *F* = 100 non corrected model rapidly departs from the reference value. Differences in flame propagation speed are then expected, even for the small resolved stretch at this instant. This is corroborated by the difference seen in Fig. 10, where the pressure peak for the non corrected TFM over-estimates the pressure peak for the first sensor, located inside the box. Conversely, the corrected TFM and unity Lewis propagation speeds are relatively close to the reference value in the range [$\bar{K} - \sigma, \bar{K} + \sigma$]. Even though the

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Probability Density Function (PDF) of the stretch *K* in the turbulent propagation phase. Estimation of the flame propagation speed as a function of stretch using the linear relation $S_l^0 - L_b K$ is added to the plot. Top: iso-surface of the progress variable $c^* = 0.8$.</figcaption></figure>

Markstein lengths are quantitatively different, and even positive in the case of the unity Lewis simulation, the difference is here seldom visible due to the low values of resolved stretch. This explains the good results obtained for both corrected and uncorrected TFM in Fig. 10.

For the second instant (Fig. 12), when the flame is outside the box, the resolved stretch is more spread out and takes larger values. In this case, we see again that the uncorrected *F* = 100 estimation of the propagation velocity is wrong. The Ma-TFM prediction is again relatively close to the reference in the range [$\bar{K} - \sigma, \bar{K} + \sigma$]. The unity Lewis estimation induces here a non-negligible difference compared to the reference, especially near the limit $K = \bar{K} + \sigma$.

To further investigate the differences between the Ma-TFM (simulation C) and unity Lewis (simulation E) models, Fig. 13 shows the shape of the gas cloud at a given instant. Fig. 13(a) represents an image acquired by a high speed video camera of the gas cloud seeded with ammonium chloride microparticles at time *t* = 155 ms [57]. At that time, the hydrogen gas cloud has left the box and burns outside. Fig. 13(b) and (c) represent the field of $Y_{H_2}$ at time *t* = 110 ms in the simulation, which takes into account the time shift $\Delta t \approx 45$ ms. The comparison is here qualitative as the experimental image does not enable to measure the actual value of hydrogen concentration. The strong difference between the Ma-TFM (b) and unity Lewis (c) simulations is however pretty clear as the burning rate in the Ma-TFM case is larger. Although very qualitative, the Ma-TFM simulation is more consistent with the experiment, where the gas cloud has expanded and hydrogen is depleted on the forefront of the gas cloud. The differences between Ma-TFM and unity Lewis can be understood from Fig. 12: as positive stretch rates are larger in the regions where the flame interacts with turbulence, the flame propagation speed is underestimated with unity Lewis due to the large positive Markstein length at *F* = 100. The present results therefore suggest a better behavior of the Ma-TFM model, even though the Markstein length is not quantitatively retrieved, compared to the unity Lewis case. Nevertheless, an important limitation of the
<figure><img src="image_16.png" /><figcaption>**Fig. 13.** Analysis of the gas ejection and combustion at a given instant ($t = 155$ ms in experiment and $t = 110$ ms in simulation to account for the time shift). (a): Experimental high speed video image of the gas seeded by microparticles of ammonium chloride. (b): $Y_{H_2}$ field for simulation C (3D instability model and Ma-TFM). (c) $Y_{H_2}$ field for simulation E (3D instability model and unity Lewis diffusion.).</figcaption></figure>

present simulations is that they assume a decoupling between turbulence and instability effects on the flame wrinkling. Additional work is needed to understand better the interactions between instabilities and turbulence.

## 6. Conclusions

Simulations of deflagrations in an industrial safety context are nototonously difficult to perform due to the large size of the domains usually considered. For the TFM model formalism considered in the present work, this translates into a very large thickening factor $F$. Additional issues then arise when considering the specific case of hydrogen: (i) due to differential diffusion, the flame is very sensitive to stretch; (ii) instabilities of the flame front lead to an increased flame propagation speed.

In the context of TFM, the flame sensitivity to stretch is artificially increased by the thickening. The recently developed Ma-TFM model [25] proposes a correction of diffusive fluxes to reach the correct Markstein length and has been selected in the present paper. It is assessed for the first time for a large thickening factor ($F = 100$). A sensitivity study shows that the Ma-TFM model becomes more sensitive as $F$ increases although its accuracy remains satisfying for the low stretch rates encountered in the present work.

Thermo-diffusive and hydrodynamic flame instabilities are also altered by the flame thickening. By computing the dispersion relationship for a hydrogen/air mixture, it is shown here that the instability scales are proportional to the thickening factor $F$. In the case of the large thickening factors encountered in industrial safety, it is thus safe to assume that these instabilities will only develop at the subgrid scale. Using this assumption, a simple model, assuming that turbulence and instabilities affect the flame independently, is considered here [31]. The

key ingredient of the model is an estimation of the flame acceleration due to instabilities, which is estimated here from (i) 2-D and (ii) 3-D DNS of laminar unstable $H_2$/air flames. We observe in particular a much higher value for the flame acceleration when the 3-D case is considered compared to the 2-D one, which is in good agreement with the literature [59].

The Ma-TFM and SGS instability models are applied for the first time together on a semi-industrial vented $H_2$ deflagration [57]. By considering the influence of instabilities on the flame propagation, an accurate prediction of the pressure peaks induced by the deflagration is achieved. The importance of using the Ma-TFM stretch correction is also highlighted. Additionally, we observe that using the flame acceleration from 3-D DNS of unstable flames gives a better agreement than the 2-D DNS values.

The necessity to use models to better predict $H_2$ flame sensibility to stretch and the effects of SGS instabilities in industrial safety is thus established in the present work in the context of thickened flames. Provided adequate models, cases with large thickening factors become feasible. The proposed methodology, however, depends on a scale separation between instabilities and turbulence. Further research is needed to gain a deeper understanding of the interactions between these two physical phenomena and to develop a more accurate LES model that can effectively capture them. There is more specifically a need to develop models able to predict the laminar-to-turbulent transition of $H_2$/air premixed deflagrations. Additionally, further work is required to compute the instability induced flame acceleration for laminar flames, as the 3-D DNS are computationally very expensive and cannot be realistically considered for a large range of conditions.

## CRediT authorship contribution statement

**Cédric Mehl:** Writing – review & editing, Writing – original draft, Validation, Software, Methodology, Data curation, Conceptualization.

**Sandy Poncet:** Writing – review & editing, Methodology. **Karine Truffin:** Writing – review & editing, Methodology, Conceptualization.

**Olivier Colin:** Writing – review & editing, Methodology, Conceptualization.

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

## Acknowledgments

The authors would like to thank Benjamin Truchot and Guillaume Lecocq from INERIS for providing the experimental results and for the fruitful discussions about the experimental results, and Antonio Attili from for providing information on the dispersion relationship calculations.

## Appendix A. Comparison of diffusion models

The diffusion model presented in Section 2.2, combining a constant Schmidt numbers assumption for species diffusion and a mixture-averaged formulation for the conductivity, is here assessed against the standard mixture-averaged model, which is here used as a reference. Fig. A.14 shows the displacement speed evolution as a function of stretch for a laminar spherical $H_2$ flame with $\phi = 0.47$, $T_0 = 300$ KL and $p = 1$ atm. The result for the constant Schmidt numbers along with a constant Prandtl number for the conductivity is also added to the plot.

We see that the constant Schmidt number model with mixture averaged conductivity is in good agreement with the full mixture-averaged model. This justifies the use of this model in this work. We see also that significant discrepancies are obtained when computing the conductivity from a constant Prandtl number. In particular, the laminar flame speed as well as the Markstein length are not well predicted.
<figure><img src="image_17.png" /><figcaption>**Fig. A.14.** Evolution of the scaled displacement speed $\dot{S}_d^b$ as a function of stretch rate $K$ for a laminar spherical $H_2$ flame ($\phi = 0.47$, $T_0 = 300$ KL, $\rho = 1$ atm). Several formulation for the diffusive fluxes are compared: (i) the mixture-averaged formulation; (ii) the constant Schmidt model with conductivity from Eq. (8); (iii) the constant Schmidt model with constant Prandtl assumption for the conductivity.</figcaption></figure>

## References

[1] Catlin C, Fairweather M, Ibrahim S. Predictions of turbulent, premixed flame propagation in explosion tubes. Combust Flame 1995;102(1-2):115-28. http://dx.doi.org/10.1016/0010-2180(94)00245-N.

[2] Fairweather M, Hargrave G, Ibrahim S, Walker D. Studies of premixed flame propagation in explosion tubes. Combust Flame 1999;116(4):504-18. http://dx.doi.org/10.1016/S0010-2180(98)00055-8.

[3] Patel S, Jarvis S, Ibrahim S, Hargrave G. An experimental and numerical investigation of premixed flame deflagration in a semiconfined explosion chamber. Proc Combust Inst 2002;29(2):1849-54. http://dx.doi.org/10.1016/S1540-7489(02)80224-3.

[4] Sathiah P, Van Haren S, Komen E, Roekaerts D. The role of CFD combustion modeling in hydrogen safety management—II: Validation based on homogeneous hydrogen-air experiments. Nucl Eng Des 2012;252:289-302. http://dx.doi.org/10.1016/j.nucengdes.2012.06.023.

[5] Ferrara G, Di Benedetto A, Salzano E, Russo G. CFD analysis of gas explosions vented through relief pipes. J Hard Mater 2006;137(2):654-65. http://dx.doi.org/10.1016/j.jhazmat.2006.03.037.

[6] Vermorel O, Quillatre P, Poinsot T. LES of explosions in venting chamber: A test case for premixed turbulent combustion models. Combust Flame 2017;183:207-23. http://dx.doi.org/10.1016/j.combustflame.2017.05.014.

[7] Molkov V, Makarov D, Schneider H. Hydrogen-air deflagrations in open atmosphere: Large eddy simulation analysis of experimental data. Int J Hydrog Energy 2007;32(13):2198-205. http://dx.doi.org/10.1016/j.ijhydene.2007.04.021.

[8] Li Y, Xie H, Bi M, Bo Y, Gao W. Effects of cloud size and built-in obstacles on hydrogen cloud explosion using large eddy simulation. J Loss Prev Process Ind 2022;77:104788. http://dx.doi.org/10.1016/j.ljp.2022.104788.

[9] Di Sarli V, Di Benedetto A, Russo G, Jarvis S, Long EJ, Hargrave GK. Large eddy simulation and PIV measurements of unsteady premixed flames accelerated by obstacles. Flow Turbul Combust 2009;83(2):227-50. http://dx.doi.org/10.1007/s10494-008-9198-3.

[10] Kirkpatrick M, Armfield S, Masri A, Ibrahim S. Large eddy simulation of a propagating turbulent premixed flame. Flow Turbul Combust (formerly Appl Sci Res) 2003;70(1-4):1-19. http://dx.doi.org/10.1023/B:APPL.0000004912.87854.35.

[11] Masri A, Ibrahim S, Cadwallader B. Measurements and large eddy simulation of propagating premixed flames. Exp Therm Fluid Sci 2006;30(7):687-702. http://dx.doi.org/10.1016/j.expthermflusci.2006.01.008.

[12] Gubba SR, Ibrahim SS, Malalasekera W, Masri AR. LES modeling of premixed deflagrating flames in a small-scale vented explosion chamber with a series of solid obstructions. Combust Sci Technol 2008;180(10-11):1936-55. http://dx.doi.org/10.1080/00102200802261852.

[13] Poinsot T, Veynante D. Theoretical and numerical combustion. third ed.. 2012.

[14] Zimont VL. Theory of turbulent combustion of a homogeneous fuel mixture at high Reynolds numbers. Combust Explos Shock Waves 1979;15:305-11.

[15] Holler T, Komen EM, Kljenak I. The role of CFD combustion modeling in hydrogen safety management - VII: Validation for hydrogen deflagration in large-scale hydrogen-air-steam experiment. Nucl Eng Des 2019;342:133-46. http://dx.doi.org/10.1016/j.nucengdes.2018.11.033.

[16] Tolias I, Venetsanos A, Markatos N, Kiranoudis C. CFD evaluation against a large scale unconfined hydrogen deflagration. Int J Hydrog Energy 2017;42(11):7731-9. http://dx.doi.org/10.1016/j.ijhydene.2016.07.052.

[17] Lipatnikov A, Chomiak J. Turbulent flame speed and thickness: phenomenology, evaluation, and application in multi-dimensional simulations. Prog Energy Combust Sci 2002;28(1):1-74. http://dx.doi.org/10.1016/S0360-1285(01)00007-7.

[18] Di Sarli V, Di Benedetto A, Russo G. Using large eddy simulation for understanding vented gas explosions in the presence of obstacles. J Hard Mater 2009;169(1-3):435-42. http://dx.doi.org/10.1016/j.jhazmat.2009.03.115.

[19] Colin O, Ducros F, Veynante D, Poinsot T. A thickened flame model for large eddy simulations of turbulent premixed combustion. Phys Fluids 2000;12(7):1843-63. http://dx.doi.org/10.1063/1.870436.

[20] Charlette F, Meneveau C, Veynante D. A power-law flame wrinkling model for LES of premixed turbulent combustion Part II: dynamic formulation. Combust Flame 2002;131(1-2):181-97. http://dx.doi.org/10.1016/S0010-2180(02)00401-7.

[21] Charlette F, Meneveau C, Veynante D. A power-law flame wrinkling model for LES of premixed turbulent combustion Part I: non-dynamic formulation and initial tests. Combust Flame 2002;131(1-2):159-80. http://dx.doi.org/10.1016/S0010-2180(02)00400-5.

[22] Veynante D, Moureau V. Analysis of dynamic models for large eddy simulations of turbulent premixed combustion. Combust Flame 2015;162(12):4622-42. http://dx.doi.org/10.1016/j.combustflame.2015.09.020.

[23] Volpiani PS, Schmitt T, Vermorel O, Quillatre P, Veynante D. Large eddy simulation of explosion deflagrating flames using a dynamic wrinkling formulation. Combust Flame 2017;186:17-31. http://dx.doi.org/10.1016/j.combustflame.2017.07.022.

[24] Quillatre P. Simulation aux grandes échelles d'explosions en domaine semi-confiné (Ph.D. thesis), Toulouse, INPT; 2014.

[25] Poncet S, Mehl C, Truffin K, Colin O. A thickened flame model adaptation to weakly stretched flames for non-unity Lewis number mixtures. Combust Flame 2024;270:113758. http://dx.doi.org/10.1016/j.combustflame.2024.113758.

[26] Popp S, Kuenne G, Janicka J, Hasse C. An extended artificial thickening approach for strained premixed flames. Combust Flame 2019;206:252-65. http://dx.doi.org/10.1016/j.combustflame.2019.04.047.

[27] Detomaso N, Hok J-J, Dounia O, Laera D, Poinsot T. A generalization of the thickened flame model for stretched flames. Combust Flame 2023;258:113080. http://dx.doi.org/10.1016/j.combustflame.2023.113080.

[28] Berger L, Kleinheinz K, Attili A, Pitsch H. Characteristic patterns of thermodynamically unstable premixed lean hydrogen flames. Proc Combust Inst 2019;37(2):1879-86. http://dx.doi.org/10.1016/j.proci.2018.06.072.

[29] Goulier J, Chaumeix N, Halter F, Meynet N, Bentaïb A. Experimental study of laminar and turbulent flame speed of a spherical flame in a fan-stirred closed vessel for hydrogen safety application. Nucl Eng Des 2017;312:214-27. http://dx.doi.org/10.1016/j.nucengdes.2016.07.007.

[30] Molkov V, Makarov D, Grigorash A. Cellular structure of explosion flames: modeling and large-eddy simulation. Combust Sci Technol 2004;176(5-6):851-65. http://dx.doi.org/10.1080/00102200490428495.

[31] Aniello A, Laera D, Berger L, Attili A, Poinsot T. Introducing thermodiffusive effects in large-eddy simulation of turbulent combustion for lean hydrogen-air flames. In: Proceedings of the summer program, center for turbulence research. 2022, p. 267-77.

[32] Hok J-J, Dounia O, Detomaso N, Jaravel T, Douasbin Q, Vermorel O. A modeling strategy for the thickened flame simulation of propagating lean hydrogen-air flames. Int J Hydrog Energy 2024;78:1133-41. http://dx.doi.org/10.1016/j.ijhydene.2024.06.298.

[33] Légier JP, Poinsot T, Veynante D. Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion. In: Proceedings of the summer program, center for turbulence research. 2000, p. 157-68.

[34] Jaravel T, Riber E, Cuenot B, Bulat G. Large eddy simulation of an industrial gas turbine combustor using reduced chemistry with accurate pollutant prediction. Proc Combust Inst 2017;36(3):3817-25. http://dx.doi.org/10.1016/j.proci.2016.07.027.

[35] Mehl C, Liu S, Colin O. A strategy to couple thickened flame model and adaptive mesh refinement for the LES of turbulent premixed combustion. Flow Turbul Combust 2021. http://dx.doi.org/10.1007/s10494-021-00261-2.

[36] Ó Conaire M, Curran HJ, Simmie JM, Pitz WJ, Westbrook CK. A comprehensive modeling study of hydrogen oxidation: A comprehensive modeling study of hydrogen oxidation. Int J Chem Kinet 2004;36(11):603-22. http://dx.doi.org/10.1002/kin.20036.

[37] Wilke CR. A viscosity equation for gas mixtures. J Chem Phys 1950;18(4):517-9. http://dx.doi.org/10.1063/1.1747673.

[38] Richards KJ, Senecal PK, Pomraning E. Converge 3.0. 2022.

[39] Hindmarsh AC, Brown PN, Grant KE, Lee SL, Serban R, Shumaker DE, et al. SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers. ACM Trans Math Software 2005;31(3):363-96. http://dx.doi.org/10.1145/1089014.1089020.

[40] Comer AL, Gallagher TP, Duraisamy K, Rankin BA. A modified thickened flame model for simulating extinction. Combust Theory Model 2022;26(7):1262-92. http://dx.doi.org/10.1080/13647830.2022.2141138.
[41] Giannakopoulos GK, Frouzakis CE, Mohan S, Tomboulides AG, Matalon M. Consumption and displacement speeds of stretched premixed flames - theory and simulations. Combust Flame 2019;208:164–81. http://dx.doi.org/10.1016/j.combustflame.2019.06.027.

[42] Frouzakis CE, Fogla N, Tomboulides AG, Altantzis C, Matalon M. Numerical study of unstable hydrogen/air flames: Shape and propagation speed. Proc Combust Inst 2015;35(1):1087–95. http://dx.doi.org/10.1016/j.proci.2014.05.132.

[43] Altantzis C, Frouzakis CE, Tomboulides AG, Matalon M, Boulouchos K. Hydrodynamic and thermodiffusive instability effects on the evolution of laminar planar lean premixed hydrogen flames. J Fluid Mech 2012;700:329–61. http://dx.doi.org/10.1017/jfm.2012.136.

[44] Bechtold J, Matalon M. Hydrodynamic and diffusion effects on the stability of spherically expanding flames. Combust Flame 1987;67(1):77–90. http://dx.doi.org/10.1016/0010-2180(87)90015-0.

[45] Sivashinsky G. Nonlinear analysis of hydrodynamic instability in laminar flames—I. derivation of basic equations. Acta Astronaut 1977;4(11–12):1177–206. http://dx.doi.org/10.1016/0094-5765(77)90096-0.

[46] Lapenna PE, Remiddi A, Molinaro D, Indelicato G, Creta F. A-posteriori analysis of a data-driven filtered wrinkled flamelet model for thermodiffusively unstable premixed flames. Combust Flame 2024;259:113126. http://dx.doi.org/10.1016/j.combustflame.2023.113126.

[47] Lapenna PE, Berger L, Attili A, Lamioni R, Fogla N, Pitsch H, et al. Data-driven subfilter modelling of thermo-diffusively unstable hydrogen–air premixed flames. Combust Theory Model 2021;1–22. http://dx.doi.org/10.1080/13647830.2021.1925350.

[48] Berger L, Attili A, Pitsch H. Intrinsic instabilities in premixed hydrogen flames: parametric variation of pressure, equivalence ratio, and temperature. part 2 – non-linear regime and flame speed enhancement. Combust Flame 2022;240:111936. http://dx.doi.org/10.1016/j.combustflame.2021.111936.

[49] Kassar SA, Berger L, Lapenna PE, Creta F, Pitsch H, Attili A. Efficient and accurate calculation of dispersion relations for intrinsically unstable premixed flames. 2023, arXiv:2306.10901 [physics].

[50] Matalon M, Cui C, Bechtold JK. Hydrodynamic theory of premixed flames: effects of stoichiometry, variable transport coefficients and arbitrary reaction orders. J Fluid Mech 2003;487:179–210. http://dx.doi.org/10.1017/S0022112003004683.

[51] Matalon M. Intrinsic flame instabilities in premixed and nonpremixed combustion. Annu Rev Fluid Mech 2007;39(1):163–91. http://dx.doi.org/10.1146/annurev.fluid.38.050304.092153.

[52] Matalon M, Matkowsky BJ. Flames as gasdynamic discontinuities. J Fluid Mech 1982;124(-1):239. http://dx.doi.org/10.1017/S0022112082002481.

[53] Berger L, Attili A, Pitsch H. Synergistic interactions of thermodiffusive instabilities and turbulence in lean hydrogen flames. Combust Flame 2022;244:112254. http://dx.doi.org/10.1016/j.combustflame.2022.112254.

[54] Goodwin DG, Moffat HK, Schoegl I, Speth RL, Weber BW. Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes. 2022, http://dx.doi.org/10.5281/ZENODO.742000.

[55] Sharpe GJ, Falle SAEG. Nonlinear cellular instabilities of planar premixed flames: numerical simulations of the reactive Navier–Stokes equations. Combust Theory Model 2006;10(3):483–514. http://dx.doi.org/10.1080/13647830500472354.

[56] Wen X, Berger L, Cai L, Parente A, Pitsch H. Thermodiffusively unstable laminar hydrogen flame in a sufficiently large 3D computational domain – Part I: Characteristic patterns. Combust Flame 2023;113278. http://dx.doi.org/10.1016/j.combustflame.2023.113278.

[57] Daubech J, Proust C, Gentilhomme D, Jamois D, Ichard M. Hydrogen-air vented explosions—new experimental data. 2013.

[58] Vyazmina E, Jallais S, Krumenacker L, Tripathi A, Mahon A, Commanay J, et al. Vented explosion of hydrogen/air mixture: An intercomparison benchmark exercise. Int J Hydrog Energy 2019;44(17):8914–26. http://dx.doi.org/10.1016/j.ijhydene.2018.07.195.

[59] Wen X, Zirwes T, Scholtissek A, Böttler H, Zhang F, Bockhorn H, et al. Flame structure analysis and composition space modeling of thermodiffusively unstable premixed hydrogen flames — Part II: Elevated pressure. Combust Flame 2022;238:111808. http://dx.doi.org/10.1016/j.combustflame.2021.111808.