# Computational Cost Optimization of Hydrogen-Air Combustion for NOx Prediction

Quentin Buisson, Julien Leparoux, Stefano Puggelli, Maxime Leroy, Clément Mirat

► To cite this version:

Quentin Buisson, Julien Leparoux, Stefano Puggelli, Maxime Leroy, Clément Mirat. Computational Cost Optimization of Hydrogen-Air Combustion for NOx Prediction. ASME Turbo Expo 2024: Turbomachinery Technical Conference and Exposition, American Society of Mechanical Engineers, Jun 2024, Londres, United Kingdom. pp.GT2024-127585, 10.1115/GT2024-127585. hal-05169040

HAL Id: hal-05169040

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

Submitted on 22 Jul 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
# COMPUTATIONAL COST OPTIMIZATION OF HYDROGEN-AIR COMBUSTION FOR NOₓ PREDICTION.

Quentin Buisson<fnref n="1" />,<fnref n="2" />*, Julien Leparoux<fnref n="3" />, Stefano Puggelli<fnref n="3" />, Maxime Leroy<fnref n="1" />,<fnref n="3" />, Clément Mirat<fnref n="1" />

## ABSTRACT

The switch from fossil fuels to hydrogen represents a break-
through path towards the decarbonization of the aeronautical
sector. However, compared to hydrocarbons, hydrogen has very
different characteristics. This challenges the current state of the
art of numerical methods that must be re-evaluated in terms of
accuracy and computational cost for aeronautical applications.
Large Eddy Simulations (LES) applied to hydrogen combustion
can easily become too costly with respect to the industrial con-
straints.

In this work, a study to reduce the time-to-solution in LES
is performed. The selected configuration is a two-staged swirled
combustor. The reference LES setup includes a detailed chemical
description by using a 21-species, 64-reactions chemical scheme
(including NOₓ) and the thickened flame model to handle turbu-
lent combustion. From this, different numerical combinations are
investigated for each phase of the simulation (from non-reactive
flow to flame stabilisation) as well in terms of chemistry descrip-
tion to minimise the LES return time. Ultimately, the switch to-
wards a lighter 9-species, 12-reactions chemical scheme and the
inclusion of NOₓ with a tabulated formalism is retained as optim-
ized setup. This choice determines all in all a 50% reduction in
terms of computational cost with respect to the reference setup,
keeping the experimental flame topology fairly well retrieved.
The general trend for NOₓ emissions at the burner outlet is as
well correctly captured, even if the chosen methodologies differ
on the absolute values.

## INTRODUCTION

As a response to the concern about its environmental impact,
the aeronautical industry is investigating new ways to reduce its
pollutant emission levels. Among the paths being explored, the
use of alternative fuels looks promising. In this context, sus-
tainable aviation fuels (SAF) and hydrogen are solutions under
investigation. Hydrogen in particular represents at the same time
an appealing and disruptive prospect, having only NOₓ as major
pollutant. However, it entails several challenges from a tech-
nological and a computational point of view. Concerning com-
bustion, thermochemical properties of hydrogen are significantly
different from kerosene. For instance, the laminar flame speed is
higher, the thermal flame front thinner, the chemistry stiffer. Hy-
drogen is also a very volatile molecule so it diffuses rapidly in air
and is subject to differential diffusion effects [1]. All in all, ac-
tual modelling practices for combustion are challenged and need
to be evaluated considering these specific properties.

Recently, several works have tried to tackle this problem. In
[2], the Thickened Flame LES model (TFLES) was compared to
the Flamelet Generated Manifold (FGM) on a hydrogen-air jet in
crossflow flame [3]. A similar study was carried out in [4] on the
swirled flame experimentally investigated in [5] increasing the
percentage of hydrogen in the fuel. Both works show that TFLES
seems capable of better handling the impact of strain and differ-
ential diffusion with respect to FGM and this has an effect on
the flame stabilization. Similarly, in [6] an atmospheric swirled
H₂/air flame was computed using TFLES. Once again, the ma-
jor flame characteristics, such as flame structures, velocity fields

<footnote n="1">¹ Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91190, Gif-sur-Yvette, France</footnote>

<footnote n="2">² Safran Aircraft Engines, Site de Villaroche, Rond-Point René Ravaud, 77550 Moissy Cramayel, France</footnote>

<footnote n="3">³ Safran Tech, Rue des Jeunes Bois, Châteaufort, 78114 Magny-Les-Hameaux, France</footnote>

<footnote>*Corresponding author: quentin.buisson@centralesupelec.fr</footnote>
and heat release rate distributions were correctly captured. LES was then used to better understand the transition from lifted to anchored flame shapes.

Nevertheless, several open points still need to be addressed concerning the application of TFLES to hydrogen and the computational cost of this type of methodology. Considering the operating conditions usually encountered in aero-engines, in order to have a sufficient resolution of the flame front and eventually of the non-premixed flame zones, the $CPU_{cost}$ can become a major limit of this approach in an industrial framework.

This consideration motivates the present study that is aimed at testing different strategies in order to optimize the $CPU_{cost}$ of TFLES for H₂ combustion. To do so, LES of the two-staged swirled hydrogen-air burner from Leroy et al. [7] are performed to evaluate different numerical setups.

More specifically, a reference setup based on a detailed transported chemistry and TFLES is initially presented, serving as reference for the rest of the study. Then, following the work of Leparoux et al. [8], a tabulated formalism coupled with TFLES is employed. In [8], it was shown that this approach, albeit able of retrieving the tendency of geometrical variations on pollutants with a very cheap $CPU_{cost}$, fails in predicting some local phenomena of hydrogen combustion such as differential diffusion and therefore, the absolute values of $NO_x$. This approach will be therefore used as a low-cost reference methodology, which can eventually be interesting in an industrial context. Finally, a novel procedure, based on the combination of a skeletal mechanism for H₂/air combustion and a tabulated model for $NO_x$ prediction, is proposed as a promising compromise between precision and computational cost. The accuracy and $CPU_{cost}$ of these three methodologies are evaluated and compared against experimental data.

