# Ammonia-air Flame-Turbulence Interaction in forced HIT environments

Solène Hoflack, Noé Monnier, Pierre Brequigny, Seif Zitouni, Christine

Mounaïm-Rousselle, Quentin Douasbin, Eleonore Riber

► To cite this version:

Solène Hoflack, Noé Monnier, Pierre Brequigny, Seif Zitouni, Christine Mounaïm-Rousselle, et al. Ammonia-air Flame-Turbulence Interaction in forced HIT environments. European Combustion Meeting (ECM), Apr 2025, Edimbourg (Ecosse), United Kingdom. hal-05250777

HAL Id: hal-05250777

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

Submitted on 11 Sep 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
Ammonia-air Flame-Turbulence Interaction in forced HIT environments

S. Hoflack<sup>a</sup>, N. Monnier<sup>b</sup>, P. Brequigny<sup>b</sup>, S. Zitouni<sup>b</sup>, C. Rousselle<sup>b,c</sup>,

Q. Douasbin<sup>a</sup>, E. Riber<sup>a</sup>

<sup>a</sup>CERFACS, 42 Avenue Gaspard Coriolis, 31057 Toulouse, France

<sup>b</sup>University of Orléans, PRISME, Château de la Source, 45072 Orléans, France

<sup>c</sup>Institut Universitaire de France, Paris, France

**Abstract**

In the context of climate change induced by anthropological greenhouse gas emissions, the need to find renewable and low-carbon energy sources is growing even more. While green hydrogen ($H_2$) is a promising carbon-free fuel, its storage and transport induce challenges, including risks of leaks and, potentially, fast deflagration or detonation. Ammonia ($NH_3$) stands out as a convenient hydrogen carrier and zero-carbon fuel. However, limited knowledge of its combustion hinders its use. Improving our understanding and modeling of ammonia combustion is then crucial. To do so, PRISME laboratory designed a spherical constant-volume chamber with integrated fans to study $NH_3$/air flames, at ambient pressure and at $T = 423$ K, under forced homogeneous isotropic turbulence (HIT). The experimental study focus on the flame propagation and the importance of the fuel mixture under turbulent conditions. In this work, two Direct Numerical Simulations (DNS) of turbulent $NH_3$/air flames are performed using the compressible Navier-Stokes solver AVBP, reproducing the experimental conditions at two turbulent intensities $u' = 0.49$ m/s and $u' = 0.65$ m/s. A turbulence forcing method, based on the addition of energy at large scales, is used to counteract the decay of turbulence. A semi-detailed Analytically Reduced Chemistry (ARC), derived from NUIG detailed chemistry, is derived to accurately reproduce laminar $NH_3$/air flames under experimental conditions. Flame topology, displacement speed, and wrinkling ratio are analysed. A good agreement is found between results from experimental 2D Mie-Scattering tomography and 2D flame iscontours from DNS. The simulations correctly capture the effect of turbulence on the flame development, as the flame distortion at larges scales. The influence of the barycentre displacement on the flame speed and the wrinkling ratio resulting from 2D DNS and experimental diagnostics is then discussed. Results from 3D DNS isosurfaces are also determined to account for the full flame dynamics. Differences between 2D and 3D results finally highlight the need for 3D diagnostics techniques to achieve greater measurements accuracy.

Keywords: Ammonia; Turbulence forcing; Turbulent spherical expanding flames; DNS simulations
# 1. Introduction

The combustion of fossil fuels account for nearly 80 % of the worldwide energy use [1] leading to green house gas emission responsible for the anthropogenic climate change.

Storing renewable energy excess using electro-fuels or e-fuels is a possible solution to decarbonize our energy sources. By doing so, the energy produced from renewable electricity is stored in the chemical bonds of e-fuels. NH₃ is one of the potential e-fuel candidates that has raised interest in the combustion research community. Its high volumetric density, low storage pressure, and long-term stability ensure easier and safer storage compared to hydrogen [2]. However, the high minimum ignition energy, low laminar burning velocity, and high auto-ignition temperature of NH₃ make its combustion challenging. One potential solution to address these issues is to blend it with a more reactive fuel, such as hydrogen. In this case, the laminar flame speed is increased, and NOₓ emissions are reduced [3]. Introducing turbulence during the combustion process can also enhance mixing and flame properties, including the rate of fuel consumption, as explained in [4]. Focus on ammonia turbulent-flame interaction has received little attention [5, 6] with greater understanding required.

