# Multifidelity & multi-objective Bayesian optimization of hydrogen-air injectors for aircraft propulsion

Philippe Farjon<sup>a,*</sup>, Nicolas Bertier<sup>a</sup>, Sylvain Dubreuil<sup>b,c</sup>, Jérôme Morio<sup>b,c</sup>

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

Check for
updates

## ARTICLE INFO

Communicated by Damiano Casalino

**Keywords:**

- Bayesian optimization
- Multifidelity
- Hydrogen
- Combustion

## ABSTRACT

The use of hydrogen as a fuel is a promising way to reduce the emissions of civil aviation but it requires the development of wholly new injectors for the combustion chamber. Thanks to the increase in available computing power, the application of optimization techniques combined with CFD computations is now possible to develop these injectors. Among the optimization approaches, Bayesian optimization is particularly relevant when the objective functions and constraints of the optimization problem are expensive to evaluate which is the case in CFD-based optimization. Besides, the use of a multifidelity strategy allows to reduce the simulation cost of the Bayesian method. Therefore, this paper investigates the application of a multifidelity and multi-objective Bayesian approach to improve the performances of a laboratory swirl injector using hydrogen and operating in conditions close to industrial targets. This optimization study combines LES simulations as high-fidelity model with 2D RANS simulations as low-fidelity.

## 1. Introduction

The mitigation of pollutants and green-house gases constitutes an essential and imperative challenge for aircraft manufacturers. To reach this goal, H₂ combustion is a promising solution to replace kerosene as it does not emit CO₂ during combustion, but thermochemical properties of H₂ and kerosene are so different that it is mandatory to rethink in depth several key elements of the plane, including the combustion chamber's injector. It's worth noting that this injector, which is dedicated to both mix the reactants and stabilize the flame, concentrates most of the design issues for the combustion chamber.

Historically, combustion chamber's injector design improvement was based on a trial and error approach made of successive experimental tests and numerical simulations without the use of an automated optimization algorithm. This approach was successful for legacy kerosene injectors and has already been applied as well for hydrogen injectors with, for example, the Micromix concept developed by the Aachen University of Applied Sciences [1] or the lean direct injection concept (LDI) proposed by NASA [2]. Nevertheless, this approach is fundamentally limited because it is both costly and tedious to explore all the given

parameter space by hand. What is more, it is sometimes difficult to envisage some configurations that are far from what is usually known, but which could nonetheless turn out to be highly effective. However, with the advances in computing science, Computational Fluid Dynamics (CFD) simulations can be considered massively in the combustor design process. Moreover, new optimization techniques have the potential to improve the exploration of the design space and to lead to injectors with better performances, more stable combustion and lower emissions.

CFD-based optimization is already quite common in aerodynamics [3] with the use of gradient-based, gradient-free or surrogate-based methods (Bayesian optimization is a sub-family of these surrogate-based approaches). It is now an essential tool for aircraft, car or wind turbine aerodynamics designs but it is not yet the case for combustion devices. Indeed, the literature is scarce for reactive flows because combustion simulations are generally more challenging and expensive than aerodynamic ones, as more physical phenomena need to be taken into account. However, the surrogate-based approach seems really promising for optimization of combustion systems ([4–6] for example) as it is a good trade-off between the cost of the method and the theoretical ability to find a global optimum. Globally, this approach is more suited

<footnote><sup>a</sup> ONERA/DMPE, Université de Toulouse, F-31055 Toulouse, France</footnote>

<footnote><sup>b</sup> ONERA/DTIS, Université de Toulouse, F-31055 Toulouse, France</footnote>

<footnote><sup>c</sup> Fédération ENAC ISAE-SUPAERO ONERA, Université de Toulouse, F-31000 Toulouse, France</footnote>

<footnote>* Corresponding author.
E-mail address: philippe.farjon@onera.fr (P. Farjon).</footnote>
than gradient-based or genetic algorithms when the evaluation of the
objective function is based on numerical simulation. Among surrogate-
based optimization algorithms, Bayesian optimization (BO) [7] appears
as a good alternative possible for CFD-based optimization with complex
cases ([4,8,9]).

Regardless of the optimization method, one of the key aspects for op-
timization studies is the choice of the CFD simulation's level of fidelity.
The chosen one should ideally reproduce the main trends of quantities
of interest during the optimization process, while preserving a compu-
tational cost compatible with the current available computing power. In
the CFD context, turbulence modeling has a huge influence on both solu-
tion quality and CPU cost. On the one hand, Reynolds-Averaged Navier-
Stokes (RANS) approaches are widely used to compute industrial scale
configurations and often lead to fairly good results for aerodynamic
problems [10]. On the other hand, large scale scalar mixing between ar-
eas of the flow made of different compositions or temperatures is a weak
point of RANS methods. For combustion problems, the use of RANS ap-
proaches often needs a fine tuning of model's constants which makes
these methods poorly predictive. Despite this significant drawback, af-
fordability of RANS methods makes them attractive and almost all op-
timization procedures are based on these approaches ([11,12,8,13,5]).
Due to the RANS lack of accuracy, the most recent and state-of-the-
art combustion chamber's simulations are carried out with Large-Eddy
Simulations (LES), an approach allowing to capture the most energetic
scales of the flow and its unsteady features ([14–17] for example, among
numerous works). However, LES methods are far more expensive than
RANS ones (several orders of magnitude) which makes their use unsuit-
able for optimization procedures for aircraft combustors design and have
never been employed yet for this purpose to the authors' knowledge.
A multifidelity RANS-LES optimization would therefore take advantage
of the predictive capability of LES while preserving an acceptable com-
putational cost thanks to considerably cheaper RANS simulations.

Moreover, it is important to emphasize that the combustion cham-
ber's injector optimization must always integrate several objectives and
constraints to fulfill its specifications. In particular, it is necessary to
satisfy safety constraints like prohibiting flashback. Another important
objective consists in reducing the pollutant emissions such as nitric ox-
ides (NOx) to respect the more and more stringent environmental norms.
These two objectives can lead to different, even opposing, design choices
like diffusion flames to prevent flashback or lean premixed flames to
minimize NOx emissions. With the use of hydrogen, this trade-off prob-
lem is further amplified by the fact that hydrogen flames are both faster
(promoting flashback) and potentially hotter (promoting NOx emis-
sions) than kerosene ones. As a consequence, it is mandatory to take
into account several objectives and constraints which may be compet-
ing in the design process to find the best trade-off possible.

In the literature of aeronautical combustion chambers optimization,
LES have never considered yet to the authors' knowledge. Moreover, the
multifidelity and multi-objective approaches have never been applied
together for such a task.

Therefore, the final objective of this work is to present, apply and as-
sess a multifidelity and multi-objective Bayesian optimization method
for a realistic H₂-powered aircraft combustor injector. To achieve this,
the MFMO-SEGO algorithm developed by Charayron et al. [18] was used
with both high accuracy LES and affordable RANS approaches that were
combined together in a multifidelity strategy [19] for the Bayesian op-
timization, with the aim of achieving a new concept leading to stable
combustion (no flashback) and low NOx emissions

This paper begins with a review of the literature of optimization
of combustors in Section 2 before a brief description of the test injec-
tor, the chosen operating point and the numerical settings used by the
automated workflow for numerical simulations in Section 3. Then, an
overview of the principle of the Bayesian optimization methods used
for this study, is provided in Section 4. The results of the multifidelity
and multi-objective study using LES simulations are given in Section 5.

Finally, general conclusions and future prospects are presented in Sec-
tion 6.

## 2. A review of previous combustion chambers optimization studies

In the literature, several multi-objective studies have already been
performed for combustion systems but in terms of CFD simulation, only
RANS or laminar approaches were used. For example, Torkzadeh et al.
[11] or Asgari et al. [12] used 2D axisymmetric RANS computations of a
combustor to build a surrogate model based on Support Vector Machine
techniques and to perform a multi-objective optimization study directly
on the surrogate. In [11], four different objectives are considered: the
combustion efficiency, the pattern factor, the pollutant emissions (NO
and CO) and entropy generation. In [12], the study is non-reactive and
aims to evaluate the effects of fuel injection position in the combustor.
To do so, the objectives are: the evaporation efficiency, the homogene-
ity index (a metric measuring the homogeneity of the fuel vapor profile
at the chamber outlet), the average value of the fuel vapor mass fraction
at the outlet and entropy generation. The cost of 3D RANS simulations is
higher than 2D ones but still affordable for an optimization study. This
level of fidelity was used for example by Duchaine et al. [4] to optimize
the combustor efficiency and the exit temperature profile. They chose
a Kriging-based metamodel combined with a lower confidence bound
type criteria. It was also considered by Jeong et al. [20] to optimize a
diesel engine whose parametrization is made of 10 design variables. The
optimization algorithm called EGOMOP (Efficient Global Optimization
for Multi-Objective Problem) is a Bayesian method allowing to solve
in this case 4 objectives linked to pollutant emissions and thermal effi-
ciency. More generally, a surrogate based approach was used by Thomas
et al. [5] to optimize a small-scale combustor. Kriging is chosen to build
a metamodel which is directly used by a Multi-Objective Genetic Al-
gorithm (MOGA). The variance of Kriging is not directly integrated in
an infill criterion but it is used to check if the solution of the optimiza-
tion problem is trustworthy or if 3D RANS simulations must be added to
the design of experiment. For this study, 15 design variables were cho-
sen and the goal was to optimize the total pressure losses, the pattern
factor and the combustion efficiency. The approach of Laraia [6] for a
Lean Premixed Prevaporized (LPP) injection system is quite similar but
an artificial neural network is built instead of Kriging to minimize the
pollutant emissions (NOx and CO).

As previously mentioned, when the simulations become too expen-
sive, multifidelity modeling [19] may reduce the overall cost of the
optimization study. This approach was put into practice for combus-
tors for the first time by Wankhede et al. [8] who combined RANS and
URANS 2D simulations with several mesh sizes and time steps to opti-
mate the shape of a backward facing step using a Kriging-based surrogate
in a Bayesian framework. More recently, Zhang et al. [9] have proposed
a multifidelity study of a fuel injector feed arm with non-reactive RANS
simulations. They show that multifidelity performs better than monofi-
delity to optimize the pressure losses of their test case. Toal et al. [13]
developed a multifidelity approach for a complex geometry really close
to industrial combustors. In their work, they built a Co-Kriging surrogate
and employed directly the predictor of the surrogate for a surrogate-
based optimization study. Table 1 sums up the presented studies.