The paper unfolds as follows: Sec. 1 firstly presents the experimental test case selected in this study as well as the operating conditions. Sec. 2 then describes the numerical setup used for both detailed and tabulated chemistry descriptions. Sec. 3 covers the cost optimization study. Results obtained with the different methodologies are compared in Sec. 4 with a focus on $NO_x$ prediction.

# 1 EXPERIMENTAL TEST CASE

**Geometrical configuration** The coaxial burner configuration experimentally studied by Leroy et al. [7] is here selected as reference test case and is sketched in Fig. 1. It is composed of two concentric injection nozzles. The inner nozzle (hereinafter referred to as channel “1”) has a diameter $d_1$ = 10 mm and delivers a hydrogen-air mixture at an equivalence ratio $\phi_1$ ranging from 2 to a non-premixed injection. The outer nozzle (hereinafter noted as channel “2”) has a diameter $d_2$ = 20 mm and injects air to dilute the mixture down to a global equivalence ratio $\phi_{tot}$ ranging from 0.35 to 0.40. The combustion chamber is 215 mm

long and has a 80 × 80 mm squared section decreasing down to 60 × 60 mm at the burner exit. Both injection nozzles are swirled. Channel 1 is swirled using swirl vanes, housed in the inner channel, with an angle of 30° giving the internal flow a geometrical swirl number $S_1$ = 0.29. Channel 2 is swirled using four tangential aerodynamic ports located in the injector housing. By varying the flow rate injected tangentially, channel 2 can yield a swirl number $S_2$ between 0 to 0.8. The experimental work includes several optical diagnostic techniques such as direct visualization and OH* chemiluminescence. Exhaust gas composition measurements are also performed at the exit of the chamber using a HORIBA VA-5111 multicomponent gas analyzer with a VS-5000 series sampling unit. According to the manufacturer, the NO sensor has a precision of ±2.5 ppm. Full diagnostics details are given in [7].

<figure><img src="image_2.png" /><figcaption>**FIGURE 1:** Schematic representation of the coaxial burner geometry adapted from Leroy et al. [9]</figcaption></figure>

**Operating conditions** Two operating cases are selected from the experimental database. Hereinafter, these cases are referred to as A and B and are summarized in Tab. 1. As shown in the table, the main difference is related to the swirl number for channel 2 that is zero for case A and 0.8 for case B. The rest of the operating conditions is identical.

# 2 NUMERICAL SETUP

Simulations were performed using the unstructured low-Mach number code YALES2 [10] with a variable density solver.

The time integration uses a fourth-order accurate TFV4A scheme [11], which is a hybrid between a centered and a Lax-Wendroff scheme. Spatial discretization relies on fourth-order accurate central finite-volume schemes. The subgrid Reynolds stresses are included using the WALE model [12]. A specific
<table><caption>TABLE 1: Numerical operating points</caption><thead><tr><th>Case</th><th>A</th><th>B</th></tr></thead><tbody><tr><td>Pressure P (Pa)</td><td>101325.0</td><td>101325.0</td></tr><tr><td>Air temperature T<sub>a</sub> (K)</td><td>298.0</td><td>298.0</td></tr><tr><td>Fuel temperature T<sub>f</sub> (K)</td><td>298.0</td><td>298.0</td></tr><tr><td>Swirl number S<sub>1</sub></td><td>0.29</td><td>0.29</td></tr><tr><td>Swirl number S<sub>2</sub></td><td>0</td><td>0.8</td></tr><tr><td>Premixed equivalence ratio &phi;<sub>prem</sub></td><td>2.0</td><td>2.0</td></tr><tr><td>Global equivalence ratio &phi;<sub>tot</sub></td><td>0.35</td><td>0.35</td></tr><tr><td>H<sub>2</sub> mass flow rate m&#775;<sub>H<sub>2</sub></sub> (g/s)</td><td>0.0839</td><td>0.0839</td></tr></tbody></table>

treatment is applied on transported scalars to prevent local extrema and sharpen discontinuities. A MUSCL scheme [13] is associated to a second-order flux Superbee limiter for scalar time integration. The time integration scheme for scalars is a RK3 scheme. The stiffness of the chemistry is handled with the CVODE stiff integrator [14].

For all the calculations reported in this study, the interactions between chemistry and turbulence are modeled using the TFLES model [15]. The model is based on an artificial thickening of the flame front while ensuring correct laminar flame burning velocity retrieval [15]. In order to integrate the effects of the flame wrinkling at the subgrid scales, the efficiency function as formulated by [16] is introduced. A flame index is used to define where the thickening is applied. In this study, Takeno index [17], denoted $\xi$, is employed to distinguish premixed ($\xi = 1$) and non-premixed ($\xi = -1$) flame zones ($\xi$ being 0 elsewhere). The thickening procedure is only applied where $\xi = 1$. The definition of the flame sensor $S$ is given in Secs. 2.1 and 2.2 for transported and tabulated chemistry respectively. No subgrid scale model is applied to non-premixed regions of the flame assuming that mesh description is sufficient to ensure the resolution of the flame behaviour. A mesh sensitivity study is carried out to ensure that the conclusions obtained do not depend on the mesh resolution (see Sec. 4.2.3). This numerical framework is used in combination with two kinetics descriptions that are detailed in the next two sub-sections.

## 2.1 Reference detailed chemistry setup