The purpose of this paper is twofold: to numerically simulate the development of the spherical NH₃/air flame using a turbulence forcing method, and to analyse experimental and numerical results to highlight the specific characteristics of NH₃/air flame-turbulence interaction. The experimental setup of Zitouni et al. [7] is first described in Section 2. Then, Section 3 focuses on the numerical setup representing the experiment of Zitouni et al. [7]. The results of the DNS are presented and discussed in Section 5 for two turbulent intensities.

# 2. Experimental framework

The experimental dataset is provided by Zitouni et al. [7]. A schematic description of the experimental set-up is given in Fig. 1. It consists of a constant-volume spherical vessel used for both premixed laminar and turbulent flame propagation. The combustion chamber has an inner diameter of 200 mm and six four-blade fans are integrated near the wall in an hexahedral configuration.

In the centre of the vessel, the fans produce a nearly homogeneous isotropic turbulence (HIT) zone. High-speed particle image velocimetry (PIV) was used to quantify the turbulent fluctuation velocity $u'$ and the integral length scale $l_t$.

The flame propagation was measured at three turbulent fluctuations, $u' = 0.49, 0.65$ and $0.82$ m/s, which corresponds to fan speeds of 3000, 4000, and 5000 rpm, correspondingly. Using two-point velocimetry, a mean integral length scale of $l_t = 2.64$ mm was measured. More details about the experimental set-up and the turbulence characterisation can

<figure><img src="image_2.png" /><figcaption>Fig. 1: Schematic view of the combustion vessel (adapted from [8]).</figcaption></figure>

be found in Galmiche et al. [8].

The vessel is filled with unburnt gases at an initial temperature of 298 K and atmospheric pressure of 0.1 MPa. Different mixtures are studied: NH₃/CH₄ and NH₃/H₂ fuel blends, and pure CH₄/air flames at an equivalence ratio of $\phi = 1$. Due to its low laminar flame speed, the temperature was increased up to 423 K for pure NH₃/air flames. Finally, laser Mie-scattering tomography is used to capture the flame propagation and analyse the flame structure.