## 3. The MICADO test case

In this section, the injector that we want to optimize is described
along with the operating conditions. We explain the different hypothesis
made to simplify the injector and the numerical settings used in the
computational workflow are presented.

### 3.1. *The baseline injector and the operating point*

The configuration considered for this study is an axial swirl injector
on the MICADO test bench [22] (see Fig. 1). This experimental rig allows
<table><caption>**Table 1**
Some studies on combustion chamber optimization.</caption><thead><tr><th>Authors</th><th>Algorithm</th><th>Multi-objective</th><th>Multifidelity</th><th>Constraints</th><th>CFD approach</th></tr></thead><tbody><tr><td>Duchaine and al. [4]</td><td>BO</td><td>✓</td><td>✗</td><td>✗</td><td>RANS</td></tr><tr><td>Toal and al. [13]</td><td>SBO</td><td>✗</td><td>✓</td><td>✓</td><td>RANS</td></tr><tr><td>Zhang and al. [9]-1</td><td>BO</td><td>✓</td><td>✗</td><td>✓</td><td>RANS</td></tr><tr><td>Zhang and al. [9]-2</td><td>BO</td><td>✗</td><td>✓</td><td>✗</td><td>RANS</td></tr><tr><td>Janiga and al. [21]-1</td><td>Simplex</td><td>✗</td><td>✗</td><td>✗</td><td>Laminar</td></tr><tr><td>Janiga and al. [21]-2</td><td>GA</td><td>✗</td><td>✗</td><td>✗</td><td>Laminar</td></tr><tr><td>Jeong and al. [20]</td><td>BO</td><td>✓</td><td>✗</td><td>✗</td><td>RANS</td></tr><tr><td>Asgari and al. [12]</td><td>SBO</td><td>✓</td><td>✗</td><td>✗</td><td>RANS</td></tr><tr><td>Torkzadeh and al. [11]</td><td>SBO</td><td>✓</td><td>✗</td><td>✗</td><td>RANS</td></tr><tr><td>Thomas and al. [5]</td><td>SBO</td><td>✓</td><td>✗</td><td>✓</td><td>RANS</td></tr><tr><td>Wankhede and al. [8]</td><td>BO</td><td>✗</td><td>✓</td><td>✗</td><td>(U)RANS</td></tr><tr><td>Laraia and al. [6]</td><td>SBO</td><td>✓</td><td>✗</td><td>✓</td><td>RANS</td></tr><tr><td>This paper</td><td>BO</td><td>✓</td><td>✓</td><td>✓</td><td>LES &amp; RANS</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Baseline version of the H₂-air injector.</figcaption></figure>

the study of combustion with representative conditions of aeronautical combustion chambers (maximum 30 bars and 900 K for the incoming fresh gases). The MICADO test rig was originally designed to be fueled with methane but will be adapted to burn hydrogen as well. The present injector is fitted with 8 axial vanes through which the air is injected. Between each vane, there is a small hole (not represented in Fig. 1) allowing to inject the fuel. The vanes allow to impose a swirling motion to the flow whose intensity is described by the swirl number *S*. It is defined as the ratio of the axial flux of tangential momentum over the characteristic radius times the axial momentum flux:

$$S = \frac{\int_{R_i}^{R_o} U_x U_\theta r^2 dr}{\int_{R_i}^{R_o} R_o U_x^2 r dr} \quad (1)$$

Under the assumption that the vanes are very thin and that the axial and azimuthal velocities are uniform, *S* can be linked to the deviation angle *φ* imposed by the vanes and to the inner (*R*<sub>i</sub>) and outer (*R*<sub>o</sub>) radius of the injector [23]:

$$S = \frac{2}{3} \left( \frac{1 - \left( \frac{R_i}{R_o} \right)^3}{1 - \left( \frac{R_i}{R_o} \right)^2} \right) \tan(\phi). \quad (2)$$

Eq. (2) allows to choose the vanes angle that will be imposed to the CAD (Computer Aided Design) model to match the desired swirl number. The two radius *R*<sub>i</sub> and *R*<sub>o</sub> and the deviation angle *φ* (in red) can be seen in Fig. 2 (with diameters instead of radius).

The swirl number is a crucial parameter because the swirling motion generated by the vanes, creates an Inner Recirculation Zone (IRZ) at the exit of the injector which is the main flame stabilization mechanism of this injector. An example of a velocity field of a swirled flow is shown in Fig. 3. If the swirl number is too low, the IRZ will be not able to stabilize efficiently the flame whose length increases. On the other hand, with a high swirl number, the size of the IRZ will become more important

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Cut view of the injector. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Axial velocity field and streamlines for a configuration of the MICADO injector. The black isoline corresponds to a null axial velocity and allows to identify the IRZ and the ORZ.</figcaption></figure>

and the flame will be very compact with good stabilization properties. However, it should be noted that a too high swirl number can create a blocking effect in the chamber because the IRZ will take too much space. As a consequence, the swirl number is an important parameter which controls the aerodynamics of the flow and the stabilization of the flame. For the real injector installed in the test rig (called the baseline in the following sections), the vane angle is 50° corresponding to around *S* = 1.0.

The fuel can be injected also through the tip of the injector to create a pilot flame to help the stabilization of the main flame. The importance of this pilot flame can be quantified with the pilot ratio *PR* which is the ratio between the hydrogen mass flow rate injected by the pilot and the global hydrogen mass flow rate. It is a crucial parameter as it controls the combustion regime of the flame. Indeed, if *PR* = 0%, all hydrogen is injected through the main inlet before the vane and premixed with air. The advantage of the premixing regime is that when the equivalence ratio *φ* is low (lean conditions), the maximum temperature of the flame remains fairly low and the production of NOx in these conditions is limited. However, in premixing conditions, the risk of flashback, i.e. the flame going upstream in the injector, increases (especially with hydrogen whose laminar flame speed is very high compared to hydrocarbons). On the opposite, with *PR* = 100%, all hydrogen is fed through the pilot
<figure><img src="image_5.png" /><figcaption>Fig. 4. CAD model of the injector.</figcaption></figure>

<table><caption>**Table 2**
Operating conditions.</caption><thead><tr><th>P [bar]</th><th>T<sub>air</sub> [K]</th><th>T<sub>H₂</sub> [K]</th><th>m&#775;<sub>air</sub> [kg.s<sup>-1</sup>]</th><th>φ</th></tr></thead><tbody><tr><td>12</td><td>750</td><td>290</td><td>0.606</td><td>0.3</td></tr></tbody></table>

inlet. In this case, hydrogen and air burn in diffusion regime, i.e. near a stoichiometric equivalence ratio. In these conditions, the flame temperature is higher than in premixing regime. The diffusion regime is prone to produce more pollutant than the premixing one but the flame does not propagate and the risk of flashback is therefore limited. For intermediate values of PR, both regimes co-exist. In the baseline, PR equals 50%. Moreover, in this paper, the fuel entering the chamber through the pilot inlet is injected with an angle of 30° to create an hollow cone of hydrogen. As S and PR are two important parameters with a great impact on the flow, they are chosen as design variables in this paper.

The injector geometry was modeled with the CAD software Freecad,<fnref n="1" /> allowing a parametric version for optimization purposes. The CAD model of the MICADO injector is depicted in Fig. 4. Hydrogen will be considered as premixed with air from inlet injection, which is a common assumption for injector simulations with so small holes (diameter smaller than 1 mm).

Moreover, a part of the air flow rate is used to cool the chamber during real experiments. To reduce the size of the mesh, the cooling system is not modeled geometrically neither but its impact will be taken into account on the effective air mass flow rate entering in the injector and through the use of a special boundary condition for the most upstream wall of the chamber.

For the walls of the chamber, the temperature is imposed to 800 K. When hydrogen is injected through the pilot, hydrogen inside the fuel injection system (not represented in the CAD model) contributes to cool the wall, so the temperature of the wall around the pilot inlet is fixed to 400 K.

The operating conditions of an injector (temperature, pressure, fuel and air mass flow rate,...) change greatly according to the flight phase (cruise, take-off, idle,...). The simulations of this study are conducted at cruise conditions representative of a small-to-medium range (SMR) aircraft. The cruise flight phase is selected because it is the longest during a flight compared for example to idle or take-off. It is worth noting that this work is thus only focused on this operating point, summarized in Table 2, because adding several operating conditions increases significantly the computational cost and is beyond the scope of this work (even if this remains a long-term goal). The presented flow mass rate is representative of a sector of the chamber (20 sectors for the whole chamber). As explained before, the cooling system of the test rig is not simulated but we suppose that 20% of the air mass flow rate is used to cool the

<figure><img src="image_6.png" /><figcaption>Fig. 5. Meshes used for the optimization process.</figcaption></figure>

walls of the chamber so only 80% of the air mass flow rate given by Table 2 will be imposed at the inlet boundary condition.

Now that the full reference injector was presented, we can describe the possible geometrical simplifications

## 3.2. Reduction of the computational domain

The full injector geometry presented in Section 3.1 allows to get reference LES solutions. These reference solutions are obtained on a mesh of nearly 50 M cells and the CPU time of this LES simulation is not compatible with an optimization study. Indeed, the restitution time is about 2-3 days to have converged fields with 2400 cores (in this paper, for all simulations, the type of processor is Intel Xeon E5-2680 v4 2.4 GHz). The full geometry is thus difficult to integrate in an optimization approach so it is necessary to find some alternatives to reduce the computational burden. We propose here two solutions: a 3D sector and a 2D axisymmetric geometry (Fig. 5).

The first alternative to the full geometry is a 45° sector with lateral periodic conditions. For LES simulations, this approach needs meshes between 1 and 10 M cells which is a significant reduction compared to the full injector. The CPU time of this kind of computation is acceptable and the restitution time is around 20 hours on 768 cores. With this geometrical solution, it becomes possible to integrate LES simulations in the optimization process. Then, the high-fidelity model chosen for the multifidelity Bayesian optimization is LES simulation of the 3D sector, which we will refer to hereafter as S3DLES to lighten the reading.

The second alternative to the full geometry is a 2D axisymmetric geometry. As turbulence is an intrinsically 3D phenomenon, this approach cannot be used for LES but is interesting for 2D RANS simulations. For this test case, a typical mesh for 2D RANS is around 30K-50K cells which is again a huge reduction compared to the previous approaches using either the full injector or the 45° sector. Such computations last about 1.5 hours on 128 cores which is fast compared to the 3D LES approaches. As