The reference detailed chemistry setup uses the reduced San Diego mechanism (UCSD) [18]. It contains 21 species and 64 reversible reactions allowing $NO_x$ prediction. Species diffusion are considered individually using a mixture averaged approach [19]. In this context, the flame sensor $S$ is defined equal to one when the source term for $H_2O \omega_{H_2O} > 10\% \omega_{H_2O,max}$ and $S = 0$ elsewhere. This methodology was applied in [8] for a similar staged

$H_2$/air flame and is considered as the reference for this study.

## 2.2 Tabulated chemistry setup

In the tabulated chemistry context, the TFLES-FPI approach is employed, similarly to what proposed in [20]. Thermochemical variables are therefore tabulated as a function of the mixture fraction $Z$ and the progress variable $Y_c = Y_{H_2O}$. A set of 1D laminar freely-propagating premixed flamelets is generated using the reduced San Diego mechanism described in Sec. 2.1. NO and $NO_2$ are as well tabulated throughout the whole computational domain using the NOMANI model [21]. This model allows to track NO formation within both the inner flame front and the post-flame region through the additional transport of $\tilde{Y}_{NO}$. The application of the NOMANI model to H$_2$-air combustion was already validated in a previous study [8]. A correction factor $\alpha_c$ proposed by Mercier et al. [22] is used to account for the differential diffusion in the direction normal to the flame front.

## 2.3 Optimized chemistry setup

The optimized setup for computational cost, which will be fully presented in Sec. 3, relies on the Boivin mechanism [23] that includes 9 species and 12 reactions. A substantial reduction on the simulation cost with respect to UCSD is therefore expected as already shown in [24].

However, the Boivin scheme does not include species and reactions for $NO_x$ formation. Therefore, to enable $NO_x$ prediction, a hybrid approach is here proposed. In particular, the NOMANI model [21] is coupled with the Boivin mechanism. To the authors knowledge, this is the first time that a skeletal mechanism is coupled to a tabulation approach to predict $NO_x$ for H$_2$-air combustion.

## 2.4 Computational domain and mesh description

The computational domain is composed of the injector and the combustion chamber described in Sec. 1. A square-shaped atmosphere downstream of the combustor outlet is also added to handle potential gas recirculation through the exhaust. The swirl in channel 1 and 2 is not imposed, but the swirlers geometries are directly included in the computations. Two mesh resolution are tested. Initially, both cases A and B are represented using the same tetrahedral mesh, counting about 46M cells. This mesh is denoted as MR0 and it is shown in Fig. 2. A uniform cell size of 1.0mm is applied in the combustion chamber, while a finer description is imposed in the injector parts. Similarly to [8], a mesh refinement is then performed. Dedicated meshes are then obtained for case A and case B, denoted MR1 (see Fig. 3), counting about 93M cells. Meshes were obtained by performing an on-the-fly mesh adaptation as described in [25].
<figure><img src="image_3.png" /><figcaption>FIGURE 2: Initial mesh MR0 for both cases A and B</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3: Refined mesh MR1 zoomed on the combustor (a) Case A (b) Case B.</figcaption></figure>

### 3 OPTIMIZATION OF THE LES RETURN TIME

Figure 4 shows the workflow of the LES performed on case A using the reference chemistry setup (i.e. the reduced San Diego mechanism) described in Sec. 2.1. The figure shows the order of magnitude of the simulation times required for each step of the calculation as percentage of the total CPU time for mesh MR0. A four-step process is employed:

1. **Initialization of the flow field:** counting 10% of the total simulation time, it is a non-reactive simulation mainly establishing the velocity flow field as well as the air-fuel mixing to reach a sufficient equivalence ratio in the domain;

2. **Ignition:** counting 5% of the total simulation time, this step allows the transition to a reactive flow field with a kernel energy deposition;

3. **Stabilization of the reactive flow:** counting 55% of the total simulation time, here, the flame and the reactive flow-field are stabilized. It is the most expensive step in the process due to the time required to converge the mass balance in the domain;

4. **Statistical accumulation:** once mass balance is reached, a sufficient accumulation time must be ensured to converge statistics. This step represents around 30% of the total simulation time.

<figure><img src="image_5.png" /><figcaption>FIGURE 4: LES workflow for the simulation of case A using the mesh MR0.</figcaption></figure>

Most of the simulation time is spent in steps 3 and 4, mainly due to scalars computations. In order to assess the impact of the different numerical methodologies on the computational cost, two metrics are used:

* **Return Time To Solution (RTTS***^2^): Expressed as following RTTS* = $\frac{t_{sim}}{WCT \cdot N_{cores}}$ where t_{sim} is the simulation time in ms, WCT is the Wall Clock Time in µs, N_{cores} is the number of cores. This metric offers an overview of the performance and cost of simulation from a user perspective.

* **CPU cost per simulation step (Cost|step):** To be comparable, computational cost was computed as follows $Cost|step = \frac{t_{sim|step} \cdot 24 \cdot 1000}{RTTS^*|step}$ where t_{sim|step} is the simulation time required for the step in ms and RTTS^*|step is the Return Time To Solution of the step in ms/24h/1000 procs.

The Flow-Through-Time (FTT) is as well introduced and defined as $\tau_c = L_c/U_{Bulk}$ where L_c is the combustion chamber length and U_{Bulk} refers to the mean axial velocity at the exit of channels. Here, the FTT is about 4.3 ms for both cases. The main objective of this section is to propose an optimized numerical setup for the detailed chemistry case that reduces both the computational cost of each step and the return time of the entire workflow. For the sake of brevity, the optimization process is presented only for case A on the mesh MR0, but the same conclusions hold for case B.

#### 3.1 Non-reactive flow initialization