In this study, only the NH₃/air mixture at stoichiometry is analysed at the two first turbulent intensities ($u' = 0.49$ m/s and $u' = 0.65$ m/s). The turbulence characteristics of these cases are summarised in Table 1, with the turbulent Reynolds number $Re_t = u'l_t/\nu$ where $\nu$ is the kinematic viscosity and $l_t$ is the integral length scale, the turbulent kinetic energy $k$, the turbulent time scale $\tau_t = l_t/u'$, the energy dissipation rate $\epsilon$, the Kolmogorov length scale $\eta_k$ and the Karlovitz number $Ka = (u'/s_L^0)^{3/2}/(l_t/\delta_L)^{1/2}$ with $s_L^0 = 1.34 \times 10^{-1}$ m/s and $\delta_L = 1.88 \times 10^{-3}$ m the laminar flame speed and flame thickness respectively.

<table><caption>Table 1: Characteristics of HIT field for the experimental cases 1 ($u' = 0.49$ m/s and $l_t = 2.64$ mm) and 2 ($u' = 0.65$ m/s and $l_t = 2.64$ mm).</caption><thead><tr><th>Variable [units]</th><th>Case 1</th><th>Case 2</th></tr></thead><tbody><tr><td>u' [m/s]</td><td>0.49</td><td>0.65</td></tr><tr><td>k [m<sup>2</sup>/s<sup>2</sup>]</td><td>0.360</td><td>0.634</td></tr><tr><td>&epsilon; [m<sup>2</sup>/s<sup>3</sup>]</td><td>44.6</td><td>104</td></tr><tr><td>&tau;<sub>t</sub> [s]</td><td>5.39 &times; 10<sup>-3</sup></td><td>4.06 &times; 10<sup>-3</sup></td></tr><tr><td>Re<sub>t</sub> [-]</td><td>42.8</td><td>56.8</td></tr><tr><td>&eta;<sub>k</sub> [m]</td><td>1.58 &times; 10<sup>-4</sup></td><td>1.28 &times; 10<sup>-4</sup></td></tr><tr><td>Ka [-]</td><td>6.17</td><td>9.43</td></tr></tbody></table>

# 3. Numerical set-up

## 3.1. Chemistry modelling
Based on the exhaustive kinetic model evaluation conducted by Girhe et al. [9], the detailed kinetic mechanism developed by Y. Zhu et al. [10] (39 species and 306 reactions) was selected, since it provides an overall good agreement with experimental data for species concentration, ignition delay time and laminar burning velocity. In this study, the chosen mechanism is denoted as "Zhu_2024" herewithin. Figure 2 illustrates the laminar flame speed versus the equivalence ratio at atmospheric pressure and for fresh gas temperatures of 298 K and 423 K (computed with the code CANTERA [11]) and compares with experimental data from [12–14].

<figure><img src="image_3.png" /><figcaption>Fig. 2: Laminar flame speed [cm/s] versus equivalence ratio [-] of NH₃/air flames at P = 1 bar, T = 298 K and T = 423 K. Experimental data from [12] (square markers), [13] (circle markers), and [14] (triangle markers). The Y. Zhu et al. mechanism [10] is represented by dotted and solid lines.</figcaption></figure>

Although the Zhu_2024 mechanism tends to overestimate the laminar flame speed for rich mixtures, it accurately reproduces the experimental data for lean and stoichiometric mixtures. The agreement is very good for the operating conditions considered in this work: $\phi = 1$ and T = 423 K.

To save computational cost and reduce temporal stiffness, this mechanism has been reduced to obtain a semi-detailed Analytically Reduced Chemistry (ARC) using the reduction code ARCANE [15]. First, only relevant species and reactions are kept to obtain a skeletal mechanism using DRGEP. In a second step, the number of transported species is further reduced using the Quasi-Steady State Approximation (QSSA) [16].

In the context of this study, two ARC mechanisms has been derived: one targeting NOₓ (12 species, 7 species with QSS assumption and 103 reactions) and another without targeting NOₓ (10 species, 9 species with QSS assumption and 101 reactions). These two ARC mechanisms have been validated and compared to the detailed Zhu_2024 mechanism for a range of temperature T = 300-500 K and equivalence ratio $\phi = 0.7-1.4$. As an example, Fig. 3 highlights the good agreement between ARC and the detailed mechanism for the final molar fraction of NH₃, H₂, H₂O, NO and N₂O species of

1D freely propagating ammonia/air flames at atmospheric pressure and a temperature of 423 K.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Final molar fraction [-] versus equivalence ratio [-] of NH₃/air flames at P = 1 bar and T = 423 K. The detailed mechanism from Zhu et al. [10] (Zhu_2024) is represented by solid lines, the ARC mechanism targeting NOₓ (ARC) by triangles, and the ARC mechanism without targeting NOₓ (ARC_NOx) by circles.</figcaption></figure>

The second overarching aspect to reproduce the experimental setup is to generate the turbulence which will interact with the flame.

## 3.2. Homogeneous Isotropic Turbulence forcing method

In the experiment of Zitouni et al. [7], a fan-generated near-HIT field is maintained at the centre of the combustion vessel throughout the flame expansion, up to a diameter of 40 mm from the centre of the vessel. To numerically reproduce such a turbulence, the objective is to impose a volumetric turbulence forcing source term in the domain to reach a statistically stationary equilibrium state. When the equilibrium state is reached, the averaged rate of energy added by forcing is equal to the averaged energy dissipation rate.

In this work, the HIT method selected is the one of Eswaran and Pope [17]. It consists in a spectral forcing method, where the forcing is applied only to a low wavenumber band, i.e. the largest scales of the turbulence. An example of application can be found in [18].

The continuity and the momentum equations can be expressed in wave number space as:

$$k_i \hat{u}_i(k, t) = 0 \qquad (1)$$

$$\frac{\partial \hat{u}_i(k, t)}{\partial t} = \hat{a}_i(k, t) + \hat{a}_i^F(k, t) \qquad (2)$$

with $k_i$ the $i^{th}$ component of the mode $k = (k_1, k_2, k_3)$, $\hat{u}_i(k, t)$ the Fourier coefficients of the velocity field, $\hat{a}_i(k, t)$ the Fourier coefficient of the rate of change of velocity with time, and $\hat{a}_i^F(k, t)$ the Fourier coefficients of the forcing vector.
A time-evolving stochastic source term $\hat{a}_i^F$ is added in Eq. (2), which can be represented as a truncated Fourier series:

$$ \hat{a}_i^F(\boldsymbol{x}, t) = \sum_{k \le k_f} e^{i\boldsymbol{k}\boldsymbol{x}} \hat{a}_i(\boldsymbol{k}, t) \quad (3) $$

with $k_f = 2\sqrt{2}k_0$ the radius of the sphere where forcing occurs, and $k_0$ the lowest wave number. In other words, $k_f$ represents the maximum wave number of the forced modes. The summation in the right-hand side of Eq. (3) is applied over all forced wave numbers in wave-space.

Then to ensure the continuity condition of Eq. (1), each complex Fourier amplitude is obtained from the divergence-free projection of a stochastic vector $\hat{b}(\boldsymbol{k}, t)$:

$$ \hat{a}_i^F(\boldsymbol{k}, t) = \hat{b}_i(\boldsymbol{k}, t) - \frac{k_i k_j \hat{b}_j(\boldsymbol{k}, t)}{k^2} \quad (4) $$

Here, $\hat{b}_i(\boldsymbol{k}, t)$ is composed of six independent random processes, called Uhlenbeck-Ornstein (UO) processes:

$$ \hat{b}_i(\boldsymbol{k}, t + dt) = \hat{b}_i(\boldsymbol{k}, t)\left(1 - \frac{dt}{\tau_f}\right) + \sigma_f \sqrt{\frac{2dt}{\tau_f}} N_t^{t+dt} \quad (5) $$

In Eq. (5), $N_t^{t+dt}$ stands for a Brownian fluctuation which follows a normal Gaussian distribution with zero mean and unit variance within the time interval $[t, t+dt]$. The parameter $\tau_f$ represents a forcing time-scale, which controls when turbulence reaches a steady state, whereas $\sigma_f$ controls the turbulent kinetic energy level, i.e. the amplitude of the forcing. $\sigma_f$ is controlled by the number of modes in the forced band $N_f$, the target energy dissipation rate $\epsilon$, the characteristic forcing time scale $\tau_f$, and the lowest wave number $k_0$. In this method, the integral length scale $l_t$ can be controlled by choosing a length $L_0$ as:

$$ l_t \sim \frac{1}{k_0} = \frac{L_0}{2\pi} \quad (6) $$

Three input parameters are finally necessary to generate a forced HIT field, i.e. $L_0$ [m], $\tau_f$ [s], and $\sigma_f$ [m/s²] using the Eswaran and Pope method, where the target energy dissipation rate $\epsilon$ and the integral length scale $l_t$ are fixed according to the experimental setup. The root mean square of the velocity fluctuations $u'$ is deduced from the scaling $\epsilon \sim u'^{3/l_t}$.

### 3.3. Numerical configuration

A single numerical set-up is used for the two cases referenced in Table 1. The configuration consists of a 3D closed spherical vessel with an adiabatic no-slip wall as boundary condition, and composed of unstructured tetrahedral cells as illustrated in Fig. 4. The

sphere has a diameter of 200 mm with a refined cubic zone at the centre of 17 mm in length, since the zone of interest is a central cube of 16 mm in length. The coarse zone allows to keep the same volume as the experiment and to maintain a constant atmospheric pressure during the interrogation window.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Numerical constant-volume spherical vessel configuration.</figcaption></figure>

The discretization in the refined zone enables to have approximately 10 points in the ammonia/air laminar flame thickness at stoichiometry and is sufficiently refined to reproduce the energy cascade generated by turbulence for both cases. In fact, the resolution criterium defined by Pope [19] is respected with $\eta_k/\Delta x > 0.5$, where $\eta_k$ is the Kolmogorov length scale, i.e. the smallest scale of the turbulence (see Table 1), and $\Delta x$ is the characteristic length of a tetrahedral cell in the refined zone. In the generated mesh, the cell size $\Delta x = 1.80 \times 10^{-4}$ m yields ratios of $\eta_k/\Delta x = 0.88$ for case 1 ($u' = 0.49$ m/s), and $\eta_k/\Delta x = 0.71$ for case 2 ($u' = 0.65$ m/s), both above the resolution criterion of 0.5.

The 3D DNS is performed with the in-house compressible Navier-Stokes solver AVBP [20] using the Lax-Wendroff scheme (2nd order accurate in both time and space).

First, the forced HIT field is generated for an air mixture and only in the refined zone for 50 turbulent time scales $\tau_t$ to let the HIT reach its statistically stationary state. Table 2 indicates the values chosen to generate the HIT field needed to match the experimental conditions (see Table 1) for the two cases.

In a second step, a NH₃/air flame is initiated with an energy deposit at the center of the domain [21]. For cases 1 and 2, an energy amount of 14 mJ deposited on a spherical kernel of diameter equal to $2.5 \times 10^{-3}$ m for 3.5 ms, is enough to ignite the flame. The ARC
<table><caption>Table 2: Input parameters for cases 1 and 2 for the HIT forcing method (see Table 1).</caption><thead><tr><th>Variable [units]</th><th>Case 1</th><th>Case 2</th></tr></thead><tbody><tr><td>L<sub>0</sub> [m]</td><td>1.66 &times; 10<sup>-2</sup></td><td>1.66 &times; 10<sup>-2</sup></td></tr><tr><td>&tau;<sub>f</sub> [s]</td><td>1.62 &times; 10<sup>-2</sup></td><td>1.22 &times; 10<sup>-2</sup></td></tr><tr><td>&sigma;<sub>f</sub> [m/s<sup>2</sup>]</td><td>6.2</td><td>10.6</td></tr></tbody></table>

1 mechanism denoted as "ARC" in Fig. 3 is used for the
2 reactive simulations.

## 4. Results and discussion

### 4.1. Homogeneous Isotropic Turbulence field

5 First, the HIT field generated with the method
6 of Eswaran and Pope [17] is analysed. In this sec-
7 tion, only case 1 is reported for clarity's sake, as case
8 2 presents very similar features. Note that results have
9 been analysed and are presented only in the region of
10 interest defined in Fig. 4, specifically the central cube
11 with an edge length of L = 16 mm. The velocity
12 magnitude field at t = 50 &tau;<sub>t</sub> is shown in Fig. 5.

<figure><img src="image_6.png" /><figcaption>Fig. 5: 2D cut of the velocity magnitude field in the X/Y plan ($z = 0$ mm) at $t = 50 \tau_t$ for case 1.</figcaption></figure>

13 Statistical stationarity is reached around $t = 5 \tau_t$.
14 To ensure statistical reliability, the turbulent quantities reported here are measured in the time interval $t \in [15 \tau_t, 50 \tau_t]$. The integral length scale is calculated from the two-point velocity correlation [19] as:

$$l_t(t) = \int_0^L \frac{\langle u(\mathbf{x}, t)u(\mathbf{x} + \mathbf{e}_1 r) \rangle}{\langle u^2 \rangle} dr \quad (7)$$

18 with $u$ the velocity along the longitudinal direction,
19 $r$ the distance between two points in the domain, and
20 $L = 16$ mm in our case. Due to the isotropic nature
21 of the considered turbulent field, the turbulent kinetic
22 energy (TKE) can be calculated as follows [19]:

$$k(t) = \frac{3}{2} u'(t)^2 \quad (8)$$

23 with $u'$ the RMS velocity fluctuations.

24 The time evolution of the volume-averaged turbulent kinetic energy (TKE) and the integral length scale

<figure><img src="image_7.png" /><figcaption>Fig. 6: Turbulent kinetic energy $k$ [$m^2/s^2$] and integral length scale $l_t$ [m] over the time $t$ [s] normalised by the turbulent time scale $\tau_t$ [s] in the zone of interest for case 1 ($u' = 0.49$ m/s).</figcaption></figure>

26 ($l_t$) in the region of interest is reported in Fig. 6. A
27 very good agreement is reached for the TKE as the
28 experimental target value is recovered with a relative
29 error of only -0.28 %. Although the integral length
30 scale reported in Fig. 6 shows slightly larger fluctua-
31 tions around its mean value than the TKE, it is also in
32 good agreement with experimental target value with a
33 relative error of -0.76 %. The relative errors on these
34 turbulent quantities are negligible compared to typical
35 experimental variability for such forced HIT configu-
36 rations.

37 Note that for case 2 (not shown here), good results
38 have also been obtained with relative errors of 0.16
39 % and 13 % on the mean TKE and the mean integral
40 length scale, respectively.

### 4.2. Flame propagation

42 The well-reproduced stationary steady-state turbulent fields obtained for both case 1 and case 2 allow to study the key aspect of this work: the flame/turbulence interaction for fully premixed ammonia/air flames.

47 Figure 7 depicts the flame shape of both cases, in the time interval $t \in [50\tau_t, 53\tau_t]$ with a sampling rate of $\Delta t = 0.5$ ms, using isocontours of progress variable based on the fuel mass fraction $c_{Y_F}$:

$$c_{Y_F} = \frac{Y_{\text{NH3}} - Y_{\text{NH3}}^u}{Y_{\text{NH3}}^b - Y_{\text{NH3}}^u} \quad (9)$$
where $Y_{NH3}$ is the fuel mass fraction and the superscripts u and b denote the unburnt and burnt gases, respectively. Here, the value of $c_{Y_F} = 0.1$ is chosen to be close to the unburnt side as in the experiment of Zitouni et al. [7]. Selecting a low value of $c_{Y_F}$ allows to study the turbulent field seen by the flame wrinkling while avoiding a bias due to the flame thermal expansion.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Flame isocontours at $c_{Y_F} = 0.1$ at a regular time interval of $\Delta t = 0.5$ ms for $t \in [50\tau_t, 53\tau_t]$ for case 1 (Top) and case 2 (Bottom).</figcaption></figure>

In agreement with the measurements, Fig. 7 shows that the flame does not remain spherical as it interacts with turbulent structures. Besides, it is noteworthy that the flame topology differs between case 1 and case 2. For case 1, with the lowest turbulent intensity, the flame has a nearly-symmetrical shape and is less stretched than for case 2.

Figure 8 shows the location of the flame barycentre over time in each direction. The barycentre coordinates $(x_0, y_0, z_0)$ at a given time $t$ is calculated as:

$$C = \frac{1}{N} \sum_{i=1}^{N} p_i \qquad (10)$$

with $C = (x_0, y_0, z_0)$ the barycentre coordinates and $p_i = (x_i, y_i, z_i)$ the coordinates of the $i^{th}$ node in the isosurface at $c_{Y_F} = 0.1$.

The flame kernels are shown to be convected away from the ignition central location by turbulent motion in both cases. The barycentre displacement is non-negligible for both cases and highligths that the flame is not exactly centered in Fig. 7. In fact, the flame

<figure><img src="image_9.png" /><figcaption>Fig. 8: Flame barycentre relative to the edge length $L$ of the region of interest for each direction $(x, y, z)$ and for the 2 cases (superscripts 1 and 2 respectively).</figcaption></figure>

in case 1 moves towards the positive $\vec{x}$, $\vec{y}$, and $\vec{z}$ directions on more than 5 % of the region of interest length. The flame displacement is even more significant for case 2, with a barycentre position which reaches more than 15 % and 10 % of the zone of interest length in the $\vec{z}$ and $\vec{x}$ directions respectively at $t = 8.25$ ms from the flame ignition (represented by the last isocontour of Fig. 7).

This phenomenon directly affects the determination of the equivalent flame displacement speed and the wrinkling ratio. These quantities are derived from two equivalent flame radii $R_p$ and $R_a$, which are established using flame isocontours, as shown in Fig. 7, resulting either from DNS in this study or from 2D Mie-Scattering tomograhy on the X/Y plan for the experimental study.

Following the work of Zitouni et al. [7], $R_p$ and $R_a$ stand for equivalent radii based on a spherical flame shape as follows:

$$R_p = \frac{P}{2\pi} \qquad (11)$$

$$R_a = \sqrt{\frac{A}{\pi}} \qquad (12)$$

hence defining a circle of radius $R_p$ for which the perimeter $P$ is equal to the length of the isocontour of the progress variable; and a circle of radius $R_a$ for which the area $A$ is equal to the burnt gases area contained inside the isocontour of the progress variable.

An equivalent flame displacement speed with respect to the burnt gases $V_T$ can, thus, be calculated as:

$$V_T = \frac{dR_a}{dt} \qquad (13)$$

The wrinkling ratio $W$ can then be defined as the surface ratio of two spheres with radius $R_p$ and $R_a$ [22]:

$$W = \frac{R_p^2}{R_a^2} \qquad (14)$$

As such, the displacement of the flame barycentre location reported in Fig. 8 will result in a lack of precision on the $R_a$ and $R_p$ calculation, as the flame kernel
<figure><img src="image_10.png" /><figcaption>Fig. 9: Equivalent flame displacement speed over the equivalent flame radius from 2D (X/X cuts) and 3D data. Comparison with experiments from Zitouni et al. [7].</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Wrinkling ratio obtained from 2D (X/X cuts) and 3D data over the equivalent flame radius. Comparison with experiments from Zitouni et al. [7].</figcaption></figure>

can be off-plan and, thus, bias their values and, subsequently, both the turbulent flame displacement speed $V_T$ and wrinkling ratio $W$.

To avoid this shortcoming due to 2D experimental diagnostics, the flame displacement speed and the wrinkling ratio can be calculated in DNS with two equivalent spheres defined from 3D isosurfaces:

$$R_v = \left( \frac{3V}{4\pi} \right)^{\frac{1}{3}} \qquad (15)$$

$$R_s = \sqrt{\frac{S}{4\pi}} \qquad (16)$$

hence defining a sphere of radius $R_v$ for which the volume $V$ is equal to the volume inside the isosurface of progress variable $c_{YF} = 0.1$; and a sphere of radius $R_s$ for which the surface $S$ is equal to the surface of the isosurface. With this, $V_T$ and $W$ can be determined by substituting $R_a$ by $R_s$ in Eq. (13) and Eq. (14), respectively.

Thereby, Fig. 9 and Fig. 10 report the equivalent flame displacement speed and the wrinkling ratio with the aforementioned methods, i.e. using 2D experimental-like diagnostics and 3D diagnostics, over the equivalent flame radius $R_s$, in the X/Y plan for the experimental data and in the X/Y, X/Z, Y/Z plans for the DNS flame isocontours, and $R_v$ for the DNS flame isosurfaces, for the two cases. We observe a good agreement with experimental data for both $V_T$ and $W$ determined from flame isocontours (2D) and isosurfaces (3D). It is also noteworthy that the results

from the 3D flame isosurfaces are subject to fewer fluctuations than for the 2D data, especially in Fig. 9. Notably, the 2D results show slight differences with the 3D results and do not encompass them.

These two types of experimental diagnostics will be further compared with the same experimental configuration using 3D Mie-Scattering measurement technique to reconstruct the flame surface as in Ahmed et al. [23]. In fact, significant discrepancies are found between 2D and 3D experimental data. For instance in Zheng et al. [24], narrower distributions and lower values of the mean curvatures were observed in 3D compared to 2D diagnostics for lean hydrogen- and methane-air spherically expanding flames.

## 5. Conclusion

In this study, two 3D DNS of NH$_3$/air flames are performed to reproduce the experimental study of spherically expanding flames reported by Zitouni et al. [7]. These flames are within a HIT environment for two different turbulent intensities: $u' = 0.49$ m/s and $u' = 0.65$ m/s.

A HIT forcing method [17] was first implemented in AVBP [20], the in-house CERFACS solver, to generate and sustain turbulence. A good agreement is obtained for the two turbulent intensities, both in terms of turbulent kinetic energy and integral length scale.

Then, a semi-detailed reduced mechanism was derived from the detailed mechanism of Zhu et al. [10]
and shows high consistency in laminar flame speed and species profiles for the experimental conditions of interest.

Finally, a direct comparison between the experimental data and the DNS results is done for the equivalent displacement speed and the wrinkling ratio. A good agreement is achieved for both cases. Comparison between 2D and 3D data from DNS suggest a bias in the 2D Mie-scattering analysis. This stresses the importance of the current experimental research effort towards the use of 3D Mie-scattering tomography, where 3D isosurfaces can directly be recovered [23].

Future studies should focus on increasing the refined zone size in DNS to get more information on the NH₃/air turbulent flame propagation.

## Declaration of competing interest

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

## Acknowledgments

This work was supported by ANR ANR-22-CE50-0022 and performed using HPC resources from GENCI-TGCC, GENCI-CINES, and GENCI-IDRIS (Grant 2025-A0172B10157).

## References

[1] M. Höök and X. Tang. Depletion of fossil fuels and anthropogenic climate change - a review. *Energy Policy*, 52:797-809, 2013.

[2] Hideaki Kobayashi, Akihiro Hayakawa, K.D. Kunkuma, A. Somarathne, and Ekenechukwu C. Okafor. Science and technology of ammonia combustion. *Proc. Combust. Inst.*, 37:109-133, 2019.

[3] W. Yang, K.K.J. Ranga Dinesh, K.H. Luo, and D. Thevenin. Direct numerical simulation of turbulent premixed ammonia and ammonia-hydrogen combustion under engine-relevant conditions. *Int. J. of Hydrogen Energy*, 47:11083-11100, 2022.

[4] A.Y. Poludnenko and E.S. Oran. The interaction of high-speed turbulence with flames: Global properties and internal flame structure. *Combust. Flame*, 157: 995-1011, 2010.

[5] L. Xu, Q. Fan, X. Liu, X. Cai, A. A. Subash, C. Brackmann, Z. Li, M. Aldén, and X. Bai. Flame/turbulence interaction in ammonia/air premixed flames at high karlovitz numbers. *Proc. Combust. Inst.*, 39: 2289-2298, 2023.

[6] T. Tian, C. Song, H. Wang, C. Xu, K. Luo, and J. Fan. The effects of turbulence on the flame structure and no formation of ammonia turbulent premixed combustion at various equivalence ratios. *Fuel*, 332, 2023.

[7] S. Zitouni, P. Brequigny, and C. Mounaïm-Rousselle. Turbulent flame speed and morphology of pure ammonia flames and blends with methane or hydrogen. *Proc. Combust. Inst.*, 39(2):2269-2278, 2022.

[8] B. Galmiche, N. Mazellier, F. Halter, and F. Foucher. Turbulence characterization of a high-pressure high-temperature fan-stirred combustion vessel using Idv, piv and tr-piv measurements. *Exp. Fluids*, 55, 2014.

[9] S. Girhe, A. Snackers, T. Lehmann, R. Langer, F. Lofredo, R. Glaznev, J. Beeckmann, and H. Pitsch. Ammonia and ammonia/hydrogen combustion: Comprehensive quantitative assessment of kinetic models and examination of critical parameters. *Combust. Flame*, 267, 2024.

[10] Y. Zhu, H. J. Curran, S. Girhe, Y. Murakami, H. Pitsch, K. Senecal, L. Yang, and C. Zhou. The combustion chemistry of ammonia and ammonia/hydrogen mixtures: A comprehensive chemical kinetic modeling study. *Combust. Flame*, 260:113-229, 2024.

[11] D. G. Goodwin, H. K. Moffat, I. Schoegl, R. L. Speth, and B. W. Weber. Cantera: An objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes. https://cantera.org/, 2023. Version 3.0.0.

[12] K. Tazikawa, A. Takahashi, K. Tokuhashi, S. Kondo, and A. Sekiya. Burning velocity measurements of nitrogen-containing compounds. *J. Harzard. Mater.*, 155:144-152, 2008.

[13] S. Zhou, B. Cui, W. Yang, H. Tan, J. Wang, H. Dai, L. Li, Z. Rahman, X. Wang, S. Deng, and X. Wang. An experimental and kinetic modeling study on nh3/air, nh3/h2/air, nh3/co/air, and nh3/ch4/air premixed laminar flames at elevated temperature. *Combust. Flame*, 248, 2023.

[14] K. P. Shrestha, C. Lhuillier, A. A. Barbosa, P. Brequigny, F. Contino, C. Mounaïm-Rousselle ad L. Seidel, and F. Mauss. An experimental and modeling study of ammonia with enriched oxygen content and ammonia/hydrogen laminar flame speed at elevated pressure and temperature. *Proc. Combust. Inst.*, 38: 2163-2174, 2020.

[15] Q. Cazères, P. Pepiot, E. Riber, and B. Cuenot. A fully automatic procedure for the analytical reduction of chemical kinetics mechanisms for computational fluid dynamics applications. *Fuel*, 303, 2021.

[16] Q. Cazères. *Analysis and reduction of chemical kinetics for combustion applications*. Ph.D. Thesis, Institut National Polytechnique de Toulouse, 2021.

[17] V. Eswaran and S. B. Pope. An examination of forcing in direct numerical simulations of turbulence. *Computers and Fluids*, 16(3):257-278, 1988.

[18] Roberto Paoli and Karim Shariff. Turbulent condensation of droplets: Direct simulation and a stochastic model. *J. Atmos. Sci.*, 2008.

[19] Stephen B. Pope. *Turbulent Flows*. Cambridge university press, 2000.

[20] https://www.cerfacs.fr/avbp7x/.

[21] G. Lacaze, E. Richardson, and T. Poinsett. Large eddy simulation of spark ignition in a turbulent methane jet. *Combust. Flame*, 156:1993-2009, 2009.

[22] B. Renou, A. Boukhalfa, D. Puechberty, and M. Trinité. Effects of stretch on the local structure of freely propagating premixed low-turbulent flames with various lewis numbers. *Proc. Combust. Inst.*, 27: 841-847, 1998.

[23] P. Ahmed, B. Thorne, M. Lawes, S. Hochgreb, G. V. Nivarti, and R. Stewart Cant. Three dimensional measurements of surface areas and burning velocities of turbulent spherical flames. *Combust. Flame*, 233, 2021.

[24] Y. Zheng, P. Ahmed, and S. Hochgreb. 3d flame surface curvature analysis from reconstructed scanning across spherical expanding flames. *Proc. Combust. Inst.*, 40, 2024.