<footnote n="1">¹ https://www.freecadweb.org/.</footnote>
<figure><img src="image_7.png" /><figcaption>Fig. 6. Automated CFD workflow.</figcaption></figure>

a consequence the 2D RANS simulation, which we will refer to hereafter as 2DRANS, is retained as the low-fidelity model.

For these two levels of geometric simplification, 2DRANS and S3DLES, the vanes of the swirler are not represented which helps the generation of the mesh particularly for the 3D sector. Nevertheless, the swirling motion imposed to the flow by the vanes is taken into account through the boundary condition of the inlet where an angular velocity is defined. The more the swirl number (or equivalently the vane angle) is important, the more the imposed angular velocity increases. With Eq. (1) and the hypothesis that $U_\chi$ is constant at the inlet and that $U_\theta$ can be written as $U_\theta = r\omega$ with $\omega$ the angular velocity, it is possible to link the swirl number $S$ to the angular velocity that will be imposed at the inlet of the numerical domain:

$$ \rho\omega = \frac{2R_0\dot{m}S}{\pi(R_0^4 - R_i^4)}, \quad (3) $$

with $\dot{m}$ the mass flow rate of fresh gases injected at the inlet. This feature allows to represent swirling flows in 2D and to have simpler meshes in 3D as the vanes can be quite complex to mesh accurately.

The informations about the meshes and the number of cores used for each computational domain are recalled in Table 3.

Now that the geometry simplifications introduced to take into account LES simulations in the multifidelity optimization framework have been described, the computational workflow as well as the numerical settings used for RANS and LES simulations will be presented in the next section.

### 3.3. The computational workflow and the numerical settings

All computations of this work are made with the ONERA in-house CFD code CEDRE [24]. This code is made of several solvers including the finite-volume solver CHARME to solve the compressible Navier-Stokes equations for reactive flows on unstructured meshes. When performing optimization studies, numerous computations are launched and it is necessary to automate as much as possible the computation process. Therefore, a workflow has been developed to run CFD computations in a Python framework. For 3D sector and 2D simulations, Gmsh [25] allows to modify the geometry and to build the meshes. The different surfaces of the injector are directly identified in Gmsh to define the boundaries conditions of the fluid domain.

Once the mesh is generated, it is transferred to CEDRE for the CFD simulation. When the simulation is over, the ONERA python package Cassiopée [26] is used for the post-processing to extract the data needed to evaluate the objective and constraint functions. The whole process is recalled in Fig. 6. All the meshes of this paper are unstructured and made of tetrahedras for the sector and triangles for 2D simulations. The meshes, whose characteristics can be found in Table 3, are refined at the tip of the injector and on the back of the chamber as it is the area where the flame is expected to be.

For both fidelity levels, the kinetics mechanisms are the San Diego one [27] for the H₂-air combustion and the Zeldovich one [28] for the thermal NO formation. As previously mentioned, the thermal NO formation is very sensitive to temperature: above 1800 K, its concentration increases following a quadratic evolution. No other pathways of NOx formations (N₂O or NNH) are included. For the LES simulation, the Smagorinsky model [29] is used as a subgrid-scale model and the interaction between chemistry and turbulence is taken into account thanks to a Thickened Flame model using the Volpiani efficiency function [30]. For the RANS simulation, the k-ω SST [31] model is chosen.

<table><caption>Table 3
Features and performances of the different CFD methodologies.</caption><thead><tr><th></th><th>Computational domain</th><th>Cells</th><th>Cores</th><th>Rest time</th></tr></thead><tbody><tr><td>2DRANS</td><td>2D axi</td><td>30K-50K</td><td>128</td><td>1.5 h</td></tr><tr><td>S3DLES</td><td>3D sector</td><td>1M-10M</td><td>768</td><td>20 h</td></tr><tr><td>F3DLES</td><td>full 3D</td><td>50M</td><td>2400</td><td>3 days</td></tr></tbody></table>

The numerical methods consist in a second order MUSCL (Monotonie Upwind Scheme for Conservation Laws) type scheme for space discretization, combined with a HLLC (Harten-Lax-van Leer-Contact) scheme [32] for hyperbolic fluxes. For time integration, an implicit scheme, associated to the GMRES (Generalized Minimal RESiduals) method to solve the linear system [33], is used.

### 3.4. Comparison of CFD approaches

A preliminary step before performing optimization studies consists in comparing the results obtained with the different CFD approaches. Indeed, as it was decided to reduce the computational domain and to use RANS simulations for the 2D injector in order to decrease the computational cost, it is necessary to check that these simplifications are not too important and allow to get results similar to a reference solution simulated with the full injector geometry. The full geometry, the sector and the 2D RANS approaches are compared for the baseline ($S = 1.0$ and $PR = 50\%$) in Figs. 7–9.

For the axial velocity, there is a good agreement between the three fields. This validates the use of a dedicated boundary condition for S3DLES and the 2DRANS simulations instead of modeling the vanes which is a very important point. For the temperature fields, 2DRANS tends to overestimate the temperature compared to the full geometry LES (F3DLES) and S3DLES. For the heat release rate, the flame brush of the full geometry is thicker than with 2DRANS and S3DLES. However, the phenomenology of the flow remains similar for all computations, with a premixed flame coming from the main inlet interacting with a diffusion flame induced by the pilot inlet. Overall, the discrepancies between each approaches remain acceptable for the baseline.

Comparing these approaches only on the baseline point is of course not enough. Indeed, to perform multifidelity optimisation, the variations of the different models for varying inputs must be well-correlated [19]. Consequently, additional comparisons were done for other values of $S$ and $PR$ (not shown here for conciseness) that confirmed that the simplification levels considered here were relevant.

## 4. Multifidelity and multi-objective Bayesian optimization

In this section, after explaining why a Bayesian optimization approach is chosen, we present Kriging and its multifidelity extension. Then, we explain the principle of the Efficient Global Optimization method (EGO) and how it can be extended to multifidelity, multi-objective and constrained problems.

### 4.1. Motivations behind the choice of Bayesian optimization

When optimizing a function defined on a continuous design space, three main techniques exist in the literature: the gradient-based method, the gradient-free one and the surrogate-based one.

For gradient-based method, the computation of the gradient is a difficult task with a CFD code. Indeed, it can be evaluated either by finite
<figure><img src="image_8.png" /><figcaption>Fig. 7. CFD approaches comparison - Units: Velocity: m.s⁻¹.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. CFD approaches comparison - Units: Temperature: K.</figcaption></figure>

differences which increases significantly the number of computations needed or by an adjoint solver. This last solution is quite common in aerodynamics [34] and allows to get the gradient at a reduced cost but the development of an adjoint solver is more difficult for reactive flows. Indeed, the source terms due to the chemical reactions are often not differentiable which decreases the accuracy of the adjoint solver. Therefore, there exists only few adjoint solvers for combustion in the literature ([35], [36]) and they are often based on simplification to derive the adjoint equations. As for gradient-free techniques such as genetic or evolutionary algorithms, they often need a prohibitive number of iterations to converge. There are few studies using gradient-based or gradient-free methods but they are limited to simple test case of combustion systems.

For instance, Janiga et al. [21] use directly a genetic algorithm or the Simplex algorithm but their simulations are 2D laminar, i.e. the turbulence is not taken into account whereas this phenomenon is crucial in a real combustor. As a consequence of the limitations of gradient-based and gradient-free methods, it was decided to choose the surrogate-based method and more particularly the Bayesian optimization one.

## 4.2. Kriging and its multifidelity extension

Kriging was firstly introduced in geostatistics by Krige [37] then mathematically formalized by Matheron [38]. It is now a popular type of surrogate model to approximate a function $y: X \subseteq \mathbb{R}^d \to \mathbb{R}$ by a Gaussian Process (GP) [39] with $d$ the dimension of the input vector $x$. The
<figure><img src="image_10.png" /><figcaption>Fig. 9. CFD approaches comparison - Units: Heat release: $W.m^{-3}$.</figcaption></figure>

construction of such a metamodel implies two steps: first, the definition of a prior GP {$Y(x), x \in \mathcal{X}$}, then the use of evaluations of y to condition this GP. The prior approximation $Y(x)$ can be divided in two parts:

$$Y(x) = \mu(x) + Z(x), \quad (4)$$

where $\mu(x)$ is a mean term of the form $\mu(x) = h'(x)\beta$ with $h'(x) = (h_1(x), ..., h_p(x))$ a set of basis function and $\beta = (\beta_1, ..., \beta_p)'$ a set of regression parameter associated to the basis function. $\mu(x)$ represents the trend of the Gaussian Process $Y(x)$ with respect to the input parameters. A common choice is to take a constant trend, i.e. $h_1(x) = 1$ and $h_p(x) = 0$ for $i > 1$. $Z(x)$ is a stationary Gaussian Process with a null mean. This GP is specified by its covariance function $k(x, x^*)$ also called the kernel of the GP, that allows to model the correlation between the different points of the design space. The choice of the kernel is essential as it controls the smoothness of the GP and thus the regularity of the y approximation. This covariance function can be written $k(x, x^*) = \sigma^2 r(x, x^*; \theta)$ with $\sigma^2$ the variance parameter and $\theta$ a vector allowing to parametrize the correlation function $r$. The different components of $\theta$ are called the hyper-parameters of the covariance function and they are estimated by Maximum Likelihood Evaluation (MLE).

With the design of experiment (DoE) $D = ((x_1, y_1), ..., (x_n, y_n))^T$ made of n evaluations of the true function and the non-conditioned GP $Y(x)$, it is possible to generate a conditioned GP $Y(x)|D$ (this notation means $Y$ conditionally to $D$). Once it is built, it allows to predict values for a new point $x^*$ that does not belong to the DoE instead of using a complex and costly numerical model.

One interesting property of Kriging compared to other kind of metamodel is that such a surrogate model provides an estimation of the uncertainty of the metamodel thanks to the Kriging variance. This feature is a huge advantage that can be used to develop optimization algorithms and to create an adaptive strategy.

Kriging was extended to multifidelity through the co-Kriging technique by [40]. It allows to build a metamodel with data from different levels of simulation fidelity. In this section, the principle of multifidelity kriging is presented for 2 levels of fidelity only but it can be extended to more levels. The chosen formulation of this paper is the Autoregressive Model of order 1 (AR1). It assumes an affine relationship between the high-fidelity and the low-fidelity functions. We use the recursive imple-