The main issue for this first step is the simulation time required to reach a sufficient air-fuel mixture in the combustor prior to the ignition step. Although it is a non-reactive flow, the reference setup of Sec. 2.1 includes the transport of all the species leading to an non-negligible CPU cost. To ensure the convergence of the flow stabilization, two metrics are considered: the mean $L^2$ norm of the velocity and the volume integral of the mixture fraction within the combustor. These two metrics lead to an
initialization time of about 5 FTT.

To optimize this step, the incompressible solver with a constant
density available in YALES2 is used instead of the reference nu-
merical setup. This avoids the resolution of species transport.
Only a transported passive scalar representing the mixture frac-
tion is included. This optimization reduces the computational
cost of this step by 97% compared to the reference case and by
8% in terms of total cost.

## 3.2 Ignition

In the reference setup, the non-reactive initialization in
Sec. 3.1 includes the transport of species. Starting from the ini-
tialization, the ignition is straightforward with the inclusion of an
ignition kernel.

In the optimized setup, the mass fraction of the species is initial-
ized in the domain at the beginning of the ignition step based on
the passive mixture fraction. The flow is then ignited by placing
a toroidal kernel of burnt gases at equilibrium located at the exit
of channel 1. The physical properties of the mixture are then cal-
culated at the first time step using the built-in chemistry solver of
YALES2. All the procedure is performed on-the-fly.

## 3.3 Stabilization of the reactive flow

This step consists in converging the mass balance in the do-
main before statistical accumulation. This is a relatively long
process, especially when the flow has stagnant regions making
the residence time distribution in the chamber non-uniform. The
mean residence time is then the most relevant parameter to estim-
ate the simulation time needed. In the present case, the mass bal-
ance in the domain is stabilized after 45 FTT for case A. The cost
reduction in the optimized setup occurs mainly on scalars com-
putation due to the lower number of species transported, leading
to a 54% reduction in the computational cost of this step and 31%
of the total cost.

## 3.4 Statistical accumulation

Similarly to the stabilization step, the statistical accumula-
tion is mainly impacted by the cost of scalars computation. Stat-
istical accumulation is performed on 16 FTT for case A. This
leads to a 54% reduction in the computational cost of this step
and 15% of the total cost.

## 3.5 Comparison of computational costs

Figure 5 summarizes the computational cost for the three
numerical setups investigated.

First, we can observe that the reference setup applied on
mesh MR0 has a cost of about 1 million CPUh. The computa-
tional cost of this methodology on the refined mesh MR1, presen-
ted in Sec. 2.4, would be at least doubled. Considering the limits

<figure><img src="image_6.png" /><figcaption>**FIGURE 5:** Computational cost comparison (in CPUh) of the reference, optimized and TFLES-FPI methodologies for a complete study of case A.</figcaption></figure>

of computational resources at the time of the study, this compu-
tation was not possible in the present investigation. However, it
is straightforward to understand that this setup cannot be used in
a design process. Instead, the TFLES-FPI setup, shows a very
affordable computational cost of about 100 kCPUh, a factor of
ten compared to the reference setup. This makes this methodo-
logy very cheap despite being not completely able of predicting
some local phenomena as described in Sec. 2.2.

Finally, the optimized setup shows a significant reduction of
the total computational cost, which is improved by a factor higher
than 2. The cost of this methodology closer to TFLES-FPI cost,
even if a factor of 4 remains. This makes it an attractive setup
from an industrial perspective. The refined mesh MR1 doubles
the total CPU cost with the optimized setup. However, it is worth
pointing out that the cost on mesh MR1 with the optimized meth-
odology is as high as the reference one with MR0. The optimized
setup can therefore be a good compromise to run simulations on
highly refined meshes. In the next section, a comparative ana-
lysis is performed to evaluate the accuracy of the different nu-
merical methodologies.

# 4 RESULTS COMPARISON AND DISCUSSION

## 4.1 Presentation of results with the reference setup

Results on cases A and B are firstly analyzed using the refer-
ence numerical setup of Sec. 2.1. Figure 6 shows the main mean
flow fields obtained using the mesh MR0.

* Case A: This case shows a jet-like behaviour with an elong-
ated flame that propagates inside the combustor. The flame shows a predominantly axial dynamic since 81.6% of the air mass flow rate is injected through the non-swirled channel 2. The effect of the swirl in channel 1 ($S_1 = 0.29$) is visible in the inner shear layer between channel 1 and 2. A double flame front configuration is retrieved. This is highlighted in the mixture fraction $\bar{Z}$ and the temperature $\bar{\omega}_T$ fields. A first conical premixed flame front is present at the exit of channel 1. A long diffusion flame front is then stabilized in the shear layer generated by the flow coming from channel 2. This second flame front is at stoichiometric conditions. High temperatures are therefore locally generated and high $NO_x$ formation is expected.

* **Case B:** This case exhibits a similar premixed flame front at channel 1 exit, but a more compact diffusion flame. The velocity field $\bar{U}$ shows that the flow topology inside the combustor is changed with the addition of swirl in the channel 2 ($S_2 = 0.8$): a large inner recirculation zone downstream of the premixed flame front is created. The premixed flame front appears very similar to case A since $S_1$ is identical between both cases. However, the non-premixed flame front is significantly shortened due to the presence of the inner recirculation zone and the high temperature zone is shrunk as shown by the temperature flow field $\bar{T}$. This suggests a reduction of $NO_x$ production with respect to case A.