mentation of Le Gratiet [41] to construct multifidelity Kriging models. The AR1 approach supposes that:

$$f_{HF}(x) = \rho f_{LF}(x) + \delta(x), \\ \text{with } f_{LF} \perp \delta, \quad (5)$$

with $f_{HF}$ the high-fidelity model, $f_{LF}$ the low-fidelity model, $\rho$ a scaling factor as both fidelities may have different scales and $\delta$ a discrepancy function that captures the differences between the high and low-fidelity model. The symbol $\perp$ means that the discrepancy function is independent of the low-fidelity model.

To ease the construction of the metamodel, the DoEs are nested. Every point evaluated with the high fidelity model is also evaluated with the low-fidelity. This property allows to express the surrogate model variance in a closed form.

In practice, the construction of a multifidelity Kriging begins with the generation of a single-fidelity Kriging model based on the low-fidelity data. Then, a Kriging model of the discrepancy function is built with the differences between the high and low fidelity evaluations as data. With these two metamodels, we have now access to a multifidelity Kriging model. Similarly to the Kriging, the advantage of such a surrogate model is the possibility to predict quickly point outside the DoE. For a new point $x^*$, this predictor gives:

$$\hat{y}_{HF}(x^*) = h'(x^*)\beta + r(x^*)'\mathbf{R}^{-1}(Y_{DOE} - \mathbf{H}\beta). \quad (6)$$

And the variance is:

$$s_{HF}^2(x^*) = \rho^2 \sigma_{LF}^2 + \sigma_{HF}^2 - r(x^*)' \mathbf{R}^{-1} r(x^*), \quad (7)$$

where $\beta = (\beta_{LF}, \beta_{HF})$, $Y_{DOE} = (y_{LF}, y_{HF})$, $y_{LF}$ the true values of the low-fidelity function, $y_{HF}$ the true values of the high-fidelity function and with $\mathbf{R}$ the matrix of correlation of the samples and $\mathbf{r}(x)$ the vector of correlation between the samples and the points to predict. Moreover, we have $h'(x) = (\rho h'_{LF}(x), h'_{HF}(x))$ and $\mathbf{H}$ and $\mathbf{R}$ are defined by:
$$H = \begin{pmatrix} \mathbf{h}'_{LF}(\mathbf{x}_I^{LF}) & 0 \\ \vdots & \vdots \\ \mathbf{h}'_{LF}(\mathbf{x}_{n_{LF}}^{LF}) & 0 \\ \rho \mathbf{h}'_{LF}(\mathbf{x}_I^{HF}) & \mathbf{h}'_{HF}(\mathbf{x}_{n_{HF}}^{HF}) \\ \vdots & \vdots \\ \rho \mathbf{h}'_{LF}(\mathbf{x}_{n_{HF}}^{HF}) & \mathbf{h}'_{HF}(\mathbf{x}_{n_{HF}}^{HF}) \end{pmatrix},$$

$$R = \begin{pmatrix} \sigma_{LF}^2 R_{LF}(X_{LF}) & \rho \sigma_{LF}^2 R_{LF}(X_{LF}, X_{HF}) \\ \rho \sigma_{LF}^2 R_{LF}(X_{HF}, X_{LF}) & \rho \sigma_{LF}^2 R_{LF}(X_{HF}) \\ & + \sigma_{HF}^2 R_{HF}(X_{HF}) \end{pmatrix}.$$

With $X_{LF} = (\mathbf{x}_I^{LF}, ..., \mathbf{x}_{n_{LF}}^{LF})$ and $X_{HF} = (\mathbf{x}_I^{HF}, ..., \mathbf{x}_{n_{HF}}^{HF})$ respectively the sampling points of the low and high-fidelity models.

One interesting result of Le Gratiet, thanks to his recursive formulation, is that it is possible to express the variance contribution of each fidelity level:

$$\sigma_{cont_{LF}}^2(\mathbf{x}) = \rho \sigma_{LF}^2(\mathbf{x}) \text{ and } \sigma_{cont_{HF}}^2(\mathbf{x}) = \sigma_{\delta}^2(\mathbf{x}). \qquad (8)$$

These two expressions will be used for the choice of the level of fidelity in the multifidelity Bayesian optimization process.

### 4.3. Principles of Bayesian optimization

Bayesian optimization is a part of surrogate-based optimization methods and was first defined by Mockus [42] and Kushner [43]. Among this class of methods, the Efficient Global Optimization (EGO) algorithm [44] is a popular method to solve unconstrained optimization problems using black box function, i.e. a function whose expression and derivatives are unknown. The main idea of EGO is the following: the objective function is assumed to be the realization of a Kriging model and with this Gaussian approximation, new points are added iteratively to the DoE to increase the quality of the metamodel and to converge toward the global optimum. To determine new points $\mathbf{x}^*$, the original version of EGO uses the Expected Improvement (EI) as an infill criterion (called also an acquisition function):

$$EI(\mathbf{x}) = \mathbb{E}[\max(y_{min} - Y(\mathbf{x})|D, 0)], \qquad (9)$$

where $y_{min}$ is the current minimum value of the objective function over the DoE. As $Y(\mathbf{x})|D$ is a Gaussian random variable, it is possible to derive an analytic formula for EI:

$$EI(\mathbf{x}) = (y_{min} - \hat{y})\Phi(\frac{y_{min} - \hat{y}}{s}) + s\Phi(\frac{y_{min} - \hat{y}}{s}), \qquad (10)$$

with $\hat{y}$ the value of the kriging predictor in $\mathbf{x}$, $s$ the kriging variance, $\phi$ and $\Phi$ are respectively the probability density function and the cumulative distribution function of a standard normal random variable. This criterion is popular as it allows a good trade-off between exploitation of the Gaussian surrogate model and the exploration of the design space. To sum it up, this algorithm consists of three main steps at each iteration:

1. Construction of the Kriging metamodel with the DoE of size n.

2. Maximization of the EI to find the new point $\mathbf{x}^*$.

3. Addition of the new point $\mathbf{x}^*$ to DoE and evaluation of it ($y(\mathbf{x}^*)$).

EGO is appropriate for unconstrained problems but most engineering design problems are subject to constraints, so it is necessary to develop this algorithm to integrate them. There are several techniques to integrate the constraints such as the expected feasible improvement (EFI) [45,46], the stepwise uncertainty reduction (SUR) [47] or the Penalty-adjusted EI [48] for example. EGO can also be extended to constrained optimization studies with the Super EGO (SEGO) approach [49,50] which is the approach retained in this work to handle constraints. In the SEGO framework, the objective function $y$ is approximated by a GP exactly as in the EGO approach. To transform the problem in a constrained one, a metamodel is generated for each constraint function $g_i, i = 1, ..., m$

(with $m$ the number of constraints). The predictor of each constraint $\hat{g}_i$ is then used during the maximization of the EI. This results in the following constrained optimization problem for the infill criterion:

$$\max_{\mathbf{x} \in \Omega_f} EI(\mathbf{x}), \qquad (11)$$

where $\Omega_f$ is the feasible domain defined by:

$$\Omega_f := \{\mathbf{x} \in \mathcal{X} \subseteq \mathbb{R}^d : \hat{g}_1 \le 0, ..., \hat{g}_m \le 0\}. \qquad (12)$$

As with the EGO approach, the best point $\mathbf{x}^*$ given by the maximization of the EI is added to the DoE and the process continues iteratively until the convergence. Even if the evaluation of $\mathbf{x}^*$ by the true function gives an unfeasible result, it is added to the DoE as it allows to gain more information about the design space to update the Kriging models.

### 4.4. Multi-objective and multifidelity extension of Bayesian optimization

For most design problems, there are several objectives functions of interest and it is important to take them all in the definition of the optimization problem. Moreover, in a lot of domains as it is the case in CFD, several numerical models of variable fidelity can be used. There is often a high-fidelity model that is very accurate but quite expensive. On the other hand, there is a low-fidelity model, which is cheap but inaccurate. It is interesting to combine both models in an optimization process to exploit advantages of each model and to limit the drawbacks of each approach. Therefore, it is essential to upgrade SEGO to tackle multifi-delity and multi-objective problems. It is done with the Multi-Fidelity and Multi-Objective Super Efficient Global Optimization (MFMO-SEGO) algorithm presented in [18]. In general, a multi-objective constrained problem can be defined as:

$$\min_{\mathbf{x} \in \mathcal{X} \subseteq \mathbb{R}^d} \mathbf{f}(\mathbf{x}) := [f_1(\mathbf{x}), ..., f_{n_{obj}}(\mathbf{x})]^T, \qquad (13) \\ \text{s.t.} \quad \mathbf{g}(\mathbf{x}) < 0,$$

with $n_{obj}$ the number of objective functions considered in the optimization problem. With Bayesian optimization technique, it is necessary to build a surrogate model for each objective and to adapt the infill criterion to the multi-objective case. There are several multi-objective infill criteria in the literature and most of them are extensions of mono-objective infill criteria. For instance, the Expected Hypervolume Improvement (EHVI) [51] is an extension from the EI criterion. In this paper, we select the extension made in [52] to multi-objective problem of the scaled Watson and Barnes 2 (WB2S) criterion [50]. WB2S is an evolution of the Watson and Barnes 2 (WB2) criterion [53] which is:

$$WB2(\mathbf{x}) = -\hat{y}(\mathbf{x}) + EI(\mathbf{x}). \qquad (14)$$

The EI criterion is often difficult to optimize and the WB2 criterion, where the mean value $\hat{y}(\mathbf{x})$ of the surrogate is subtracted from the EI, allows to make the optimization easier as it is less multimodal. However, after some iterations, the EI can become very small compared to the predictor of the surrogate and this criterion loses its capacity of exploration of the design space. Therefore, the WB2S criterion was proposed in [50] to solve this issue. It consists in adding a scaling to the WB2 criterion:

$$WB2S(\mathbf{x}) = sEI(\mathbf{x}) - \hat{y}(\mathbf{x}), \qquad (15)$$

with $s$ a non-negative scale. WB2 and WB2S criteria were originally made for single-objective optimization but can be extended to multi-objective cases. In this paper, the multi-objective WB2S (MOO-WB2S) combined with the EHVI as sub-criterion can be written:

$$MOO - WB2S(\mathbf{x}) = \gamma EHVI(\mathbf{x}) - \chi(\hat{\mu}_f(\mathbf{x})), \qquad (16)$$

where $\gamma \in \mathbb{R}$ is a constant parameter which acts as a scaling factor and $\hat{\mu}_f : \mathbb{R}^d \rightarrow \mathbb{R}^{n_{obj}}$ is the GP prediction of all the objectives gathered in $\mathbf{f}$. $\chi : \mathbb{R}^{n_{obj}} \rightarrow \mathbb{R}$ is a scalarization operator. There are several choices for this operator as a maximization function or a sum function for example.
In this paper, a sum function is chosen ($\chi(\hat{\mu}_f(x)) = \sum_{i=1}^{n_{obj}} \hat{\mu}_{f_i}(x)$). In [52], this scalarization operator was compared with a maximization function and the results were quite similar for both operators, showing a limited influence of the scalarization operator on the results of the optimization study.

To perform multifidelity Bayesian optimization, a crucial aspect is that at each iteration, the most promising point and the appropriate level of fidelity at which it must be evaluated, are determined. The MFMO-SEGO algorithm follows the same approach as the Multi-Fidelity SEGO (MFSEGO) algorithm [54] where this aspect is divided in two steps. First, the infill criterion is optimized to find the most promising point. Then, to choose the fidelity level at which the point will be evaluated, the variance contribution of each fidelity is used. To be more explicit, with the expression of the variance contribution, it is possible to express the variance reduction of the high-fidelity model $\sigma_{red}^2(l, x^*)$ with $l \in [0, 1]$ (0 for the low-fidelity model and 1 for the high-fidelity one). This variance reduction can be written:

$$ \sigma_{red}^2(0, \mathbf{x}^*) = \sigma_{cont\_LF}^2(\mathbf{x}^*), \\ \text{and } \sigma_{red}^2(1, \mathbf{x}^*) = \sigma_{cont\_LF}^2(\mathbf{x}^*) + \sigma_{cont\_HF}^2(\mathbf{x}^*). \qquad (17) $$

When the point is evaluated at the high-fidelity level, it is also evaluated at the low-fidelity level to keep the nested DoE structure. In the MFSEGO approach, the variance reduction is penalized by the computational cost to find the most appropriate level of fidelity. With $c_0$ and $c_1$, the costs associated with the low-fidelity and high-fidelity model respectively, the enrichment level criterion can be written:

$$ l^* = \underset{l \in [0,1]}{\operatorname{argmax}} \frac{\sigma_{red}^2(l, \mathbf{x}^*)}{\left(\sum_{i=0}^{l} c_i\right)^2}. \qquad (18) $$

In the MFMO-SEGO framework, the enrichment level criterion defined in (18) is computed for each objective and then the maximum is selected:

$$ l^* = \max\{l_1^*, \dots, l_o^*, \dots, l_{n_{obj}}^*\}. \qquad (19) $$

Note that modification may be sub-optimal as the maximization step may lead to evaluate all objectives with the high-fidelity level whereas it would be useful only for one objective.

To sum it up, the MFMO-SEGO algorithm consists of four main steps at each iteration:

1. Construction of a Multifidelity Kriging metamodel for each objective and constraint.

2. Maximization of the MOO-WB2S criterion under constraints to find $x^*$ that solves Eq. (16).

3. Choice of the level of fidelity according to Eq. (19).

4. Add the new point $x^*$ to DoE and compute $f(x^*)$ and $g(x^*)$ at the retained fidelity.

It is worth noting for step 4, that when the high-fidelity level is chosen, the low-fidelity model is evaluated as well, in order to keep the nested structure of the DoE.

## 5. Results

The MFMO-SEGO methodology described in Section 4.4 is applied to the MICADO test case presented in Section 3.

### 5.1. Definition of the optimization problem

The problem presented here is a constrained optimization counting two design variables, with two fidelity levels and two objectives. As the main goal is to assess the MFMO-SEGO algorithm for the MICADO test case and as the LES simulations are quite costly, it is decided to begin with a rather simple optimization problem with a low-dimension. In a

future work, this number of parameters will be increased to come closer to industrial problems. As explained in Section 3, the swirl number S and the pilot ratio PR are two major parameters with a huge impact on the flow and on the pollutants formation. Therefore, they are chosen as design parameters for this study. The ranges of these parameters are given in Table 4. These parameters are linearly scaled to [0, 1] for the optimization process:

$$ x_{scaled} = \frac{x - x_{min}}{x_{max} - x_{min}}, \qquad (20) $$

with $x_{min}$ and $x_{max}$ respectively the lower and upper bound of the parameter.

The two objectives to minimize are the following: the NO mass flow rate $\dot{m}_{NO}$ at the outlet of the chamber and the flame length (FL). NO is the main pollutant species produced by the combustion of hydrogen with air, therefore it is essential to limit its production. It is also important to limit the size of the combustor to facilitate its integration in the whole engine and to reduce its mass. To reach this goal, it is necessary to have a compact flame and therefore to minimize the flame length FL. $\dot{m}_{NO}$ values are very small for cruise conditions (between $10^{-8}$ kg.s$^{-1}$ and $10^{-4}$ kg.s$^{-1}$) and preliminary tests have shown that a log transform of the NO mass flow rate makes the optimization easier. For similar reasons, the flame length FL is normalized by the length of the chamber CL. Two constraints are added for this problem: the absence of flashback $g_{FB}(x)$ and the absence of unburnt hydrogen at the outlet of the combustor $g_{UH}(x)$. For the flashback constraint, the heat release rate in the injector ($HR_{injector}$) is measured and then divided by the total heat release rate $HR_{total}$ (injector + chamber). If this ratio is superior to a fixed threshold (1% in this work), it means that there is a flashback. For the unburnt hydrogen, the $H_2$ mass flow rate at the outlet of the combustor is divided by the total $H_2$ mass flow rate injected. Similarly to the flashback constraint, if this ratio of $H_2$ mass flow rate reaches a value superior to 1%, it is considered that too much hydrogen is unburnt. These two constraints enforce realistic solutions. The choice of the objectives and constraints to define the optimization problem is an important task and can be quite difficult. Minimizing the pollutant production, preventing the flashback and limiting the unburnt hydrogen to ensure a high combustion efficiency are primary goals taken into account by a lot of studies in the literature and by the industrial designers. These three aspects are selected in this work too. As for the flame length, it is a more unusual and original metric. However, the minimization of this metric is of great interest to modify directly the shape of the flame which is not the case with other usual functions used in the literature. In the numerous studies reviewed in Section 2, two functions appear particularly important for the manufacturer: the pressure loss and the pattern factor (a metric measuring the homogeneity of the exit temperature profile). These two functions are not retained in this work. In the MICADO test case, the swirl vanes are one of the main sources of pressure loss. However, as the vanes are not represented in 2D or with the sector modelisations, it is not relevant to take into account the pressure loss in this optimization problem. Therefore, the absence of the vanes is an advantage to simplify the computational domain but it limits the choice of the objective function. Concerning the pattern factor, it is a function of interest to ensure the durability of the turbine located after the combustor in a real engine. Nevertheless, the MICADO test case is a laboratory injector without any turbine after the combustion chamber and is long enough to ensure an homogeneous temperature profile at the exit of the domain in most cases. Therefore, the pattern factor is not chosen either.

The two fidelity levels chosen for this problem are the 3D sector LES as high-fidelity model and the RANS 2D as low-fidelity model.

The optimization problem mentioned in the following as Case 1, can be written as:
<table><caption>**Table 4**
Ranges of variations of the parameters.</caption><thead><tr><td>Parameter</td><td>S [-]</td><td>PR [%]</td></tr></thead><tbody><tr><td>Lower bound</td><td>0.4</td><td>0</td></tr><tr><td>Upper bound</td><td>1.2</td><td>100</td></tr></tbody></table>

<table><caption>**Table 5**
Initial DoE (parameters are linearly scaled to [0, 1]).</caption><thead><tr><td>Calculation number</td><td>S [-]</td><td>PR [-]</td><td>Fidelity level</td></tr></thead><tbody><tr><td>1</td><td>0.12</td><td>0.82</td><td>LF</td></tr><tr><td>2</td><td>0.09</td><td>0.47</td><td>LF</td></tr><tr><td>3</td><td>0.51</td><td>0.94</td><td>LF</td></tr><tr><td>4</td><td>0.41</td><td>0.35</td><td>LF</td></tr><tr><td>5</td><td>0.94</td><td>0.27</td><td>LF</td></tr><tr><td>6</td><td>0.73</td><td>0.46</td><td>LF/HF</td></tr><tr><td>7</td><td>0.81</td><td>0.82</td><td>LF/HF</td></tr><tr><td>8</td><td>0.11</td><td>0.22</td><td>LF/HF</td></tr><tr><td>9</td><td>0.33</td><td>0.77</td><td>LF/HF</td></tr><tr><td>10</td><td>0.49</td><td>0.11</td><td>LF/HF</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Initial DoE (x0: swirl number S, x1: Pilot ratio PR).</figcaption></figure>

$$ 
\begin{aligned}
\mathbf{x}^* &= \underset{x \in [0;1]^2}{\operatorname{argmin}}[\log(\dot{m}_{NO}^{HF}(\mathbf{x})), \frac{FL^{HF}(\mathbf{x})}{CL}], \\
\text{s.t.} \quad &g_{FB}^{HF}(\mathbf{x}) \le 0 \text{ and } g_{UH}^{HF}(\mathbf{x}) \le 0, \\
\text{with} \quad &g_{FB}(\mathbf{x}) = \frac{HR R_{injector}}{HR R_{total}} - 0.01 \qquad (21) \\
\text{and} \quad &g_{UH}(\mathbf{x}) = \frac{\dot{m}_{H2}^{outlet}(\mathbf{x})}{\dot{m}_{H2}^{total}(\mathbf{x})} - 0.01.
\end{aligned}
$$

This optimization study will be compared to three other optimization studies: a multifidelity single objective constrained study for the flame length (Case 2), a multifidelity single-objective constrained study for the NO mass flow rate (Case 3) and a monofidelity multi-objective constrained study (Case 4). For this last case, the retained fidelity is the high-fidelity one. For all cases, the constraints are the same. We also performed an additional test where NSGA-II [55] was applied directly to multifidelity surrogates built from a nested LHS of size equivalent to the total cost of the 4 previous cases. This last test is not shown here for conciseness but this non-adaptive approach is less efficient than a Bayesian strategy.