Experimental images of inverse Abel transform of $OH^*$ chemiluminescence are also presented in Fig. 6. The images show a good correlation with $\bar{\omega}_T$ fields. Premixed flame heights are found in good agreements with experiments. Flame regions with high $OH^*$ signal intensity in experiments correlate well with high $\bar{\omega}_T$ values in simulations. Therefore, for both cases, flame topologies obtained with the reference setup qualitatively agree with the experiments. The optimized and TFLES-FPI setups, not shown here for the sake of brevity, retrieve the same flame shape as the reference case. Ultimately, the three numerical setups are qualitatively able to replicate the flame topology. The comparison in the next section is then focused on the capabilities of the optimized methodology in predicting $NO_x$ emission levels compared to the reference and TFLES-FPI setups.

## 4.2 Numerical setups comparison

To be comparable with experimental values, $NO_x$ are post-processed from numerical simulations following the definition of $X_{NO_x,dry}$ measured by Leroy et al. [9]. Expression to compute the molar fraction $X_{NO_x}$ from the mass fraction reads:

$$X_{NO_x} = \left( \frac{Y_{NO}}{M_{NO}} + \frac{Y_{NO_2}}{M_{NO_2}} \right) \bar{M}, \quad (1)$$

where $\bar{M}$ is the average molar mass of the flue gases, $M_{NO}$

<figure><img src="image_7.png" /><figcaption>**FIGURE 6:** Average velocity magnitude $\bar{U}$, average temperature $\bar{T}$, average mixture fraction $\bar{Z}$ and average heat release rate $\bar{\omega}_T$ fields retrieved using the reference setup compared with experimental images of inverse Abel transform of $OH^*$ chemiluminescence for (a) case A (b) case B.</figcaption></figure>

<table><caption>**TABLE 2:** $X_{NO_x}$ (in ppm) levels retrieved experimentally, using the reference setup, the optimized setup and the TFLES-FPI setup.</caption><thead><tr><th>X<sub>NO<sub>x</sub></sub><br>(ppm)</th><th>Exp.</th><th>Ref.</th><th>Optim.</th><th>Optim.</th><th>TFLES-FPI</th></tr></thead><tbody><tr><td>Case A</td><td>17.3 ±2.5</td><td>22.0</td><td>21.2</td><td>27.6</td><td>14.3</td></tr><tr><td>Case B</td><td>10.2 ±2.5</td><td>11.9</td><td>9.9</td><td>12.3</td><td>5.1</td></tr></tbody></table>

refers to the molar mass of NO and $M_{NO_2}$ refers to the molar mass of $NO_2$. In simulations, according to experiments, $X_{NO_x}$ was averaged in time at the exhaust gas analyzer position located at the exit plane of the combustion chamber.

Experimental and numerical values are presented in Tab. 2. Experimentally, case A exhibits a higher value of $NO_x$ than case B. The absence of central recirculation zone promotes $NO_x$ production because of the long diffusion flame. This hierarchy is
well retrieved by LES with the three methodologies. However, discrepancies exist:

1. The reference setup overestimates $X_{NO_x}$ by 27% for case A. This overestimation is consistent with the findings on a similar burner studied in [8]. For case B, $X_{NO_x}$ is within the uncertainty window. This quantitative comparison shows that this methodology, albeit intensive from a computational point of view, allows to quantitatively replicate the experimental evidences when the diffusion flame region is small.

2. The optimized setup achieves similar, but smaller, $X_{NO_x}$ values compared to the reference setup. This indicates that this setup does not significantly alters the NO$_x$ predictions compared to the reference setup. However, it still provides overestimated values for case A.

3. The optimized setup using the refined mesh MR1 shows significant increase in $X_{NO_x}$ compared to MR0. This results in a large overestimation for case A while it remains satisfactory for case B. A possible explanation refers to the diffusion flame region where no sub-grid scale model is applied and can lead to a non-linear phenomenon with the cell size refinement. Another explanation is about the assumption of adiabatic walls in the simulation.

4. The TFLES-FPI setup underestimates the experimental levels with error as high as 50% for case B. Nonetheless, the tendency between the geometries is well retrieved. This can already be an interesting information from an industrial point of view.

To go further, in-depth comparisons are performed, firstly between the optimized and reference setups and then between the optimized and TFLES-FPI setups to better understand the differences on NO$_x$ between these methodologies.

### 4.2.1 Comparison of reference and optimized cases

Figure 7 shows contours of the mass fraction $Y_{NO}$ obtained with both methods overlaid on premixed and non-premixed flame regions. The flame regime is identified by using the Takeno index as described in Sec. 2.

For case A, albeit some differences appear in the premixed flame region, the main interest is in the diffusion flame zone where NO production is mainly located. There, contours of the $Y_{NO}$ mass fraction are found very similar between reference and optimized setups. One notable difference is located on the axis of the combustion chamber, upstream of the diffusion flame region, where a finger-like shape is present in the reference setup and is not recovered in the optimized setup. This difference suggests probably that the two methods handle differently the effects of preferential diffusion that can be important downstream of the premixed flame front as already shown in [8].

For case B, both methods correctly reproduce NO production zones in the diffusion flame region: one attached to the in-

<figure><img src="image_8.png" /><figcaption>**FIGURE 7:** Contours of the mass fraction $Y_{NO}$ (solid black lines for the reference setup and solid blue lines for the optimized setup) overlaid on premixed flame region (red dots) and diffusion flame region (orange dots). (Top) case A (Bottom) case B. Jector outlet in the shear layer between the two channels and the other at the diffusion flame front.</figcaption></figure>

To go further in the comparison, Fig. 8 shows the joint Probability Density Function (PDF) of the NO source term $\overline{\omega}_{NO}$ against the mean mixture fraction $\bar{Z}$ and the progress variable $\overline{Y_C/Y_{C,eq}}$, for case A. The reference setup is plotted in solid black lines while the optimized setup is in solid blue lines. The same plots are shown for case B in Fig. 9. PDFs are computed using a Kernel Density Estimation method (KDE) from mesh nodes on the meridian plane.