## Database Pareto front at the end of the MFMO SEGO optimization

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Case 1 optimization: Multi-objective and multifidelity (obj1: NO mass flow rate, obj2: flame length).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Case 4 optimization: Multi-objective with only high-fidelity simulations (obj1: NO mass flow rate, obj2: flame length).</figcaption></figure>

### 5.2. Parameters of the optimization algorithm

For each case, we made one optimization run as the computations are quite costly for the LES 3D sector. For the multifidelity cases, before the beginning of the iterative process, we build an initial nested DoE of 5 HF points and 10 LF points. The size of the initial HF DoE follows a rule of thumb found by the authors in previous works that this size should equal $2 * n_{dim} + 1$ with $n_{dim}$ the dimension of the optimization problem. This rule of thumb allows that the most important part of the computation budget is used during the iterative process to adapt the metamodel. The initial DoE, shown in Fig. 10 and detailed in Table 5, is used for cases 1 to 3. This DoE is built with the Latin Hypercube Sampling (LHS) technique
<figure><img src="image_14.png" /><figcaption>Fig. 13. Case 1 evaluations: Multi-objective and multifidelity (objective 1: NO mass flow rate, objective 2: flame length).</figcaption></figure>

[56]. To construct this DoE, we use the SMT<fnref n="2" /> opensource toolbox [57, 58] where we choose the Enhanced Stochastic Evolutionary (ESE) [59] algorithm to create the LHS. The method implemented in SMT to build nested LHS [54] follows the one proposed by Le Gratiet in [41]. Firstly, a LHS $D_{HF}$ is built for the high-fidelity model. Secondly, a LHS $D_{LF}$ is built for the low-fidelity model. Then, for each point $x_i$ in $D_{HF}$, the closest point in $D_{LF}$ is replaced by $x_i$. This ensures the nested structure of the high and low fidelity design of experiments. For the monofidelity case, the initial DoE is made of the 5 HF points of the nested DoE. Once the initial DoE is generated and evaluated, the optimization process is performed with the Super Efficient Global Optimization with Mixture of Experts (SEGOMOE) package developed jointly at ONERA and ISAE-SUPAERO [50]. To build the different surrogate models needed for the optimization, the implementation in the SEGOMOE library is identical to the SMT one. For this paper, the surrogate models were built with a constant trend (i.e. $\mu(x) = \beta_1 * h_1(x)$ with $h_1(x) = 1$). For the kernel, we selected the square exponential or Gaussian kernel that can be written:

$$k(\mathbf{x}, \mathbf{x}') = \sigma^2 \prod_{l=1}^{n_{\text{dim}}} \exp(-\theta_l (\mathbf{x}_l - \mathbf{x}'_l)^2). \qquad (22)$$

As presented in Section 4.2, the infill criterion used for the multi-objective studies is the *MOO-WB2S* (Case 1 and 4) and the classical WB2S for single objective ones (Case 2 and 3). This criterion is optimized with the SLSQP technique [60] which is a quasi-Newton method. As the criterion can be multimodal, its optimization is based on a multistart approach to avoid to find only a local optimum. For the two dimension optimization problems of this work, the multistart approach consists in generating a LHS of 20 points. Then, these points are used as initial point of the SLSQP optimization.

For the multifidelity cases, the cost associated with the low-fidelity model is fixed to 1, and the cost for the high-fidelity model is 15. These costs were chosen by comparing the restitution time of the two CFD models.

<table><caption>Table 6
Repartition of the fidelities.</caption><thead><tr><th></th><th>HF</th><th>LF</th><th>Equivalent HF</th></tr></thead><tbody><tr><td>Case 1</td><td>13</td><td>29</td><td>14.93</td></tr><tr><td>Case 2</td><td>10</td><td>72</td><td>14.8</td></tr><tr><td>Case 3</td><td>10</td><td>71</td><td>14.73</td></tr><tr><td>Case 4</td><td>15</td><td>0</td><td>15</td></tr></tbody></table>

<table><caption>Table 7
Optimum values.</caption><thead><tr><th></th><th>Optimal parameters [%;]</th><th>Optimum value</th></tr></thead><tbody><tr><td>Baseline</td><td>[0.97; 0.5]</td><td>[-4.20; 0.31]</td></tr><tr><td>Case 1</td><td>[1.14; 0.0]</td><td>[-7.22; 0.123]</td></tr><tr><td>Case 2</td><td>[1.108; 9.4]</td><td>[-4.65; 0.189]</td></tr><tr><td>Case 3</td><td>[0.982; 1e<sup>-11</sup>]</td><td>[-7.16; 0.131]</td></tr><tr><td>Case 4</td><td>[0.844; 1e<sup>-13</sup>]</td><td>[-7.15; 0.143]</td></tr></tbody></table>

For all cases, we have a budget of 15 equivalent HF computations, initial DoE included (if we perform 15 LF computations, it equals 1 HF computations in terms of budget).

## 5.3. Results of the optimization studies

Table 6 summarizes the number of HF and LF simulations performed for each case and Table 7 shows the different optimum values for each case study.

For Cases 1 and 4 that are multi-objective, as the number of iterations is quite low due to the cost of the simulations, there is only one point in the Pareto Front (Fig. 11 and Fig. 12). For Fig. 13–16, the blue zone between iteration 0 and 9 represents the initial LF DoE and the orange one represents the initial HF DoE. The red line corresponds to the current best HF evaluation. The Case 1, i.e. the multifidelity and multi-objective study, gives the best result at the end of the computation budget with this initial DoE. Cases 3 and 4 give an optimum with slightly lower performances. The optimum found by the Case 2 gives a more important value for FL (Fig. 14) and the NO mass flow rate is quite high too. As it is not an objective of this case, it explains why the NO mass flow rate is so high compared to other cases but it proves the interest of

<footnote n="2">² https://github.com/SMTorg/smt.</footnote>
<figure><img src="image_15.png" /><figcaption>Fig. 14. Case 2 evaluations: Mult fidelity single-objective for flame length (objective 2: flame length).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 15. Case 3 evaluations: Mult fidelity single-objective for NO mass flow rate (objective 1: NO mass flow rate).</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. 16. Case 4 evaluations: Multi-objective with only high-fidelity simulations (objective 1: NO mass flow rate, objective 2: flame length).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. CFD fields of the optimum (HF) - Units: Temperature: K, Velocity: $m.s^{-1}$, Heat release: $W.m^{-3}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 18. Iso-surface of instantaneous heat release rate for a value of $10^{10} W.m^{-3}$ (blue) with the instantaneous temperature field. Results of the 45° injector duplicated on a 360° geometry.</figcaption></figure>

performing a multi-objective optimization to get an efficient injector. The solution of case 3 (Fig. 15) is quite close to the solution of the case 1 so it could be possible to perform only a single objective multifidelity study and to find a low-NOx injector with a compact flame but there is no guarantee. Performing a multi-objective study remains a safer choice to get an injector with improved performances.

The results of the different cases suggest that the best parameters are a swirl number superior to 0.8 (which is rather a high value) and a pilot hydrogen repartition of 0%, i.e. no hydrogen injected through the pilot. No pilot flame means an injector working in premixing regime.

With Case 1 and Case 4, it is possible to evaluate the interest of using a multifidelity approach for a given budget. The optimum of Case 4 is found after adding 6 HF computations to the initial DoE, i.e. 11 HF computations (Fig. 16). For Case 1, 14 LF and 3 HF are needed with the initial DoE to find the optimum, i.e. 9.6 equivalent HF simulations (Fig. 13). Therefore, Case 1 allows to find a better optimum than Case 4 in less iterations which shows the interest of multifidelity for this case and this initial DoE. This result will need a more in-depth validation with the use of other initial DoE for example, but the multifidelity seems really interesting as the use of LF simulations allows to explore the design space at a reduced cost and the first HF simulations added to the DoE gives already a significant improvement of the injector performances for both objectives whereas it is not the case for the monofidelity case. The difference in terms of equivalent HF computations between case 1 and
<figure><img src="image_20.png" /><figcaption>(a) Case 1 final DoE: Multi-objective multifidelity</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Case 2 final DoE: Multifidelity FL</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(c) Case 3 final DoE: Multifidelity NO</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(d) Case 4 final DoE: Multi-objective monofidelity</figcaption></figure>

<caption>**Fig. 19.** Final DoE for the 4 cases. x0: swirl number S, x1: pilot repartition PR. HF points are in orange, LF points in blue and the points of the initial DoE are surrounded by black circles.</caption>

4 to reach the optimum is 1.4 and as the restitution time of a LES sector is around 20 hours, there is a gain of around one day in terms of restitution time between Case 1 and Case 4.

The different DoEs at the end of the optimization are presented in Fig. 19. For Case 1, most points are added in the “corners” of the design space. In general, the variance of a surrogate is high near the edges of the design space. That’s why the algorithm add new points in the “corners” for exploration. The best point for Case 1 is found in a zone with a high swirl and a zero PR. Other HF points are added in this zone for exploitation. The analysis of Fig. 13 and Fig. 19a shows that after adding some LF points for exploration, the first HF point at iteration 26 allows to reduce significantly the NO mass flow rate and to decrease the flame length too. After this first HF point, other HF simulations are added by the algorithm for exploitation. After this step, the last HF points allows more exploration of the domain. A similar behavior is visible for Case 2 and 3 where most points are added toward the edges of the domain. For Case 2, the last HF point of the initial DoE gives a value of 0.22 which is already quite low. The first iterations (until the second added HF evaluation) select points close to the last HF point of the initial DoE. Therefore, the beginning of the optimization is an exploitation phase. This phase does not allow to reduce the flame length. Then, as explained for Case 1, there is an exploration phase with points in the lower right corner and in

the upper left corner. After this exploration, the algorithm detects that the region with a high swirl ($S \ge 0.8$) and a low PR ($PR \le 0.2$) is the zone where the flame is the more compact and two HF evaluations are added including the current optimum. Contrary to Case 1, no HF point is added for a null value of PR and a swirl superior to 0.9. In Case 1, this zone allows to reach a smaller value of flame length than for Case 2. In Case 2, the variance reduction of adding a HF simulation seems not sufficient compared to the cost of the simulation and only LF points are computed in this promising zone. However, with a higher budget, it is likely that HF points would be evaluated in this region. For Case 3, the adaptive phase begins with an exploration phase of the lower left corner, i.e. the zone with a very low swirl (lower than 0.1) and PR close to 0.0. This exploration allows to find very promising solutions with a significant reduction of the NO mass flow rate compared to the initial DoE. After detecting this promising zone, there is an exploitation phase followed by an exploration phase of higher swirl value ($S \ge 0.2$ with PR still closed to zero). The exploration continues similarly to Case 1 and 2 with the edges of the domain. This exploration is mainly done with LF simulations which allows a wide exploration at a reduced cost. At the end of the budget, a last HF point is computed for a normalized swirl number of 0.73 which allows to reduce again the NO mass flow rate. This Case and the previous ones show that when PR = 0.0%, the NO
<figure><img src="image_24.png" /><figcaption>Fig. 20. Correlation for both metrics of interest after the end of Cases 1-3.</figcaption></figure>

mass flow rate is low and its variations with respect to the swirl number are weak. Therefore, the performance of the best injector of Case 3 is very satisfying even if the NO mass flow rate is higher than for Case 1. Indeed, computing new HF points with a higher swirl and PR = 0.0% would offer only a limited gain of NO mass flow rate for a higher computational cost. For Case 4, in the initial DoE, the best point for both objectives is found for S = 0.49 and PR = 0.11. The first iterations after the initial DoE add new simulations close to this point, it is therefore a phase of exploitation. In Fig. 19d, it is possible to identify a cluster of points due to this exploitation of the surrogates. Most points are gathered in a zone delimited by $0.5 \le S \le 0.6$ and $0.0 \le PR \le 0.11$. The PR parameter is progressively reduced until the edge of the domain is reached. However, with this low budget, there is a lack of exploration which explains why Case 4 finds an optimum with lower performances compared to Case 1.

## 5.4. Comparison of the fidelities post-optimization

For Cases 1 to 3, both S3DLES and 2DRANS simulations are used. Each time a HF simulation is performed, a LF simulation at the same design point is done to keep the nested property of the DoE. The values for both fidelities are compared for the NO mass flow rate in Fig. 20a and for the flame length in Fig. 20b. It is worth noting that most points are close to the red line $y = x$ showing a pretty good correlation between the fidelities with less than 2% of relative error for most design points. However, some of them display huge discrepancies with relative errors superior to 20%. For the flame length, there is one point reaching 99% of relative error because the RANS simulation predicts a flashback whereas it is not the case with the LES simulation. These results combined with the two $r^2$ criterion [61] visible on Fig. 20 underline that for both metrics of interest, the level of correlation is satisfying but not perfect. It shows that it is not possible to trust completely the LF fidelity and that is relevant to combine HF and LF simulations in a multifidelity strategy.

## 5.5. CFD results of the optimum

As mentioned in Section 5.3, the optimum found by Case 1 has the best performances among all the computations made for this paper. The Swirl number is 1.14 which is a strong value. With the use of Eq. (2), it gives an angle of 54.3 deg. The H₂ pilot ratio is zero which means that no hydrogen is injected through the pilot. All the hydrogen is perfectly premixed with air and injected through the main inlet.

The temperature, heat release rate and axial velocity mean fields of the 3D sector in the central plane are shown in Fig. 17. In Fig. 18, the instantaneous heat release rate is shown along with the temperature field. The heat release rate allows to visualize the position of the flame. For this set of parameters, one gets a premixed flame attached to the lips of the injector. As the combustion regime is premixed and that there is no pilot flame, it burns at a low equivalence ratio $\varphi$ of around 0.4 which allows a rather “low” temperature. Indeed, the maximum temperature in this computation is around 1800 K whereas when there is a pilot flame, a part of hydrogen burns in a diffusion regime at stoichiometry and the maximum temperature can rise to 2500 K. As the production of NO increases with the temperature, it explains why this optimum produces few pollutants. The velocity field shows that thanks to the swirling effect, there is an important IRZ that helps to the stabilization of the flame.

# 6. Conclusion and perspectives

In this work, we applied the MFMO-SEGO algorithm to the CFD optimization of a laboratory swirl injector using hydrogen: the MICADO test case. To perform the optimization with this multi-objective, multifidelity and constrained approach, 2D RANS simulations were combined with 3D LES sector computations. The use of LES remains quite rare in the literature due to the computational cost but it allows to integrate simulations with more accuracy as more physical phenomena can be simulated. The MICADO test case was parametrized with two parameters and the optimization study with MFMO-SEGO was compared to studies with monofidelity or with single-objective. The MFMO-SEGO case outperforms the other ones. Indeed, it converges faster than the monofidelity approach and the multi-objective feature is necessary to ensure efficient solutions. In future work, it is planned to perform other optimization studies to understand more in-depth the influence of the initial DoE and of the numerical parameters of the algorithm such as the kernel. Moreover, we plan to continue to improve the fidelity of the CFD simulation for combustors optimization by integrating full injector 3D LES computations. We also aim to increase the dimension of the optimization problem to come closer to the industrial constraints and also to combine cruise simulations with take-off simulations to get a more efficient injector on a larger range of operating conditions. Indeed, with the cruise flight phase which is the longest one, the take-off phase is where the conditions are the most extreme. Therefore, combining these two phases in the optimization problem will allow to have a better design for the injector.
CRediT authorship contribution statement

**Philippe Farjon:** Writing – original draft, Visualization, Software, Methodology, Investigation, Formal analysis. **Nicolas Bertier:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition. **Sylvain Dubreuil:** Writing – review & editing, Validation, Supervision, Methodology. **Jérôme Morio:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition.

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.

Data availability

Data will be made available on request.

Acknowledgements

This work was funded by la Région Occitanie, "RHYO/Défi clé H2" and the French Government (DGAC, French Civil Aviation Authority) as part of the PHYDROGENE project ("France 2030" plan).

Appendix A. Supplementary material

Supplementary material related to this article can be found online at
https://doi.org/10.1016/j.ast.2024.109383.

References

[1] H.H. Funke, N. Beckmann, S. Abanteriba, An overview on dry low NOx micromix combustor development for hydrogen-rich gas turbine applications, Int. J. Hydrog. Energy 44 (13) (2019) 6978–6990, https://doi.org/10.1016/j.ijhydene.2019.01.161.

[2] C. Marek, T. Smith, K. Kundu, Low emission hydrogen combustors for gas turbines using lean direct injection, in: 41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, 2005, pp. 1–27.

[3] S.N. Skinner, H. Zare-Behtash, State-of-the-art in aerodynamic shape optimisation methods, Appl. Soft Comput. 62 (2018) 933–962, https://doi.org/10.1016/j.ASOC.2017.09.030.

[4] F. Duchaine, T. Morel, L. Gicquel, Computational-fluid-dynamics-based Kriging optimization tool for aeronautical combustion chambers, AIAA J. 47 (3) (2009) 631–645, https://doi.org/10.2514/1.37808.

[5] N.R. Thomas, M.P. Rumpfkeil, A.M. Briones, T.J. Erdmann, B.A. Rankin, Multiple-objective optimization of a small-scale, cavity-stabilized combustor, in: AIAA Scitech 2019 Forum (January), 2019, pp. 1–20.

[6] M. Laraia, M. Manna, S. Colantuoni, P. Di Martino, A multi-objective design optimization strategy as applied to pre-mixed pre-vaporized injection systems for low emission combustors, Combust. Theory Model. 14 (2) (2010) 203–233, https://doi.org/10.1080/13647831003746187.

[7] P.I. Frazier, A tutorial on Bayesian optimization, arXiv:1807.02811, 2018.

[8] M.J. Wankhede, N.W. Bressloff, A.J. Keane, Combustor design optimization using co-Kriging of steady and unsteady turbulent combustion, J. Eng. Gas Turbines Power 133 (12) (09 2011) 121504, https://doi.org/10.1115/1.4004155.

[9] X. Zhang, D.J.J. Toal, A.J. Keane, F. Witham, J. Gregory, M. Ravikanti, E. Aurifeille, S. Stow, M. Rogers, M. Zedda, Isothermal combustor prediffuser and fuel injector feed arm design optimization using the Prometheus design system, J. Eng. Gas Turbines Power 138 (6) (2015) 061504, https://doi.org/10.1115/1.4031711.

[10] P. Catalano, M. Amato, An evaluation of rans turbulence modelling for aerodynamic applications, Aerosp. Sci. Technol. 7 (7) (2003) 493–509, https://doi.org/10.1016/S1270-9638(03)00061-0.

[11] M. Torkzadeh, F. Bolourchifard, E. Amani, An investigation of air-swirl design criteria for gas turbine combustors through a multi-objective CFD optimization, Fuel 186 (2016) 734–749, https://doi.org/10.1016/j.fuel.2016.09.022.

[12] B. Asgari, E. Amani, A multi-objective CFD optimization of liquid fuel spray injection in dry-low-emission gas-turbine combustors, Appl. Energy 203 (2017) 696–710, https://doi.org/10.1016/j.apenergy.2017.06.080.

[13] D.J. Toal, X. Zhang, A.J. Keane, C.Y. Lee, M. Zedda, The potential of a multifidelity approach to gas turbine combustor design optimization, J. Eng. Gas Turbines Power 143 (5) (2021) 1–9, https://doi.org/10.1115/1.4048654.

[14] P. Wolf, G. Staffelbach, L.Y. Gicquel, J.-D. Müller, T. Poinsot, Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers, Combust. Flame 159 (11) (2012) 3398–3413, https://doi.org/10.1016/j.combustflame.2012.06.016.

[15] V. Moureau, P. Domingo, L. Vervisch, From large-eddy simulation to direct numerical simulation of a lean premixed swirl flame: filtered laminar flame-pdf modeling, Combust. Flame 158 (7) (2011) 1340–1357, https://doi.org/10.1016/j.combustflame.2010.12.004.

[16] L. Selle, G. Lartigue, T. Poinsot, R. Koch, K.-U. Schildmacher, W. Krebs, B. Prade, P. Kaufmann, D. Veynante, Compressible large eddy simulation of turbulent combustion in complex geometry on unstructured meshes, Combust. Flame 137 (4) (2004) 489–505, https://doi.org/10.1016/j.combustflame.2004.03.008.

[17] P. Benard, G. Lartigue, V. Moureau, R. Mercier, Large-eddy simulation of the lean-premixed preccinsta burner with wall heat loss, Proc. Combust. Inst. 37 (4) (2019) 5233–5243, https://doi.org/10.1016/j.proci.2018.07.026.

[18] R. Charayron, T. Lefebvre, N. Bartoli, J. Morlier, Towards a multi-fidelity & multi-objective bayesian optimization efficient algorithm, Aerosp. Sci. Technol. 142 (PB) (2023) 108673, https://doi.org/10.1016/j.ast.2023.108673.