First, looking at case A (see Fig. 8), PDFs of $\overline{\omega}_{NO}$ overlap well in general between setups. Two distinct zones of NO production are located along the $\bar{Z}$ axis in lean and rich mixtures, respectively. While a good agreement appears in the rich region, higher $\overline{\omega}_{NO}$ values are attained with the reference setup. This behaviour is more pronounced for lean mixtures: the reference setup presents a higher NO production with $\overline{\omega}_{NO,max} \simeq 0.006s^{-1}$ whereas $\overline{\omega}_{NO,max} \simeq 0.004s^{-1}$ with the optimized methodology. However, $\overline{\omega}_{NO}$ is slightly higher around the stoichiometry with the optimized setup than with the reference one. NO production zones are also shown to overlap well along the $\overline{Y_C/Y_{C,eq}}$ axis. Notable differences are still found for $\overline{Y_C/Y_{C,eq}} \in [0.9; 1.0]$ where $\overline{\omega}_{NO}$ is significantly higher with the reference setup than with the optimized setup.

Similar conclusions can be drawn from Fig. 9 for case B. The most important interesting information is that $\overline{\omega}_{NO}$ distributions differs from those shown in Fig. 8. The distribution along $\bar{Z}$ still shows a production zone for lean regions ($\bar{Z} \in [0.0; 0.017]$) and another for rich regions ($\bar{Z} \in [0.04; 0.06]$). As for case A, $\overline{\omega}_{NO}$ distributions overlap well for both zones. Differences can
<figure><img src="image_9.png" /><figcaption>**FIGURE 8:** Case A: Joint Probability Density Function of the mean $Y_{NO}$ source term density plotted against (Left) mean mixture fraction $\bar{Z}$ (Right) mean progress variable $\bar{Y}_C$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 9:** Case B: Joint Probability Density Function of the mean $Y_{NO}$ source term density plotted against (Left) mean mixture fraction $\bar{Z}$ (Right) mean progress variable $\bar{Y}_C$.</figcaption></figure>

be highlighted, especially in the lean region where $\dot{\omega}_{NO,max} \simeq 0.003s^{-1}$ with the reference methodology while $\dot{\omega}_{NO,max} \simeq 0.0015s^{-1}$ with the optimized methodology. Still, the range of mixture fraction where distributions do not overlap is quite narrow thus showing the accuracy of the optimized methodology. The PDFs plotted along $Y_C/Y_{C,eq}$ lead to the same conclusions. Ultimately, this analysis demonstrates that the optimized methodology is able to correctly replicate the quantitative evolution of $NO_x$: the zone of production in the space of mixture fraction and progress variable are retrieved even if some local differences exist on the absolute maximum values. This can be possibly ascribed to a slightly different response of the NOMANI model in capturing the dependence of $NO_x$ production to strain and to the intermediate species like O and OH due to the tabulation [26].

### 4.2.2 Comparison of optimized and TFLES-FPI cases
Figure 10 shows contours of the $NO_x$ mass fraction with the two setups for case A and B.

For case A, both methods globally find similar NO production zones: one attached to the injector outlet and another at the tip of the diffusion flame front. Nevertheless, the contours differ significantly. The production zone found at diffusion flame tip appears more confined around the combustor axis with TFLES-FPI than for the optimized setup. Furthermore, the production zone is less extended in the axial direction.

Differences are even more pronounced for case B since TFLES-FPI misses the production zone attached to the injector
<figure><img src="image_11.png" /><figcaption>**FIGURE 10:** Contours of the mass fraction $Y_{NO}$ (solid blue lines for the optimized setup and solid green line for TFLES-FPI) overlaid on premixed flame region (red dots) and diffusion flame region (orange dots). (Top) case A (Bottom) case B.</figcaption></figure>

outlet. It appears that NO is produced just in presence of the diffusion flame.

Figure 11 shows the PDFs of $\overline{\omega}_{NO}$ plotted against the mean mixture fraction $\bar{Z}$ for case A and case B. For case A, the shape of the PDFs obtained differs between the optimized and TFLES-FPI setups. The TFLES-FPI setup presents a continuous NO production ranging from $\bar{Z} = 0$ to $\bar{Z} = 0.58$ and its PDF is characterised by a plateau at $\overline{\omega}_{NO} \approx 0.0015s^{-1}$ for $\bar{Z} \in [0.01; 0.045]$. On the contrary, the optimized setup presents two distinct regions of significantly higher NO production.

More significant differences can be highlighted on case B. $\overline{\omega}_{NO}$ levels are up to 2.4 times lower with the TFLES-FPI setup. This demonstrates again that the optimized setup is capable of retrieving sharper $\overline{\omega}_{NO}$ distribution over the computational domain while TFLES-FPI shows a uniform distribution. These observations help explaining why $X_{NO_x}$ presented in Tab. 2 is much lower with TFLES-FPI.

### 4.2.3 Influence of the mesh resolution

Finally, a mesh influence is performed to evaluate the effect of the mesh resolution on the previous conclusions. Mesh MR1 is therefore tested (see Fig. 3). Only results obtained using the optimized methodology are shown since the reference methodology was not affordable in terms of computational costs. Figure 12 and Fig. 13 show the PDFs of $\overline{\omega}_{NO}$ plotted against the mean mixture fraction and the progress variable for MR0 and MR1 meshes for case A and case B, respectively.

For case A, Fig. 12 exhibits a broadened distribution of $\overline{\omega}_{NO}$ with the MR1 mesh compared to the MR0 mesh. On all values of $\bar{Z}$ between 0.0 and 0.06, a larger range of $\overline{\omega}_{NO}$ is obtained with