[19] M. Giselle Fernández-Godino, C. Park, N.H. Kim, R.T. Haftka, Issues in deciding whether to use multifidelity surrogates, AIAA J. 57 (5) (2019) 2039–2054, https://doi.org/10.2514/1.J057750.

[20] S. Jeong, Y. Minemura, S. Obayashi, Optimization of combustion chamber for diesel engine using Kriging model, J. Fluid Sci. Technol. 1 (2) (2006) 138–146, https://doi.org/10.1299/jfst.1.138.

[21] G. Janiga, D. Thévenin, Reducing the CO emissions in a laminar burner using different numerical optimization methods, Proc. Inst. Mech. Eng. A, J. Power Energy 221 (5) (2007) 647–655, https://doi.org/10.1243/09576509JPE387.

[22] A. Cochet, V. Bodoc, C. Brossard, O. Dessornes, C. Guin, R. Lecourt, M. Orain, A. Vincent-Randonnier, ONERA test facilities for combustion in aero gas turbine engines, and associated optical diagnostics, Aerospace Lab. 11 (2016) 1–16, https://doi.org/10.12762/2016.AL11-01.

[23] Y. Huang, V. Yang, Dynamics and stability of lean-premixed swirl-stabilized combustion, Prog. Energy Combust. Sci. 35 (4) (2009) 293–364, https://doi.org/10.1016/j.pecs.2009.01.002.

[24] A. Refloch, B. Courbet, A. Murrone, C. Laurent, J. Troyes, G. Chaineray, J.B. Dargaud, F. Vuillot, CEDRE software, J. Aerospace Lab 2 (2011) 1–10.

[25] C. Geuzaine, J.F. Remacle, Gmsh: a 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Methods Eng. 79 (11) (2009) 1309–1331, https://doi.org/10.1002/nme.2579.

[26] C. Benoit, S. Péron, S. Landier, Cassiopee: a CFD pre- and post-processing tool, Aerosp. Sci. Technol. 45 (2015) 272–283, https://doi.org/10.1016/j.ast.2015.05.023.

[27] P. Saxena, F.A. Williams, Testing a small detailed chemical-kinetic mechanism for the combustion of hydrogen and carbon monoxide, Combust. Flame 145 (1–2) (2006) 316–323, https://doi.org/10.1016/j.combustflame.2005.10.004.

[28] Y.B. Zeldovich, The oxidation of nitrogen in combustion and explosions, Acta Physicochem. USSR 21 (1946) 577–628.

[29] J. Smagorinsky, General circulation experiments with the primitive equations, Mon. Weather Rev. 91 (3) (1963) 99, https://doi.org/10.1175/1520-0493(1963)091<0099>GCEWTP>2.3.CO;2.

[30] P.S. Volpiani, T. Schmitt, O. Vermorel, P. Quillatre, D. Veynante, Large eddy simulation of explosion deflagrating flames using a dynamic wrinkling formulation, Combust. Flame 186 (2017) 17–31, https://doi.org/10.1016/j.combustflame.2017.07.022.

[31] F.R. Menter, R.B. Langtry, S.R. Likki, Y.B. Suzen, P.G. Huang, S. Völker, A correlation-based transition model using local variables—Part I: Model formulation, J. Turbomach. (2004) 413–422, https://doi.org/10.1115/1.2184352.

[32] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, Berlin Heidelberg, 2009.

[33] Y. Saad, M.H. Schultz, Gmres: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986) 856–869, https://doi.org/10.1137/0907058.

[34] J. Peter, R. Dwight, Numerical sensitivity analysis for aerodynamic optimization: a survey of approaches, Comput. Fluids 39 (2010) 373–391, https://doi.org/10.1016/j.comfluid.2009.09.013.

[35] M. Lemke, J. Reiss, J. Sesterhenn, Adjoint based optimisation of reactive compressible flows, Combust. Flame 161 (10) (2014) 2552–2564, https://doi.org/10.1016/j.combustflame.2014.03.020.

[36] A. Kord, J. Capecelatro, A discrete-adjoint framework for optimizing high-fidelity simulations of turbulent reacting flows, Proc. Combust. Inst. 39 (4) (2023) 5375–5384, https://doi.org/10.1016/j.proci.2022.06.021.

[37] D.G. Krige, A statistical approach to some basic mine valuation problems in the Witwatersrand, J. Chem. Metall. Min. Soc. S. Afr. 52 (6) (1951) 119–139, https://doi.org/10.2307/3006914.

[38] G. Matheron, Principles of geostatistics, Soc. Econ. Geol. 58 (8) (1962) 1246–1266, https://doi.org/10.2113/gsecongeo.58.8.1246.

[39] C.E. Rasmussen, C.K.I. Williams, Gaussian Processes for Machine Learning, in: Adaptive Computation and Machine Learning, MIT Press, 2006.

[40] M.C. Kennedy, A. O'Hagan, Predicting the output from a complex computer code when fast approximations are available, Biometrika 87 (1) (2000) 1–13, https://doi.org/10.1093/biomet/87.1.1.
[41] L. Le Gratiet, Multi-fidelity Gaussian process regression for computer experiments, Theses, Université Paris-Diderot - Paris VII, Oct. 2013.

[42] J. Močkus, On Bayesian methods for seeking the extremum, in: G.I. Marchuk (Ed.), Optimization Techniques IFIP Technical Conference, Novosibirsk, July 1-7, 1974, Springer Berlin Heidelberg, Berlin, Heidelberg, 1975, pp. 400-404.

[43] H.J. Kushner, A versatile stochastic model of a function of unknown and time varying form, J. Math. Anal. Appl. 5 (1) (1962) 150-167, https://doi.org/10.1016/0022-247X(62)90011-2.

[44] D.R. Jones, M. Schonlau, W.J. Welch, Efficient global optimization of expensive black-box functions, J. Glob. Optim. 13 (4) (1998) 455-492, https://doi.org/10.1023/A:1008306431147.

[45] M. Schonlau, W.J. Welch, D.R. Jones, Global versus local search in constrained optimization of computer models, Lect. Notes Monogr. Ser. 34 (1998) 11-25, http://www.jstor.org/stable/4356058.

[46] B. Talgorn, S. Le Digabel, M. Kokkolaras, Statistical surrogate formulations for simulation-based design optimization, J. Mech. Des. 137 (2) (2015) 021405, https://doi.org/10.1115/1.4028756, https://asmedigitalcollection.asme.org/mechanicaldesign/article-pdf/137/2/021405/6225912/md_137_02_021405.pdf.

[47] V. Picheny, A stepwise uncertainty reduction approach to constrained global optimization, in: S. Kaski, J. Corander (Eds.), Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, in: Proceedings of Machine Learning Research, vol. 33, PMLR, Reykjavik, Iceland, 2014, pp. 787-795, https://proceedings.mlr.press/v33/picheny14.html.

[48] M.J. Sasena, P. Papalambros, P. Goovaerts, Exploration of metamodeling sampling criteria for constrained global optimization, Eng. Optim. 34 (3) (2002) 263-278, https://doi.org/10.1080/03052150211751.

[49] M.J. Sasena, Flexibility and Efficiency Enhancements for Constrained Global Design Optimization with Kriging Approximations, Ph.D. thesis, University of Michigan, 2002.

[50] N. Bartoli, T. Lefebvre, S. Dubreuil, R. Olivanti, R. Priem, N. Bons, J.R. Martins, J. Morlier, Adaptive modeling strategy for constrained global optimization with application to aerodynamic wing design, Aerosp. Sci. Technol. 90 (2019) 85-102, https://doi.org/10.1016/j.ast.2019.03.041.

[51] M.T. Emmerich, K.C. Giannakoglou, B. Naujoks, Single- and multiobjective evolutionary optimization assisted by gaussian random field metamodels, IEEE Trans. Evol. Comput. 10 (4) (2006) 421-439, https://doi.org/10.1109/TEVC.2005.859463.

[52] R. Grapin, Y. Diouane, J. Morlier, N. Bartoli, T. Lefebvre, P. Saves, J. Bussemaker, Regularized infill criteria for multi-objective Bayesian optimization with application to aircraft design, in: AIAA AVIATION 2022, Chicago, United States, 2022.

[53] A.G. Watson, R.J. Barnes, Infill sampling criteria to locate extremes, Math. Geol. 27 (5) (1995) 589-608, https://doi.org/10.1007/BF02093902.

[54] M. Meliani, N. Bartoli, T. Lefebvre, M.A. Bouhlel, J. Martins, J. Morlier, Multi-fidelity efficient global optimization: methodology and application to airfoil shape design, in: AIAA Aviation 2019 Forum, 2019.

[55] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Trans. Evol. Comput. 6 (2) (2002) 182-197, https://doi.org/10.1109/4235.996017.

[56] F.A. Viana, A tutorial on Latin hypercube design of experiments, Qual. Reliab. Eng. Int. 32 (5) (2016) 1975-1985, https://doi.org/10.1002/qre.1924.

[57] M.A. Bouhlel, J.T. Hwang, N. Bartoli, R. Lafage, J. Morlier, J.R.R.A. Martins, A python surrogate modeling framework with derivatives, Adv. Eng. Softw. (2019) 102662, https://doi.org/10.1016/j.advengsoft.2019.03.005.

[58] P. Saves, R. Lafage, N. Bartoli, Y. Diouane, J. Bussemaker, T. Lefebvre, J.T. Hwang, J. Morlier, J.R. Martins, Smt 2.0: a surrogate modeling toolbox with a focus on hierarchical and mixed variables gaussian processes, Adv. Eng. Softw. 188 (2024) 103571, https://doi.org/10.1016/j.advengsoft.2023.103571, https://www.sciencedirect.com/science/article/pii/S096599782300162X.

[59] R. Jin, W. Chen, A. Sudjianto, An efficient algorithm for constructing optimal design of computer experiments, J. Stat. Plan. Inference 134 (1) (2005) 268-287, https://doi.org/10.1016/j.jspi.2004.02.014.

[60] D. Kraft, Algorithm 733: Tomp-fortran modules for optimal control calculations, ACM Trans. Math. Softw. 20 (3) (1994) 262-281, https://doi.org/10.1145/192115.192124.

[61] D.J. Toal, Some considerations regarding the use of multi-fidelity Kriging in the construction of surrogate models, Struct. Multidiscip. Optim. 51 (6) (2015) 1223-1245, https://doi.org/10.1007/s00158-014-1209-5.