large increase of $\overline{\omega}_{NO}$ for $\bar{Z} < 0.01$. The refinement from MR0 to MR1 enlarges the range of $\bar{Z}$ where high NO production occurs. This analysis shows that the switch from MR0 to MR1 increases $X_{NO_x}$ by notably enhancing its production rate in the diffusion flame. Another significant difference between MR0 and MR1 is that $\overline{\omega}_{NO}$ is now non-zero on the whole range of $\overline{Y_C}/Y_{C,eq}$ from 0.0 to 1.0 even though it remains small with $\overline{\omega}_{NO} < 0.0005s^{-1}$ for $\overline{Y_C}/Y_{C,eq} \in [0.15; 0.8]$. This means that the switch to MR1 also increases $X_{NO_x}$ by triggering the production of NO over the whole range of $Y_C$.

Interestingly, Fig. 13 shows that PDFs obtained with MR1 on case B have shapes very similar to those of case A. The analysis of results obtained on MR0 in Sec. 4.2 showed that NO production regions are located only in lean and rich zones of the mixture. Instead, the mesh refinement allows NO production across all $\bar{Z}$ values from 0.0 to 0.06. For all $\bar{Z}$ values, $\overline{\omega}_{NO}$ is significantly higher with MR1 than with MR0. For example, $\overline{\omega}_{NO,max}$ is still obtained for lean mixture but its value doubles from $0.0015 s^{-1}$ with MR0 to around $0.003 s^{-1}$ with MR1. Looking at the distribution along $\overline{Y_C}/Y_{C,eq}$, the same behaviour is retrieved. The overall distribution of NO production zones across the flame is thus drastically modified. These observations help to explain why $X_{NO_x}$ presented in Tab. 2 are significantly higher with the MR1 than with MR0.

## 5 CONCLUSIONS AND PERSPECTIVES

In the present work, the optimisation of a numerical setup using a detailed chemistry description for hydrogen-air combustion is performed and validated on a two-staged, swirled hydrogen-air burner.

The cost optimization is performed by simplifying the numerical setup using the 12-step skeletal mechanism of Boivin and using the NOMANI model to retrieve $NO_x$ emissions. This optimization determines a reduction of 51% in computational cost.

The resulting optimized methodology is compared to the reference methodology in terms of $NO_x$ emissions. Results show an overall good agreement between these methodologies on NO production zones for both simulated cases with slight discrepancies around the diffusion flame front. Quantitative analysis of the mean NO source term $\overline{\omega}_{NO}$ aimed at understanding those differences are then performed. Results showed that both methodologies give similar results on the major part of the mixture fraction and progress variable spaces. Overall, $NO_x$ levels achieved for both methodologies are close to each other. The experimental value is well retrieved for a compact flame (ie. small diffusion flame regions) whereas levels are overestimated for flame with a longer diffusion flame front. The mesh refinement study showed an important effect mainly caused by a wider spatial distribution of NO production, especially in the diffusion flame leading to large overestimation of $NO_x$ levels for case A. An explanation for this phenomenon is that the assumption of fully adiabatic walls
<figure><img src="image_12.png" /><figcaption>**FIGURE 11:** Comparison between the optimized setup and TFLES-FPI on the joint Probability Density Function of the mean $Y_{NO}$ source term density plotted as function of the mean mixture fraction (Left) case A (Right) case B.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**FIGURE 12:** Case A: Comparison between MR0 and MR1 on the joint Probability Density Function of the mean $Y_{NO}$ source term density plotted against (Left) mean mixture fraction $\bar{Z}$ (Right) mean progress variable $\bar{Y}_C$.</figcaption></figure>

in the combustion chamber may be too strong. Having wall heat transfer could lead to a longer diffusion flame front, thus lowering local temperature and $NO_x$ production. Another explanation is that no sub-grid scale model is applied to the diffusion flame region which can lead to a non-linear phenomenon with the cell size refinement. These aspects will be investigated in the next future starting from the basic conclusions obtained in the present work.

The TFLES-FPI methodology is also tested showing a significant gap with the experimental values and the reference methodology. The production of $NO_x$ is strongly underestimated both in terms of spatial distribution and absolute levels. However, the approach remains able of capturing the hierarchy of pollutants levels with a computational cost really interesting from an indus-

trial point of view.

Overall, the cost optimization performed in this study enables an important decrease in the $CPU_{cost}$ and satisfactory results in terms of pollutants. The proposed methodology can be employed to have efficiently an accurate estimation of $NO_x$ for a given case. TFLES-FPI instead can just be employed to screen from an industrial point of view several configurations and obtain a first hierarchy between concepts.

## ACKNOWLEDGMENT

This work was granted access to the HPC resources of TGCC under the allocation A0142B13460 made by GENCI.
<figure><img src="image_14.png" /><figcaption>FIGURE 13: Case B: Comparison between MR0 and MR1 on the joint Probability Density Function of the mean $Y_{NO}$ source term density plotted against (Left) mean mixture fraction $\bar{Z}$ (Right) mean progress variable $\bar{Y}_C$.</figcaption></figure>

References

[1] D. Cecere, E. Giacomazzi, and A. Ingenito. A review on hydrogen industrial aerospace applications. *International Journal of Hydrogen Energy*, 39(20):10731–10747, 2014.

[2] R. Meloni, S. Orsino, N. Ansari, R. Yadav, D. Bessette, S. Castellani, P. Nassini, A. Andreini, and I. Boxx. Partially premixed hydrogen-methane flame simulations at relevant gas turbine conditions with a thickened flame model enhancement. *Proceedings of ASME Turbo Expo*, 2023. Paper No. GT2023-102427.

[3] P. Saini, I. Chterev, J. Pareja, M. Aigner, and I. Boxx. Effect of Pressure on Hydrogen Enriched Natural Gas Jet Flames in Crossflow. *Flow, Turbulence and Combustion*, 2020.

[4] S. Castellani, P. Nassini, A. Andreini, R. Meloni, E. Pucci, A. Valera-Medina, S. Morris, B. Goktepe, and S. Mashruk. Numerical modelling of swirl stabilised lean-premixed H2-CH4 flames with the artificially thickened flame model. *Proceedings of ASME Turbo Expo*, 2023. Paper No. GT2023-101994.

[5] D.G. Pugh, P.J. Bowen, R. Marsh, A.P. Crayford, J. Runyon, S. Morris, A. Valera-Medina, and A. Giles. Dissociative influence of H2O vapour/spray on lean blowoff and nox reduction for heavily carbonaceous syngas swirling flames. *Combustion and Flame*, 177:37–48, 2017. ISSN 0010-2180.

[6] A. Aniello, D. Laera, S. Marragou, H. Magnes, L. Selle, T. Schuller, and T. Poinsot. Experimental and numerical investigation of two flame stabilization regimes observed in a dual swirl h2-air coaxial injector. *Combustion and Flame*, 249:112595, 2023. ISSN 0010-2180.

[7] M. Leroy, C. Mirat, A. Renaud, S. Pugelli, S. Zurbach, and R. Vicquelin. Structure and NOX emissions of stratified

hydrogen-air flames stabilized on a coaxial injector. *Proceedings of ASME Turbo Expo*, 2023. Paper No. GT2023-102876.

[8] J. Leparoux, R. Mercier, S. Pugelli, M. Cailler, and V. Moureau. Numerical investigation of a hydrogen-air flame for nox prediction. *ASME*, 2023.

[9] M. Leroy, C. Mirat, A. Renaud, S. Pugelli, S. Zurbach, and R. Vicquelin. Structure and NOx emissions of stratified hydrogen-air flames stabilized on a coaxial injector. *Proceedings of the 11th European Combustion Meeting*, 2023.

[10] V. Moureau, P. Domingo, and L. Vervisch. Design of a massively parallel cfd code for complex geometries. *Comptes Rendus Mécaniques*, 2011.

[11] M. Kraushaar. Application of the compressible and low-Mach number approaches to Large Eddy Simulation of turbulent flow in aero-engines. PhD thesis, INPT, 2011.

[12] F. Ducros, F. Nicoud, and T. Poinsot. Wall-adapting local eddy-viscosity models for complex geometries. *Conference of Numerical Methods in Fluid Dynamics*, 1998.

[13] B. Van Leer. Towards the ultimate conservative difference scheme v. a second order sequel to godunov's method. *Journal of computational physics*, 1979.

[14] Lawrence Livermore National Laboratory cvode. https://computing.llnl.gov/projects/sundials/cvode. Accessed: 2023-07-19.

[15] O. Colin, F. Ducros, D. Veynante, and T. Poinsot. A thickened flame model for large eddy simulation of turbulent flows. *Physics of Fluids*, 2000.

[16] F. Charlette, C. Meneveau, and D. Veynante. A power-law flame wrinkling model for LES of premixed turbulent combustion part i: non-dynamic formulation and initial tests. *Combustion and Flame*, 131(1):159–180, 2002.
[17] H. Yamashita, M. Shimada, and T. Takeno. A numerical study on flame stability at the transition point of jet diffusion flames. *26th Symposium (International) on Combustion*, 1996.

[18] University of California San Diego san diego chemical-kinetic mechanisms for combustion applications. http://combustion.ucsd.edu. Accessed: 2023-07-19.

[19] J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird. *Molecular Theory of Gases and Liquids*. John Wiley and Sons, Inc., New York, NY, 2 edition, 1964.

[20] G. Kuenne, F. Seffrin, G. Fuest, T. Stahler, A. Ketelheun, D. Geyer, J. Janicka, and A. Dreizler. Experimental and numerical analysis of a lean premixed stratified burner using 1D Raman/Rayleigh scattering and Large Eddy Simulation. *Combustion and Flame*, 159(8):2669–2689, 2012.

[21] F. Pecquery, V. Moureau, G. Lartigue, L. Vervisch, and A. Roux. Modeling nitrogen oxide emissions in turbulent flames with air dilution: Application to les of a non-premixed jet flame. *Combsution and Flame*, 2014.

[22] R. Mercier, P. Auzillon, V. Moureau, Darabiha N., O. Gicquel, D. Veynante, and B. Fiorina. Les modeling of the impact of heat losses and differential diffusion on turbulent stratified flame propagation: Application to the tu darmstadt stratified flame. *Flow, Turbulence and Combsution 3*, 2014.

[23] P. Boivin. *Reduced-Kinetic Mechanisms for Hydrogen and Syngas Combustion Including Autoignition*. PhD thesis, Universidad Carlos III de Madrid, 2011.

[24] N. C. W. Treleaven, S. Puggelli, R. Mercier, J. Leparoux, X. Sun, and B. Sethi. High Altitude Relight Performance of Hydrogen-Air Micromix Combustion Systems. Volume 3B: Combustion, Fuels, and Emissions:V03BT04A063, 06 2023. doi: 10.1115/GT2023-103970.

[25] S. Puggelli, P. Leparoux, C. Brunet, R. Mercier, L. Liberatori, S. Zurbach, G. Cabot, and F. Grisch. Application of an automatic mesh convergence procedure for the Large Eddy Simulation of a multi-point injection system. *Proceedings ASME Turbo Expo, Paper No. GT2022-82272*, 2022.

[26] A. Porcarelli, B. Kruljevic, and I. Langella. Suppression of NOx emissions by intensive strain in lean premixed hydrogen flamelets. *International Journal of Hydrogen Energy*, 49:413–431, 2024.