Total of combustion energy Traditional biomass Coal Gas Oil Biofuels

Keywords: Quantification d Quantifying uncertainties in Large Eddy Simulations of auto-ignition-stabilized H 2 Uncertainty Quantification, Turbulent combustion, Hydrogen, Digital twin

Quantification d'incertitudes dans les simulations aux grandes échelles de flammes d'hydrogène stabilisées par auto-allumage

Dans cette thèse, nous propageons les incertitudes provenant du mécanisme de cinétique chimique, de mesures expérimentales et du modèle de turbulence à travers la simulation aux grandes échelles (LES) d'une flamme turbulente. Nous choisissons de travailler sur la flamme de Cabra H 2 car son mécanisme de stabilisation principal est l'auto-allumage, ce qui la rend extrêmement sensible aux paramètres qui affectent les taux de réaction. Nous simulons la flamme de Cabra en utilisant une approche basée sur la tabulation de flammelettes instationnaires et sur des PDF présumées. La propagation d'incertitudes repose sur la réduction a priori de la dimension incertaine en utilisant une configuration simplifiée : un réacteur homogène isobare qui s'auto-allume. Une étude Monte-Carlo sur ce réacteur dévoile un comportement très non-linéaire dans la plage d'incertitudes considérée.

Ce comportement a été identifié à la transition entre les auto-allumages "faible" et "fort". Cette étude a aussi permis de réduire les incertitudes venant du mécanisme chimique et de la température du co-flow à un espace de dimension 2.

Dans cet espace réduit, une approche basée sur une surface de réponse est réalisée en utilisant seulement un faible nombre d'échantillons LES. Nous étudions également l'influence du coefficient de dissipation de la variance de fraction de mélange. Dans un premier temps, seule l'incertitude de la hauteur de stabilisation de la flamme est quantifiée.

Dans un second temps, nous quantifions les incertitudes de toutes les quantités mesurées expérimentalement à toutes les positions dans le domaine de calcul. Cela est rendu possible par une analyse en composantes principales (PCA) et trois surfaces de réponses pour prédire les coefficients des plus importants modes de la PCA. source of energy in our society. To minimize its pollutant emissions in a context where operating conditions and computational models are uncertain, we need to quantify the uncertainties of our predictive simulation results.

In this thesis, we propagate kinetic, experimental, and model uncertainties in the Large Eddy Simulation of a turbulent flame. We chose to work with the H 2 Cabra flame because its main stabilization mechanism is auto-ignition, which makes it extremely sensitive to parameters that impact reaction rates. We simulate the Cabra flame using an unsteady flamelet and presumed PDF approach. The uncertainty propagation relies on an a priori uncertain dimension reduction using a simplified configuration: an auto-igniting homogeneous reactor at constant pressure. An a priori Monte-Carlo study allows highlighting an extremely non-linear phenomenon in the uncertainty range: the transition between weak and strong auto-ignition. It also allows reducing the uncertainty coming from the co-flow temperature and kinetic mechanism to a mere two dimensions.

In this reduced space, a surrogate-modelbased approach is viable using a limited number of LES observations. We later investigate the sensitivity of the simulation to the dissipation coefficient of mixture fraction variance. We first quantify the uncertainty of the Cabra flame lift-off height. Then, we turn our gaze towards the quantification of uncertainties of every experimentally measured quantity at every location in the computational domain. This is done using a Principal Component Analysis on the training samples and three surrogate models to predict the coefficients of the most contributing modes brought forth by the PCA.

iii

Table of Contents

Introduction

The primary motivation: Uncertainty Quantification Uncertainty Quantification (UQ) is a framework that allows taking into account unknown (epistemic uncertainties) and/or randomly fluctuating (aleatoric uncertainties) variables in a given model, in order to determine the uncertainty of a quantity of interest. These uncertainties can stem from many sources: measurement errors, inadequate computational models, unknown parameter values... From this wide concept, diverse procedures can be derived. Among them, forward and inverse uncertainty propagations are becoming more and more commonplace:

• Forward uncertainty propagation is a procedure that uses the uncertain inputs of a model to evaluate the probability distribution of an output quantity of interest in the model. This can be useful, for example, to evaluate risks, compute safety margins, or estimate performance variability.

• Conversely, given a model and a set of uncertain outputs, inverse uncertainty propagation infers the probability distributions of the input variables of the model. This allows calibrating computational model parameters given uncertain measurement results.

One of the major interests of the uncertainty quantification framework, in particular its non-intrusive branch, is that "the model" can be practically anything that takes a set of variables as inputs and produces another set of output variables. "The model" can be a computer simulation, an analytical model, an experiment, etc. Therefore, uncertainty quantification can be used for extremely various problems. Hence, the methods presented and developed in this thesis are not limited to our particular use case. Three, very diverse, examples of fields that use uncertainty quantification are given below:

• Geosciences can use uncertainty quantification to predict the probability distributions of natural phenomena. One such example is the study by Li et al. [START_REF] Li | Quantifying initial and wind forcing uncertainties in the Gulf of Mexico[END_REF] which evaluated the impact of initial conditions and wind-forcing uncertainties on the currents of the Gulf of Mexico. Another example is the seminal work of Krige [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF] who introduced the Kriging method as a prospection tool for the Witwatersrand gold mines in South Africa.

• Modern Portfolio Theory, first introduced by Markowitz in 1952 [START_REF] Markowitz | Portfolio Selection[END_REF], is the foundation of the modern financial sector. By applying the probability theory to the financial markets, he developed a new notion of risk and proposed to evaluate investment portfolios by evaluating both the expectation and variance of return. His work paved the way for previously-unseen profits by banks and investment funds. However, over-reliance on his risk-management method and its successors, coupled with over-simplifications of the market mechanisms notoriously led to periodic financial crises in the 20th and 21st centuries (1987: Black Monday, 1990: LTCM collapse, 1994: Bond market crisis, 2007: Subprime mortgage crisis...) [START_REF] Kay | Radical Uncertainty, Decision-Making Beyond the Numbers[END_REF].

• Finally, engineering sciences are fond of uncertainty quantification. In particular, the combustion community, which initiated this thesis, is in demand for uncertainty quantification for purposes of robustifying efficient combustors, better predicting pollutant emissions, etc.

A more practical motivation: Combustion

Modern human civilization is addicted to combustion. From the myth of Prometheus to the industrial revolution and the wonders of the Saturn V moon rocket, humans have used fire as a means to exploit the chemical energy contained in their environment. A major tipping point was the 17th and 18th centuries when combustion moved from being a source of thermal energy (for cooking, metallurgy, or simply as a source of heat) to being a source of mechanical energy thanks to a flourishing of innovations in the steam engine. By the 19th century, combustion had become ubiquitous in human activities: Draught animals were replaced by engines; wood and peat were gradually replaced by coal, gas, and oil. Mechanical energy, made cheaper and way more plentiful by this new technology, was consumed at an ever-increasing rate [START_REF] Stanley | The Coal Question[END_REF]. Among other examples, transportation habits have dramatically evolved since then. Where an intercontinental journey might have taken a few weeks aboard a sailing ship, the same journey is now done in a few hours aboard any commercial airplane. As a result, intercontinental flights have become commonplace for citizens of developed countries. This rebound effect probably reached its climax with the Concorde, a supersonic airliner which was routinely used for back-and-forth trips between Europe and America for business meetings.

For several decades now, the combustion of fossil fuels has been known to degrade our environment and its ecosystems by upsetting the global climate. Yet, despite a clear scientific consensus and collective awareness, mankind continues to prospect and exploit new fossil fuel deposits. At this point, it seems rather clear that our economic organization traps us in a prisoner's dilemma: the first institution, either company or country, trying to dispense with combustion would probably do without potential profits and spending power respectively. This is highlighted by several studies linking the per-capita Gross Domestic Product (GDP) with the per-capita greenhouse gas emissions [START_REF] Tukker | Environmental Impacts of Products: A Detailed Review of Studies[END_REF][START_REF] Hertwich | Carbon Footprint of Nations: A Global, Trade-Linked Analysis[END_REF]. Therefore, it seems that mankind will likely continue to rely on combustion for as long as it finds fuel to burn. This tendency is suggested by Figure 1, which plots the evolution and shares, by source, of the world's total primary energy supply between 1800 and 2021 [88]. Since 1997 and the signature of the Kyoto Protocol, fossil fuel consumption has not shown any sign of decrease. On the contrary, mankind, as a whole, has burnt more and more fossil fuel since then, with only small drops in 2008 and 2020 caused by the subprime crisis, the Covid-19 pandemic, and the respective global recessions that followed.

The primary means of action to reduce greenhouse emissions is therefore to dramatically reduce the global energy supply and the share of combustion in this supply. Hence, it is mostly a drastic cultural and political change that is needed. Greenhouse emissions Figure 1: Global primary energy consumption, in TWh, by source, between 1800 and 2021. Source: OWID [88] and global warming are not problems that engineering sciences alone can solve.

For the activities that will be judged absolutely necessary for society, and for which combustion is definitely irreplaceable, the environmental cost must be minimized. In this prospect, some potential alternatives to fossil fuels are presently gaining traction. Among them, biofuels and dihydrogen are particularly in demand these days and draw in human resources and corporate and state funding.

Despite significant drawbacks at the production stage of these alternative fuels (soil depletion [START_REF] Stagnari | Sustainable Soil Management[END_REF] and overall debatable carbon balance [START_REF] Xu | Life-cycle energy use and greenhouse gas emissions of palm fatty acid distillate derived renewable diesel[END_REF], and electric power requirements respectively), they are seen by decision-makers as valid substitutes that may reduce pollution without thwarting too much the energetic bulimia of our productivist society. This belief may lead to yet another rebound effect and should be considered with great caution.

Another area of work is the increase of fuel efficiency of (present and irreplaceable future) combustion devices and the reduction of their toxic byproducts like CO, NO x , soot, etc. The motivation of this thesis stems from this context: Cautious uncertainty quantification is required to ensure that industrial systems stay in their safe and regulatory operating ranges, the alternative being the implementation of conservative safety margins which typically decrease performance.

As engineering has come to rely more and more on numerical simulation, Computational Fluid Dynamics (CFD) is now as commonplace in design offices as in research labs. Still, these simulations are subject to uncertainties, coming from the operating conditions (temperature, mixture composition in the case of biofuels, etc.), the combustion or turbulent models, or even from the chemical kinetic model itself.

Kempf et al. [START_REF] Kempf | Error analysis of large-eddy simulation of the turbulent non-premixed sydney bluff-body flame[END_REF], for example, investigated the influence of the Smagorinsky constant, a key parameter of the turbulence model in Large Eddy Simulation (LES). Figure 2 presents the qualitative influence of its value on the dynamics of the Sydney flame, a nonpremixed bluff-body configuration. This example exhibits that the choice of turbulence model constants is paramount when trying to correctly predict a turbulent flame in LES. Fiorina et al. [START_REF] Fiorina | Challenging modeling strategies for LES of nonadiabatic turbulent stratified combustion[END_REF] published another interesting study, where several research groups were asked to simulate the same flame. It corresponds to the unconfined Turbulent Stratified Flame from Darmstadt. Each group used its own numerical setup: each defined its own computational domain, mesh type and mesh size, kinetic mechanism and turbulent combustion model. Figure 3 presents the radial profiles of the Favre-averaged gas temperature at an axial location of 15 mm from the burner found by the different groups (plain and dashed lines denote different wall boundary conditions). Although all results are in global agreement with each other and with the experiment, some discrepancies remain between the prediction of the different groups, due to their modelization choices.

Figure 3: Radial profiles of the gas temperature in the Turbulent Stratified Flame, as simulated by each research group. Adapted from Fiorina et al. [START_REF] Fiorina | Challenging modeling strategies for LES of nonadiabatic turbulent stratified combustion[END_REF] The two previous examples show that even high-fidelity simulations like Large Eddy Simulation are still plagued by epistemic uncertainties resulting from the modelization choice. In other cases, experimental measurement uncertainties also make measurements and simulations seemingly disagree. This is the case of the H 2 Cabra flame [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]: a turbulent, non-premixed, auto-ignition-stabilized flame. Figure 4 shows a collection of simulations and experimental measurements of the Cabra flame lift-off height. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF] are the black symbols, simulations [START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF] are the colored lines. Note the measurement uncertainty of the co-flow temperature (black error bar) Again, in this case, discrepancies appear between the simulation results of those authors. Perhaps more importantly, the experimental measurements show both a strong sensitivity of the lift-off height to the co-flow temperature and a significant uncertainty in the measurement of this temperature. In this configuration, proper uncertainty quantification is much-needed to evaluate the agreement between experiments and simulations and highlight the underlying phenomena.

Uncertainty quantification in inexpensive combustion simulations

In the literature, most of the intersection between combustion and uncertainty quantification revolves around the uncertainties in kinetic mechanisms. A kinetic mechanism is a set of reactions and equations to determine their rate constants as a function of the ambient conditions. In the vast majority of cases, rate constants are found using an Arrhenius law:

k j = A j T β j exp -E a j RT
Where k j is the rate constant of reaction j, A j is its pre-exponential factor, T is the gas temperature, β j is a dimensionless constant, E a j is the energy activation of reaction j, and R is the universal gas constant. A j , β j , and E a j are either inferred from experiments [START_REF] Hong | A new shock tube study of the H+O2→OH+O reaction rate using tunable diode laser absorption of H2O near 2.5microm[END_REF] or with collision theory computations [START_REF] Karkach | Ab initio analysis of the transition states on the lowest triplet H2O2 potential surface[END_REF] which are currently affordable only for extremely simple molecules. In the former case, measurement uncertainties and dependencies between all the coefficients lead to uncertainties of the coefficients in the Arrhenius laws of the kinetic mechanism.

Theoretically, a joint probability distribution should be considered for all the coefficients of all the reactions in the mechanism. In practice, those coefficients are often assumed independent, and the pre-exponential constants are the only ones considered uncertain. Their uncertainty is typically represented by a log-normal distribution parametrized by an "Uncertainty Factor" (UF), such that:

ξ j = log(A j /A 0 j ) 1 3 log(UF j ) ∼ N (0, 1)
Forward uncertainty propagation methods are readily available to propagate the variability of some simulation inputs to a given quantity of interest. They have been applied to simple combustion problems such as homogeneous reactors and one-dimensional laminar flames [START_REF] Phenix | Incorporation of parametric uncertainty into complex kinetic mechanisms: Application to hydrogen oxidation in supercritical water[END_REF][START_REF] Reagan | Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection[END_REF][START_REF] Reagan | Spectral stochastic uncertainty quantification in chemical systems[END_REF][START_REF] Gy | Uncertainty analysis of updated hydrogen and carbon monoxide oxidation mechanisms[END_REF][START_REF] Klippenstein | Uncertainty driven theoretical kinetics studies for CH3OH ignition: HO2+CH3OH and O2+CH3OH[END_REF]. For example, Phenix et al. [START_REF] Phenix | Incorporation of parametric uncertainty into complex kinetic mechanisms: Application to hydrogen oxidation in supercritical water[END_REF] propagated the kinetic uncertainties from an 8-reaction reduced hydrogen combustion mechanism through the simulation of a plug-flow reactor in a supercritical state. Figure 5 reports the results of a Monte-Carlo study of this configuration. Significant uncertainties are observed. In particular, Phenix et al. show that they are much larger than the difference between simulations performed with and without real-gas corrections in the equation of state.

In most practical cases, non-intrusive methods are used, as they do not require any modification of the CFD solver. It must be noted, though, that intrusive methods have also been investigated [START_REF] Reagan | Spectral stochastic uncertainty quantification in chemical systems[END_REF]. In the non-intrusive studies listed above, random sampling was used to propagate uncertainties. This is not a realistic approach for more expensive simulations because the simplest laboratory-scale configurations can cost several hundred to several thousand hours of CPU time to simulate and several thousand runs may be needed to converge meaningful statistics. Surrogate-based approaches are more suited to this kind of problem because a response surface is built using relatively few observations, corresponding to equally few deterministic simulation runs. This response surface is then resampled to collect statistics. This approach has been implemented by Sheen and co-workers [START_REF] Sheen | Spectral uncertainty quantification, propagation and optimization of a detailed kinetic model for ethylene combustion[END_REF][START_REF] Sheen | The method of uncertainty quantification and minimization using polynomial chaos expansions[END_REF][START_REF] Wang | Combustion kinetic model uncertainty quantification, propagation and minimization[END_REF] or more recently Yousefian et al. [START_REF] Yousefian | Uncertainty Quantification of NOx Emission Due to Operating Conditions and Chemical Kinetic Parameters in a Premixed Burner[END_REF] for simple combustion problems. However, surrogate models still require many training samples when the [START_REF] Phenix | Incorporation of parametric uncertainty into complex kinetic mechanisms: Application to hydrogen oxidation in supercritical water[END_REF] uncertain dimension is high, which can be problematic for costlier numerical setups.

Uncertainty quantification in turbulent combustion simulations

More recent works have propagated uncertainties through turbulent combustion simulations. Notably Mueller and co-workers [START_REF] Michael | Chemical kinetic uncertainty quantification for Large Eddy Simulation of turbulent nonpremixed combustion[END_REF] used a two-step uncertainty propagation to build an uncertain tabulation for the chemistry. This allowed them to deal with only a reduced set of dependent random variables to propagate through the Large Eddy Simulation (LES) of the Sandia flame. A few years later, Mueller and Raman [START_REF] Mueller | Model form uncertainty quantification in turbulent combustion simulations: Peer models[END_REF] propagated model form and parameter of the mixture fraction variance dissipation model in that same simulation. Figure 6 shows the distribution of time-averaged temperature and CO mass fraction, estimated using a collocation method, along several radial profiles. Again, notable uncertainties are found by these authors on the quantities of interest, especially on the most downstream profile.

Concomitantly, Khalil et al. [START_REF] Khalil | Uncertainty quantification in LES of a turbulent bluff-body stabilized flame[END_REF] propagated the uncertainties of the turbulent model (Smagorinsky coefficient and turbulent Prandtl and Schmidt numbers) through the LES of a bluff-body-stabilized flame. They did not need to reduce the uncertain dimension as only three random variables were considered, and uncertain chemistry was disregarded. A similar study was performed by Lucor et al. [START_REF] Lucor | Sensitivity analysis of large-eddy simulations to subgrid-scale-model parametric uncertainty using polynomial chaos[END_REF] in the context of non-reactive LES. All these uncertain turbulent flame studies cited above were performed with either [START_REF] Mueller | Model form uncertainty quantification in turbulent combustion simulations: Peer models[END_REF] a limited number of uncertain variables [START_REF] Lucor | Sensitivity analysis of large-eddy simulations to subgrid-scale-model parametric uncertainty using polynomial chaos[END_REF][START_REF] Khalil | Uncertainty quantification in LES of a turbulent bluff-body stabilized flame[END_REF][START_REF] Mueller | Model form uncertainty quantification in turbulent combustion simulations: Peer models[END_REF] or an a priori dimension reduction [START_REF] Michael | Chemical kinetic uncertainty quantification for Large Eddy Simulation of turbulent nonpremixed combustion[END_REF]. In a radically different approach, other studies made an a posteriori dimension reduction of the uncertain space using the concept of "active subspace" [START_REF] Constantine | Active Subspace Methods in Theory and Practice: Applications to Kriging Surfaces[END_REF]. In this framework, the physical model is first sampled randomly on the entire uncertain space, and the observations are used to infer a smaller uncertain space on which the surrogate is built. Constantine and co-workers [START_REF] Constantine | Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet[END_REF] used it to quantify the uncertainty of the integrated exit pressure of a scramjet by propagating experimental and boundary conditions uncertainties. Ji et al. [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] used an active subspace to estimate the distribution of the H 2 Cabra flame lift-off height considering only the kinetic mechanism to be uncertain. This approach is only affordable for numerical setups that can be run several times the number of uncertain dimensions, i.e. it is limited to highly symmetric setups or RANS simulations and/or to cases where the uncertain dimension is low.

Objective of the thesis

The literature for uncertainty quantification in combustion contains few examples of CFD, and even fewer examples of Large Eddy simulations. LES is more precise than the Reynolds-Averaged Navier-Stokes approach, and it is also much more expensive due to its spatial and temporal resolution. However, LES is not immune to uncertainties, in particular when the physical configuration features high sensitivities to certain parameters.

This thesis aims to expand on the studies cited above and synthesize high-fidelity simulations and high-dimensional uncertainties. In other words, we aim to propagate high-dimensional uncertainties through expensive computer simulations, that cost several tens of thousand to several million of CPU hours for each run. In particular, we aim to propagate kinetic, turbulence model, and measurement uncertainties in the Large Eddy Simulation of a turbulent flame. We also choose to work on a flame stabilized by auto-ignition because their lift-off height is first-order-determined by the chemical kinetics of the mixture. Hence, our preferred choice of configuration was the Cabra flame [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF], which fulfills our criteria and has already been extensively studied, both numerically and experimentally. Furthermore, we choose to work on the H 2 variant of the Cabra flame because the contribution of auto-ignition is clearer in this flame than in the CH 4 variant [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Domingo | Auto-ignition and flame propagation effects in LES of burnt gases diluted turbulent combustion[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | A numerical study of auto-ignition in turbulent lifted flames issuing into a vitiated co-flow[END_REF][START_REF] Chun | A 3D DNS study of the stabilization of a turbulent lifted hydrogen/air jet flame in an autoignitive heated coflow[END_REF] and because uncertain H 2 /O 2 combustion mechanisms are readily available.

In short: We aim to propagate several sources of uncertainty through an LES of the H 2 Cabra flame. We will take into account the co-flow temperature (uncertain experimental measurement), the detailed kinetic mechanism of H 2 oxidation, and a closure coefficient of the turbulent combustion model.

Structure of the manuscript

• Chapter 1 presents a set of uncertainty quantification methods that will be useful in the following: random sampling, surrogate modeling, sensitivity analysis, dimension reduction. We hope that this chapter may serve as a point of entry to uncertainty quantification for future researchers. This chapter also summarizes the methods implemented in the in-house code SURE during this thesis.

• Chapter 2 is a preliminary study on auto-igniting homogeneous reactors and 1D laminar flames. Its purpose is to highlight the uncertainty of the auto-ignition delay in typical combustion configurations. It also serves as a first benchmark of the methods presented in Chapter 1.

• Chapter 3 presents this thesis' main case of study: the Cabra flame and our LES of this flame. It details the key equations and the choice of turbulence and combustion models which rely on tabulated chemistry and presumed filtered-density function, which were implemented in the AVBP LES solver. Finally, it presents the deterministic simulation results and their validation.

• Chapter 4 studies the uncertainty of the lift-off height of the simulated Cabra flame. Using a reduction of the initial uncertain space, a surrogate-based approach is leveraged to estimate the probability distribution of the lift-off height. 39 LES were run to train the surrogate and 20 more were run to serve as validation samples. This study allows the synthesis between seemingly disagreeing experimental observations and simulation results from several research groups.

• Chapter 5 extends the methodology introduced in Chapter 4 in order to quantify the uncertainties of every averaged LES field at every location in the computational domain. These local uncertainties are obtained by applying Principal Component Analysis to a training set of averaged flames and building surrogate models to estimate the weights of the most important modes. From these local field predictions, the uncertainty of global quantities, like the lift-off height, can be rediscovered.

Chapter 1

Handling uncertainties

This chapter gives an overview of the methods of the uncertainty quantification (UQ) methods employed during this thesis. It provides the necessary elements to understand and implement these methods, as well as references to further investigate the points of detail. As such, it is written in the form of a beginner's guide towards uncertainty propagation, in the hope that it may provide a helpful entry point for future Ph.D. students.

The framework of this thesis is forward uncertainty propagation, which means that we have information on the uncertainty of uncertain inputs, and we seek to determine statistics (probability distribution, sensitivity, or at least expectation and variance) of a Quantity of Interest (QoI).

Unless specified otherwise, the methods presented in this chapter have been implemented in the in-house UQ code SURE. This chapter begins with a quick overview of the different classes of methods that will be developed later. The second part focuses on how to sample an uncertain space. The third part is about surrogate modelling, or how to emulate a physical model given only a few observations. The fourth part deals with sensitivity analysis. Part 5 presents a few ways to reduce the uncertain dimension. Finally, part 6 introduces the error metrics used in this thesis to verify the quality of the methods presented in this chapter.

Methods overview

The methods employed in this thesis have been hand-picked to answer the specific needs of propagating kinetic and experimental uncertainties in a Large Eddy Simulation (LES):

• The simulation software are already coded with deterministic applications in mind, which means going for non-intrusive methods only.

• One observation in LES costs at least several tens of thousands of CPU-hours. Therefore, sample-efficient methods are wished for.

• The dimensionality of kinetic is high (cf Chapter 2), which means we need methods able to handle several tens of uncertain inputs, and/or ways to reduce the uncertain dimension 1.1.1 Deterministic methods

Linearization technique

A local gradient of the model is evaluated at the nominal input point. This allows performing a local sensitivity analysis by comparing the components of the gradient. Furthermore, if the input distributions are peaked enough, one can propagate uncertainties through a local linear surrogate of the model and retrieve an approximate output distribution, either by resampling or by analytical developments. On the one hand, this method is quite cheap, because the minimum number of evaluations of the model is d + 1 in a non-intrusive framework. Even better, when an adjoint solver is used, gradients are already part of the simulation results. For "simple" analytical models, the gradient can even be analytically retrieved. On the other hand, this analysis is local. It is not valid if the input distributions are wide enough that the uncertain inputs frequently exit the neighborhood of the point of gradient evaluation.

Perturbation techniques and Neumann expansions

These techniques obey the same idea as the linearization technique but function with higher-order derivatives, making them more precise. Still, they remain local.

Moments methods

The idea behind these methods is not to solve the model for a given set of inputs, but rather to solve a system of equations on the moments of the inputs to retrieve moments of the output. These methods suffer from a problem of closure since the equations on the moments are often not closed and closure terms have to be modeled. Furthermore, they enter the domain of intrusive methods and do not allow the use of a legacy deterministic solver. Therefore, they are out of the scope of this thesis.

Simulation methods

Simulation methods rely on the evaluation of the stochastic process on a large number of samples, either directly or indirectly, to approximate the Probability Distribution Function (PDF) of the uncertain output.

Sampling techniques

Sampling techniques, as their name suggests, rely on sampling extensively the uncertain input space according to its joint probability distribution, evaluating the physical model, and storing the quantity of interest for each sample. The probability distribution and relevant statistics of the QoI can then be retrieved from the data set. These methods are very robust and do not require any particular property from the physical model (non-continuous outputs are fine, for example), and they allow the use of legacy/deterministic solvers. However, they suffer from slow convergence rates, meaning that a considerable number of observations may be necessary to retrieve significant statistics. This is not a problem for cheap models featuring simple or reduced physics, but it will become an issue for more complex and expensive models.

Surrogate modeling

Surrogate modeling alleviates the cost of sampling techniques by building a surrogate, also called "response surface" or "emulator", using a few observations of the physical model. The surrogate is in turn extensively sampled instead of the physical model. As such, it is way cheaper to implement than a naive sampling technique, when the evaluation of the physical model is expensive. This cost may lay, for example, in the monetary costs of an experiment or in the computational power required to run a simulation. The surrogate, by construction, is very cheap to evaluate, which makes resampling affordable whatever the number of samples needed to obtain significant statistics. However, the precision of the resulting PDF relies entirely on the accuracy of the surrogate. Monte-Carlo simply consists in randomly sampling the input space and computing the corresponding outputs. An example of Monte-Carlo sampling on the unit hypercube, in 2D, is given in Figure 1.1. Let F (S ) be our quantity of interest, where S is the output of the physical model. The function F is defined depending on what you need to compute. For example, if you want to evaluate:

• the expectation of the output: F (S ) = S

• the second moment of the output:

F (S ) = S 2
• the probability of the value exceeding a given threshold:

F (S ) = 1 if S > Threshold 0 otherwise
Convergence of a Monte-Carlo estimator F (S ) is a function of the simulation S which is itself a deterministic function of the uncertain inputs. Therefore, F is a function of the uncertain inputs ξ and we will denote indifferently F (S ) = F (S (ξ)) = F (ξ) in the following.

The random vector ξ is defined on the set Ξ of joint probability π ξ . Let F N be the Monte-Carlo estimator of E[F ] defined by:

F N = 1 N N i=1 F (S (ξ (i) )) (1.1)
Where N is the number of input samples ξ (i) drawn according to the probability distribution π ξ .

The central limit theorem gives, in the asymptotic regime, the following property:

F N = E[F ] + ε N , with ε N ∼ N 0, V[F ] N (1.2)
Where V[F ] is the variance of F . Substituting the estimator V N for V[F ], we can show that once in the asymptotic regime, a Monte-Carlo estimator has a universal convergence rate of 1/2, meaning that the slope of the upper bound for 50% confidence on the residual ε N has a -1/2 slope in a log-log diagram, as shown in Figure 1.2. We can also say that it converges in O(1/ √ N ).

Relative error

Relative error of the MC estimator Upper bound for 50% confidence When trying to converge a Monte-Carlo estimate, it is good practice to make sure you at least double the number of samples when the error estimator is too high. This is because F N is random and fluctuates. Thus, adding just a few more samples will not necessarily help with the convergence and may actually worsen the estimator (see figure 1.2). Finally, it is important to note that the error estimate ε N is an absolute value. Only the relative error ε N / F N should be used to characterize the estimator's level of convergence.

Sampling randomly from an arbitrary distribution

Random numbers generators (or rather pseudo-random generators [START_REF] Ecuyer | History of uniform random number generation[END_REF]) are typically only capable of sampling random variables X that follow a uniform distribution over the segment [0, 1] (X ∼ U ([0, 1])). However, uncertain variables that do not follow a U ([0, 1]) distribution are plentiful. Thus, we need a way to sample from an arbitrary PDF. Two cases can be identified:

• When we know (or can approximate) its cumulative distribution function (CDF) and we can inverse it

• When we only have access to the value of the probability density at any value taken by the variable These two cases can be treated respectively by a measure change and Markov-chain Monte-Carlo sampling. These methods will be detailed in the following sections.

Measure change

The measure change is an analytical transformation that generates samples of a random variable using its CDF.

Let ξ be a one-dimensional random variable of cumulative distribution function F ξ . If F ξ is strictly monotonic, it is invertible and the random variable F -1 ξ (u), where u follows a uniform distribution on [0, 1], follows the distribution of ξ. This property is illustrated in Figure 1.3. In short:

u ∼ U ([0, 1]) ⇐⇒ F -1 ξ (u) ∼ ξ (1.3) F -1
ξ is called the quantile function of ξ. Some examples are given here for usual distributions:

• Normal distribution:

F -1 ξ (u) = µ + σ √ 2 erf -1 (2u -1) (1.4)
Where µ is the expectation of the distribution and σ its standard deviation

• Log-normal distribution:

F -1 ξ (u) = exp µ + σ √ 2 erf -1 (2u -1) (1.5)
Where µ and σ are defined by the associated normal distribution It follows that, for a given quantity of interest Q,

E[Q] = Ξ Q(ξ)π(ξ)dξ = 1 0 Q F -1 ξ (u) π(u)du = 1 0 Q • F -1 ξ (u)du (1.6)
A Monte-Carlo estimator of Q will then be:

E[Q] ≈ Q N = 1 N N i=1 Q • F -1 ξ (u (i) ) (1.7)
Or more generally, for any weighted estimator:

E[Q] ≈ Q N = N i=1 w (i) Q • F -1 ξ (u (i) ) (1.8)
where u (i) are samples of the random variable u and w (i) the weights associated with the integration point u (i) . This is an artificial way of sampling from the distribution of ξ because what is ultimately sampled is still u ∼ U ([0, 1]). However, it allows to efficiently evaluate an expectation estimator, and it is compatible with every Monte-Carlo convergence acceleration technique presented in the following. In the case of a random vector made of a number of independent random variables verifying the conditions for the measure change, the method can be generalized and is then called the Rosenblatt inverse transform [START_REF] Rosenblatt | Remarks on a Multivariate Transformation[END_REF].

This procedure works only in the case where the probability distribution is simple/canon enough so that a quantile function, or at least a good approximation of it, is available. In the multi-dimensional case, it works only if the components of the random vector are independent.

Markov-Chain Monte-Carlo

Markov Chain Monte-Carlo (MCMC) [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] is a method that allows drawing samples directly from an arbitrary distribution, of which we only need to be able to evaluate the probability distribution at any point. This means we do not need an analytical form for the Probability Distribution Function (PDF), Cumulative Distribution Function (CDF), or quantile function. MCMC is, therefore, the default choice for sampling non-standard distributions or random vectors whose components are not independent.

MCMC is based on the concept of Markovian processes, or Markov chains, which are memory-less random walks. In practice, the Metropolis-Hastings algorithm [START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] and Gibbs sampling [START_REF] Geman | Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images[END_REF] are widely used, but other, more advanced algorithms of Markov-Chain Quasi Monte-Carlo [START_REF] Pierre | A Randomized Quasi-Monte Carlo Simulation Method for Markov Chains[END_REF] can be used to combine the ability of sampling from any distribution with the advantages of Quasi-Monte-Carlo (presented in Section 1.2.3.2).

Improving the convergence of the Monte-Carlo estimator

Monte-Carlo estimators do not converge quickly. For this reason, it is common to try and accelerate this convergence, which means minimizing the residual ε N for a given computational effort. To achieve this you can either:

• Decrease the variance of the QoI (Linked to the intercept of the convergence curve)

• Increase the rate of convergence (The slope of the convergence curve)

Reducing the variance of the quantity of interest

Separation of the principal part

The idea here is to find a "principal part" H of F , that is, a function with known expectation E[H] that approximates the behavior of F . We then define the error Θ as: Θ(ξ) = F (ξ) -H(ξ). The expectation being a linear operator, we thus have:

E[F ] = E[H] + E[Θ].
We can then define a new estimator:

F N = Θ N + E[H]
(1.9)

This new estimator converges as σ(Θ)

√ N = σ(F -H) √ N ≪ σ(F ) √ N This idea is used in:
• Surrogate Variance Reduction / Reduced Model Variance Reduction (exactly the method presented)

• Multifidelity (degrade physics) and Multilevel (degrade discretization) Monte-Carlo methods.

The separation of the principal part is illustrated in Figure 1. [START_REF] Bachoc | Cross Validation and Maximum Likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF].

In practice, it can be simply implemented by spending a few samples to get the rough shape of F and construct a surrogate (for example a linear fit) to get a suitable H which is easily samplable/integrable. The better the surrogate, the faster the convergence.

Importance sampling

The idea of Importance sampling is to draw more samples in the areas where F has large values. That way, the integral computation is more efficient. This is done by using a variable change and an alternative distribution π * to draw the samples.

E π [F ] = Ξ F (ξ)π(ξ)dξ = Ξ F (ξ) π(ξ) π * (ξ) π * (ξ)dξ = Ξ G(ξ)π * (ξ)dξ = E π * [G]
(1.10)

Where G(ξ) = F (ξ) π(ξ) π * (ξ)
. Of course, π * Must be non-zero on Ξ. The only difficulty resides in finding a probability distribution π * such that

V π * [G] < V π [F ].
Analytically, the π * that minimizes V π * [G] is defined as:

π * (ξ) = |F (ξ)| π(ξ) Ξ F (η)π(η)dη (1.11)
The denominator of this expression is not known. Indeed, it is the absolute value of E π [F ], the final result of the study. However, equation 1.11 states that any distribution proportionnal to |F (ξ)| π(ξ) will be suitable. Fortunately, the Metropolis-Hastings algorithm, mentioned in Section 1.2.2.2, allows sampling from any arbitrary distribution known within a multiplication factor.

Increasing the convergence rate: Low discrepancy sampling

In this section, we will restrict ourselves to the sampling of uniform random variables on the hypercube [0, 1] d , without loss of generality.

Latin Hypercube Sampling

The idea behind Latin Hypercube Sampling (LHS) [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code[END_REF] is to control the sampling so that it covers the input space more evenly than Monte-Carlo, while keeping it random. In 1D, the segment [0, 1] is partitioned into N non-overlapping intervals having the same probability (equal lengths for uniform distributions). In d-D, the space is partitioned into N rows in each dimension, and the algorithm draws one and only one sample per row in each dimension. This is illustrated in Figure 1.5. In practice, what can be done is an LHS is 1D, then one random permutation of this initial sampling is done for each remaining dimension. For example, to draw 4 LHS points in the uniform unit cube (3D), we first draw 4 LHS points in the segment [0, 1]: {0.2, 0.37, 0.51, 0.82}. We then perform two random permutations of this list (two remaining dimensions), which gives us: {0.37, 0.51, 0.2, 0.82} and {0.82, 0.37, 0.2, 0.51}. Then, we read the three coordinates of the first 3D LHS point as the first elements of the three lists: (0.2, 0.37, 0.82). The second point is (0.37, 0.51, 0.37), and so on. This ensures better coverage of the uncertain space than brute-force Monte-Carlo. In terms of even coverage, an improvement on this method is orthogonal sampling [START_REF] Tang | Orthogonal Array-Based Latin Hypercubes[END_REF]. In orthogonal sampling, the input space is first divided into N equiprobable subspaces, and one sample is randomly drawn from each subspace in such a way that the total set of samples follows the rules of LHS.

Mc Kay et al. [START_REF] Mckay | A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code[END_REF] showed that LHS could not have a worse convergence rate than brute-force Monte-Carlo, and that its convergence is O(1/N ) if F (ξ) is monotonous in each of its arguments.

However, LHS presents one major issue. A fixed number of samples is decided a priori. In the case where the convergence of the quantity of interest's estimator is not satisfactory, increasing the number of samples while retaining the properties of LHS is not trivial. Tong [START_REF] Tong | Refinement strategies for stratified sampling methods[END_REF] gives useful pointers in that regard.

The implementation of LHS used in this thesis is the one proposed by the OpenTURNS library [START_REF] Baudin | Open-TURNS: An industrial software for uncertainty quantification in simulation[END_REF].

Quasi Monte-Carlo

Quasi-Monte-Carlo sampling (QMC) makes use of deterministic space-filling (or "lowdiscrepancy") sequences. They are designed to fill the space as uniformly as possible but lose in the process the randomness of standard Monte-Carlo sampling.

Their advantage is that the sampling is sequential, and the sequences are not periodic, so samples can always easily be added to the sample set if the estimator is not converged enough. This is a very interesting property compared to LHS. The easy refinement of the sampling comes at a cost, though, and this cost is the ). This is asymptotically always better than a brute-force Monte-Carlo estimator, but it may not be as good as LHS depending on the case. However, they are impaired by some important drawbacks:

• The sampling sports patterns. If the model is oscillatory, some features may be difficult to retrieve

• The efficiency of QMC tends to decrease with the increase of dimensionality, while the convergence rate of Monte-Carlo and LHS is independent of the dimension. In large dimensions, some parts of the sampling space will not be covered at all if the number of samples is not sufficient, because the exploration of the space is stratified. An optimal number of samples, for example N = 2 n for the Sobol' sequence, is required to maximize the uniformity of sampling.

• As QMC is purely deterministic, Monte-Carlo-style error estimates cannot be computed for this kind of sampling, because they rely on the intrinsic randomness of the sampling and the Central limit theorem. This drawback can be alleviated by the use of Randomized Quasi-Monte-Carlo sampling, which is presented in the next section.

The implementation of QMC used in this thesis is the one proposed by the OpenTURNS library [START_REF] Baudin | Open-TURNS: An industrial software for uncertainty quantification in simulation[END_REF].

Randomized Quasi-Monte-Carlo

The randomization of QMC allows the determination of the estimator's error in the same way as Monte-Carlo [START_REF] Tuffin | Randomization of Quasi-Monte Carlo Methods for Error Estimation: Survey and Normal Approximation[END_REF]. The randomization is not a trivial task if the space-filling properties of QMC are to be kept. Dumont's thesis [START_REF] Dumont | Méthodes numériques et modèle réduit de chimie tabulée pour la propagation d'incertitudes de cinétique chimique[END_REF] provides an interesting review of the different randomization methods.

Surrogate modeling

Surrogate modeling is the science of reducing the cost of sampling expensive functions by building "surrogates", also called "emulators" or "response surfaces". These surrogates are functions that should be cheap to evaluate and that give a good enough approximation of their underlying "true" physical model. At this point, you might be wondering one of the following three things:

1. How do I recognize a "good enough" approximation of my function? 2. How cheap is "cheap to evaluate"? 3. Theory is all very well, but where do I start in practice? First, a "good enough" approximation depends on the problem you want to tackle. It is an approximation that allows you to retrieve the quantities that you seek with sufficient accuracy. In practice, you might find that it relates to recovering the probability distribution's shape and moments of the QoI you want to emulate, but watch out, the level of accuracy is highly case-dependent. This is best illustrated with two -perhaps caricatural -scenarii: Scenario 1: Global performance analysis Say you are designing a brand new run-of-the-river hydroelectricity plant to power a city in a remote area. It happens that the only available river is capricious, meaning that its flow rate is far from constant. This flow rate does not appear to be predictable for now, but you can model it as a random variable of a given probability distribution based on previous observations. You then aim to propagate this randomness -this uncertainty -through a model of the plant to assess the uncertainty of your power output. Alas, the decision has been made to design a field of immersed stream generators to reduce environmental impact. The precise computation of the power output of such a field requires advanced -and expensive -CFD simulations due to the interactions between the stream generators. Thus, you need a surrogate model of the CFD simulation to obtain meaningful statistics for your power output at a reasonable cost. Now, what are you really trying to figure out with this study? The expectation of the power output, to be sure, and probably its standard deviation to design adequate buffers or annex generators. This means that while a very coarse surrogate only allows us to capture the first two moments of the power output distribution, it should be enough to suit your needs.

Scenario 2: Rare events prediction Still in the energy business, you are now working in a nuclear power plant as a safety engineer. Your goal is now to assess the probability of the plant blowing up due to, for instance, excessive pressure in the primary coolant circuit. Let's assume that the QoI is a function of the mechanical strength of the circuit and of the stress applied by the coolant's pressure. Uncertain inputs, on the other hand, refer to the composition, crystal structure and mechanical fatigue of the pipes, coolant temperature and pressure, and so on. Coupled simulations where heat transfer, fluid mechanics, and solid mechanics interact are very expensive. Therefore, you need a surrogate again to propagate uncertainties in this problem. Now, a nuclear plant blowing up is -hopefully -a rare event, so the events we seek to predict are not in the belly of a bell curve, but rather on the tails of the distribution -regions of low probability. Thus, the surrogate must be accurate enough so that the shape of the probability distributions, and in particular the tails, are well predicted.

Second, as a rule of thumb, a cheap function relies on an analytical formula, even a complicated one, ranging from simple piece-wise constant expressions to complex neural networks. The key here is that no complex iterative process should be involved in the evaluation of the surrogate, as these are generally expensive.

Last but not least, you can start building surrogates by using the two relatively common surrogate types used in this thesis. In the following, you will find a brief presentation of these surrogate models and how to build them.

Polynomial Chaos Expansion

Despite its daunting denomination, a Polynomial Chaos Expansion (PCE) is nothing more than a polynomial approximation of a stochastic process, which is itself a function of random variables. It must be noted that the random input variables must be independent to build a valid PCE. Polynomials are very practical functions to approximate continuous physical phenomena. It is in fact possible to demonstrate that PCEs have a spectral convergence towards any continuous stochastic process [START_REF] Le | Spectral Methods for Uncertainty Quantification[END_REF]. That is, the more polynomials in the expansion, the better the approximation of the process, provided that the weights are exactly computed.

Formally, a PCE is defined as follows: Let (F : ξ ∈ Ξ) ∈ R d → R be the physical model where ξ = (ξ 1 , . . . , ξ d ) are independent, real-valued, random variables of joint-PDF π ξ . Let α = {α i } ∈ N d be a multi-index and {ψ α } be the set of d-variate orthogonal polynomials:

ψ α (ξ) = d i=1 ψ (i) α i (ξ i ) (1.12)
If F ∈ L 2 (Ξ, π ξ ), its PCE exists and is written:

F (ξ) = lim N →∞ ∥α∥≤N ψ α (ξ)F α , ∥α∥ = d i=1 |α i | (1.13)
In essence, a Polynomial Chaos Expansion is therefore the combination of a polynomial basis {ψ α } and a set of weights {F α }. Here is how it shakes down:

1.3.1.

The basis

Building the polynomial basis

The polynomial basis must be an orthogonal family of polynomials. Being orthogonal means that the dot product between two basis elements is zero. You may remember that the dot product between the polynomials ψ 1 and ψ 2 is written as follows:

⟨ψ 1 , ψ 2 ⟩ = Ω ψ 1 (ξ) • ψ 2 (ξ) • π(ξ)dξ (1.14)
The nullity of this expression obviously depends on π(ξ), the probability distribution of the random vector ξ. Consequently, a different polynomial basis must be built for each distribution.

Case 1: Canonical distributions

For canonical distributions, canonical polynomial bases are available (see Table 1.1). These bases can be tensorized for multivariate PCEs. In the case where two random input variables are canonical but do not have the same probability distribution, one can (and must) compose the bases. For example, if ξ 1 is normal and ξ 2 is uniform, the polynomial basis will be the tensorization of Hermite and Legendre polynomials. In the case of arbitrary distributions, a custom polynomial basis can be built iteratively using the Gram-Schmidt orthonormalization process.

distribution Polynomials Support ξ ψ k (ξ) Continuous RV Gaussian Hermite R γ Laguerre R + β Jacobi [a, b] Uniform Legendre [a,

Truncating the polynomial basis

The polynomial bases constructed previously are infinite families of polynomials. The complete expansion -i.e. the infinite sum of correctly-weighted basis polynomials -is a perfect representation of the stochastic process. However, in practice, you must truncate the basis. This means keeping a certain number of relevant polynomials while discarding all the other ones. The most common ways of truncating the basis are:

1. Total order truncation 2. Hyperbolic truncation

Adaptive basis construction

Total order truncation This is the simplest of the truncation methods. All polynomials in the basis which are under a given total order are retained. This allows describing high-order interactions. However, this comes at the cost of rapidly increasing the cardinal of the basis. The number of basis polynomials M for a PCE of total order p+1 involving d uncertain variables is given by Equation 1.15.

M = (d + p)! d!p! (1.15)
As illustrated in Figure 1.7, the number of polynomials (and therefore the cost of finding the weights) reaches quickly unreasonable amounts for uncertain spaces of high dimensionality. This is called the "Curse of Dimensionality". Total-order truncation is therefore affordable only when the dimensionality of the problem is not too high.

Hyperbolic truncation

According to the "sparsity-of-effects principle", or "hierarchical ordering principle" [START_REF] Hamada | Experiments: Planning, Analysis, and Optimization[END_REF], systems are usually governed by main effects (involving a single variable) and low-order interactions. Therefore, neglecting high-order interactions should not result in a significant loss of representativeness for the PCE, in the majority of cases. A systematic way to do so is to perform a so-called hyperbolic truncation, as introduced by Blatman and Sudret in [START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF].

A quasi-norm, so called q-quasi-norm ∥•∥ q is defined as:

∀ 0 < q < 1, ∥α∥ q = d i=1 α i 1/q
and ∥α∥ 0 = d i=1

1 {α i >0} (1.16)
The truncated PCE of total order p and q-norm is then defined as:

F (ξ) ≈ Fp,q (ξ) = ∥α∥ q ≤p ψ α (ξ)F α (1.17)
The 1-quasi-norm is actually the absolute norm of N d and corresponds to the total order truncation. For 0 < q < 1, the truncation follows a hyperbolic-like pattern illustrated in Figure 1.8, which cuts presumably-unneeded high-order interactions from the PCE. By weighting the different directions in the definition of the quasi-norm, you can also get an anisotropic hyperbolic truncation, which can help describe in more detail the evolution of the QoI in certain directions relative to others.

Adaptive basis

To further complexify the construction of the basis and retain only the useful terms of the polynomial basis, Blatman [START_REF] Blatman | Chaos polynomial creux et adaptatif pour la propagation d'incertitudes et l'analysede sensibilité[END_REF] developed adaptive methods to determine conjointly the optimal basis polynomial and weights of the PCE based on Least-angle Regression [START_REF] Efron | Least angle regression[END_REF][START_REF] Blatman | Adaptive sparse polynomial chaos expansion based on least angle regression[END_REF].

The weights

There are several ways of computing the PCE weights, both intrusive and nonintrusive. In the following, we will discuss only the non-intrusive methods, which are projection and regression.

a. Projection

Because the {ψ α } form an orthogonal basis of the Hilbertian space, the weights F α can be found by simply projecting F on the basis:

F α = ⟨F, ψ α ⟩ ⟨ψ α , ψ α ⟩ = 1 ⟨ψ α , ψ α ⟩ Ξ F (x)ψ α (x)π ξ (x)dx (1.18)
The integral can then be approximated with a finite number of samples using a proper quadrature rule which is a proper numerical way to compute the integral: The choice of the quadrature rule (ξ (i) , w (i) ) is crucial as it conditions the number of samples needed to evaluate a given weight with a given precision. Most of these rules can be split into three main families:

F α ≈ 1 ⟨ψ α , ψ α ⟩ N i=1 w (i) F (ξ (i) )ψ α (ξ (i) ) (1.19)
• Brute-force integration

• Integration of polynomials that interpolate the integrand

• Subtler and more powerful change of variables-based approaches An extensive discussion and benchmarks of numerous quadrature methods can be found in Dumont's thesis [START_REF] Dumont | Méthodes numériques et modèle réduit de chimie tabulée pour la propagation d'incertitudes de cinétique chimique[END_REF]. This section will only give a quick overview of these families, as well as the formal definition of the method primarily used during this thesis.

Brute force integration

This approach relies on Monte-Carlo estimates of the integral obtained using the methods described in Section 1.2. This approach is obviously not very efficient in this context because a large number of observations (i.e. a large training set) will be needed to compute precisely the integrals, which defeats the purpose of surrogate modeling.

Integration of polynomials that interpolate the integrand

This family of rules relies on the exact (analytical) integration of polynomials that interpolate the integrand. It comprises the Rectangle method (polynomials of order 0), the Trapezoidal rule (order 1), and the Newton-Cotes formulas (general case). Obviously, the closer the interpolating polynomial is to the integrand, the more precise the evaluation of the integral. However, high-order polynomials tend to oscillate out of hand and do not necessarily improve the precision of the estimation. Hence the popularity of the Composite Simpson rule, which is a special case of the Newton-Cotes rule, where the integrand is approached by a piecewise second-order polynomial on a regular grid.

Change of variable-based approaches

This family includes the well-known Gauss quadrature rules, which are designed to return the exact integral of a polynomial of order 2n + 1 using only n quadrature points, i.e. n observations, which must be taken at very specific locations. The Gauss quadrature rules have the best convergence rate when the integrand is well approximated by polynomials. However, it is not possible to re-use observations when a more precise estimation of the integral is desired. For any refinement of the numerical integration, previous samples have to be put aside and entirely new observations have to be made. This is an important drawback because the number of quadrature points needed to get a "good enough" projection is not known a priori, in the general case. To tackle this problem, nested quadrature rules need to be used. "Nested" means that there are quadrature "levels" such that for a given level k, the quadrature level k + 1 reuses all the observations of level k, and a reduced number of additional samples are evaluated. This allows to refining the estimation of the integral by observing a reduced amount of new samples for each new quadrature level. This also allows sparse tensorization, which will be discussed in the following. The Clenshaw-Curtis [START_REF] Clenshaw | A method for numerical integration on an automatic computer[END_REF] rule and the Fejér quadrature of the second kind [START_REF] Fejér | Mechanische Quadraturen mit positiven Cotesschen Zahlen[END_REF] are nested rules that feature a similar convergence rate as Gauss quadrature rules [START_REF] Trefethen | Is Gauss Quadrature Better than Clenshaw-Curtis?[END_REF]. Clenshaw-Curtis is a closed quadrature rule, meaning that it includes the boundaries of the evaluation interval, while the second Fejér rule is an open rule, making it suitable to integrate a quantity over an open interval (typically R for normal random variables).

The evaluation points and the associated weights are presented for the two methods, without loss of generality, for an integration over the uniform interval [-1, 1]. Clenshaw-Curtis quadrature points and weights are given respectively by Equations 1.20 and 1.21:

x i = cos iπ N , i ∈ 0, N (1.20)              w i = c i N   1 - E[N/2] j=1 b j 4j 2 -1 cos 2jiπ N   b j = 1, j = N/2 2, j < N/2 , c i = 1, i = 0 or i = N 2, else (1.21)
The Fejér quadrature of the second kind uses the same evaluation points, excluding the boundaries (-1 and 1), with different weights. They are given respectively in Equations 1.22 and 1.23.

x i = cos iπ N , i ∈ 1, N -1 (1.22) w i = 4 N sin iπ N E[N/2] j=1 sin (2j-1)iπ N 2j -1 (1.23)
In the Clenshaw-Curtis rule, a quadrature of level k which uses N k = 2n + 1 evaluations re-uses all the evaluations of the quadrature with N k-1 = n + 1 points. For the second Fejér rule of level k, the quadrature with N k = 2n -1 points re-uses the evaluation of the quadrature with N k-1 = n -1 points.

For multi-dimensional integration, the quadrature rules need to be tensorized. Sometimes, the tensorized quadrature rule is called the "cubature rule". The full tensorization follows the following principles:

• For each quadrature point in a given dimension, the other dimensions must be sampled on each quadrature point. As a consequence using a quadrature rule which needs n evaluations in 1D, and tensorizing it in a space of dimension d requires making n d evaluations. This exponential growth of the number of required evaluations is also a curse of dimensionality.

• The weights of a multi-dimensional quadrature rule are computed as the product of the 1D quadrature weights used in each dimension

To alleviate the curse of dimensionality in the integration, a sparse tensorization method may be used, such as Smolyak's tensorization [START_REF] Smolyak | Quadrature and interpolation formulas for tensor products of certain classes of functions[END_REF]. More recently, Gernster and Griebel [START_REF] Gerstner | Numerical integration using sparse grids[END_REF] have proposed an adaptive sparse tensorization method to improve on Smolyak's idea and prioritize the precision of the cubature on certain dimensions relative to others. We will not dive further into the mathematical details of sparse tensorization. Instead, we refer to Dumont [START_REF] Dumont | Méthodes numériques et modèle réduit de chimie tabulée pour la propagation d'incertitudes de cinétique chimique[END_REF] for details on these methods and benchmarks of their respective performance. We showcase the dramatic decrease in evaluation cost, for a given level of quadrature, with Smolyak's method compared to the full tensorization in Figure 1.9.

Sparse cubature rules have some negative weights. Therefore, the positivity of the approximation of the dot product (defined in equation 1.14) is not guaranteed. While this may seem problematic from a mathematical point of view, we did not, in practice, encounter any unexpected behavior using Smolyak's tensorization. However, Constantine et al. [START_REF] Constantine | Sparse pseudospectral approximation method[END_REF] and Conrad et al. [START_REF] Conrad | Adaptive Smolyak Pseudospectral Approximations[END_REF] explored this issue. They proposed alternative sparse grid algorithms for orthogonal projection to ensure, in particular, the absence of internal aliasing (ensure that ⟨ψ i , ψ j ⟩ = δ i,j ).

b. Regression

The other family of methods is based on the minimization of the surrogate's error. We search the optimal PCE F of F such that: The choice of the distance is critical to get a good surrogate F . The basic choice is the Ordinary Least-Square method, which transforms the problem into:

F = argmin F * ∈Ψ E distance(F * (ξ), F (ξ)) , Ψ = Span({ψ α }) (1.
F = argmin F * ∈Ψ 1 N N i=1 ∥F * (ξ (i) ) -F (ξ (i) ))∥ 2 , Ψ = Span({ψ α }) (1.25)
Where N is the number of observations F (ξ (i) ) used to calibrate the PCE.

The least-square method can lead to an overfit of the surrogate if the number of samples is not large enough compared to the number of polynomials in the truncated basis. This phenomenon is illustrated in Section 1.3.1.3. A good rule of thumb is to use about three to five times as many samples as coefficients to calibrate.

To improve the quality of the PCE, you can also improve the distance function:

• By adding ε • p k=1 |F k | 2 , ε > 0 : This is called the "Tikhonov regularization" or "L 2 regularization" [START_REF] Tikhonov | Use of the regularization method in non-linear problems[END_REF]. It is a precision/variance trade-off, which means that it trades some discrepancy with the observations for lower PCE weights. This reduces the variance of the PCE and its tendency to overfit.

• By adding γ • p k=1 |F k |, γ > 0 : This is called "Compressed sensing" or "L 1 regularization" [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]. This promotes sparsity in the weights by pulling their values to zero.

• By adding both of the previous regularizations. This is called the "Elastic net" [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF].

The hyperparameters ε and γ can be adjusted using Leave-One-Out or Cross-validation algorithms [START_REF] Allen | The Relationship Between Variable Selection and Data Agumentation and a Method for Prediction[END_REF]. Finally, more advanced algorithms like Least Absolute Shrinkage and Selection Operator (LASSO) [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF] or Least-Angle Regression (LARS) [START_REF] Efron | Least angle regression[END_REF] can be used for the L 1 regularization.

This family of methods is much more computation-intensive than the projection method because there is a minimization problem to solve instead of simple arithmetics to evaluate. This minimization problem can even become overwhelming as the cost increases exponen-tially with both the dimensionality and the highest polynomial order. Still, for a similar sampling cost, regression-based PCE can yield better precision than their projection-based counterparts if the pitfall of overfitting is avoided. This performance topic is addressed in Chapter 2.

1.3.1.3 About that coveted spectral convergence... The spectral convergence of the PCE is no unicorn in real life. Whatever the evaluation process, one can not simply enrich the polynomial basis and hope the approximation will improve. The spectral convergence is achieved if and only if the coefficients given to each basis function are exactly computed. In practice, for non-trivial stochastic processes, this comes at a cost:

• For regression-based PCE determination, increasing the number of polynomialsincreasing the number of coefficients to determine -will lead to overfitting if new samples are not computed to train the surrogate. That is, the surrogate will eventually pass through every sample point but will oscillate between them. In this situation, the surrogate becomes too specialized to efficiently generalize the dataset and emulate the "true" physical model. This is illustrated in the left part of Figure 1.10.

• For projection-based PCE determination, the effect is a bit subtler. A quadrature method will always give a result for the coefficient of a given polynomial, so you can always compute as many weights as you want given a dataset. The trick is that the higher the degree of the polynomial, the more points are needed to accurately compute its coefficient. A PCE of high degree built with an insufficient number of points will sport ill-determined high-degree terms that will make the surrogate oscillate out of hand. This does not qualify as overfit, though, because the surrogate will not even pass through the sample points. This is illustrated in the right part of Figure 1.10.

-2 On the right-hand side, a PCE is trained by projection using a level 3 Fejer quadrature of the second kind (N = 7).

0

Kriging

Kriging, also known as Gaussian process regression, is another surrogate modeling method first introduced by Krige [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF]. Consider a set of training data

D = {(ξ (i) , y (i) ), i ∈ 1, N } (1.26)
Such that ∀i, y (i) = F (ξ (i) ).

The Ordinary Kriging approximates F with the mean of the stochastic process Y defined as:

Y (ξ) = µ + Z(ξ) (1.27)
Where µ is the sum of the global mean contribution and local variations are represented by the stationary Gaussian process Z(ξ). The mean contribution µ can be replaced by a deterministic function of the coordinates µ(ξ) when F has a principal part. This is the principle behind the Universal Kriging introduced by Matheron [START_REF] Matheron | Le krigeage universel[END_REF]. The following developments will be restricted to Ordinary Kriging. However, the Universal Kriging case is covered if F is replaced by

F ′ (ξ) = F (ξ) -µ(ξ).
The Kriging Y is defined, at each point ξ by its predictor ỹ(ξ) and its variance σ(ξ). ỹ(ξ) and σ(ξ) are determined using the vector y of the observations {y (i) } and the autocorrelation matrix R:

R ij (θ) = R(ξ (i) , ξ (j) , θ) (1.28)
where R is a correlation function, or "kernel", and θ is a set of correlation parameters, "hyperparameters", or "length-scales". Common one-dimensional kernel choices are given in Table 1.2 and illustrated in Figure 1.11. Table 1.2: Common choices of one-dimensional kernels. d ij = |ξ (i)ξ (j) | is the distance between the samples i and j.

Name

R(ξ (i) , ξ (j) , θ) Linear max(0, 1 - Multi-dimensional kernels, in dimension d, can be obtained by tensorizing 1D kernels:

d ij /θ) Exponential exp(-d ij /θ) Squared-Exponential exp(-d 2 ij /θ) Matérn 3/2 (1 + k 3/2 )exp(-k 3/2 ) k 3/2 = √ 3d ij /θ Matérn 5/2 (1 + k 5/2 + k 2 5/2 /3)exp(-k 5/2 ) k 5/2 = √ 5d ij /θ
R ij = d k=1 R(ξ (i) k , ξ (j) k , θ k ) (1.29)
The kriging predictor ỹ at any point ξ is defined as the expectation of the realizations of the gaussian process presented in Equation 1.27:

F (ξ) ≈ µ Ŷ (ξ) = ỹ (1.30)
And is computed with:

ỹ = μ + r T R -1 (y -1μ) (1.31)
Where y is the vector of the observations, and 1 is a vector of size N containing 1 for each coordinate. The prior estimation μ of the global mean µ is found with:

μ = (1 T R -1 1) -1 1 T R -1 y (1.32)
And the vector r is defined by:

r i = R(ξ, ξ (i) , θ) ∀i ∈ 1, N (1.33)
The variance of the gaussian process σ can also be computed for any given point ξ, with:

σ2 = σ 2 Ŷ (ξ) = σ2 1 -r T R -1 r + u 2 (1 T R -1 1) -1 (1.

34)

With:

u = 1 T R -1 r -1 and σ2 = 1 N (y -1μ) T R -1 (y -1μ) (1.35)
The fact that the variance of the gaussian process can be computed at any point means that credible intervals can be predicted along the expectation of the QoI.

The choice of an appropriate kernel and its proper calibration are extremely important, as illustrated in Figure 1.12. Ideally, in the case where a large number of observations is available, a variogram of the data can be plotted. A variogram, first defined in [START_REF] Matheron | Principles of geostatistics[END_REF], is a plot of the squared difference between the values of each couple of observations as a function of the distance between these observations. A kernel that fits the data can then be chosen and calibrated on this variogram. In the absence of a sufficient amount of data to plot a meaningful variogram, the hyperparameters can be determined using a Maximum-Likelihood approach [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF][START_REF] Bachoc | Cross Validation and Maximum Likelihood estimations of hyperparameters of Gaussian processes with model misspecification[END_REF]. As for the choice of a kernel, Santner et al. [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] recommend using the Matérn 5/2 kernel as a default choice. Furthermore, Schobi et al. [START_REF] Schöbi | Combining Polynomial Chaos Expansions and Kriging[END_REF] state that the influence of the kernel is low in small dimensional problems.

Kriging can also be trained on noisy data, by adding a regularization constant to the diagonal of the correlation matrix R. This regularization constant is called the 'nugget' variance, and translates a lack of correlation with distance between neighboring samples. This lack of correlation itself can result from uncertain realizations (noisy measurements, variance introduced by hidden variables...) or from a very rough underlying function. This nugget variance is the intercept of the variogram of the data, i.e. its value for a zero distance. Even in the case of non-noisy data, introducing a small regularization constant can also be used to facilitate the numerical inversion of R.

Finally, Kriging estimates can be used to efficiently resample the uncertain space to sharpen the surrogate model where it is needed. A state-of-the-art review of such techniques can be found in [START_REF] Fuhg | State-of-the-Art and Comparative Review of Adaptive Sampling Methods for Kriging[END_REF].

-2 On the left-hand side, a Regression-Kriging with a linear trend and a Matérn 5/2 kernel of hyperparameter θ = 0.5.

0
On the right-hand side, an Ordinary-Kriging (no trend) and a Matérn 5/2 kernel of hyperparameter θ = 0.2.

Preconditionning

In certain cases, the physical model does not have an agreeable-enough behavior for the accurate calibration of a surrogate model. For example, the model may not be very smooth or comprise a plateau. It may have saturation behaviors. In some cases, the emulator may have positivity constraints. In these cases, the naive application of a surrogate modeling method may not lead to a satisfactory result. To alleviate some of these issues, it may be interesting to apply a pre-conditionning to the observations of the physical model.

A pre-conditionning Φ is an invertible transformation of the observations such that:

F ′ (ξ) = Φ • F (ξ) (1.36)
And F ′ , the pre-conditionned model, has more interesting properties than F . The best pre-conditionning is case dependent but here are a few pointers:

• To ensure the positivity of a surrogate, it is common practice to use a logarithmic pre-conditionning: Φ = log

• If F ′ has a tighter spectrum than F , it will be more accurately represented by a low-order PCE

• The substraction of the physical model's principal part (cf sec. 1.2.3) can help with the quality of a Kriging surrogate. This is the idea behind Regression Kriging [START_REF] Matheron | Le krigeage universel[END_REF] and PC-Kriging [START_REF] Schobi | Polynomial-Chaos-Based Kriging[END_REF] 1

.4 Sensitivity analysis

Sensitivity analysis is a step sideways from the previous sections. While Monte-Carlo algorithms and surrogate modeling are primarily aimed at estimating a Quantity of Interest or its probability distribution given a set uncertain input, sensitivity analysis aims at determining the relative impact of the uncertain parameters on the QoI.

Local sensitivity analysis

Local sensitivity analysis, as its name suggests, investigates the immediate vicinity of the nominal conditions. The simplest method of this family is gradient analysis, which is very close to the linearization method discussed in Section 1.1.1.1. The gradient of the model F with regard to the components of ξ is computed at a nominal point ξ * . The components of the gradient may then be multiplied by the variability of their respective components of ξ and finally normalized to obtain local sensitivity coefficients.

The advantage of this method is that it is quite cheap because it requires a minimum of d + 1 observations of the model in the non-intrusive framework.

Its downsides are:

• The choice of the nominal point to compute the gradient is arbitrary: It may be more relevant, for a given case, to use either the expectation (ξ = E[ξ]), the median value ( ξ), or the location of the most probable value of the uncertain input.

• Local derivatives may not fairly represent the variability of the model on the whole input space. This is the same issue as with the other local UQ methods discussed in Sections 1.1.1.1 and 1.1.1.2.

• Interaction effects cannot be retrieved.

Some of these issues may be mitigated by using higher-order derivatives, but the locality of the method remains an important limitation. For this reason, advanced UQ studies tend to use global sensitivity analysis.

Global sensitivity analysis

Global sensitivity analysis aims at filling the gaps of local sensitivity analysis. Ergo, its desirable properties are:

• Taking into account the behavior of the QoI on the whole range of possible values of the uncertain inputs

• Taking into account the probability distribution of the uncertain inputs

• If possible, representing the effect of interactions between random inputs

A powerful metric for sensitivity analysis is the set of Sobol' indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF]. Each index represents the portion of the variance of the QoI explained by the associated variable or group of variables, which is precisely what we need to discriminate uncertainties of consequence. In most practical cases, the Sobol' indices are not accessible analytically, and two main methods co-exist to approximate them numerically.

Brute-force approach

The first approach, originally proposed by Sobol' [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF] is a specialized Monte-Carlo algorithm that allows for a precise and controlled approximation of any Sobol' index. The downside is that a specific design of experiment is necessary to evaluate each Sobol' index. The number of evaluations needed to evaluate n S Sobol' indices is N (n S + 1) using this method, where N is the number of evaluations needed to estimate the QoI at the desired level of convergence. This very large number of evaluations makes this approach impractical when several indices need to be computed and/or when observations are expensive. Saltelli et al. [START_REF] Saltelli | A Quantitative Model-Independent Method for Global Sensitivity Analysis of Model Output[END_REF] proposed an alternate method called FAST to evaluate all the first-order Sobol' indices using only the N evaluations used to estimate the QoI. While this approach is much cheaper for the first-order Sobol' indices (the variance explained by one, and only one, variable at a time), it does not allow to compute higher-order indices, meaning the influence of interacting effects of two or more variables.

Indirect approach via the PCE

The other approach, proposed by Schöbi and Sudret [START_REF] Schöbi | PCE-based imprecise Sobol'indices[END_REF], uses a PCE of the quantity of interest. Algebraic operations are performed on the weights of the PCE to find all the Sobol' indices at once. This is an obvious advantage over the other method. However, the precision of the estimation of the Sobol' indices depends on the fidelity of the PCE. This entails that the precision of this method is less controlled. Furthermore, when the uncertain dimension is high, the maximum order and the number of interactions considered must be limited so that the determination of the PCE is tractable. Still, this method gives access to sound approximations of the Sobol' indices and allows ranking the effects of all the input variables, which is the key feature of sensitivity analysis when trying to reduce the uncertain dimension. Ergo, this second approach is the one that we selected in this thesis.

To simplify notations, let us assume that the components ξ (i) of ξ are independent uniform random variables on the [0, 1] interval. That way, no pdfs will show in the following formulas. F is then defined on the unit hypercube U d . Let us also assume that F is square-integrable, that is:

F ∈ L 2 (U d ) ⇐⇒ U d F (ξ) 2 dξ < ∞ (1.37)
The Sobol' indices of the QoI appear in its Sobol'-Hoeffding (S-H) decomposition, which writes:

F (ξ) = F (ξ 1 , . . . , ξ d ) = F ∅ + d i=1 F i (ξ i ) + d i=1 d j=1 F i,j (ξ i , ξ j ) + . . . + F 1,...,d (ξ 1 , . . . , ξ d ) (1.38)
Or in a more compact form:

F (ξ) = ι⊆D F ι (ξ ι ) (1.39)
With F the QoI, ξ the vector of its uncertain inputs and D the set of indices of the components of ξ, ι a subset of D, ξ ι the corresponding input vector, and F ι the associated S-H functional.

The Sobol' indices are defined for each group of variables ι as:

S ι = V[F ι ] V[F ] (1.40)
Where V[X] is the variance of the functional X.

To approximate the Sobol' indices, the S-H functionals can be approximated by their corresponding truncated Polynomial Chaos (PC), that is:

F ι (ξ ι ) ≈ α∈Aι ψ α (ξ ι )F α (1.41)
Where ψ α is the α-th basis polynomial, F α its weight in the PC expansion of F , and A ι the set of indices of PC basis polynomials that involve the group of variables ι in the truncated PCE.

Replacing Equation 1.41 in Equation 1.40, we get:

S ι (F ) = V[F ι ] V[F ] = ⟨F ι , F ι ⟩ ⟨F, F ⟩ ≈ ⟨ α∈Aι ψ α F α , α∈Aι ψ α F α ⟩ ⟨ α∈A ψ α F α , α∈A ψ α F α ⟩ (1.42)
Where A is the set of indices of all PC basis polynomials in the truncated PCE.

Finally, the orthonormality of the PC basis gives:

S ι (F ) ≈ α∈Aι F 2 α α∈A F 2 α (1.43)
This is a simple algebraic formula which uses only the coefficient of a PCE of the QoI, which makes the evaluation of any Sobol' index very cheap once a PCE has been computed.

Uncertain dimension reduction

When the dimensionality of an uncertain problem is too high, some propagation methods (other than Monte-Carlo whose convergence is dimension-independent) may not be tractable. In particular, as a rule of thumb, surrogate modeling with PCE or Kriging becomes too expensive past about ten dimensions. As a perspective, when dealing with chemical kinetic uncertainties, the number of uncertain parameters ranges from a few tens, for simple fuels, to several hundred when complex hydrocarbons are involved. In this context, reducing the uncertain dimension is crucial to enabling uncertainty propagation methods. Two categories of methods can be distinguished: supervised and unsupervised methods.

Supervised methods

Supervised methods are applied to problems where the dataset comprises both the coordinates {ξ (i) } and values {y i (i)} of the observations. In the context of uncertain dimension reduction, it means determining which components or combination of components of ξ are important to explain the variations of y.

Additive problems

An additive problem is the simplest possible case for uncertain dimension reduction. It means that the uncertain inputs have an independent impact on the QoI. This translates into cross-sensitivity indices being null or negligible. In terms of Sobol' indices, this can be formalized by:

S ι (F ) ≪ 1 , ∀ι | Card(ι) > 1 (1.44)
If this condition is verified, the interactions between uncertain variables can be neglected and an uncertainty propagation study can be performed for each dimension independently, which greatly reduces the cost.

Dominant variables

An uncertainty propagation problem with dominant variables is another relatively simple case to deal with. As its name suggests, it is a problem where a few uncertain variables explain most of the variance of the QoI. The starting point of dealing with such a problem is to rank the effects of variables and their interactions according to a relevant metric, for example, their Sobol' indices. Then, it may be possible to retain only the most prominent uncertain variables and interactions in the study. Several courses of action can be considered for this segregation:

• Choosing an effect threshold (for example: "A variable or interaction should explain at least 1% of the variance to be kept in the study").

• Choosing a total portion of explained variance in the study: Sobol' indices are summed in decreasing order until this chosen portion of variance is reached. All other variables and interactions are discarded.

• Spotting a significant drop in the values of the Sobol indices sorted by decreasing order. Variables and interactions prior to the drop are kept. All the others are discarded.

Of course, this approach requires a priori knowledge of the behavior of the physical model. This can be achieved, depending on the case, by relying on previous knowledge of the model, by running an a priori study with degraded physics or resolution, by examining an approximate analytical model, etc.

Active subspace

An active subspace (AS) is a set of directions of the total uncertain space Ξ on which the QoI varies most [START_REF] Constantine | Active Subspace Methods in Theory and Practice: Applications to Kriging Surfaces[END_REF]. This subspace may be aligned with the directions defined by the original uncertain variables. In this case, the "Dominant variables" approach of Section 1.5.1.2 applies. If the subspace is not aligned, a rotation can be applied to the uncertain space by performing a base change. The goal of this base change is to sort the directions in which the QoI F encounters the greatest gradients. This rotation of the basis is found with:

C F = Ξ (∇ ξ F (x)) (∇ ξ F (x)) T π ξ (x)dx = W ΛW T (1.45)
Where W = (w 1 , ..., w d ) is a unitary matrix defining the rotated basis and Λ the diagonal matrix of eigenvalues sorted in decreasing order:

Λ = Diag(λ 1 , ..., λ d ) | ∀i ∈ 1, d -1 , λ i > λ i-1 (1.46)
The set {w 1 , ..., w d } defines a new basis of the uncertain space. This basis can be truncated at the first rank r such that λ r ≫ λ r+1 . Span(w 1 , ..., w r ) is called the active subspace.

The expectations for an AS are that it should explain most of the variability of F and that a surrogate model built on this active subspace should be a good approximation of F .

In practice, the necessary number of observations M to accurately determine an active subspace is

M = α βlog(d) (1.47)
in the case of an adjoint simulation where the gradient is obtained "for free" with each sample, and

M = α β • d • log(d) (1.48)
when the gradient information is obtained by finite differences, where d is the initial dimension of the problem, α an over-sampling factor, and β the largest dimension of active subspace acceptable.

In these situations, the cost of finding an active subspace is untenable when observations are expensive. Constantine et al. propose an alternate solution [START_REF] Constantine | Exploiting active subspaces to quantify uncertainty in the numerical simulation of the HyShot II scramjet[END_REF] if certain conditions are met [START_REF] Cook | Regression graphics: Ideas for studying regressions through graphics[END_REF]:

• F is monotonous in each direction ξ i • An active direction (i.e.

an AS of dimension 1) can be discovered

Then, an active direction can be determined with a linear regression of F :

F (ξ) = b T ξ + b 0 + ε F (ξ) (1.49)
The active direction w 1 is expressed as:

w 1 = b T /∥b∥ 2 (1.50)
And the cost of finding such a direction is only:

M = α • d (1.51)
Finally, if this method was actually appropriate, a single direction w 1 on which most of the variability of the QoI F is found. A reasonable low-dimensional approximation can be written as:

F (ξ) ≈ G(w T 1 ξ) (1.52)
A surrogate-based study is then very easy to conduct in this new 1D space.

Finally, Constantine and Diaz [START_REF] Constantine | Global sensitivity metrics from active subspaces[END_REF] have proposed to derive global sensitivity indices from the determination of the active subspace. When the AS is one-dimensional, they recommend using the components of the eigenvector w 1 . These components provide insight into both the magnitude and sign of the effect of each original variable on the QoI. When the AS in multi-dimensional and the brute-force approach is used to determine the eigenvectors {w i } and associated eigenvalues {λ i }, they recommend using their "activity scores" α i defined as:

α i = r j=1 λ j w 2 i,j , i ∈ 1, r (1.53)

Unsupervised methods

Unsupervised methods used datasets where only observations are available. They try to identify patterns in the dataset to build a more concise representation of the data.

Principal Component Analysis

Principal Component Analysis (PCA) [START_REF] Pearson | On lines and planes of closest fit to systems of points in space[END_REF], also called discrete Karhunen-Loève Transform (KLT), Proper Orthogonal Decomposition (POD), or Singular Value Decomposition (SVD) is a statistical tool to reduce, or compress, the dimensionality of a dataset while controlling the variance left on the table by the compression.

In concrete terms, PCA is the process of finding the principal components, or modes, of a dataset. The i-th principal component is the direction of the line that best fits the dataset, with the constraint of being in the subspace orthogonal to the first i -1 principal components. This is illustrated in Figure 1. [START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF].

In practice, we used the SVD formalism implemented in the numpy.linalg library [START_REF] Harris | Array programming with NumPy[END_REF]. For a dataset matrix M of shape (n,d), where n is the size -or the number of samples of the dataset -and d is the dimension -or the number of features -we can write:

M = U * Σ * V T (1.54)
Where * denotes the matrix multiplication, Σ is the diagonal matrix of the singular values of M, and U and V are the singular matrices, containing orthonormal vectors in their rows and columns respectively.

The principal components are the rows of P , such that:

P = Σ 1/2 * V T (1.55)
The portion of variance V i explained by the mode i is found with:

V i = Σ 2 i,i n (1.56)
The new basis of modes can then be truncated to compress the data optimally -in the sense of retained variance. The possible truncature strategies are the same as those mentioned in Section 1.5.1.2. This process is useful to compress correlated inputs. An excellent example of this is [START_REF] Li | Quantifying initial and wind forcing uncertainties in the Gulf of Mexico[END_REF] where the authors reduce an uncertain field of data (hundreds or thousands of correlated uncertain initial conditions) into a few uncertain coefficients associated with PCA modes.

Following the same idea, instead of propagating uncertainties for a very large number of correlated QoIs (such as fields of data in a fluid flow, for example), it is possible to compress these uncertain outputs and perform the propagation of uncertainties on the few coefficients associated to the modes of the output.

Finally, PCA can be used to "clean" data by projecting noisy observations on the truncated basis of modes, which allows representing most of the variance of the data. This cleaning must be exercized with caution because it may erase useful data contained in the original, "noisy" observations. PCA is very sensitive to the amplitude of variations in each direction. Therefore, a proper normalization of the data is required to get the best of a PCA.

The formalism of PCA and Active Subspaces are strikingly similar, at least at first glance. In both cases, particularly interesting directions of the initial uncertain space are derived, ranked, and selected using the factorization of a certain matrix. However, the nature of the matrix that is factorized is extremely different between those cases. As a result, this allows both methods to pursue completely different objectives:

• In the case of the AS: The matrix C F , defined in Equation 1.45 is a square and symetric matrix of size d×d, where d is the dimension of the input space of the scalar function F . The components of C F represent the expected squared derivative of F in each direction in the original basis. Its diagonalization highlights new directions of the input space in which the scalar F varies most.

• In the case of the PCA: The matrix M, defined in Equation 1.54, is rectangular of size n × d, and is the juxtaposition of n realizations of a random vector of R d . Its factorization highlights a subspace of R d where the n realizations can be projected with a minimal loss of variance. Each basis vector of this subspace is also called a "mode" of the dataset.

V. Compression

Error metrics

The methods presented in Sections 1.3 and 1.5 have in common that they induce a loss of information compared to the reference uncertain problem:

• A surrogate model blurs the features of the physical model it emulates to avoid overfitting. It may also introduce artifacts, like the oscillations of high-order polynomials. Furthermore, any prediction from a surrogate model is only "the best guess" that we could make given the training observations and prior knowledge. A new observation of the physical model at a new coordinate could always yield a different, surprising value.

• A dimension reduction neglects the influence of some variables or combination of variables.

A visual verification: the summary plot

A summary plot of a method is a comparison of observations and predictions for a set of samples. It is mainly applicable to surrogate models, but it can also be used in the case of some dimension reduction methods. Figure 1.14 shows two examples of summary plots. For each sample, the original model observation is plotted in x-axis, while the prediction is reported in y-axis. If the method which is investigated does not lose or modify any information, all the points collapse on the first bisector (x = y). In practice, most predictions are not perfect and the points are not exactly on the first bisector. This is the case of the plot on the left side of Figure 1.14, were a "good" surrogate is evaluated. On the right side, however, a very poor method (i.e. with a lot of information loss or modification) is evaluated. While the summary plot gives only visual, qualitative information on the agreement between predictions and observations, it is enough to rapidly rule out blatantly substandard surrogates or dimension reductions.

Distance of probability distribution functions

For more quantitative evaluations, a reference metric is the distance between the predicted PDF and the PDF of observations. Obviously, this metric is only appropriate when observations are sufficiently cheap to allow an extensive sampling of the uncertain space.

The metric used in this thesis is the Wasserstein distance between the reference distribution for a given case -i.e. the result of a brute-force Monte-Carlo study on the original model, accounting for all the uncertainties -and the distribution predicted by the method currently studied. The absolute value of the distance does not mean much in itself. However, comparing the distances obtained with different methods allow ranking them by precision: the lower the distance, the more precise the method.

The Wasserstein distance, also known as the "earth-mover's distance" or "Kantorovitch distance", was first introduced by Kantorovitch in [START_REF] Kantorovich | Mathematical Methods of Organizing and Planning Production[END_REF]. It can be understood as the cost of changing the value of predictions to match the distribution of observations. The higher the number of observations to change and the higher the change amplitude, the higher the distance. Its precise mathematical definition can be found in [START_REF] Vaserstein | Markov Processes over Denumerable Products of Spaces, Describing Large Systems of Automata[END_REF]. Finally, this choice of metric is an arbitrary choice from our part, and other error metrics like the energy distance give similar results.

Probability distributions can also be compared by evaluating and comparing their successive moments: expectation, variance, skewness, kurtosis, etc.

Performance estimation with validation sets of low cardinality

When observations are more expensive, a reference distribution may not be available. Therefore, another metric needs to be used in order to check a method's quality. In this thesis, we used two such metrics: the Normalised Root Mean Squared Error (NRMSE) and the Normalised Maximum Absolute Error (NMAE) between a set of validation observation and the predictions made by the method at the same coordinates. These metrics are defined as:

NRMSE = ∥Y obs -Y pred ∥ 2 ∥Y obs ∥ 2 (1.57) NMAE = ∥Y obs -Y pred ∥ ∞ ∥Y obs ∥ ∞ (1.58)
Where Y obs is the vector of the validation observations and Y pred is the vector of the predictions at the same coordinates. The NRMSE, whose definition is given in Equation 1.57, measures the global fit of the predictions with the observations. The NMAE, defined by Equation 1.58, measures the maximum discrepancy between the predictions and the observations. Ergo, it is more of a local indicator of the maximum error that may cripple the predictions. On the one hand, the NRMSE is an efficient -though global -measure of the surrogate fit, even for validation sets of low cardinality, but it may overlook significant local discrepancies between the predictions and the observations. On the other hand, the NMAE will catch these local errors of prediction, provided that the validation set covers the relevant areas. Still, the NMAE is blind to the fact that the prediction may be excellent over most of the uncertain space in the event that it catches one important discrepancy. This may be especially penalising for polynomial surrogates, which tend to oscillate on the fringes of the uncertain space while having a good fit in high-probability regions.

Both metrics have their pros and cons, but their combination gives a rather complete overview of the quality of a surrogate's predictions, both in terms of global fit and local errors.

Cross-validation

When a separate validation dataset is not available, or when trying to optimize a set of hyperparameters with regards to a cerain error metric, cross-validation [START_REF] Stone | Cross-Validatory Choice and Assessment of Statistical Predictions[END_REF][START_REF] Stone | An Asymptotic Equivalence of Choice of Model by Cross-Validation and Akaike's Criterion[END_REF] can be used. Cross-validation is a resampling method, where the available dataset is split randomly in a training and a validation dataset. The model is trained on the former and the error metric is evaluated on the latter. To limit the dependence of the error evaluation on the split of the dataset, it is good practice to repeat this process a number of times and average the results.

A commonly used variation of this method is k-fold cross validation, where the dataset is randomly split in k mutually exclusive subsets of equal size. The model is trained k times using k -1 subsets and evaluated on the remaining subset. Finally, Leave-One-Out (LOO) is an extreme form of k-fold cross-validation, where k = n, the number of samples in the dataset.

Chapter 2

Uncertainty propagation in H 2 /O 2 auto-ignition simulations

The goal of this chapter is twofold. First, the importance of quantifying uncertainties in realistic hydrogen combustion simulations will be substantiated. Second, an efficient framework will be investigated to tackle kinetic uncertainty propagation in CPU-heavy combustion simulations. The case study will be a 1D strained and auto-igniting nonpremixed flame. This configuration emulates the physical complexity of more expensive combustion simulations while still being affordable to sample with a Monte-Carlo method. This property is essential as the Monte-Carlo study will provide a reference distribution against which different uncertainty propagation methods can be measured, which allows for validating the developed framework. With the help of the framework developed in this chapter, uncertainty quantification in CPU-expensive simulations will be done in Chapters 4 and 5.

The crux of the framework is to build a surrogate model, or response surface, which is an analytical formula that approximates sufficiently well the behavior of the costly simulation for a negligible fraction of the computational effort. This surrogate, in turn, lends itself perfectly to extensive Monte-Carlo resampling from which any statistic of interest can be retrieved.

The challenge is to build a satisfying surrogate model using as few training samples as possible because each training sample needs a whole simulation to be run. While this may not be particularly difficult in low-dimensional uncertain spaces, it certainly becomes a hurdle for high-dimensional spaces. De facto, this is the case when considering an uncertain kinetic mechanism, even as "simple" a mechanism as one for hydrogen combustion. To tackle this problem, we perform an a priori uncertain dimension reduction using a simplified physical problem that features the same dominant physical phenomena as the original one.

The general framework is then a two-stage process. First, we define the simplified physical problem, which can undergo an extensive Monte-Carlo simulation, and identify the dominant uncertainties among all the input uncertainties. Only uncertainties leading to significant amounts of variance in the quantity of interest are retained. This truncation of this uncertain space results in a much-reduced set of uncertain variables, making surrogate building manageable with few training samples in the second stage. Four surrogate modeling techniques (projection-PCE, regression-PCE, Kriging, and PC-Kriging) are compared to propagate the uncertainty in the 1D igniting flame case using a minimal number of simulation runs.

The following sections detail the corresponding framework and their results on the 1D strained auto-igniting non-premixed flame. Section 2.1 presents the strained auto-igniting non-premixed 1D flame configuration along with its simplified physical proxy, an autoigniting homogeneous reactor. A qualitative comparison of the results from several kinetic mechanisms for H 2 combustion is given in Section 2.2. The considered high-dimensional uncertain space is presented in Section 2.3. The methods chosen for the study are very briefly presented in Section 2.4, which refers to Chapter 1 for all the mathematical details of the methods. An extensive Monte-Carlo study is performed for the simplified setup, and an analysis of the resulting uncertainties and underlying physical phenomena is presented in Section 2.5. In particular, quantitative sensitivity analysis is performed in Section 2.5.2 to reduce the uncertain dimension space. Finally, uncertain dimension reduction and surrogate modeling techniques are applied to the target configuration in Section 2.6, which presents and analyses the propagated uncertainties.

Physical configurations and numerical setups

This section presents the two configurations explored in this chapter and the numerical setups which simulate them. The complex configuration, a strained auto-igniting nonpremixed 1D flame, is presented in Section 2.1.1. Quantifying uncertainties in this configuration using as few realizations of the simulation as possible is the goal of this chapter. Then, a simplification of this setup is presented in Section 2.1.2. This simplified configuration, an auto-igniting homogeneous reactor, will be helpful to analyse the physics and uncertainties that dominate hydrogen auto-ignition in a vitiated mixture. In the following, this analysis will allow sample-efficient uncertainty propagation methods to be used on the complex case.

The target configuration: a strained and auto-igniting nonpremixed 1D flame

This section presents a setup that emulates complex combustors. The interest is double: (i) we will verify that a non-trivial uncertain behavior can be found in realistic setups; (ii) we will benchmark several UQ methods to handle this uncertain problem using only a few samples, while still being able to afford the complete reference distribution, using Monte-Carlo sampling.

Table 2.1: Nominal boundary conditions for the auto-igniting 1D flame: temperature and species molar fractions. The pressure is set to 101 325 Pa.

Fuel side

Oxidizer side

T fuel [K] 305 T ox [K] 1020 X H 2 [-] 0.25 X H 2 O [-] 0.1 X N 2 [-] 0.75 X O 2 [-] 0.15 X N 2 [-] 0.75
The setup that will be studied in the following is a 1D auto-igniting non-premixed H 2 /O 2 flame. Details of the boundary conditions are given in Table 2.1 and are similar to conditions found in the Cabra H 2 flame [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF], only with a slightly colder oxidizer side.

These conditions mimic what is typically found in the later stages of staged industrial burners and combustion chambers. The initial conditions are given by the pure mixing solution from the boundary conditions, which corresponds, in the mixture fraction space, to a linear evolution of the mass fraction composition and enthalpy from one boundary to the other. A unity Lewis number is assumed so that this 1D flame is representative of turbulent combustion [START_REF] Peters | Cambridge Monographs on Mechanics[END_REF]. The corresponding evolution of species mass fractions Y k (with unity Lewis number) in the mixture fraction space is then driven by Equation 2.1:

ρ ∂Y k ∂t = ρχ 2 ∂ 2 Y k ∂z 2 + ρ ωk (2.1)
With:

χ(z) = a π exp -2 erf -1 (2z -1) 2 (2.2)
Where ρ is the density of the mixture, z is the mixture fraction, ωk is the reactive source term of species k, and a is the strain rate of the flame. A similar equation can be written for the evolution of the temperature. In the early stages of auto-ignition, if the transport term -proportional to the strain rate a -is too large compared to the source term ρ ωk for the intermediate species like O, H, and OH, the flame never ignites. A low strain rate is chosen to ensure the ignition of all sampled flames, even for low initial temperatures or slow realizations of the kinetic mechanism: a = 10 s -1 . This allows a more straightforward comparison between surrogate modeling methods, but similar results can be obtained with higher strain rates by conditioning the surrogate models to the observations where the flame actually ignites.

The 1D flame equations are solved on a uniform 400-point mesh using the 2 nd -order stabilized explicit Runge-Kutta method (ROCK2) with error control [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF]. The autoignition delay is determined when the mass fraction of the hydroxyl radical reaches the threshold Y OH = 600 ppm in the computational domain. The computation time per sample can rise up to a dozen minutes, making it just affordable in Monte-Carlo. This setup, while featuring species and heat transport and being the backbone of flamelet-based tabulated chemistry for 3D simulations [START_REF] Pitsch | Unsteady flamelet modeling of differential diffusion in turbulent jet diffusion flames[END_REF][START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF], is driven by auto-ignition. Its temporal behavior is presented in Figure 2.1. Starting from the pure mixing line, the gas temperature first starts to rise at the most-reactive mixture Z MR = 0.04 [START_REF] Mastorakos | Ignition of turbulent non-premixed flames[END_REF], which is a very lean and hot mixture. Then, heat and reactive species are transported towards the slower-reacting mixtures, progressively igniting the whole flame.

The simplified configuration: Auto-igniting homogeneous reactor

This section presents an auto-igniting constant-pressure homogeneous reactor, which will serve as a physics-reduced simplification of the setup presented in section 2.1.1. Its initial composition and temperature are chosen as the conditions initially reigning at the most reactive mixture Z MR of the non-premixed flame. Thus, it represents the chemistrydominated part of the evolution of the 1D flame, where ignition first occurs. As shown in Table 2.2, the initial temperature is already quite hot and the initial composition includes both N 2 and H 2 O. This case is much simpler, as it features a purely chemical behavior and omits the "transport" part of the equations that govern the evolution of the physical quantities (Equation 2.3). Therefore, it is a suitable reduced model of the costly setup when studying only chemical kinetics.

ρ ∂Y k ∂t = ρ ωk (2.3)
Table 2.2: Nominal initial conditions for the auto-igniting homogeneous reactor (Pressure, Temperature and species molar fractions)

P [Pa] 101325 T [K] 987 X O 2 [-] 0.1424 X H 2 [-] 0.0126 X H 2 O [-] 0.0950 X N 2 [-] 0.7500
Similarly to the 1D flame setup, auto-ignition is detected when Y OH = 600 ppm. The temporal integration is performed using the stiffly-accurate implicit Runge-Kutta method of fifth-order accuracy (Radau IIA) with error control [START_REF] Wanner | Solving Ordinary Differential Equations II. Stiff and differential-algebraic problems[END_REF]. The computational cost per simulation is about 0.1 CPU-second. This makes this setup suitable for a priori studies with extensive Monte-Carlo sampling. As shown in section 2.5, this minimalistic setup (H 2 /O 2 chemistry only) is already enough to showcase a non-trivial uncertain behavior.

A qualitative comparison of deterministic H 2 /O 2 mechanisms

This section provides an (admittedly non-exhaustive) overview of the variety of ignition delay times (IDT) that can be found in the case of the auto-igniting homogeneous reactor presented in Section 2.1.2. [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF], and Konnov's 2019 [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] deterministic mechanisms, at different initial temperatures.

On the one hand, this Figure shows that the IDT varies widely with the initial temperature of the reactor, with two different regimes clearly emerging: "weak" or "mild" combustion at low initial temperatures and "strong" or "sharp" combustion for high initial temperatures [START_REF] Saytzev | Study of combustion of an adiabatically-heated gas mixture[END_REF][START_REF] Yetter | A combined stability-sensitivity analysis of weak and strong reactions of hydrogen/oxygen mixtures[END_REF]. In industrial systems and, in particular, in the context of staged combustion, this initial temperature may not be exactly known, leading to potentially large uncertainties in the behavior of the system.

On the other hand, while the different kinetic mechanisms agree for the most part in both regimes, the transition zone between them is subject to strong disagreement between them. For example, two orders of magnitude separate GRI-mech 3.0 [START_REF] Gregory | GRI-Mech 3.0[END_REF] and [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF] [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF] at initial temperatures of about 1000K. Perhaps more interestingly, Figure 2.2, also compares Konnov's 2008 [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF] and Konnov's 2019 [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] mechanisms, the latter being an update of the former. These two closely related mechanisms still present an order of magnitude of difference in the transition zone. This discrepancy entails that the IDT is very sensitive to the parameters of the kinetic mechanism. Furthermore, Arrhenius coefficients that parametrize the kinetic mechanisms are uncertain because they are calibrated from sets of measurements of finite precision. Uncertainty factors are reported in some kinetic mechanisms, which means it is both possible and interesting to properly propagate kinetic uncertainties through an auto-ignition simulation, to assess the precision of numerical auto-ignition predictions.

Uncertainties in hydrogen flame ignition

While hydrogen combustion features a relatively simple kinetic mechanism, the simulation of ignition is still subject to uncertainties that can introduce a significant variance in the results. In the configuration studied here, a 1D non-premixed autoigniting flame emulating the staged combustion, only initial oxidizer temperature and kinetic constants are delved into.

Initial temperature

The scope of this study is on staged combustion, and more specifically on the second combustion stage of that process. In such circumstances, the initial temperature of the vitiated gases can be uncertain. This uncertainty can come from thermal losses in the combustor or mixing inhomogeneities, for example. Furthermore, precisely measuring the temperature inside the combustor can be challenging. Laser-based measurements have a typical uncertainty of ±3% [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF], and thermocouples are difficult to calibrate in radiative environments such as heated pipes or combustors, also leading to measurement uncertainties. For these reasons, we set a representative uncertainty of ±30 K for the vitiated gases. We also choose a uniform distribution for this temperature, as described in Equation 2.5, because it is the least-informative distribution for a given range of values. For the auto-igniting 1D flame, this corresponds to an initial oxidizer temperature ranging from 990K to 1050K. The fuel temperature, on the other hand, is considered to be precisely known. For the homogeneous reactor at Z MR , this translates into an initial temperature ranging from 959K to 1016K.

T ∼ U (T min , T max ) (2.4)
To represent this uncertainty in a normalized way, a unit uniform random variable ξ T is introduced:

ξ T = T cof -T Min T Max -T Min ∼ U (0, 1) (2.5)

Kinetic mechanism

In 2008, Konnov first published a detailed combustion mechanism for H 2 containing uncertainty data for all reactions [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF]. It models the Arrhenius pre-exponential factors as log-normal random variables, while the other kinetic parameters remain constant:

k j = A j T β j exp -E a j RT (2.6) ξ j = log(A j /A 0 j ) 1 3 log(UF j ) ∼ N (0, 1) (2.7)
With k j being the rate constant of reaction j, A j the pre-exponential factor, A 0 j its nominal value, and UF j its uncertainty factor. In this study, we choose to work on an updated version of this mechanism, published in 2019 [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF], which brings state-of-the-art nominal Arrhenius coefficients and estimations of their uncertainty factor. Konnov's 2019 mechanism contains excited species and reactions that involve them. While these are useful for predicting species observed in experimental diagnostics, this extra complexity only leads to more expensive simulations without changing the evolution of predominant species significantly [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF]. Therefore, we boiled-down the kinetic mechanism to retain only non-excited species and the relevant reactions. In Figure 2.3, we plot the probability distribution of the auto-ignition delay time of the homogeneous reactor introduced in Section 2.1.2 for several temperatures in the operating range, using the original mechanism from Konnov [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] (in red) and its reduction to non-excited species and relevant reactions (in turquoise color). This boiled-down mechanism is detailed in Annex A. In total, 32 uncertain variables are accounted for in the study: {ξ 1 , . . . , ξ 31 , ξ T }. This relatively high dimensionality means that only Monte-Carlo-like methods will be tractable to tackle this uncertain problem directly.

Numerical methods for uncertainty propagation

To propagate these uncertainties through the costly simulation, specific uncertainty methods will be needed. The brute-force approach of uncertainty propagation is called Monte-Carlo sampling and is presented in Section 1.2. Although robust and suitable for high dimensional uncertain spaces, the Monte-Carlo approach has a slow convergence rate, meaning that many samples are needed to obtain meaningful statistics of the quantity of interest. This leads to tractability issues when the computation time per sample is significant. A way to alleviate this problem is to build a surrogate of the costly simulation. A surrogate model, or response surface, emulates the costly simulation for a fraction of the cost, making an extensive sampling affordable. Two kinds of methods are explored in this study: Polynomial Chaos expansion, with its two variants, projection-PCE and regression-PCE, introduced in Section 1.3.1, and Kriging, presented in Section 1.3.2. Alas, surrogate models need a rising number of training samples when the uncertain dimension increases. This makes the construction of surrogate models impractical, or even intractable when the uncertain dimension is too high and the computation time per sample is significant. To allay this problem, we propose to reduce the uncertain dimension a priori by analyzing a cheaper simulation representative of the physical phenomena at play. The mathematical tools that we used are presented in Section 1.5.

Analysis of uncertainty in the homogeneous reactor's results and dimension reduction

This section determines and analyses the uncertainty of the simplified case's autoignition delay time. This case is the homogeneous reactor defined in Section 2.1.2, which features hydrogen auto-ignition in a vitiated mixture. A brute-force uncertainty propagation is performed in Section 2.5.1. An analysis of the resulting uncertainty is then carried out, with a search of the dominant input uncertainties (Section 2.5.2) and an interpretation of the shape of the auto-ignition delay time's probability distribution (Section 2.5.3).

Brute force uncertainty propagation

The homogeneous reactor's behavior is piloted by the physical phenomenon of autoignition. Therefore, a relevant quantity to study is the auto-ignition delay time of the reactor. However, this delay can easily span several orders of magnitude for slightly different initial conditions, as shown in Figure 2.2. This makes the analysis and uncertainty quantification difficult. Therefore, a better Quantity of Interest (QoI) is the logarithm of the auto-ignition delay: QoI = log 10 (τ ). This new QoI is contained in one order of magnitude only, and it has the added benefit of being defined over the entire real axis, whereas the auto-ignition delay must be strictly positive.

The simulations of auto-igniting homogeneous reactors are extremely cheap, which makes a brute-force Monte-Carlo approach affordable. The 32-dimension space Span(ξ 1 , . . . , ξ 31 , ξ T ), defined in Section 2.3 is therefore sampled extensively. The temporal evolution of the temperature of a handful of these samples is presented in Figure 2.4 to illustrate the wide variety of thermal runaway starting points and slopes. 2 19 QMC samples are drawn using the Sobol' sequence, and the resulting probability distribution function of the QoI is shown in Figure 2.5. Table 2.3 gathers major statistics of this distribution. First of all, it should be noted that the PDF of the auto-ignition delay time τ = 10 QoI spans two orders of magnitude. This uncertainty is far from negligible, and the uncertainty propagation framework will be helpful in this case. Second, the PDF of the QoI is bimodal. This shape highlights strong non-linearities in the underlying physical phenomena. Therefore, this non-trivial statistical behavior must be represented by its PDF because usual metrics like mean and standard deviation are ill-suited to describe it.

Furthermore, such a behavior cannot be retrieved using local (linearized) sensitivity analysis. This justifies the use of a full-blown global uncertainty propagation approach for combustion simulations.

Quantitative sensitivity analysis and uncertain dimension reduction

From these QMC samples, Sobol' sensitivity indices are obtained using the indirect method described in section 1.4.2. As a reminder, the Sobol' indices are the portion of QoI variance explained by each variable or group of variables. The PCE used to retrieve the Sobol' indices is determined using least square regression with 10-fold cross-validation using 2 19 samples. The PC basis is truncated at a maximum order of 5 and first-order interactions (p=5, q=0.5) to contain the computational cost of determining the PCE weights.

This PCE is verified visually using the summary plot of Figure 2.6. The fit is not perfect because we constrained the highest order of polynomials and limited ourselves to first-order interactions due to tractability reasons. However, the behavior of the QoI is mostly retrieved, and we will proceed with confidence in the highest Sobol' indices we can extract from it. These highest indices are plotted in decreasing order in Figure 2.7. • The initial temperature T : As already exposed in Section 2.2, the QoI is highly sensitive to the initial temperature. The large uncertainty range associated with this variable ( T ∼ U (959 K, 1016 K) ) makes it the most impacting uncertainty in this configuration. Its index stands alone in this study, dwarfing all effects from other uncertain parameters.

• The pre-exponential constant A 12 : The 12 th reaction, H + O 2 = OH + O, has been known for its very high sensitivity for a long time [START_REF] Baulch | Evaluated kinetic data for high temperature reactions Volume I, Homogeneous gas phase reactions of the H2-O2 systems[END_REF]. Because of this, it has since been extensively studied, and as a result, its uncertainty factor (UF 12 = 1.1) is the lowest of Konnov's kinetic mechanism. However, it remains the most impacting uncertain kinetic parameter as already underlined for laminar flame speed [START_REF] Wang | Combustion kinetic model uncertainty quantification, propagation and minimization[END_REF][START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF].

Interpretation of the bimodal distribution

To understand the bimodal distribution observed for this uncertain problem in Figure 2.5, we propose Figure 2.8 which represents the evolution of the auto-ignition delay as a function of the initial temperature, while also taking into account the uncertainties of the kinetic mechanism. This allows for visually isolating the effect of the initial temperature, which dominates the uncertain problem, from all other uncertain effects.

This figure showcases a "step" in the auto-ignition delay between the initial temperatures of 970 K and 1000 K. This is the transition between the high-and the low-temperature regimes (or "strong" and "weak" regimes respectively). This transition can be observed at slightly different temperatures, depending on the initial composition. For numerical studies like this one, its position also obviously depends on the kinetic mechanism in use.

Considering the probability distribution of the initial temperature -a uniform distribution between 959 and 1016 K -the samples populate both regimes, in which the QoI is much more probable than in the transition zone between the two. Hence, the two bellies of the distribution. A similar conclusion can be reached using Figure 2.3, which shares the same axes but aims at determining the probability distribution of auto-ignition delay at several deterministic initial temperatures.

Analysis of the auto-ignition delay time's uncertainty in the complex configuration

Section 2.5 proved the necessity of proper uncertainty propagation for the simulation of the auto-ignition of the mixtures typically found in the second stage of staged This study is allowed by the still manageable cost of such a simulation. Other uncertainty propagation methods are explored in the following because the Monte-Carlo approach is not affordable for typical 3D simulations of turbulent flames. To tackle this problem, the two-stage framework consists first of an uncertainty dimension reduction based on a simplified configuration. The reduction proposed in the previous section is assessed on the actual 1D problem in Section 2.6.2. This allows surrogate models to be built in a second stage using a limited number of training samples. A comparison of surrogate modeling methods is given in Section 2.6.3. Finally, error metrics for the surrogates, usable with a modest set of validation samples, are presented in Section 2.6.4.

Reference distribution from a Monte-Carlo study

The 1D auto-igniting non-premixed flame is still accessible in Monte-Carlo despite being several orders of magnitude more expensive to simulate than a homogeneous reactor. 2 13 QMC samples are drawn in the 32-dimension uncertain space as a reference for the methods investigated in the following. The resulting probability distribution of the QoI is shown in Figure 2.9. Significant statistics are presented in Table 2.4.

A bimodal distribution is found, much like in the homogeneous reactor case, although this distribution is skewed by heat and species diffusion. This confirms that complex, uncertain kinetic behaviors can be observed in more physically complete setups. In turn, this entails that proper uncertainty quantification can be needed when such situations arise.

In the case of the design of a complex burner, though, a Monte Carlo experiment is not tractable because the cost of a 3D turbulent combustion simulation typically ranges from a few hundred CPU hours (in RANS) to tens or hundreds of thousand CPU hours (in LES). Therefore, only a few observations can be obtained, on which surrogate models can be built. However, surrogate modeling is problematic in high-dimensional uncertain spaces. This problem will be addressed in the next section.

Impact of the uncertain dimension reduction

Building a surrogate in high dimensions requires many samples, preferably twice as many as uncertain dimensions, to capture at least first-order effects in all directions. In the case of simulations featuring H 2 /O 2 combustion, this means at least 64 simulations of the physical model. For complex systems that may be expensive to simulate, especially in the context of design and optimization where several iterations of the system will be compared, such a high number of simulations to build a surrogate can be unaffordable. Therefore, we propose to perform a priori uncertain dimension reduction to drastically reduce the number of samples required to build the surrogate. The most forward way to do it is to conduct a global sensitivity analysis, i.e. sorting the variables whose uncertainty lead to the greatest uncertainty in the QoI. The researcher can then truncate the uncertain space at the level of uncertainty that is chosen to be retained.

The problem, however, is that a global sensitivity analysis requires either sufficiently abundant Monte-Carlo samples or a good surrogate of the phenomenon, which we cannot afford in the general case. Then, the idea is to use a simplified setup that requires much less computational effort to simulate. It should feature approximately the same physical phenomena and be sensitive to the same uncertain variables. In these conditions, the sensitivity analysis results can be -at least qualitatively -conveyed onto the complex system. In our case, we need a simple setup that is piloted by auto-ignition in the same conditions as the 1D flame. The homogeneous reactor presented in Section 2.1.2, being chosen at the most reactive mixture of the 1D flame, matches these requirements. In this study, we choose to retain only the uncertain parameters whose influence on the QoI's uncertainty is greater than 1%, i.e. those whose Sobol' index is greater than 10 -2 . This leaves only the initial temperature T and the pre-exponential constant A 12 . A new uncertain space is then assembled, spanned by these two variables alone.

The reduced uncertain space, Span(ξ 12 , ξ T ), is then resampled in a separate Monte-Carlo simulation, and the simplified physical model is re-evaluated to check the validity of the uncertain space reduction. As can be seen on the left of Figure 2.10, the two pdfs mostly collapse, indicating that most of the original uncertainty of the QoI is retained in the reduced uncertain space. The more complex physical model is also resampled on the reduced uncertain space as a means of step-by-step validation. Of course, in the general case, neither this resampling nor the original "reference" sampling can be afforded and one will have to rely on the error metrics described in section 1.6. The results of this resampling are shown on the right of Figure 2.10. Again, the pdfs are very similar, indicating that the sensitivity analysis and the resulting uncertain dimension reduction could indeed be extrapolated to the more complex setup.

Comparison of the surrogate modeling methods

The performance of the four surrogate modeling methods presented in Section 2.4 (Regression PCE, Projection PCE, Kriging, and PC-Kriging) are now compared on this particular setup. Surrogate models will be trained on the reduced uncertain space Span(ξ 12 , ξ T ), using different sizes of training sets: The projection-PCE will use 17, 49, and 129 samples, corresponding to the number of points needed in the second, third, and fourth levels of the sparse Fejer cubature. The regression-PCE, Kriging, and PC-Kriging will be trained on a Quasi Monte-Carlo training set, using 10, 15, and 40 observations. To avoid overfitting, the PCEs are limited to a total order of 4. As shown in Figure 2.11 on the example of the regression PCE, the quality of the surrogate increases with the number of samples in the training set, which translates into a better fit between the PDF predicted by the surrogate and the reference PDF in the reduced uncertain space. This better fit between PDFs can be quantitatively measured by their Wasserstein distance. The evolution of this distance for each type of surrogate and each size of training set is plotted in Figure 2.12.

This plot highlights that projection-PCE is not performing as well as Kriging and regression-PCE in this configuration. The projection-PCE has a spectral convergence towards the "real" evolution of the QoI, which means that given enough exactly computed polynomials, the real evolution of the QoI can be retrieved with as much accuracy as desired. However, in the case of expensive simulations, the weights of higher-order polynomials can not be adequately computed because their evaluation requires a large number of samples. Therefore, with few available samples, projection-PCE is stuck with low-order effects and interactions. This is not a problem per se, as highlighted by the performance of the regression-PCE which uses the same limited polynomial basis. However, the weights computed using projection do not optimize the global fit of the surrogate to the observations but are merely projections on basis polynomials.

On the other hand, both Kriging and regression-PCE make good use of the available samples to mimic the behavior of the QoI. They have a higher convergence rate, and the fit is already more than acceptable for low training set sizes (15 samples). However, with these methods, one must take care of the pitfall of overfitting. For the PCE, a rule of thumb is to avoid using a number of polynomials higher than half of the cardinality of the training set. However, this may not be enough to avoid overfitting, and we recommend training PCE surrogates with several sizes of polynomial basis and comparing them using the available quality metrics. For Kriging, overfitting is avoided by correctly choosing the Kernel and length scales. In our experience, Matérn 5/2 kernel and MLE determination of length scales work in most cases.

It is also important to note that Kriging and regression-PCE can use the same training samples, i.e. random or quasi-random samples, while projection-PCE needs to use its own set of training samples. Kriging and regression-PCE also yield similar performances while the projection-PCE surrogate lags behind.

Projection-PCE also suffers from a very specific issue. For some sample points, it may not be possible to properly define an observed QoI. In the case of a 1D auto-igniting flame, if the kinetics are too slow, heat and intermediate species are diffused too fast in comparison, and the flame never ignites. This can typically happen at the edges of the uncertain space, which are reached by the quadrature scheme. Discontinuities of that sort can be handled with regression methods by building a surrogate conditioned by auto-ignition and ignoring ill-defined observations in the training process. However, with projection methods, this is not possible because a value for the QoI must be specified at each quadra-ture point. The chosen value impacts the shape of the response surface, and therefore the fidelity of the response surface. In this study, for the sake of comparing different surrogate modeling methods, this problem was avoided by setting a very low strain rate for the 1D flames, as specified in section 2.1.1, ensuring that all sampled flamelets actually ignite.

Finally, in this case, PC-Kriging using a 2nd order regression-PCE as trend outperforms the other surrogate model techniques by combining the smooth global trends of PCE with the possibility of local variations offered by Kriging.

Therefore, to maximize both the performance and the robustness of the surrogate modeling step, we recommend sampling the reduced uncertain space with a space-filling sampling method, like a Sobol' sequence, then building regression-PCE, Kriging, and PC-Kriging surrogates over these training points. The best surrogate model can then be selected a posteriori depending on which performs better in each specific situation.

Error measurement with a modest validation set

In the general case, when the reference distribution is not available, the alternative metrics presented in Section 1.6 can be used. Using 20 validation samples drawn with Latin Hypercube Sampling (LHS), we can compute the Normalized Root Mean Squared Error (NRMSE) and Normalized Maximum Absolute Error (NMAE) of each surrogate model. Figure 2.13 presents their evolution with the size of the training set. On the one hand, the evolution of the NRMSE is in qualitative agreement with the reference metric, making it suitable to compare the global quality of surrogates when the reference distribution can not be obtained. The NMAE, on the other hand, gives a measure of the maximum error of the surrogate. In particular, it can translate the increasing oscillations of the regression-PCE at the edges of the uncertain space when the maximum order of polynomial increases. The NRMSE and NMAE also have the advantage of being more quantitatively intuitive measures than the Wasserstein distance. Overall, their combination is suitable to evaluate the quality of a surrogate method, both in terms of global fit and potential localized error.

Takeaway points

On the physical interpretation side:

• This study has shown that typical uncertainties in the initial temperature of a mixture can significantly impact H 2 /O 2 combustion when the auto-ignition phenomenon is involved.

• Simulations of such configurations are also impacted by the uncertainties in the kinetic mechanism itself.

• These effects were shown to happen not only in the relatively simple case of homogeneous reactor auto-ignition but also in the more complex case of auto-igniting non-premixed 1D flames.

• Significant non-linearities in the transition zone between high-and low-temperature chemistry can lead to heavily skewed or even bimodal probability distributions of the auto-ignition delay.

• A two-dimensional space, Span(ξ 12 , ξ T ), was shown to explain most of the variability of both auto-igniting homogeneous reactors and non-premixed 1D flames in the investigated conditions.

On the methodological side:

• In configurations where auto-ignition is involved, a proper uncertainty propagation framework must be used to tackle these uncertainties in combustion simulations to provide researchers and designers with reliable, uncertainty-enriched simulation results

• Propagating high dimensional uncertainties in combustion simulations is challenging because the Monte-Carlo approach is typically not affordable, and surrogate modeling methods require large training samples due to the curse of dimensionality.

• We put forward a framework based on a simplified configuration that shares the same underlying physical phenomena. Extensive sampling and sensitivity analysis can be performed on this less costly simulation, allowing a priori uncertain dimension reduction for the more complex configuration.

• We recommend sampling the reduced uncertain space with a space-filling method that allows building both Kriging surrogates and regression-PCE and picking whichever surrogate method works best. The comparison between the two methods can be done using the NRMSE and NMAE error metrics over a set of validation samples.

• With these tools, we are ready to tackle uncertainty propagation in expensive simulations such as 3D Large Eddy Simulations of reacting flows.

Chapter 3

The H 2 Cabra flame

This chapter presents the H 2 Cabra flame [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] and its Large Eddy Simulation. Section 1 presents the experimental setup and the flame's characteristics of interest. Section 2 presents the numerical setup used in this thesis to simulate the Cabra flame. Section 3 presents the results of the nominal simulation. Section 4 studies the convergence of the results with regards to the mesh size and length of the averaging time period. Sections 5 and 6 explore qualitatively the sensitivity of the Cabra flame to various parameters.

The experimental setup

The Cabra flame is a lifted non-premixed flame in a vitiated co-flow. As shown in Figure 3.1, a jet of H 2 diluted in N 2 exits from the central nozzle at ambient temperature. It exits in a co-flow of hot, burned gases composed of H 2 O, O 2 , and N 2 . This hot co-flow is obtained by surrounding the base of the central nozzle with a myriad of small premixed H 2 -air burners in a lean mixture.

Choice of operating conditions

Many combustor designs use the recirculation of hot combustion products. This can be achieved for example by using a swirled injector which creates an inner and an outer recirculation zone around the cone of injection. With lower local velocities and bountiful activation energy coming from the high temperature of the burnt gases, this recirculation helps the stabilization of flames in highly turbulent flows. The modelization of such combustors is both difficult and computationally intensive due to the coupling between a complex, turbulent flow and the chemical kinetics. A corollary of this state of affairs is that simulations of such combustors are difficult to validate, because the sources of error are plenty, and modelization errors can sometimes compensate, leading to a "good-looking" solution for the wrong reasons. For this reason, Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] introduced the Cabra flame configuration, which partly decouples the complex flow configuration of industrial burners from the chemical kinetics by employing a jet in co-flow configuration. The fuel is injected in a cold, diluted central jet while the oxidizer is present in the hot, vitiated co-flow at atmospheric pressure. This allows studying the kinetics of a non-premixed flame in a vitiated environment in a rather simple flow configuration. On the one hand, from an experimental point of view, this is practical because it allows easy optical access to the flame, which in turn allows many different diagnostics to be performed. On the other hand, from the numerical point of view, this configuration also eliminates sources of error by using a simple and well-known flow configuration and reducing the sensitivity to wall boundary conditions.

Lastly, the flame is lifted in its nominal operating conditions, meaning that the flame front is not attached to the lips of the burner. It rather hovers a few jet diameters above them. This property is desirable in many practical applications because a lifted flame does not heat the injector as much as an attached flame.

Two versions of this configuration were proposed by Cabra et al.: A methane flame [START_REF] Cabra | turbulent Jet Flames Into a Vitiated Coflow[END_REF] and a dihydrogen flame [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Cabra | turbulent Jet Flames Into a Vitiated Coflow[END_REF]. In this thesis, we studied only the H 2 flame. Quantitative details of the operating conditions are given in Table 3.1.

Flame structure and behavior

Cabra et al. noted in [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] that the flame spontaneously ignites downstream soon after the fuel inlet is opened. This event happens at H/D jet = 660, which is far above the steady state flame front which hovers at about 10 diameters of the fuel inlet. This first observation indicates that the Cabra burner can -and does -auto-ignite. This initial auto-ignition is followed by "a sequence of loud auto-ignition events, each successively closer to the nozzle exit" [START_REF] Cabra | turbulent Jet Flames Into a Vitiated Coflow[END_REF], until the lifted flame stabilizes at its steady-state position.

Cabra et al. also asserted that the stabilization of the flame in its steady-state is mostly driven by auto-ignition, based on intermediate species budgets from a RANS- [START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF], another RANS-PDF simulation by Gordon et al. [START_REF] Gordon | A numerical study of auto-ignition in turbulent lifted flames issuing into a vitiated co-flow[END_REF], and a Direct Numerical Simulation by Yoo et al. [START_REF] Chun | A 3D DNS study of the stabilization of a turbulent lifted hydrogen/air jet flame in an autoignitive heated coflow[END_REF].

Finally, the structure of the steady-state flame can be described as follows:

• The fuel and oxidizer mix in the peripheral regions of the central jet.

• Small auto-igniting pockets form in places where a highly-reactive, hot, and lean mixture is found.

• These pockets swell and fuse while being convected downstream, ultimately forming a plume of burnt gases

The Quantity of Interest: the flame lift-off height

A relevant quantity of interest to describe the Cabra flame is its lift-off height. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] defines this lift-off height as the axial location where the Favre-filtered hydroxyl mass fraction Y OH first reaches 600 ppm.

Sensitivity to the co-flow temperature

The auto-igniting nature of the Cabra flame entails that the flame position and shape are affected by how reactive the initial mixture is.

The co-flow is formed by a hot, oxygenated gas (1045 K in the nominal conditions). It is the thermal energy brought by this gas that jumpstarts the combustion reaction with the cold fuel jet. The initial reactivity of an auto-igniting reactor varies exponentially with its initial temperature. Therefore, even small variations in the co-flow temperature should lead to significant modifications of the flame position and/or shape. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF]. It shows the normalized lift-off height plotted against the measured co-flow temperature. As expected, the sensitivity of the lift-off height to the co-flow temperature is significant in the vicinity of the nominal operating conditions. Also interesting is the fact that these three authors seem to be in severe disagreement unless the uncertainty on the measurement of the co-flow temperature is taken into account. This uncertainty, evaluated as ±3% in [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF], substantiates the agreement between the four experiments.

However, it also highlights the difficulty to validate a deterministic simulation, or worse, to discriminate between two modelization choices, based on these experimental results.

Numerical simulation of the cabra flame

Since the Cabra flame is mainly piloted by the local and transient phenomenon of autoignition, and since we wish to simulate it in the most precise possible way, the framework of choice is Large Eddy Simulation (LES).

LES is a framework of turbulent simulation that resolves the large scales of fluid motion, while smaller scales are filtered and their effect are modelled. Cost-and precisionwise, it falls between Direct Numerical Simulation (DNS) which resolves all scales of turbulence and Reynolds-Averaged Navier-Stokes simulations (RANS) which resolves only the mean fields. The large scales are geometry-dependent and contain most of the kinetic energy of the fluid. The smaller scales are more geometry-independent and contain less energy, which allows them to be modeled more accurately than the large scales.

To perform said LES of a reactive flow, a number of physical modelization choices need to be made, along with very practical discretization choices. To be precise, we need to choose:

Favre-filtered Navier Stokes equations

A Large Eddy Simulation solves a set of Favre-filtered Navier-Stokes equations. In the following, for any variable ϕ, Favre-filtered quantities are noted ϕ, and Reynolds-filtered quantities are noted ϕ, such that ρ ϕ = ρϕ.

The equations that will need to be solved, regardless of the modelization choices, are the following:

Mass balance ρ ∂ρ ∂t + ∇ • ρũ = 0 (3.1)
Where ρ is the gas density and u is its velocity vector.

Momentum balance ρũ ∂ρũ ∂t + ∇ • (ρũũ) = -∇P + ∇ • τ + ∇ • τ u (3.2)
Where P is the gas pressure, τ is the laminar stress tensor and τ u is the subgrid stress tensor.

The laminar and subgrid stress tensors are defined respectively as:

τ ij = µ ∂ ũi ∂x j + ∂ ũj ∂x i - 2 3 µ ∂ ũk ∂x k δ ij ; τ u ij = -ρ( u i u j -ũi ũj ) (3.3)
With µ the laminar dynamic viscosity of the fluid.

Total energy balance ẽt

∂ρ ẽt ∂t + ∇ • (ρũ ẽt ) = -∇ • (P ũ) + ∇ • (τ u) -∇ • q -∇ • q t (3.4)
e t is the gas total energy, q is the laminar heat flux and q t is the subgrid heat flux.

The viscous flux τ u can also be written τ ũ. The laminar and subgrid heat fluxes are respectively defined as:

q i = λ Cp ∂ h ∂x i ; q t i = (ρe t + P )u i -(ρ ẽt + P )ũ i (3.5)
Where λ is thermal conductivity of the gas, Cp is its laminar heat capacity at constant pressure and h is its enthalpy.

Physical choices

Turbulence model

The subgrid stress tensor τ u is modeled using the eddy-viscosity assumption:

τ u ij - 1 3 δ ij τ u ll = 2ρν t Sij - δ ij 3 Skk (3.6)
Where:

Sij = 1 2 ∂ ũi ∂x j + ∂ ũj ∂x i (3.7)
The subgrid-scale viscosity ν t is modeled in this study using the so-called Sigma model [START_REF] Nicoud | Using singular values to build a subgrid-scale model for large eddy simulations[END_REF], which is written:

ν t = (C σ ∆) 2 • σ 3 (σ 1 -σ 2 )(σ 2 -σ 3 ) σ 2 1 (3.8)
where ∆ is the filter characteristic length, C σ = 1.5 is the model constant and σ 1 ≥ σ 2 ≥ σ 3 ≥ 0 are the singular values of the resolved velocity gradient tensor. This turbulence model is the recommended model for 3D simulation in the AVBP software [START_REF] Schonfeld | Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP[END_REF]. It was developed to alleviate some drawbacks of the previous static turbulence models, notably by converging to zero subgrid-scale effects in laminar flows and having the correct asymptotic behavior near solid boundaries.

The subgrid heat flux q t is modeled using the Standard Gradient Diffusion Hypothesis (SGDH):

q t = - λ t Cp ∇ h (3.9)
Where the turbulent thermal conductivity λ t is given by:

λ t = ρν t Cp P r t (3.

10)

Where P r t is the turbulent Prandtl number whose value is set at P r t = 0.6, the default in the AVBP software.

Kinetic mechanism

The ultimate goal of this thesis is to study the uncertainties in the Cabra flame, and in particular the uncertainties introduced by the kinetic mechanism. Therefore, for the LES of the Cabra flame, we chose to use the boiled-down Konnov 2019 [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF], which was introduced in chapter 2, section 2.3.2.

Combustion and flame-turbulence interaction model

Since the flame is a diffusion flame and is auto-igniting [START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF], we leverage the Unsteady Flamelet Progress Variable (UFPV) approach with a presumed PDF closure model, proposed by Ihme and See [START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF]. The UFPV approach consists of the tabulation of nonpremixed auto-igniting flamelets. It is important to note that this approach will predict accurate source terms solely for the auto-ignition of a mixture in the flamelet regime.

The flamelet model

The flamelet model, on which the UFPV tabulation is based, assumes that the flame front is thin compared to its wrinkling scales and other flow scales. This allows to assimilate the flame front to a locally-1D laminar flame, or flamelet, whose structure depends only on time and on the coordinate normal to the flame front. This structure can be pre-computed and stored in a table, from which thermophysical properties, composition and reaction source terms can be read at each time step and at every point of a 3D simulation. In our case, a non-premixed jet flame, the flamelet of choice is a diffusion flamelet, which describes the structure of a counter-flow, non-premixed flame. The relevance of this choice is illustrated in Figure 3.3, which presents, in a simplified manner, the mixing in a shear layer behind a splitter plate, which can also be seen as a cross-section of one side of a non-premixed jet. This figure shows that the shear layer can very locally be seen as a counter-flow flame, which calls for a non-premixed flamelet modelization. Furthermore, it also illustrates that the properties of the flow, in this case the mixture fraction, can vary significantly across very short distances in a turbulent flow. These short distances may be smaller than the mesh size, which means that these variations must be modeled and taken into account when building the table. This is where the presumed PDF closure model comes into action. We also chose unsteady flamelets to correctly represent the auto-ignition of the mixture. Finally, a unity Lewis number is chosen for the simulation of the flamelets, as it is representative of turbulent combustion [START_REF] Peters | Cambridge Monographs on Mechanics[END_REF]. As a reminder, these flamelets are the ones studied in Chapter 2.

The UFPV table

The flamelets whose solutions are used to build the table should span the domain where they can auto-ignite. This domain can be determined by plotting the so-called S-curve.

The S-curve, shown in Fig. 3.4, displays the maximum temperature of the flamelet as a function of the strain rate it is subject to, all other parameters being equal. It contains two critical points:

• The first critical point corresponds to the quenching strain rate a q beyond which heat and radical species diffusion become too strong for the flamelet to sustain combustion.

• The second critical point corresponds to the ignition strain rate a i beyond which heat and radical species diffusion become too strong for the cold mixture to auto-ignite. This second critical point is the upper bound of our table. The table should include flamelets with a strain rate going from zero up to the second critical point of the S-curve. Alas, the transient simulation of the auto-igniting flamelet is not stable for too small strain rates. Thus, an empirical value of a = 5 s -1 was chosen as a lower bound of the table, because it was the smallest strain rate to guarantee numerical stability.

A set of 10 flamelets are computed, each one with a different strain rate, ranging from this empirical lower bound up to the strain rate of the second critical point of the S-curve. The evolution of the gas temperature in such a flamelet is plotted in Figure 3.5.

From these flamelets, we could pre-compute the auto-ignition of non-premixed laminar flames, and tabulate, for each quantity ϕ, the function F ϕ :

ϕ = F ϕ (z, Λ, χ z,st ) (3.11)
Where z is the mixture fraction, χ z,st is the stoichiometric scalar dissipation rate and Λ is a mixture fraction-independent reaction progress parameter [START_REF] Ihme | Prediction of local extinction and reignition effects in non-premixed turbulent combustion using a flamelet/progress variable approach[END_REF]. In this work, Λ is defined as:

Λ = Y C (z Λ ) , z Λ = 0.1 (3.12)
Where Y C is the progress variable defined in Equation 3.14. The choice of z Λ is arbitrary.

In this work, we chose to place it near the most reactive mixture z mr = 0.04. This definition of Λ allows a unique mapping of instantaneous flamelet solutions in the (t, a)space (time and strain rate) to the (Λ, χ z,st )-space (progress parameter and stoichiometric scalar dissipation rate). The progress variable is typically defined as a linear combination of the mass fractions Y k of the reactive species:

Y c = N species k=1 α k Y k (3.13)
The "correct" way to define the weights α k is still an open problem for complex fuels and reactions [START_REF] Niu | An optimization-based approach to detailed chemistry tabulation: Automated progress variable definition[END_REF], but it is much more manageable for simpler reactions such as H 2 -O 2 combustion. The progress variable, in our case, must be defined such that:

1. It is an increasing function of time during ignition 2. Its derivative is non-zero at t = 0. This condition ensures that auto-ignition can occur in the simulation To enforce the first condition, it is common to give a strictly positive weight to at least one product of the reaction: in this case, H 2 O. However, this is not enough to enforce the second condition, because the final products of the reaction are not immediately produced. So-called intermediate species are first produced, and are in turn consumed to produce the final products. Therefore, an additional negative weight on at least one reactant ensures a strict, non-zero increase of the progress variable from the very start of the reaction. In our case, we chose to define Y c as:

Y c = Y H 2 O -0.3778Y O 2 (3.14)
The weight of O 2 could have any negative value, but we chose this one to set Y c = 0 for Z = 0 and Z = 1. It also leads, that, for all Z values, the value Y C,u of the progress varibale in the unburnt conditions is null: ∀Z, Y C,u (Z) = 0. This is a purely esthetical consideration and this arbitrary choice has no impact on the transport of the progress variable or the construction of the table.

For turbulent flames in LES, a relation equivalent to Equation 3.11 must be defined for filtered quantities. A presumed joint PDF P for the mixture fraction, reaction progress parameter, and stoichiometric scalar dissipation rate is defined, so that:

φ = F ϕ (z * , Λ * , χ * z,st ) P (z * , Λ * , χ * z,st )dz * dΛ * dχ * z,st (3.15)
Λ is statistically independent of z and χ z,st by definition, and we assume that z and χ z,st are also independent. This gives:

P (z * , Λ * , χ * z,st ) = P (z * ) P (Λ * ) P (χ * z,st ) (3.16)
Furthermore, we assumed a beta-PDF for z [START_REF] Jimnez | A priori testing of subgrid models for chemically reacting non-premixed turbulent shear flows[END_REF][START_REF] Wall | An evaluation of the assumed beta probability density function subgrid-scale model for large eddy simulation of nonpremixed, turbulent combustion with heat release[END_REF], a Dirac distribution for χ z,st [START_REF] Pitsch | An unsteady/flamelet progress variable method for LES of nonpremixed turbulent combustion[END_REF], and we chose a Dirac distribution for Λ. This last choice is a significant simplification compared to Ihme and See [START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF], who chose a statistically most likely distribution for the progress parameter. This greatly simplified the implementation of the model in the AVBP solver, reduced the computer memory taken by the table, and reduced the computational time required to run a simulation. Most importantly, as shown in Section 3.3.2, this choice did not prevent the simulation from accurately reproducing the experimental data. With these assumptions, we get:

P (z * , Λ * , χ * z,st ) = β(z * , z, z v )δ(Λ * -Λ)δ(χ * z,st -χ z,st ) (3.17)
Joining Equations 3.15 and 3.17 gives:

φ = G ϕ (z, z v , Λ, χ z,st ) (3.18)
However, the definition of Λ makes it difficult to transport. The progress variable Y C is easier to transport, and its normalization C is used as a coordinate of the table, using the unique transformation C such that:

C = C ( Λ, z, z v , χ z,st ) (3.19)
With C defined as:

C = Y C -Y C,u Y C,eq -Y C,u (3.20)
Where Y C,eq and Y C,u are respectively the values taken by the progress variable at equilibrium and in unburnt conditions. We can also replace the subgrid variance of mixture fraction z v with the mixture fraction segregation S z , defined in Eq. 3.25. Finally, this gives, for every quantity ϕ the relation:

φ = T ϕ (z, S z , C, χ z,st ) (3.21)
Where the function T ϕ is precomputed and tabulated on the mesh specified in Table 3.2. Special attention is given to the discretization in C to ensure a correct prediction of the source term at the very start of the auto-ignition. The C-mesh is refined towards both the lower and upper bounds. For C near 0, the mesh needs to be very fine, because the reaction is exponential in its first instants. Therefore, any error committed early on would have important repercussions on the predicted auto-ignition delay-time. Conversely, for C near 1, we must correctly refine the reaction termination, especially in the richer mixtures, where every quantity still varies while they are already nearly settled at z Λ . Figure 3.6 provides a verification that the C-mesh is adequate, by comparing the autoignition of a homogeneous reactor simulated with both detailed chemistry and tabulated chemistry using the same C-mesh as the UFPV table. A selection of slices of this table for the hydroxyl mass fraction Y OH are plotted in Figure 3.7 as a function of the Favre-filtered mixture fraction z and of the Favre-filtered normalized progress variable C, for a few strain rate and mixture fraction segregation values. Note the significant shape change and overall amplitude decrease with increasing mixture fraction segregation. This species, which is a marker of the reaction, has a very similar behavior as the source term of progress variable. Slices of this source term, and other selected species mass fractions and thermophysical properties can be found in Annex B.

At each time step and at every point, a set of coordinates must be computed

As shown in Equation 3.21, the UFPV table used in this work involves four coordinates:

• The Favre-filtered mixture fraction balance z

• The mixture fraction segregation S z

• The Favre-filtered scalar dissipation rate at stoichiometry χ z,st

• The normalized Favre-filtered progress variable C z, S z , and C are found by solving transport equations, while χ z,st is computed using an algebraic equation. These equations are detailed in the following:

∂ρz ∂t + ∇ • (ρũz) = ∇ • ρD∇z -∇ • τ z (3.22)
Where D is the mass diffusivity, determined using the unity Lewis number assumption: Le = Sc/P r = 1. The subgrid flux of mixture fraction τ z is defined as:

τ z = ρuz -ρũz (3.23)
This subgrid flux is modeled using a turbulent diffusivity D t = Sc t ν t , where Sc t is the turbulent Schmidt number. Although the determination of the optimal value for the turbulent Schmidt number remains an open question [START_REF] Longo | CFD dispersion study based on a variable Schmidt formulation for flows around different configurations of ground-mounted buildings[END_REF], we used the default value proposed by the AVBP software: Sc t = 0.6, which is well within the range of acceptable values found in the literature. Finally, τ z is modeled using the SGDH:

τ z = -ρD t ∇z (3.24) Mixture fraction segregation S z S z = z v z(1 -z) (3.25)
Where z v = z2 -z2 is the subgrid variance of the mixture fraction and respects the following transport equation:

∂ρz v ∂t + ∇ • (ρũz v ) = ∇ • (ρD∇z v ) -∇ • τ zv -2τ z • ∇z -s χz (3.26)
The subgrid flux of mixture fraction variance τ zv , is defined as:

τ zvi = ρ u i z 2 -ũi z 2 -2zτ z i (3.27)
And it is modeled using the SGDH:

τ zv = -ρD t ∇z v (3.28)
Finally, the subgrid mixture fraction scalar dissipation rate s χz is modeled using a linear relaxation assumption:

s χz = C D ρ z v ∆ 2 /ν T = C χz C ε C u ρ z v ∆ 2 /ν T (3.29)
Caution must be exercised when dealing with this equation, for the definition of C D is not standard in the literature. Other works, such as [START_REF] Domingo | Large-eddy simulation of a lifted methane jet flame in a vitiated coflow[END_REF] introduce a factor 2 before s χz in Equation 3.26, and a factor 0.5 in the definition of C D . For the sake of simplicity, we kept Ihme's definition [START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF] which removes these unnecessary factors. With this definition, we used a nominal value of C D = 20.0 in the following. The value chosen for this coefficient is discussed in Section 3.6.2.

Favre-filtered stoichiometric scalar dissipation rate χ z,st

The stoichiometric scalar dissipation rate χ z,st is computed from the following algebraic equation:

χ z,st = χ z F (z st ) F (z) (3.30)
With the local scalar dissipation rate χ z defined as:

χ z = 2D |∇z| 2 (3.31)
And:

F (z) = 1 π exp -2 erfc -1 (2z) 2 (3.32)
This function can also be evaluated with the following equivalent expression, as seen in [117]:

F (z) = 1 π exp -2 erf -1 (2z -1) 2 (3.33)
The stoichiometric scalar dissipation rate is equivalent to the strain rate a, which follows the relation [START_REF] Peters | Cambridge Monographs on Mechanics[END_REF]:

χ z,st = aF (z st ) (3.34)
In LES, χ z is not accessible because it is not a filtered quantity. Instead, we can compute the Favre-filtered scalar dissipation rate χ z with:

χ z = 2D |∇z| 2 + s χz ρ (3.35)
The Favre-filtered stoichiometric scalar dissipation rate is then computed with:

χ z,st = χ z F (z, S z ) (3.36)
Where F is defined as:

F (z, S z ) = 1 0 F (z * ) F (z st ) P (z * )dz * (3.37) Favre-filtered normalized progress variable C C = Ỹc -Y C,u Y C,eq -Y C,u (3.38)
Where Ỹc is the Favre-filtered progress variable and respects the following balance equation:

∂ρ Ỹc ∂t + ∇ • (ρũ Ỹc ) = ∇ • (ρD∇ Ỹc ) -∇ • τ Yc + ρ ωYc (3.39)
Where ωYc is the production term of the progress variable, and is read from the table. The subgrid flux of progress variable, τ Yc , is defined as:

τ Yci = ρ u i Y c -ũi Y c (3.40)
It is then modeled with the SGDH:

τ Yc = -ρD t ∇ Y c (3.41)
Finally, Y C,eq and Y C,u are the values taken by the progress variable respectively at equilibrium and in unburnt conditions. They are functions of z, S z , and χ z,st only. Hence, they can be tabulated along with other values of interest and retrieved for the computation of C.

Numerical choices

Computational domain and spatial mesh

The computational domain is a cylinder in which the fuel inlet tube protrudes. It is presented in Figure 3.8. It is large enough that the jet is unaffected by the domain boundaries in the region of interest. Figure 3.8 (a) presents an overview of the domain and the gas flows involved.

The mesh is composed of tetrahedral elements defining 5.5 million nodes. An axial slice of this mesh is presented in Figure 3.8 (b), with a zoom on the fuel inlet zone in Figure 3.8 (c). The domain is 100 diameters long and 65 diameters wide. The cone, where the jet develops and the flame is located, is finely discretized.

Stability issues were encountered when trying to simulate the inlet tube itself, as acoustic waves were growing uncontrollably in this tube. Therefore, it was decided to remove the tube from the computational domain. Instead, the fuel inlet boundary is placed on a surface leading directly to the main volume of the computational domain. It is highlighted by the red disk on Figure 3.8 (d). The lips of the inlet tube are highlighted in blue in this same figure.

Boundary conditions

Proper boundary conditions must be chosen for the different surfaces enclosing the computational domain, which are:

• Fuel inlet • Co-flow inlet • Inlet tube walls • Side walls • Outlet
For all fluid boundaries, i.e fuel inlet, co-flow inlet, and outlet, Navier-Stokes Characteristics Boundary Conditions (NSCBC) [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] are used. These conditions do not enforce any one quantity at the boundary, but rather control the waves crossing the boundaries, such that:

• The mean quantities are kept at their target values

• Acoustic waves can be generated or let through without reflection

• When appropriate, turbulence can be injected.

Co-flow

Co-flow

Fuel inlet

The fuel inlet is probably the most critical boundary condition, because it conditions the shape of the jet, especially regarding the length of the potential core. A non-reflective (NSCBC) inlet condition is used, with turbulence injection. The injection turbulence method applied is described in [START_REF] Smirnov | Random Flow Generation Technique for Large Eddy Simulations and Particle-Dynamics Modeling[END_REF]. The quantities controlled are the velocity, temperature, and composition.

To emulate the effect of the inlet tube on the velocity and turbulence profile of the fuel inlet, a RANS simulation was performed on a tube of the same diameter, with a length L = 200 D to ensure the convergence of the various profiles of the flow. It was run using the Ansys Fluent software, with a k -ω turbulence model and resolved flow at the wall (y + = 1).

The resulting profiles for the axial velocity U and root mean square velocity U rms are presented in Figures 3.9 However, using the U rms profile resulted in a much too short potential core. This may be explained by the fact that the profile is under-resolved in our LES, in particular near the wall. The turbulent intensity is thus interpolated only from the parts of the profile where it is high, ignoring the area close to the wall where the turbulent intensity plummets. This may lead to a premature destabilization of the jet and a shortening of the potential core. This is an important problem because it perturbs the mixing process, which will in turn modify the flame lift-off height. Because the most important feature to reproduce is the potential core length, and because the intensity profile only affects this critical length, an arbitrary profile can be used for U rms as long as the potential core length fits the experiment. Therefore, an empirical, uniform value of U rms = 0.5 m/s was used. This arbitrary choice of a uniform U rms profile does not affect the flow beyond the potential core, because the flow becomes independent of the turbulence injection conditions once it reaches the developed jet zone.

Co-flow inlet

For the co-flow inlet, a non-reflective condition was again used. This time, however, due to the low mean velocity (3.5 m/s) compared to the jet velocity (107 m/s), a flat velocity profile without turbulence injection was used. Furthermore, we did not simulate all the small laminar flames present in the experiment setup, but rather injected directly the equilibrium (post-combustion) mixture, as described in Section 3.1.

Inlet tube walls

The inlet tube walls are defined as adiabatic walls with a wall law. It sets the normal heat flux and velocity to zero and computes the wall shear stress according to the logarithmic law of the wall.

Side walls

The side walls of the domain are defined as adiabatic slip walls. They are placed far enough from the jet that they do not influence the central jet.

Outlet

The co-flow outlet is a partially reflective (relaxed) outlet condition with controlled pressure, which acts as a high-frequency filter for the acoustic waves crossing the boundary.

Numerical methods

The numerical scheme applied is the Two-step Taylor-Galerkin C (TTGC) [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF], which is third-order in space and time. It emphasizes low-dissipation and dispersion. Although about twice more expensive as a Lax-Wendroff scheme [START_REF] Ni | A Multiple-Grid Scheme for Solving the Euler Equations[END_REF], it was chosen to minimize simulation error and to allow us to focus on other sources of uncertainties later on. We used the implementation of this numerical scheme provided by the AVBP solver [START_REF] Schonfeld | Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP[END_REF]. We implemented in this solver the UFPV model described in Section 3.2.2.3, using the TTC formalism for compressible equations [START_REF] Auzillon | A Filtered Tabulated Chemistry Model for Large Eddy Simulation of Reactive Flows[END_REF].

Results of the nominal simulation

Qualitative behavior of the simulated flame

The flame does ignite spontaneously, as described by Cabra et al. [START_REF] Cabra | turbulent Jet Flames Into a Vitiated Coflow[END_REF]. Furthermore, once the steady-state is reached, the correct behavior, explained in section 3.1.2 prevails: mixing happens in the periphery of the jet. Then some pockets of mixture auto-ignite, swell, fuse, and ultimately form a continuous flame front. Nota bene: reaching the steadystate does not mean that quantities become constant at every point of the flame. Rather, it means that, on average, the flame occupied the same position. A y-normal cut of an instantaneous solution is presented in Figure 3.11. It presents, on the upper part, the field of Favre-filtered mixture fraction z, and on the lower part, the field of Favre-filtered normalized progress variable C.

Comparison of the mean quantities with the experimental results

In this section, experimental measurements from Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] are compared to the simulated average fields of the Cabra flame. Measurement uncertainties in this setup, using the laser Raman-Rayleigh-LIF method are estimated in [START_REF] Barlow | Scalar profiles and NO formation in laminar opposed-flow partially premixed methane/air flames[END_REF].

For brevity's sake, not every available profile is shown and commented on. For a more exhaustive display of the comparison between experimental and simulation profiles, please refer to Annex C. For clarity's sake, in this section, unless specified otherwise, the averaged Favre-filtered scalar { φ} will simply be denoted ϕ.

A synthetic quantity of the mixing of the jet is the mixture fraction of the gas z. The axial profile of this scalar is plotted in Figure 3.12, and several radial profiles are presented in Figure 3 LES Exp Figure 3.12: Axial profile of the mixture fraction Z in the LES, compared to the measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] The axial profile of z shows the correct length of the potential core, as was expected from the work on the turbulent injection. However the initial slope of decrease (from x/D ≈ 6 onwards is too steep by about 30%. While this may be concerning, this issue has already been seen in works that computed similar jet flames in LES with AVBP, for example [START_REF] Schulz | A criterion to distinguish autoignition and propagation applied to a lifted methane-air jet flame[END_REF] and [117] for the methane Cabra flame. However, radial profiles of z, plotted in Figure 3.13, show a satisfactory agreement between experimental and simulated data. In particular, the evolution of the width of the jet is very well reproduced.

On the combustion side, the radial profiles of the reaction marker OH are plotted in Figure 3.14. A good agreement is found in the vicinity of the flame foot, i.e. for heights around x = 10 D. However, the mass fraction of OH is over-predicted higher up in the flame. This is visible in the profile taken at x = 26 D, and is best illustrated by Figure 3.15 which compares the entire fields of Y OH in the simulation and experiment.

Finally, as shown in Annex C, Temperature profiles and major species' profiles show a correct flame temperature and consumption rate of reactants, meaning that most of the reaction has already occurred when the discrepancies in hydroxyl mass fraction appear. Overall, the agreement between experimental and simulation data is very satisfactory.

Finally, the simulation of the Cabra flame costs about 100k CPU hours on the Irene KNL supercomputer with Intel Xeon Phi 7250 1.4 GHz processors.

Convergence study

Convergence is paramount to putting any trust in the simulation's result. We can distinguish two kinds of convergence in this case: 2. Mesh convergence: The mesh must be fine enough so that its influence on the simulation is negligible. The repartition between resolved and subgrid quantities will of course evolve when the mesh is refined, but the total quantities must remain stable.

Temporal mean convergence

The quantities of interest are averaged over 7.5 ms of simulated time in the nominal simulation. This corresponds to roughly 7 flow-through times in the region covered by experimental measurements, i.e from the fuel inlet to 26 diameters in the axial direction.

The convergence of the averaged quantities can be assessed by comparing profiles of the simulated flame, averaged respectively on the nominal simulation time and double this time. The simulation will be considered as converged if no discrepancies are spotted. Figures 3.16 and 3.17 display such plots for Y OH and z respectively. Other comparisons of profiles of species and quantities of interest are plotted in Annex D at several axial locations.

Additionally, since we are interested in the lift-off height of the Cabra flame, it is also interesting to quantify the convergence of this global property. To that end, a bootstrap estimate of the lift-off height is performed.

The lift-off height H is extracted from the averaged solution of the flame via the function H , as follows:

H = H (E[ψ]) (3.42)
Where ψ represents the fields of the simulated flame. The ergodicity of the flow is leveraged to replace the statistical expectation of the fields by its temporal average. However, this temporal average can not be perfectly converged, leading in turn to an incomplete convergence of the predicted lift-off height. To characterize the convergence error on the lift-off height, a bootstrap procedure is proposed.

The successive instantaneous solutions of the flame (i.e. snapshots of the flame at each time step of the simulation) are dependent, making the straightforward bootstrap procedure inapplicable. Instead, a stationary block bootstrap scheme [START_REF] Politis | The stationary bootstrap[END_REF] is implemented: instantaneous solutions of the flame are averaged in N non-overlapping blocks of data.

Each of these blocks produces a partially converged Cabra flame ψ i . This is illustrated in figure 3.18. The resulting partially converged flames are much less dependent than the original instantaneous solutions. Therefore, they are assumed independent in the following.

N bootstrap resamples are then drawn with replacement from this set of partially converged flames, and averaged to form a "simulated bootstrap sample" ψk . Such a simulated bootstrap sample may be, for example (with N = 10):

ψk = 1 N i∈I k ψ i , I k = {1, 1, 3, 5, 6, 6, 6, 6, 8, 9} (3.43)
The lift-off height Ĥk = H ( ψk ) is computed for this simulated bootstrap sample. Finally, this "draw, average, and post-process" procedure is repeated a large number of times B, giving birth to the empirical bootstrap distribution of the simulated Cabra flame's lift-off height Ĥ. We used N = 20 blocks, each averaging 500 instantaneous solutions spanning 0.4 ms of simulated time, and a resampling size of B = 1000. This resulted in the bootstrap mean and standard deviation being respectively:

Ĥ = 11.75D ; σ 2 ( Ĥ) = 0.074D ≈ 0.6% of Ĥ (3.44)
Therefore, with the time accumulation used for the nominal simulation, we are confident that the lift-off estimate will not be off by more than a percent, which is negligible before the variations we expect to witness in the following.

Finally, we can conclude that the average of the simulated flame is converged enough in our case.

Mesh convergence

The mesh convergence can be performed in two ways:

1. Verifying that most of the biggest turbulence scales are well resolved, meaning that the variance of any quantity of interest is mostly resolved rather than mostly modeled.

2. Verifying that the averaged solution of the flame does not change when the mesh is refined.

The first point can be addressed by comparing the portion of resolved (RV) and modeled (MV) mixture fraction variance. The total variance of mixture fraction Z var is written: To address the second point, a new mesh was generated on the same computational domain, with a uniform factor of refinement over the mesh used for the nominal simulation. This new mesh featured 7.5 million nodes, i.e. about 1.3 times more nodes than the original mesh. Figure 3.21 features a selection of radial profiles comparing the simulated hydroxyl mass fractions on the original and finer mesh. More comparisons of profiles can be found in Annex D. The profiles overlap for the most part, which indicates that the simulation is converged mesh-wise.

Z var = {z v } Modeled variance + {z 2 } -{z} 2
Finally, the lift-off height computed on the finer mesh is only 0.7% higher than the original, which is considered converged in our application. To verify whether the LES of the Cabra flame reproduces this pattern, a set of four additional simulations were performed at varying temperatures, spanning the range of uncertainty reported by Cabra. One new UFPV table was generated for each of these conditions. Figure 3.22 shows the mean hydroxyl mass fraction Y OH for each of these simulations. This figure does highlight a noticeable sensitivity of the flame to the co-flow temperature.

Figure 3.23 reprises the experimental data from Figure 3.2 to compare with the liftoff height predicted by the simulation. The nominal simulation, which features a coflow temperature of 1045 K, is in good agreement with Cabra's experiment. The trend given by [START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF] and [START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF] is also correctly reproduced. The exponential increase of the liftoff height with decreasing co-flow temperatures was unfortunately not observed in this study. However, we believe that simulations with co-flow temperatures of about 1000 K would show this increase, because the exponential increase of the auto-ignition delay-time with decreasing oxidizer temperature was observed with this kinetic mechanism during the uncertainty propagation studies presented in chapter 2, both for homogeneous autoigniting reactors and non-premixed auto-igniting flamelets.

For reference, we present in Figure 3.24 a comparison of this LES with other numerical works from the literature. Again, our simulation fares well against other state-of-theart simulations. Nota Bene: In this case, the results' discrepancies are not associated [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF].

with co-flow temperature uncertainties but with different modeling choices, and modeling uncertainties (chemical kinetics, turbulence, turbulent combustion).

Finally, it must be noted that the measurement uncertainty on the co-flow temperature forbids any categorical statement about whether our simulation is better or worse than any other simulation that falls within the uncertainty range reported by Cabra. 3.6 Sensitivity to the numerical setup

Sensitivity to the chemical kinetics

In chapter 2, we showed that the auto-ignition delay times of homogeneous reactors and non-premixed flamelets were very sensitive to the kinetic mechanism employed to model the chemical reactions. Different kinetic mechanisms lead to significantly different results. Kinetic mechanisms which reported uncertainties, such as [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF], allowed the determination of wide probability distributions of the auto-ignition delay. Consequently, we can expect a noticeable sensitivity of the lift-off height -and therefore of the flame shape -to the chemical kinetics.

Sensitivity to the mixture fraction variance dissipation coefficient

The flame-turbulence interaction model used in this LES contains the coefficient C D , also written [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF]:

C D = C χz C ε C u (3.46)
C χz is the mechanical-to-scalar time scale ratio:

C χz = τ sgs u τ sgs z = {k sgs } {ε sgs } { χ z sgs } { z ′′2 } (3.47)
And C ε and C u are defined as:

C ε = {ε sgs } ∆ {k sgs } 3/2 , C u = {ν t } ∆{k sgs } 1/2 (3.48)
Where ∆ is the LES filter size, k sgs is the subgrid-scale turbulent kinetic energy and ε sgs is the subgrid-scale turbulent energy dissipation rate.

Although there is a consensus on the proper value to choose for other model coefficients like the Sigma model constant C σ , this mixture fraction variance dissipation coefficient takes significantly different values in previous works. Domingo et al. [START_REF] Domingo | Large-eddy simulation of a lifted methane jet flame in a vitiated coflow[END_REF] or Vicquelin [117] set its value at C D = 2.0, for example, while Ihme and See [START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF] use C D = 4.0. Ihme's thesis [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF] (pp. [START_REF] Gordon | A numerical study of auto-ignition in turbulent lifted flames issuing into a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Harris | Array programming with NumPy[END_REF][START_REF] Hastings | Monte Carlo sampling methods using Markov chains and their applications[END_REF] gives more details on how to best determine this coefficient for a given simulation. Ihme shows that the value of C D is dependent on the Taylor-scale Reynolds number, on the filter width, and on the Kolmogorov length scale. This makes the ideal C D value dependent on the flow configuration and the local mesh size. Therefore, any attempt to define a global C D value for a given simulation must be a compromise designed to minimize the discrepancy between the experiment and the simulation.

• A value of C χz = 2.0 is proposed by Ihme for typical LES filter size, which is in agreement with other propositions in the literature ( [START_REF] Peters | Cambridge Monographs on Mechanics[END_REF][START_REF] Fox | Computational Models for Turbulent Reacting Flows[END_REF]).

• Ihme also shows that the ratio This sensitivity was expected because the equilibrium value of z v is proportional to C -1 D . This is deduced from equations 3.26 and 3.29. Therefore, as a general trend, the higher C D , the lower the subgrid mixture fraction variance. Additionally, less subgrid variance entails more reaction source term. This can be seen, indirectly, in Figure 3.7, or directly in Annex B. Finally, a higher reaction source term leads to a shorter lift-off height. This is illustrated in Figure 3

C ε /C u is in the range π 2 ≤ C ε /C u ≤

Takeaway points

• A Large Eddy Simulation of the Cabra flame was performed.

• The UFPV framework was implemented in the AVBP solver to tackle the modelization of the turbulent combustion in the non-premixed auto-igniting flamelet regime.

• Although some discrepancies between the experimental measurements and the simulation were found, the LES of the Cabra flame was deemed satisfactory overall, in particular in its reproduction of the lift-off height.

• In light of the uncertainties on the co-flow temperature and the chemical kinetics, and the substantially different values of mixture fraction variance dissipation coefficient found in the literature, the deterministic approach to simulation is not sufficient to ascertain the agreement between the experiment and simulation of the Cabra flame.

• Additional work is therefore required, and Chapter 4 tackles this issue with the framework of uncertainty propagation.

Chapter 4

Quantifying the uncertainty of the simulated Cabra flame's lift-off height

This chapter investigates the uncertainties introduced in the LES of the Cabra flame and their influence on its lift-off height. It does not aim to be an exhaustive uncertainty propagation, meaning that we do not claim to have taken into account every input uncertainty. The first part of this chapter, (Sections 1-4) details the propagation of uncertainties directly linked to chemical kinetics. The second part introduces a new uncertain parameter from the turbulence-combustion model and performs again the uncertainty propagation with the addition of this new coefficient.

Section 1 introduces the kinetic and temperature propagated in the first study. Section 2 reduces the input uncertain dimension to a manageable number. Section 3 details the computer experiment and section 4 presents the results of the uncertainty propagation. Finally, section 5 repeats the same procedure with the addition of the uncertain dissipation coefficient of the mixture fraction variance.

The input uncertainties 4.1.1 Kinetic mechanism

The combustion model used in the LES of the Cabra flame relies on a kinetic mechanism, that is, a set of elementary chemical reactions and laws to determine their associated reaction rates. Chapter 2 outlined the uncertainties associated with these reaction rates and the effects that they have on the auto-ignition delay time of a reactive mixture. As explained in Chapter 3, the driving process of flame stabilization in the H 2 Cabra setup is auto-ignition. Therefore, uncertainties in the kinetic mechanism must have an impact on the flame lift-off height. Similarly to the study presented in Chapter 2, we use the boiled-down Konnov (2019) mechanism [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] detailed in Annex A. This mechanism models the reaction rates of each reaction j with an Arrhenius law:

k j = A j T β j exp -Ea j RT (4.1)
Where T is the gas temperature, Ea j is the activation energy of reaction j, β j is a corrective dimensionless number and R is the universal gas constant. A j is the preexponential factor of the Arrhenius law and is modeled as a log-normal random variable:

ξ j = log(A j /A 0 j ) 1 3 log(U F j ) ∼ N (0, 1) (4.2)
Where A 0 j is the nominal value of A j and U F j is its "uncertainty factor". ξ j is a unit normal random variable.

The boiled-down Konnov mechanism has 31 reactions and therefore 31 uncertain preexponential factors are considered in the following.

The results of the following study are of course sensitive to the choice of kinetic mechanism and considered uncertainties. In particular, the fact that only the preexponential constants are assumed uncertain is a modelization choice from the authors of the mechanisms and their sources. We choose to use Konnov's (2019) mechanism [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] because it is state-of-the-art at the time of this thesis, but the methodology would be compatible with other mechanisms should better inputs become available in the future.

Co-flow temperature

As stated in the original publication by Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF], the co-flow temperature T cof can be measured with an uncertainty of about ±3%, which corresponds to ±30K. This uncertainty is dependent on the measurement method used (the laser Raman-Rayleigh-LIF method in the case of Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] and Wu et al. [START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF]), but Gordon et al. [START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF] report the same uncertainty using a K-type thermocouple. Figure 3.22, in Chapter 3, shows that such uncertainty on the co-flow temperature introduces significant variability in the lift-off height. Therefore, it must be taken into account to compare simulations with the experimental results. As far as we are aware, no more information is available about this uncertainty. The leastinformative probability distribution -the uniform distribution between T Min = 1015 K and T Max = 1075 K -is consequently chosen for the co-flow temperature:

T cof ∼ U (T Min , T Max ) (4.3)
To represent this uncertainty in a normalized way, a unit uniform random variable ξ T is introduced:

ξ T = T cof -T Min T Max -T Min ∼ U (0, 1) (4.4)
In total, 32 uncertain variables are accounted for in the study: {ξ 1 , . . . , ξ 31 , ξ T }.

Reduction of the input uncertain dimension

Large Eddy simulations of the Cabra flame are expensive (about 100k CPU hours on the Irene KNL supercomputer) which means that a brute-force, Monte-Carlo approach is unaffordable in this case. A surrogate-based approach will therefore be preferred. However, Chapter 1 pointed out that building a surrogate model in high dimensions requires a large number of observations, which, again, are not affordable.

Therefore, an a priori uncertain dimension reduction is mandatory. In this study, we chose to use the "dominant variables" approach (cf Section 1.5.1.2).To this end, a global sensitivity analysis has to be performed to sort out important variables from less important ones. That is, to highlight the variables whose uncertainties lead to the greatest uncertainty in the quantity of interest. Again, global sensitivity analysis cannot be implemented directly on the 3D flame lift-off height for tractability reasons. Instead, a representative physical problem of much reduced computational cost must be used.

A representative reduced physical problem

Auto-ignition being the dominant stabilization mechanism for the Cabra flame in most of the operating conditions [START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF], it can be assumed that the lift-off height of the flame is first-order dependent on the auto-ignition delay (IDT) of the most reactive mixture in the jet [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF]. In Figure 4.1, we compare the normalized auto-ignition delay time of the most reactive mixture to the normalized lift-off height for different co-flow temperatures, with different nominal kinetic mechanisms. The auto-ignition is detected in the same way as lift-off is detected in the 3D flame, i.e. when the hydroxyl mass fraction first reaches Y OH = 600 ppm. The nominal initial conditions of the reactor are given in Table 4.1. Figure 4.1 shows a very similar response of the auto-ignition delay of 0D reactors and the Cabra lift-off height in actual LES to both temperature and kinetic mechanism changes. 0D auto-ignition simulations are, therefore, a representative physical problem for most of the Cabra flame's operating conditions. Furthermore, each 0D simulation costs about 0.1 CPU-second, making them suitable for a Monte-Carlo-based sensitivity analysis. A Monte-Carlo study is conducted to account for the uncertainties of the co-flow temperature and the kinetic mechanism on the 0D auto-ignition delay. 2 17 ≈ 100k samples Table 4.1: Initial conditions for the 0D reactor at the most reactive mixture z mr = 0.04 [START_REF] Mastorakos | Ignition of turbulent non-premixed flames[END_REF] in the nominal case (corresponding to T cof = 1045 K)

P (Pa) 101325 Y O 2 0.1669 Y H 2 O 0.0626 T (K) 1011.3 Y H 2 0.0009 Y N 2
0.7696 T : Temperature; P : Pressure; Y : Mass fraction are drawn from the 32D Sobol' sequence and a homogeneous reactor autoignition is simulated for each sampled condition. Figure 4.2 shows the temporal evolution of the mixture temperature for the first 20 realizations of the homogeneous reactor. Table 4.2 records some statistics of this study. The uncertainty on the QoI is considerable, as shown by the 95% confidence interval, which spans an entire order of magnitude. Moreover, the important discrepancy between the mean of the distribution and the nominal case suggests substantial non-linearities in the underlying physical phenomenon.

Global sensitivity analysis

The global sensitivity analysis method we chose is Sobol' sensitivity analysis [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF], approximated using the PCE approach proposed by Schöbi and Sudret [START_REF] Schöbi | PCE-based imprecise Sobol'indices[END_REF] (cf Section 1.4.2). The target quantity of the surrogate is the logarithm of the IDT log 10 (τ ), to enforce the positivity of the emulated auto-ignition delay time.

The PCE is determined using least square regression with an elastic net and 10-fold cross-validation based on 2 17 ≈ 100k Quasi Monte-Carlo samples. To attenuate the computational cost of determining the PCE weights in high dimensions, the PC basis is truncated at a maximum order of 3 and first-order interactions (p = 3, q = 0.75, see Figure 1.8). A summary plot of the PCE is given in Figure 4.3 as a means of visual validation of the surrogate. The large uncertainty range associated with this variable makes it the most impacting uncertainty in this configuration. Its index stands alone in this study, dwarfing all effects from other uncertain parameters.

• The pre-exponential constant A 12 : The 12 th reaction, H + O 2 = OH + O, has been known for its very high sensitivity for a long time [START_REF] Baulch | Evaluated kinetic data for high temperature reactions Volume I, Homogeneous gas phase reactions of the H2-O2 systems[END_REF]. Because of this, it has since been extensively studied and, as a result, its uncertainty factor (1.1) is the lowest of Konnov's kinetic mechanism. However, it remains the most impacting uncertain kinetic parameter.

Dimension reduction

Even if the PCE is truncated, the resulting Sobol indices are sufficient to rank variables and get an adequate estimate of the portion of variance explained by each variable or group of variables. In this study, we choose to truncate the uncertain space to ensure that 95% of the variance of the QoI is retrieved, based on the prediction of the PCE. Only two uncertain variables are needed to explain that portion of the variance of the QoI: ξ T and ξ 12 , corresponding to the uncertainty on the co-flow temperature and 12 th reaction rate respectively.

The reduced uncertain space spanned by these two variables is resampled to validate the truncation. The resulting probability distribution of the QoI is compared to its original distribution in Figure 4.5. The shape of the distribution is well captured. The peak has both the correct position and the correct amplitude and the kurtosis of the distribution is also well reproduced.

Therefore, the uncertainties of the co-flow temperature and of the 12 th reaction preexponential factor are enough to explain most of the variability of the auto-ignition delaytime of the homogeneous reactor corresponding to the most reactive mixture in the Cabra flame configuration.

Verification of the dimension reduction on a slightly more complex problem

The dimension reduction proposed in the previous section is valid for a homogeneous reactor. In this section, we aim, to verify if it still holds in a configuration more representative of the real Cabra flame. Therefore, we test it on a strained and auto-igniting non-premixed laminar 1D flame. As shown in Chapter 3, this flamelet is the basis of the combustion model used in the LES. Therefore, if the dimension reduction holds well for this configuration, we trust that it will also hold for the uncertain behavior of the lift-off height in the Cabra flame's LES.

A pair of Monte-Carlo experiments are conducted in the original, 32D space and in the reduced, 2D space, for reference and evaluation of the reduced uncertain space respectively.

The conditions are very similar to the flamelets presented in chapter 2. The only difference is the boundary temperature of the oxidizer size, which is slightly higher in the case of the Cabra flame. Boundary conditions in the nominal case are summarized in Table 4.3. Table 4.3: Nominal boundary conditions for the auto-igniting 1D flame: temperature and species molar fractions. The pressure is set to 101 325 Pa.

Fuel side

Oxidizer side

T fuel [K] 305 T ox [K] 1045 X H 2 [-] 0.25 X H 2 O [-] 0.1 X N 2 [-] 0.75 X O 2 [-] 0.15 X N 2 [-] 0.75
For the same reasons as in the study presented in Chapter 2, the strain rate a is set to a = 10 s -1 to ensure the ignition of the flamelet even for slow realizations of the uncertain mechanism or low initial temperatures.

The flamelet equations (Equation 2.1) are solved in the mixture fraction space on a uniform 400-point mesh using the 2 nd -order stabilized explicit Runge-Kutta method (ROCK2) with error control [START_REF] Abdulle | Second order Chebyshev methods based on orthogonal polynomials[END_REF]. The auto-ignition delay is determined when the mass fraction of the hydroxyl radical first reaches the threshold Y OH = 600 ppm in the computational domain. The quantity of interest, again, is the logarithm of the auto-ignition delay time. The computation time per sample can rise up to a dozen minutes. The cost of these simulations is still moderate but no longer negligible like it was for the homogeneous reactor. For this reason, the PDFs are approximated using only 2 11 = 2048 Quasi-Monte-Carlo samples drawn from the 32D Sobol' sequence for the reference PDF and 2 11 more Quasi-Monte-Carlo samples from the 2D Sobol' sequence for the PDF in the reduced uncertain space.

The resulting PDFs are presented in Figure 4.6. Again, the two PDFs collapse for the most part, which indicates that the dimension reduction holds for strained autoigniting non-premixed flamelets. Therefore, we are confident that the uncertain behavior of the Cabra flame will also mostly take place in the same reduced uncertain space: Span(ξ T , ξ A 12 ).

- A single LES of the Cabra flame costs about 100k CPU hours on the Irene KNL supercomputer that we used. Our CPU time allotment was enough for 30 simulations. We chose to reserve a third of these simulations for validation, while the two other thirds were assigned to the surrogate training.

We use Quasi Monte Carlo sampling (a Sobol' sequence), to draw the 20 training samples. This choice is motivated by:

• The superiority of regression surrogate methods over projection methods observed in the preliminary study in Chapter 2 which disqualifies cubature sampling in this case

• The low-discrepancy property of the sequence, although we recognize that a set of 20 samples does not maximize this property

• The convenience of increasing the cardinality of the training set, if deemed necessary and if more resources were available later on

The 10 validation samples were drawn using Latin Hypercube sampling, because:

• It has a good space coverage

• Its randomness protects against potential sampling bias Each of these samples is a Large Eddy Simulation of the Cabra flame, with a setup identical as the one presented in Chapter 3. Instantaneous solutions of these simulations are shown in Figure 4.7 for the nominal case and for the samples where the highest and lowest lift-off heights are found. fraction first reaches 600 ppm. A logarithmic preconditioning is then applied to these observations of the lift-off height in order to guarantee its positivity in the following.

Finally, we chose to use a Regression Kriging as a surrogate model because it gave the best results among tested methods in this case. The chosen surrogate's properties are given in Table 4.4 and the response surface itself is plotted in Figure 4.9. The other surrogate types that were investigated are detailed, plotted, and evaluated in Annex E. .9 suggests a good agreement between the data and the surrogate. This agreement is quantified in the following section.

ξ1 2 [ -] -2 0 2 ξ T [-] 0.0 0.5 1.0 H/D [-]

Surrogate validation

The first step of validating the surrogate is to produce a summary plot. While the model plot of Figure 4.9 informs on the overall shape of the response surface, it does not really frame the agreement between the surrogate and the observations. The summary plot, shown in Figure 4.10 drops any notion of surrogate shape to focus solely on the comparison between the observations and their approximations by the surrogate. The summary plot in Figure 4.10 shows that the Kriging surrogate reproduces correctly its training set. More importantly, the validation samples which, by definition, were not included in the training of the surrogate, are also correctly reproduced. This indicates that there is little to no overfitting in this surrogate.

To quantify the error of the Kriging surrogate relative to the validation samples, we use the error metrics introduced in Section 1.6: The Normalized Root Mean Squared Error (NRMSE), and the Normalized Maximum Absolute Error (NMAE):

NRMSE = ∥y validation -y predicted ∥ 2 ∥y validation ∥ 2 = 1.4% (4.5) NMAE = ∥y validation -y predicted ∥ ∞ ∥y validation ∥ ∞ = 2.1% (4.6)
These error results are extremely satisfying, especially compared to those of the best PCE (presented in Annex E) which were about 4 times higher for both metrics. Finally, the Regression Kriging whose response surface is plotted in Figure 4.9 will be used in the following.

4.4 Uncertainty propagation results for the Cabra flame's lift-off height in Span(ξ T , ξ 12 )

Uncertain behavior analysis

As a whole, the lift-off height of the Cabra flame is mostly sensitive to the co-flow temperature, but the 12 th reaction of the mechanism (H + O 2 = OH + O) also has a sensible influence, especially at low temperatures. This behavior is consistent with the sensitivity analysis performed a priori. The surrogate also shows a sudden increase in slope, both ξ T -and ξ 12 -wise when the temperature decreases. This is consistent with the results of Figure 2.3 (from Chapter 2, reproduced below), which plots the probability distribution of the auto-ignition delay time of a homogeneous reactor at z mr for different initial temperatures. This sudden increase in sensitivity to both variables at the low end of the temperature range translates the transition from two auto-ignition modes allowed by the kinetic mechanism:

• High-temperature auto-ignition, very fast, found in the "plateau" region which begins at about ξ T ≈ 0.5

• Low-temperature auto-ignition, much slower, which would be found for even lower co-flow temperatures

We did not investigate the Cabra flame in the conditions of weak auto-ignition [START_REF] Saytzev | Study of combustion of an adiabatically-heated gas mixture[END_REF][START_REF] Yetter | A combined stability-sensitivity analysis of weak and strong reactions of hydrogen/oxygen mixtures[END_REF] in this study. However, this combustion mode is encountered in Chapter 2 for homogeneous reactors and non-premixed auto-igniting flamelets. The increase in sensitivity to ξ 12 only reflects the higher uncertainty of the mechanism in the transition zone between weak and strong auto-ignition.

The surrogate model, duly validated, can then be resampled extensively to obtain the probability distribution of the lift-off height in the uncertain space spanned by ξ T and ξ 12 . 2 19 ≈ 500k samples are drawn from the 2D Sobol' sequence and the surrogate is evaluated for each sample. The resulting PDF is presented in Figure 4.11. The peak of the distribution, located around H/D ≈ 8.5, corresponds to the plateau of high-temperature auto-ignition, while the near-uniform component between H/D = 10 and H/D = 20 corresponds to the transition zone, for ξ T < 0.5, or T cof < 1045 K.

Comparison with Ji et al. (2019)

Ji et al. [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] proposed an uncertainty propagation study for the lift-off height of the simulated H 2 Cabra flame in 2019. Their simulation framework is Reynolds Averaged Navier-Stokes simulations, with a k-ε turbulence model, the Transported Probability Density Function approach to model combustion, and Interaction-by-Exchange-with-the-Mean approach for micro-mixing. They used the kinetic mechanism from Li et al. ( 2004) [START_REF] Li | An updated comprehensive kinetic model of hydrogen combustion[END_REF] with uncertainty factors borrowed from Konnov (2008)'s uncertain mechanism [START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF]. They did not propagate the uncertainties on the temperature measurement. However, they did qualitatively investigate the influence of co-flow temperature on the flame's liftoff height.

Using an active direction approach to reduce the kinetic uncertain dimension and a second-order PCE as a response surface, they obtained the PDF that we plotted in Figure 4.12 in the dotted line. It is visually different from the one plotted in their article [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] (Fig. 4b.) because they used a logarithmic scale in abscissa. Otherwise, the data points are strictly identical. We plotted on the same graph, in plain line, a PDF of the lift-off height obtained by resampling a slice of our surrogate model (trained with UF 12 = 1.1 as given by the Konnov 2019 mechanism) at ξ T = 0.5 ⇔ T cof = 1045 K in order to place ourselves in the same conditions as Ji et al.. The PDF is much narrower than Ji's proposed PDF and the expected value is also much higher: H/D expected = 11.8.

The divergence in expected value between [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] and our study can be explained by the differences between Li (2004)'s and Konnov (2019)'s nominal kinetic mechanisms. In particular, for the most sensitive reaction of both mechanisms, H + O 2 = OH + O, Li et al. propose a rate constant of 3.55×10 15 while Konnov proposes a rate constant of 1.04×10 14 .

We have three leads to explain the difference in PDF width showcased in Figure 4.12:

• Konnov (2008)'s uncertainty factor for H + O 2 = OH + O was 1.5, while its updated (2019) version only has an uncertainty factor of 1.1 for this reaction. This uncertainty reduction, stemming from a decade of further research, must have a significant influence on this problem. To evaluate this influence, we resampled our surrogate model for an uncertainty factor of 1.5 on A 12 . The response surface was not trained for such a wide range of possible values, but the linear trend of the Kriging guarantees an asymptotic behavior which is at least sound compared to a PCE. The result of this resampling is also plotted in Figure 4.12 in dashed line. This manipulation widens the PDF about three times and brings out the signature of the transition between low-and high-temperature combustion (or "weak" and "strong" combustion respectively): the onset of a bipolar distribution with the right-side bump in the PDF. However, this lead is not enough to explain the whole difference between Ji's and ours PDF.

• Our study reduced the uncertain dimension by discarding all variables that were not sensitive enough. Most of the uncertain reactions were discarded because of their insignificance before the uncertainty on the co-flow temperature. In particular, the next reaction in terms of sensitivity, H + HO 2 = 2 OH has a first-order Sobol' index about 2 times lower than H + O 2 = OH + O in 0D homogeneous reactor simulations in the conditions of the Cabra flame (cf Figure 4.4). Taking this uncertain parameter into account may help to further explain the difference in PDF width between [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] and our study.

• Finally, our most speculative lead is that Ji et al.'s choice of surrogate type may have artificially widened their proposed PDF. Fig. 4a. of their article [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] suggests that their second-order polynomial fit lacks curvature, which would result in an over-representation of high lift-off heights in the resampling.

As a conclusion, further work is needed to benchmark simulations and uncertainty propagation methods from different teams with modern uncertainty data. 4.13 shows a good agreement between the experimental measurements and our simulation: the predicted lift-off heights are well within the range of measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF]. Its evolution with regards to the co-flow temperature is also well reproduced. It also confirms that, in the H 2 Cabra flame configuration, taking into account the uncertain kinetic mechanism is only a refinement compared to the uncertainties introduced by uncertain measurements of the co-flow temperature. The reduction of kinetic uncertainties over the last decade, compiled in [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF], allows for much less uncertain simulations nowadays. [START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF], and a violin plot of the uncertainty quantification study from [START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF] However, there is another major axis that introduces uncertainties in numerical simulations of the H 2 Cabra flame: the turbulence and turbulence-combustion interaction models. The next section investigates this new aspect of the problem. Section 3.6.2 showed that the flame shape is extremely sensitive to the dissipation coefficient of the mixture fraction variance C D . Moreover, while the theoretical value of C D should lie between 2π 2 and 40 for typical LES filter sizes, much lower values can be found in the literature [START_REF] Domingo | Large-eddy simulation of a lifted methane jet flame in a vitiated coflow[END_REF]117,[START_REF] Ihme | Prediction of autoignition in a lifted methane/air flame using an unsteady flamelet/progress variable model[END_REF]. This disagreement between the values proposed by different authors and values obtained by theoretical developments can be considered as a parametric uncertainty in the turbulence-combustion interaction model. Without any more information on the probability distribution of C D , we assume that it is equiprobable over the entire range of possible values. Therefore, we model C D as a uniform random variable:

Comparison with deterministic simulation results and experimental measurements

ξ C D = C D -C D,Min C D,Max -C D,Min ∼ U (0, 1) (4.7)
Where C D,Min = 2.0 is the lowest value we found in the literature [START_REF] Domingo | Large-eddy simulation of a lifted methane jet flame in a vitiated coflow[END_REF] and C D,Max = 40.0 is the highest possible value predicted by analytical developments [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF].

Therefore, the uncertainties must now be propagated from the 3D space: Span(ξ A 12 , ξ T , ξ C D ). To this end, 29 more LES are run for this study: 19 more training samples are drawn from the 3D Sobol' sequence and 10 more validation samples are drawn with LHS.

The 30 LES drawn from Span(ξ A 12 , ξ T ) are also recycled to train and evaluate the surrogate model. In total, 39 training samples and 20 validation samples are used in the following.

As a side note, the first sample sample from the 2D and 3D Sobol' sequence are located at the same coordinates: they correspond to the nominal case. This is why only 19 samples were drawn from the 3D sequence, starting from the second element of the sequence.

Surrogate modeling

Several surrogate modeling methods were tested to emulate the behavior of the liftoff height in Span(ξ A 12 , ξ T , ξ C D ). They are presented and evaluated in Annex E. The best candidate among them was a PC-Kriging, whose parameters are given in Table 4.5. However, the performance of this surrogate was unsatisfying: The resulting response surface is much more satisfactory. Its root mean square is slightly lower and its maximum error is significantly reduced. A slice of the surrogate in the hyperplane ξ 12 = 0 is presented in Training and validation observations are absent from this plot because most of them don't belong in the hyperplane ξ 12 = 0 and it would make no sense to represent them.

As expected from the qualitative study of Section 3.6.2, low dissipation coefficients of mixture fraction variance lead to large simulated flame lift-off. Furthermore, no significant influence of this parameter is found for values of ξ C D > 0.5, which loosely corresponds to the range of values that should be used for typical LES sizes (2π 2 < C D < 40 [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF]). This is interesting for future deterministic LES of this type of flame because it suggests a "perfect" value of C D does not need to be determined prior to the simulation. At least in this configuration, the lift-off height of the flame is only weakly dependent on this parameter, as long as it is in the range recommended by Ihme [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF]. However, the flame shape is more and more sensitive to C D with decreasing values. Increased caution should [START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF] As expected, the resulting probability distributions are much wider than those obtained when considering only the chemical kinetics uncertain, as seen in Figure 4.14. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF]. For enhanced readability, measurement uncertainties of the co flow temperature are only plotted for the experiment from Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]. The error bars are the same for the other experiments.

The abrupt bottom side of the violin plots is due to the same artifact of the surrogate model that caused the minor bulge in Figures 4. [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] Taking into account the uncertainty of the turbulence-combustion interaction model can explain some of the variety of simulation results found in the literature. Moreover, the state-of-the-art Konnov (2019) mechanism [START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] is much less uncertain than other, more ancient mechanisms that were used when the other simulations were run and Section 4.4.3 already showed that the uncertainties in kinetic mechanisms of that era were already enough to explain most of the discrepancies between the numerical results of the literature. This explains why our study, as expected, using modern uncertainty factors, does not cover the whole range of previous simulation results. Finally, our uncertainty propagation study also explains most of the experimental measurements, considering that their uncertainty on the co-flow temperature is ±3%.

Only the smallest of lift-off heights observed in the experiments and predicted by other simulations and, most importantly, observed by Gordon et al. are seemingly not explained by our study. However, looking more closely at the foot of the flame for our observation with the shortest flame lift-off, we identified an early zone of progress variable production at the foot of the jet. This is illustrated in Figure 4. 21. This early ignition zone leads to a small -but visible -increase in hydroxyl mixture fraction near the jet exit. This suggests that smaller lift-off heights are accessible with our simulation setup, but that they would need slightly higher co-flow temperatures to appear. We are confident that future investigations at higher temperatures would confirm

Takeaway points

On the methodological side:

• A procedure to take into account highly dimensional kinetic and experimental uncertainties in LES is proposed in this chapter.

• Uncertainties not directly related to chemical kinetics can also be taken into account.

We demonstrate this by taking into account an uncertain variable in the turbulent combustion model.

On the physical interpretation side:

• Taking into account state-of-the-art kinetic and temperature uncertainties shows that deterministic simulations of the Cabra flame are not suited to predict its behavior. It may not be a very useful configuration to compare quantitatively turbulence or combustion models in a deterministic framework.

Chapter 5

Uncertainty Quantification in averaged LES fields

In this chapter, we aim to quantify the uncertainties of the averaged LES fields of the Cabra flame. That is, we aim to evaluate the PDF of every time-averaged quantity at every point in the computational domain. This will allow comparing the experimental and numerical ranges of uncertainty on local quantities. Furthermore, this work is a proof of concept for the prediction of uncertain averaged fields using LES for turbulent reacting flows.

We propagate the uncertainties coming from the kinetic mechanism, the co-flow temperature, and the dissipation coefficient of the mixture fraction variance. The fields considered in this study are the fields that were measured experimentally by Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF].

Sections 1 and 2 deal with uncertain dimension reduction of the input and output space respectively. The third section presents and validates the surrogate models used to predict uncertain fields. Section 4 validates the method as a whole and presents the results of the uncertainty propagation on the averaged LES fields of the Cabra flame.

Reducing the input uncertain dimension

In this study, we consider the uncertainties from the kinetic mechanism, the co-flow temperature, and the mixture fraction variance dissipation coefficient. We employ the same reduction of the input uncertain space as the one used in Section 4.5: The dominant variable approach on the auto-ignition delay time of a homogeneous reactor corresponding to the most reactive mixture of the Cabra flame eliminates a large number of weakly sensitive reaction rates. Finally, the reduced uncertain space is Span(ξ

Reducing the output uncertain dimension

The computational mesh of the LES of the Cabra flame comprises 5.5 million nodes. Considering that Cabra [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] This brute-force approach is quite expensive, and clearly suboptimal:

• There is a spatial auto-correlation of each field and a cross-correlation between fields which means that there would be some redundant work in constructing every surrogate independently.

• With only 39 training observations, it is highly probable that the surrogates would not be perfectly determined. Therefore, inconsistencies may appear, such as spatial discontinuities or the sum of mass fractions locally not equal to 1.

These reasons call for a drastic dimension reduction of the uncertain outputs.

Preliminary work on the spatial mesh

The averaged fields are interpolated on a structured, axisymmetric, hexahedral mesh which discretizes the region covered by experimental measurements (x ∈ [0, 26D] and r ∈ [0, 3D]). This 3D mesh has 315 × 250 × 20 ≈ 1.5M nodes. An azimuthal mean is then applied to the fields to enforce the axisymmetric nature of the fields and further reduce the mesh size to n nodes = 315 × 250 ≈ 79k.

This preliminary work brings the uncertain output dimension to about 4.8 million, which is still extremely large. Moreover, it does not address the two problems cited above. Further dimension reduction is therefore needed.

Uncertain dimension reduction using Principal Component Analysis

The ideal dimension reduction tool would have the following properties:

• Reduce the dimension of the problem to a tractable number (d out < 10)

• Allow the prediction of the averaged fields from the reduced variables without too much information loss

• Preserve the auto-correlation and cross-correlations of the fields Principal Component Analysis (PCA) fills these specifications and is used in this work.

For each observation, the averaged fields are centered and normalized using the "standard score" normalization [START_REF] Kreyszig | Advanced Engineering Mathematics[END_REF]. Then, the normalized averaged fields are queued up in a unique "features' vector" of size n features = n nodes × n fields ≈ 4.7 × 10 5 .

A PCA is applied to the set of observations. The modes obtained with this PCA represent the variations of the six different fields that were initially considered. This way, their cross-correlation is conserved.

The first eigenvalues of the PCA are plotted in Figure 5.2. The basis of PC modes is truncated such that only modes that correspond to more than 1% of the variance of the observations are retained. This truncature is materialized by the red line in These modes are enough to account for most of the possible variations of the Cabra flame in the uncertain space. A set of three weights -which can also be considered as coordinates in the basis formed by the three first modes of the PCA -are therefore enough to characterize a flame shape. The output uncertain dimension is reduced to three, compared to about 33 million initially. Equation 5.1 summarizes the prediction process of all the fields:

Y prediction (ξ) = Y + N w 1 (ξ) • M 1 + w 2 (ξ) • M 2 + w 3 (ξ) • M 3 (5.1)
Where Y prediction (ξ) is the vector of predicted field values, for each considered field and at each spatial location. It can be written as: Finally, N is a (n features , n features ) diagonal matrix containing the n fields normalization coefficients applied to each field before computing the PCA. This allows rescaling the weighted modes and obtain dimensionalized quantities: As with any type of dimension reduction, we can compute error metrics to measure the loss of information induced by the compression. For each field, we concatenate the field value at each point and for each validation sample in a single vector. This vector is created for the LES averaged fields and their projections on the first three PCA modes.

Y prediction (ξ) = Y OH,
The results are given in Table 5.1. These error scores are satisfying for the most part, in particular regarding the extremely low mean error (NRMSE). Only the Y OH field's projection is worse, with the validation samples VS12 and VS16 contributing the most to the NRMSE and VS12 being responsible for the disappointing NMAE. To better understand where these errors come from, we propose Figure 5.5. Each subfigure of Figure 5.5 presents the fields of hydroxyl mass fraction Y OH computed by the LES (left), compared with their projection on the first three modes of the PCA (right) for every validation sample. This visual representation of the projection suggests that no significant amount of information was lost in the compression. In fact, in most cases, the flame shapes are very well preserved by the compression. In particular, the flame lift-off height and axial evolution of the hydroxyl mass fraction are well reproduced. Furthermore, the explanation of the relatively high error scores observed with samples VS12 and VS12 appears clearly. Sample VS12 has a higher lift-off flame than any of the training samples. As a reminder, the PCA only generalizes its training dataset and cannot accomplish miracles for samples outside its training range. The case of VS16 is a bit subtler: The lift-off height (≈ 17D) is well within the range of previously-seen values. However, the foot of the flame -i.e. the area where Y OH values are non-negligible but lower than the auto-ignition threshold -is much longer than any of the training samples.

In both of these problematic cases, using more PCA modes did not improve significantly the agreement between the averaged LES fields and their projection. In terms of error metrics, using 5 modes, for example, the NMAE for the Y OH field stayed around 22% and the NRMSE decreased slightly to reach 5.5%. To improve the quality of the projection, a more extensive sampling of the input uncertain space would be neededmeaning more LES should be run.

Overall, the compression satisfyingly preserves flame shapes given our training sample budget.

Surrogate modeling

Preliminary remarks

The dimension reduction of inputs and outputs defined in the previous sections results in only three input variables (ξ A 12 , ξ T cof , ξ C D ) and three output variables (the weights of the three most important modes of the PCA, see Equation 5.1). The input variables are independent by nature. The three coefficients of the PCA modes may be statistically dependent. This is not a problem for us because we determine them as a function of the initial, independent uncertain parameters.

Choice of surrogates

Unlike in Section 4.5.2, we did not find any further simplification of the problem of emulating the coefficients of the PCA modes. Therefore, three surrogates must be fitted in the 3D input space Span(ξ A 12 , ξ T , ξ C D ). We select PC-Kriging (Kriging using a Polynomial Chaos Expansion as a trend) to emulate the PCA coefficients. PCEs of total order 2 are used for the first two coefficients, and a PCE of total order 1 (a linear trend) is used for the third coefficient. The hyperparameters of the Kriging models are fitted by minimizing the NRMSE in a Leave-One-Out cross-validation procedure (Section 1.6.4).

Verification of Surrogates

Unfortunately, 39 observations, 20 of which are in the 2D plane Span(ξ A 12 , ξ T ) (see Section 4.5.2), do not fill the 3D input space very densely. Therefore, it is doubly important to check the quality of the surrogate models. Figure 5.6 and Table 5.2 give a summary of the quality of the surrogates used to emulate the three PCA modes.

.6 suggests that the surrogate for the coefficient of mode 1 is relatively satisfying, although one of the validation samples is poorly predicted. This one poor prediction mostly drives the relatively high values of NRMSE and NMAE computed for this surrogate. The surrogates of the coefficients of modes 2 and 3 are a bit worse, both in terms of error metrics and visual agreement on the summary plot. Nevertheless, both of them correctly reproduce the trend of the coefficients they emulate, and they both avoid overfitting. Therefore, we consider these surrogate models satisfying despite their imperfections. Moreover, Figure 5.2 showed a steep decrease in the variance explained by the successive modes. This means that the greater the mode index, the less a slightly off coefficient prediction matters in the reproduction of the fields.

Validation of the field prediction, as a whole

Both the projection phase and the surrogate modeling of the coefficients have been verified. However, the combination of these two methods to predict a flame shape requires the validation of the method as a whole. Similarly to Section 5.2.2, we compute error scores for the prediction of the validation samples. The results are given, per field, in Table 5.3. The error scores are only slightly higher than those obtained in Section 5.2.2 to evaluate the projection step. This means that while the surrogate models that predict the coefficients of the modes are not perfect, most of the prediction error comes from a lack of richness of the PCA modes, in the sense that they restrict the possibilities of flame shape variation. Ultimately, as seen in Section 5.2.2, this comes from a lack of variety in the training observations. To enrich the training dataset, more LES should be run, each one costing 100k hours of CPU time. Due to our finite computational budget, we did not include any more observations in the study. Besides, the fidelity of the fields prediction is already more than satisfying in the probable regions of the uncertain space.

Finally, we provide a visual verification of the predictions with Figures 5.7 and 5.8 which compare the surrogate predictions and LES observations for the hydroxyl mass fraction and temperature fields respectively. The comparisons of the other fields can be found in Annex H. Chapters 3 and 4 highlighted the importance and the variability of the Cabra flame's lift-off height and its variations in the uncertain space. Therefore, a proper prediction of the time-averaged fields should also lead to a correct prediction of the lift-off height. Figure 5.9 shows a summary plot of the lift-off height, as observed in the LES (x-axis) and predicted using the method presented in this chapter (y-axis). The agreement between the observations and the predictions is almost perfect. This further validates our field prediction method.

It is also possible to resample extensively the 3D uncertain space Span(ξ A 12 , ξ T , ξ C D ), predict the time-averaged fields, and determine a PDF of the lift-off height. This is shown in Figure 5.10. Figure 5.10 shows a satisfying agreement between the PDF determined with the direct surrogate (in dotted line, from Chapter 4) and the PDF determined by post-processing the predicted fields. The general shape of the distribution is mostly respected. Perhaps more importantly, the key properties reported in Table 5.4 are extremely well reproduced.

Finally, Figure 5.10 shows a fundamental difference between both approaches. The lift-off height distribution determined from the predicted fields is bounded, as opposed to the non-bounded result of the direct approach from Chapter 4. This difference is related to an issue already reported in Section 5.2.2: The predicted fields are linear combinations of the PCA modes, which are themselves linear combinations of the fields observed in the training set. Therefore, while lift-off heights inside the training set range are extremely well predicted, as shown in Figure 5.9, our method cannot predict lift-off heights outside

In this section, we complete the work of Section 3.3.2 by plotting the nominal, 50%, and 95% confidence intervals of the uncertain simulation for all the profiles investigated by Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]. Figure 5.11 shows the axial profile of dioxygen mass fraction, Figure 5.12 shows several radial profiles of water vapor mass fraction, and Figure 5.13 shows several radial profiles of hydroxyl mass fraction. In these figures, we also represent the measurement uncertainties as reported by Cabra et al. All the other profiles are given in Annex I. The profiles shown in Section 3.3.2 showcased some minor disagreements between the experimental measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] and our deterministic LES. However, for the most part, these measurements fall well within our confidence interval when taking into account the uncertainties coming from the kinetic mechanism, the co-flow temperature, and the turbulence-combustion interaction model. Therefore, our modelization choices are suitable to explain most of the experimental measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF].

Two discrepancies of note remain:

• The initial slope of the axial mixture fraction decrease, already outlined in Section 3.3.2, is still off (see Figure I.6). This issue is also noticeable in other works that used the same LES solver in a similar configuration (the CH 4 Cabra flame) [117,[START_REF] Schulz | A criterion to distinguish autoignition and propagation applied to a lifted methane-air jet flame[END_REF]. There may be an issue with the numerical setup that our numerous tests could not ferret out, or a problem with the experimental measurements.

• The hydroxyl mass fraction confidence interval is still 50 to 100% higher than the measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] in the far-downstream locations. These high values are directly linked to the mixture fraction values found in the UFPV table. In the nominal case, Figure B.4 shows maximum OH mass fraction values higher than 1000 ppm at C = 1 (end of reaction), increasing with the strain rate. The relatively high values of OH found downstream in the LES are therefore not surprising. However, we did not have the time to investigate further the discrepancy with experimental measurements in the Cabra flame.

The 95% confidence intervals displayed in these profiles are extremely large. Figure 5.13, in particular, is striking: At the flame foot, between X = 9D and X = 14D, the hydroxyl mass fraction can range from 0 to over 1000 ppm, depending on whether auto-ignition has already occurred. This puts into perspective the "high fidelity" characterization of Large Eddy Simulation if uncertainties are not taken into account.

In light of this study, we conclude that the Cabra flame is an excellent benchmark case for uncertainty quantification in reactive flow simulations, which would warrant further work from other groups, but a poor benchmark case for turbulence or combustion models in a deterministic framework.

Limitations of the study

This study presented in this chapter is subject to a few limitations, due to work hypotheses that would warrant further verifications and/or corrections. They are presented in the following in increasing order of importance.

Sampling uniformity

The sampling of the 3D uncertain space Span(ξ A 12 , ξ T , ξ C D ) was not performed according to its probability distribution: of the 39 training samples and 20 validation samples, 20 and 10 respectively were actually sampled from the 2D uncertain space Span(ξ A 12 , ξ T ). This is due to our desire to recycle the LES results obtained for the study of Section 4.3. The idea was to enrich our knowledge of the 3D uncertain space with already-available observations in one of its subspaces.

To the author's knowledge, this inclusion of additional observations does not constitute an issue for the determination of the various surrogate models. On the contrary, it allows more accurate prediction on the plane Span(ξ A 12 , ξ T ) and does not degrade the predictions elsewhere.

On the contrary, over-representing an area of the uncertain space (compared to its associated probability distribution) may tamper with the optimality of variance explained by the PCA. In this case, however, we found that the inclusion of the observations from the 2D study did not fundamentally alter the shape of the three principal components of the Cabra flame. Rather, it resulted in "smoother" modes by removing some minor artifacts.

Overall, in this case, the addition of more samples in the training set of the PCA improved the quality of the modes. Although, we would recommend, when adding new observations in such a training set, to sample according to the uncertain space's probability distribution when possible.

The dimension reduction of input uncertainties

In this chapter, we used the same dimension reduction of the kinetic uncertainties as in Chapter 4 without justifying further. We showed in Chapter 4 that using a dominant variable approach allowed neglecting most of the uncertain Arrhenius constants. However, this dimension reduction was shown valid only for the lift-off height determination. In the case of the whole flame shape, a more complete study of the sensitivity of each quantity (mass fractions, gas temperature, reaction source term...) in homogeneous reactors as a function of time should be performed to guarantee the relevance of the dimension reduction. Such a study can be found in Dumont's thesis [START_REF] Dumont | Méthodes numériques et modèle réduit de chimie tabulée pour la propagation d'incertitudes de cinétique chimique[END_REF]. However, given the predominance of the sensitivity of the lift-off height (or autoignition delay in the case of the homogeneous reactor) to the co-flow temperature (or initial temperature), compared to kinetic parameters, we are confident that the predicted PDF for each quantity at each location in the Cabra flame would not change much if more kinetic parameters had to be taken into account.

Choice of the probability distribution of C D Our choice of the probability distribution of C D , both in terms of PDF shape and amplitude, is debatable. We chose what appeared to be the most sensible option at the time:

• A uniform distribution that reflects our lack of a priori knowledge on the "best" value to use

• A range that covered both the empirical range of values found in the literature and the theoretical range of values proposed by Ihme [START_REF] Ihme | Pollutant formation and noise emission in turbulent non-premixed flames[END_REF] for typical LES sizes. Confusingly, these two ranges were distinct, which lead us to consider a wider interval that included both empirical and theoretical ranges.

This wide distribution can be considered as a prior estimation of the parameter C D that could be narrowed down in the future using bayesian optimization. We show in Figure 5.14 the impact of choosing the most extreme cases of probability distributions for the Y H 2 O radial profile at X = 11D. On the left is the uncertain radial profile obtained with the uniform distribution of C D considered in the whole chapter. On the right is the same uncertain profile obtained with a Dirac distribution of C D at its nominal value (C D ∼ δ(20.0)).

Takeaway points

On the methodological side:

• In this chapter, we propose a method to quantify the uncertainties of averaged LES fields.

• A combination of PCA, to reduce drastically the number of uncertain outputs, and surrogate modeling allows predicting every desired field at every location of the uncertain space.

• While our choices of uncertain inputs may be debatable, the method proposed in this chapter can be applied again, should further research reveal a need to change our uncertain inputs.

On the physical interpretation side:

• In this chapter, we predicted LES averaged fields using dimension reduction and surrogate modeling. This allowed the propagation of high-dimensional uncertainties and the determination of the probability distribution of every relevant field at any location in the computational domain.

• Taking into account state-of-the-art kinetic, temperature, and turbulence model uncertainties shows, again, that deterministic simulations of the Cabra flame are not suited to predict its behavior. This point is redundant with the previous chapter, but the confidence intervals of the averaged LES fields are even more eloquent.

Conclusion Main achievements and results

• A set of uncertainty quantification methods, suitable for forward uncertainty quantification in the context of very expensive observations, has been selected and implemented in the in-house code SURE, which was entirely developed in this thesis. Its modularity and non-intrusive nature allow for combining different methods and interfacing with any "observations generator". For now, it has been interfaced with the in-house 0D-1D code AGATH, and with the LES code AVBP.

• The extreme uncertainty of H 2 /O 2 auto-ignition simulations has been demonstrated in configurations where the transition between strong and weak combustion is involved. This uncertainty stems from the kinetic mechanism itself -which remains quite uncertain even when using state-of-the-art mechanisms -and from the uncertainty on the initial temperature of the mixture: uncertainties as small as ±3% in the vicinity of 1000K can lead to orders of magnitude of uncertainty in the auto-ignition delay. The resulting uncertainty of the auto-Ignition Delay Time is non-trivial: the probability distribution of the IDT can even be bimodal in the worst-case scenario.

• The deterministic Cabra flame was satisfyingly simulated using LES and UFPV tabulation.

• Thanks to an a priori study on a simplified configuration, which dramatically reduced the input uncertain dimension, we made affordable the surrogate-based uncertainty quantification of the Cabra flame lift-off height stemming from the kinetic mechanism and co-flow temperature. We also investigated the uncertainties introduced by the combustion-turbulence interaction model.

• Finally, we introduced a much richer uncertainty propagation study: by predicting entire averaged fields at any location in the uncertain space without resorting to LES other than for a few training observations, we were able to determine the probability distribution of every measured quantity and compare them with experimental measurements. This study also proved capable of rediscovering the results of the direct uncertainty quantification of the lift-off height.

• This work shows that, while recent works in chemistry greatly reduced the uncertainty of H 2 /O 2 auto-ignition, the simulated Cabra flame remains very uncertain, mostly due to the experimental uncertainty on the co-flow temperature. This uncertainty affects first and foremost the lift-off height of the flame. In turn, every measured fields are eventually impacted.

A few words of nuance about the uncertainties uncovered in this thesis

The specificities of the Cabra flame

The large confidence intervals discovered in this study are likely not the norm in more conventional flame configurations. The Cabra flame is cursed -or gifted, depending on the point of view -with a very high sensitivity to parameters that happen to be slightly uncertain (the kinetic mechanism and the co-flow temperature) or that are lacking a scientific consensus, as far as the author is aware (the dissipation coefficient of mixture fraction variance).

The high sensitivity of the Cabra flame to these parameters is due to the driving mechanism of flame stabilization: auto-ignition delay. In more conventional configurations, the driving stabilization mechanism tends to be either mixing (for nonpremixed flames) or the presence of low velocity and strain-rate areas behind a bluff body or in the recirculation of a swirler. In these configurations, the uncertain parameters investigated in this study have a much lesser sensitivity.

The choice of uncertain inputs

• We worked with a state-of-the-art kinetic mechanism, which assumed every parameter to be certain, except for the pre-exponential constants of the Arrhenius laws. This is a choice of the chemistry community to represent kinetic uncertainties that way. However, to the author's knowledge, there is no reason other than convention to assume the other parameters to be certain. Should further research from the chemistry community come up with an uncertain mechanism featuring uncertain activation energies as well as pre-exponential constants, the results of the uncertainty of the Cabra flame lift-off height may be different. In any case, the methodology proposed in this thesis can work with any parametric uncertainty in the kinetic mechanism.

• Section 5.6 already pointed out the arbitrariness of the choice of the C D probability distribution. We chose to use the least informative distribution which covered both the empirical values found in the literature and those deduced from analytical developments. We are not claiming to hold any absolute truth about this probability distribution and another choice may have been wiser. Then again, our methodology can accommodate other distributions for C D .

• Finally, other sources of uncertainty may have been overlooked.

Perspectives Further improving the prediction of the averaged fields

This section proposes a few ideas to improve the study from Chapter 5.

• The statistical dependence of the mode coefficients may even have been an asset if we used a joint surrogate model of mode coefficients that takes into account the statistical dependence of its outputs, like co-kriging [START_REF] Mc Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF][START_REF] Le | Recursive co-kriging model for design of computer experiments with multiple levels of fidelity[END_REF]. Joint surrogate models were not implemented during this thesis due to a lack of time but would be an interesting improvement in future works.

• We saw that an important limitation of our results was our inability to predict rare events, in particular very high lift-off heights. We also saw that this inability was due to the lack of such examples in the training set, not due to an insufficient number of modes retained in the PCA. Training observations being expensive LES runs, adding more of them was not an option during this thesis. This problem is inherent to the PCA algorithm, and we leave finding a satisfying solution to future work.

Bayesian inference of the distributions of the co-flow temperature and dissipation coefficient of the mixture fraction variance: narrowing down our uncertain inputs

Our ability to predict the averaged fields of the Cabra flame nearly anywhere in the uncertain space indicates that we could also infer conjointly the probability distribution of the experimental co-flow temperature and of the "ideal" value for C D using the experimental multi-scalar measurements from Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]. We did not have the time to investigate this during this thesis. However, we believe that it would be very interesting for the combustion community:

• Reducing the uncertainty of the co-flow temperature in Cabra's experiment would allow numerical scientists to produce more relevant deterministic simulations, which would in turn allow more relevant combustion model comparisons on this flame.

• Fine-tuning the value of C D would bring fresh arguments to the scientific debate over the choice of dissipation coefficient value.

Applying the methods to other use cases

The methods presented in this thesis could be applied to other combustor configurations. Given their non-intrusive nature, they could also be applied in other domains. In the chemistry sector, they could for example be applied to improve [START_REF] Verleysen | How does a resilient, flexible ammonia process look? Robust design optimization of a Haber-Bosch process with optimal dynamic control powered by wind[END_REF] the (currently) extremely polluting Haber-Bosch process that keeps afloat our intensive agriculture model. Finally, they may also be used in any scientific domain where expensive observations meet high-dimensional uncertainties and where cheaper sub-problems can be used to reduce the uncertain dimension a priori. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

  Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 1: Handling uncertainties . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Methods overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 2 :

 2 Figure 2: Influence of the value of the Smagorinsky constant on the dynamics of the LES of the Sydney bluff-body flame. Adapted from Kempf et al. [50]

Figure 4 :

 4 Figure 4: Lift-off height of the Cabra flame. Measurements[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF] are the black symbols, simulations[START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF] are the colored lines. Note the measurement uncertainty of the co-flow temperature (black error bar)

Figure 5 :

 5 Figure 5: Temporal evolution of the probability distribution of species concentration in the plug-flow reactor, using 15000 Monte-Carlo samples. Solid lines represent median values and dashed values encompass the 95% confidence interval. Adapted from Phenix et al. [81]

Figure 6 :

 6 Figure 6: Radial profiles of time-averaged Temperature and CO mass fraction. Solid lines are the LES mean, dashed lines are ±3σ, and symbols are experimental measurements. Adapted from Mueller and Raman[START_REF] Mueller | Model form uncertainty quantification in turbulent combustion simulations: Peer models[END_REF]

Figure 1 . 1 :

 11 Figure 1.1: 800 points of a Monte-Carlo sampling of a 2D uniform random vector on the unit hypercube

Figure 1 . 2 :

 12 Figure 1.2: Illustration of the convergence of a Monte-Carlo estimator

Figure 1 . 3 :

 13 Figure 1.3: Illustration of the measure change method for a reduced and centered normal distribution

Figure 1 . 4 :

 14 Figure 1.4: Illustration of the separation of the principal part

Figure 1 . 5 :

 15 Figure 1.5: LHS sampling of a 2D uniform random vector on the unit hypercube

Figure 1 . 6 :

 16 Figure 1.6: Examples of the 800 first points of low-discrepancy sequences of a 2D uniform random vector on the unit hypercube

Figure 1 . 7 :

 17 Figure 1.7: Evolution of the number of basis polynomials M with the number of uncertain variables d and the total order of the basis p

Figure 1 . 8 :

 18 Figure 1.8: Illustration of the hyperbolic truncation [95]: Retained basis terms in the PCE when varying the parameter q and the total degree p = 3, 4, 5, 6 for d = 2 are highlighted by solid black points. The x and y axes correspond to the partial degree of the polynomial in ξ 1 and ξ 2 respectively.

Figure 1 . 9 :

 19 Figure 1.9: Example of tensorization of level 4 quadrature rules on the hypercube [0, 1] 2 . On the left are the evaluation points needed for the full tensorization. On the right are the points needed by Smolyak's method. Up is Clenshaw-Curtis' rule, and bottom is Fejer's second rule.

Figure 1 . 11 :

 111 Figure 1.11: Examples of kernel functions. From left to right and from up to down: Exponential, Squared-Exponential, Matérn 3/2 and Matérn 5/2, for d ij ∈ [0, 2]

Figure 1 . 13 :

 113 Figure 1.13: Illustration of a PCA in a 2D dataset, compressed to 1 dimension

Figure 1 . 14 :

 114 Figure 1.14: Illustration of the summary plot of a satisfying UQ method (left) and of a poor one (right)

Figure 2 . 1 :

 21 Figure 2.1: Instantaneous solutions of an example of a 1D auto-igniting flame with the Konnov 2019 mechanism, plotted in mixture fraction space. Solutions are plotted darker with increasing time from t = 0 ms up to 60 ms every 2 ms.

Figure 2 . 2 :

 22 Figure 2.2: Comparison of the auto-ignition delay time of the reactor using GRI-Mech 3.0 [101], Hong et al.'s 2011 [40], Konnov's 2008[START_REF] Konnov | Remaining uncertainties in the kinetic mechanism of hydrogen combustion[END_REF], and Konnov's 2019[START_REF] Konnov | Yet another kinetic mechanism for hydrogen combustion[END_REF] deterministic mechanisms, at different initial temperatures.

Figure 2 . 3 : 2 Figure 2 .

 2322 Figure 2.3: Violin plots of the logarithm of the homogeneous reactor's auto-ignition delay time at different initial temperatures, using both the original and boiled-down versions of Konnov's 2019 kinetic mechanism, with the initial composition from Sec. 2.1.2

Figure 2 . 4 :

 24 Figure 2.4: Temporal evolution of the temperature for the first 20 QMC samples

Figure 2 . 5 :

 25 Figure 2.5: PDF of the logarithm of the homogeneous reactor's auto-ignition delay time using 2 19 QMC samples

Figure 2 . 6 :

 26 Figure 2.6: Summary plot of the PCE used to compute the Sobol' indices

Figure 2 . 7 :

 27 Figure 2.7: Most significant Sobol' indices for the homogeneous reactor, plotted in decreasing order

Figure 2 . 8 :

 28 Figure 2.8: Logarithm of the auto-ignition delay time of the homogeneous reactor as a function of the temperature. The points are QMC samples from Span(ξ 1 , . . . , ξ 31 , ξ T ). Note the consistence of the results with the per-temperature investigation shown in Figure 2.3

Figure 2 . 9 :

 29 Figure 2.9: PDF of the logarithm of the 1D flame's auto-ignition delay time

Figure 2 . 10 :

 210 Figure 2.10: Comparison of Monte-Carlo experiments in the full and in the reduced uncertain space Span(ξ 12 , ξ T ). Left: On the homogeneous reactor. Right: On the 1D autoigniting flame.

Figure 2 . 11 :

 211 Figure 2.11: Convergence of the pdf obtained by resampling the regression-PCE, for several sizes of training sets

Figure 2 . 12 :

 212 Figure 2.12: Evolution of the Wasserstein distance between the reference pdf (in the 2D uncertain space) and the resampling for each surrogate modeling method as a function of the size of the training set

Figure 2 . 13 :

 213 Figure 2.13: Evolution of the alternative metrics for each modeling method as a function of the size of the training set, using 20 validation LHS samples. Left: NRMSE. Right: NMAE

Figure 3 . 1 :

 31 Figure 3.1: The experimental Cabra flame setup[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

  Figure 3.2 presents the results of experiments by Cabra et al., Gordon et al., and Wu et al.

Figure 3 . 2 :

 32 Figure 3.2: Experimental measurements of the normalized lift-off height of the Cabra flame as a function of the measured co-flow temperature.[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF]. The red error bar denotes the experimental uncertainty on the measurement of the co-flow temperature as reported in[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure 3 . 3 :

 33 Figure 3.3: Mixing in a shear layer behind a splitter plate[START_REF] Peters | Cambridge Monographs on Mechanics[END_REF]

Figure 3 . 4 :

 34 Figure 3.4: S-Curve computed in the operating conditions of the Cabra H 2 flame. The first and second critical points are circled in red. The extinct branch is shown in dotted line.

Figure 3 . 5 :

 35 Figure 3.5: Evolution of the gas temperature in a flamelet of strain rate a = 10 s -1 , plotted in mixture fraction space. Solutions are plotted every 10 ms, darker with increasing time.

Figure 3 . 6 :

 36 Figure 3.6: Comparison of detailed and tabulated chemistry on the auto-ignition of a homogeneous reactor at z = z mr = 0.04 The evolution of temperature is plotted in plain line for the detailed chemistry and in dotted line for the tabulated chemistry.

a = 5 s - 1 a = 470 s - 1 SFigure 3 . 7 :

 1137 Figure 3.7: Slices of the filtered mass fraction of OH Y OH as a function of the filtered mixture fraction z and of the filtered normalized progress variable C, for different strain rates a and different mixture fraction segregation values S z . Note the evolution of the vertical scale

Figure 3 . 8 :Figure 3 . 8 :

 3838 Figure 3.8: Mesh used in the LES of the Cabra flame. (a): 3D view and fluid boundary conditions. (b): Axial slice in the y-normal plane

Figure 3 . 9 :

 39 Figure 3.9: Velocity profile obtained from the RANS simulation

Figure 3 . 10 :

 310 Figure 3.10: Root mean square velocity profile obtained from the RANS simulation (in black), and empirical value found to fit the potential core length with the experiment (in red)

Figure 3 .

 3 Figure 3.11: y-normal cut of the Favre-filtered instantaneous fields of mixture fraction and normalized progress variable

1 .Figure 3 . 13 :Figure 3 . 14 :

 1313314 Figure 3.13: Radial profiles of the mixture fraction Z in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure 3 . 15 :

 315 Figure 3.15: Hydroxyl mass fraction mean field. Left: Experimental measurements from [11]. Right: This study.The white dots on the left-hand side denote the locations of multi-scalar measurements, reported in Figures 3.12, 3.13 and 3.14, as well as in Annex C for every other available scalar of interest

, 3 .

 3 Figure 3.15: Hydroxyl mass fraction mean field. Left: Experimental measurements from [11]. Right: This study.The white dots on the left-hand side denote the locations of multi-scalar measurements, reported in Figures 3.12, 3.13 and 3.14, as well as in Annex C for every other available scalar of interest

Figure 3 . 16 :Figure 3 . 17 :

 316317 Figure 3.16: Temporal convergence of the Hydroxyl mass fraction mean field.

Figure 3 . 18 :

 318 Figure 3.18: Illustration of the partial averaging of instantaneous solutions to produce independent blocks of data suitable for the bootstrap procedure

Figure 3 . 19 :Figure 3 . 20 :

 319320 Figure 3.19: Comparison of resolved and modeled mixture fraction variance on the flame axis

Figure 3 . 21 :

 321 Figure 3.21: Comparison of two hydroxyl mass fraction profiles for the original and finer mesh

Figure 3 .

 3 Figure 3.2 shows that the experimental flame is sensitive to the co-flow temperature, with the lift-off height exponentially increasing when the co-flow temperature decreases. To verify whether the LES of the Cabra flame reproduces this pattern, a set of four additional simulations were performed at varying temperatures, spanning the range of uncertainty reported by Cabra. One new UFPV table was generated for each of these conditions. Figure3.22 shows the mean hydroxyl mass fraction Y OH for each of these simulations. This figure does highlight a noticeable sensitivity of the flame to the co-flow temperature.Figure3.23 reprises the experimental data from Figure3.2 to compare with the liftoff height predicted by the simulation. The nominal simulation, which features a coflow temperature of 1045 K, is in good agreement with Cabra's experiment. The trend given by[START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF] and[START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF] is also correctly reproduced. The exponential increase of the liftoff height with decreasing co-flow temperatures was unfortunately not observed in this study. However, we believe that simulations with co-flow temperatures of about 1000 K would show this increase, because the exponential increase of the auto-ignition delay-time with decreasing oxidizer temperature was observed with this kinetic mechanism during the uncertainty propagation studies presented in chapter 2, both for homogeneous autoigniting reactors and non-premixed auto-igniting flamelets.

TFigure 3 . 22 :

 322 Figure 3.22: Hydroxyl mass fraction mean fields plotted for five simulations with varying coflow temperature (1015, 1030, 1045, 1060 and 1075 K)

Figure 3 . 23 :

 323 Figure 3.23: Lift-off height as a function of the co-flow temperature. Numerical simulations are compared with results from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF].

Figure 3 . 24 :

 324 Figure3.24: Lift-off height as a function of the co-flow temperature. Numerical simulations from this thesis are compared with other numerical results from[START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF].

  [START_REF] Constantine | Sparse pseudospectral approximation method[END_REF] for typical LES filter size, with an asymptotic value of π 2 for large Reynolds numbers.In this study, a value of C D = 20 (C χz = 2.0 and C ε /C u = 10) was found to lead to a satisfying agreement between the simulation and Cabra's experimental data. However, given the wide variety of values found in the literature, between theoretical arguments and actual values implemented, we deemed it important to check the sensitivity of the macroscopical behavior of the flame to this coefficient.

Figure 3 .Figure 3 . 25 :

 3325 Figure 3.25 presents the mean field of OH mass fraction Y OH computed with different mixture fraction variance dissipation coefficient values. We witnessed a significant sensitivity of the flame shape, and in particular of the flame lift-off, to the value of C D .

ωFigure 3 . 26 :

 326 Figure 3.26: From left to right: Subgrid mixture fraction variance, Progress variable source term and Hydroxyl mass fraction mean fields plotted for two different simulations of the nominal Cabra flame, for the extreme values of C D (2.0 and 40.0)

Figure 4 . 1 :

 41 Figure 4.1: Comparison of the evolution of IDT (0D) to lift-off height (LES) for 2 kinetic mechanisms. All auto-ignition delays are log-preconditioned and normalized by the IDT at a co-flow of 1045K with the Konnov (2008) mechanism. All lift-off heights are normalized by the lift-off height at a co-flow of 1045K with the Konnov (2008) mechanism.

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the temperature, as a function of time, of the first 20 realizations of the homogeneous reactor

Figure 4 . 3 :

 43 Figure 4.3: Summary plot of the PCE used to compute the Sobol indices

Figure 4 . 4 :

 44 Figure 4.4: First-and second-order Sobol indices of the auto-ignition delay

Figure 4 . 5 :

 45 Figure 4.5: Comparison of auto-ignition delay probability distribution of the homogeneous reactor, between the reference sampling of the 32D uncertain space and the sampling of the reduced uncertain space spanned by ξ T and ξ 12

0 log 10

 10

Figure 4 . 6 : 12 4. 3

 46123 Figure 4.6: Comparison of auto-ignition delay probability distribution of the strained flame, between the reference sampling of the 32D uncertain space and the reduced uncertain space spanned by ξ T and ξ 12

Figure 4 .

 4 Figure 4.7: y-normal slices of the instantaneous fields of water vapor mass fraction for the nominal (left), lowest (center), and highest (right) observed flames in the training set

Figure 4 . 8 :

 48 Figure 4.8: Mean fields of hydroxyl mass fraction for every training sample

Figure 4 . 9 :

 49 Figure 4.9: Response surface of the Regression Kriging, plotted alongside the training and validation samples

Figure 4 . 10 :

 410 Figure 4.10: Summary plot of the Regression Kriging

Figure 2 . 3 (

 23 Figure 2.3 (reproduced): Violin plots of the auto-ignition delay τ in the conditions of the Cabra flame for different co-flow temperatures for the original and boiled-down Konnov (2019) [55] mechanisms

Figure 4 . 11 :

 411 Figure 4.11: Probability distribution of the lift-off height in the uncertain space spanned by ξ T and ξ 12 using 2 19 surrogate evaluations

  /D) (scaled) ξ T = 0.5 , UF 12 = 1.1 ξ T = 0.5, UF 12 = 1.5 Ji et al. (2019)

Figure 4 . 12 :

 412 Figure 4.12: Probability distribution of the lift-off height as estimated by Ji et al. [45] and in our work for the fixed, nominal co-flow temperature of T cof = 1045 K and uncertain kinetic mechanism

Section 4 . 4 .

 44 2 showed the PDF of lift-off height conditioned to T cof = 1045 K. This procedure can be generalized, to several co-flow temperatures in the range of measurement uncertainties. The resulting violin plots are shown inFigures 4.13

  and 4.14, compared respectively to experimental measurements and to other numerical works.

Figure

  Figure 4.13 shows a good agreement between the experimental measurements and our simulation: the predicted lift-off heights are well within the range of measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF]. Its evolution with regards to the co-flow temperature is also well reproduced. It also confirms that, in the H 2 Cabra flame configuration, taking into account the uncertain kinetic mechanism is only a refinement compared to the uncertainties introduced by uncertain measurements of the co-flow temperature.

Figure 4 .

 4 Figure 4.14, on the other hand, shows that state-of-the-art kinetic uncertainties are not enough to explain the spread of numerical results. Section 4.4.2 already highlighted that older mechanisms were more uncertain. According to the PDF proposed by Ji et al., the relatively high uncertainties from [54] are representative of the spread of kinetic mechanisms from that era and suffice to explain the different experimental results from the literature.

  Cabra et al. (2002) Wu et al. (2006) Gordon et al. (2005) a. Gordon et al. (2005) b.This study (LES/UFPV)

Figure 4 . 13 :

 413 Figure 4.13: Violin plots of the lift-off height of the Cabra flame for several co-flow temperatures, compared to the measurements from [11, 120, 34]. For enhanced readability, measurement uncertainties of the co-flow temperature are only plotted for the experiment from Cabra et al. [11]. The error bars are the same for the other experiments.

  Stankovic et al. (LES/CMC) Cao et al. (RANS/PDF) Han et al. (LES/PDF) Han et al. (LES/PDF/ISAT) Naud et al. (RANS/UFPV) Patwardhan et al. (RANS/CMC) This study (LES/UFPV) Ji et al. (2019)

Figure 4 . 14 :

 414 Figure 4.14: Violin plots of the lift-off height of the Cabra flame for several co-flow temperatures, compared to deterministic simulations from[START_REF] Stanković | LES-CMC Simulations of Different Autoignition Regimes of Hydrogen in a Hot Turbulent Air Co-flow[END_REF][START_REF] Richard Cao | Turbulent lifted flames in a vitiated coflow investigated using joint PDF calculations[END_REF][START_REF] Han | LES/PDF modeling of autoignition in a lifted turbulent flame: Analysis of flame sensitivity to differential diffusion and scalar mixing time-scale[END_REF][START_REF] Naud | RANS modelling of a lifted H2/N2 flame using an unsteady flamelet progress variable approach with presumed PDF[END_REF][START_REF] Patwardhan | CMC simulations of lifted turbulent jet flame in a vitiated coflow[END_REF], and a violin plot of the uncertainty quantification study from[START_REF] Ji | Quantifying kinetic uncertainty in turbulent combustion simulations using active subspaces[END_REF]

4. 5

 5 Investigating the uncertainty in the turbulence-combustion model 4.5.1 The uncertain dissipation coefficient of mixture fraction variance

Figure 4

 4 response surface at ξ 12 = 0

Figure 4 . 16 :

 416 Figure 4.16: Response surface of the aggregated surrogate in the hyperplane ξ 12 = 0

Figure 4 . 18 :Figure 4 . 19 :

 418419 Figure 4.18: Probability distribution of the lift-off height in the uncertain space spanned by ξ C D , ξ T , and ξ 12 using 2 19 surrogate evaluations. Note that the x-scale is different from the PDF in the 2D uncertain space shown in Figure 4.11

Figure 4 . 20 :

 420 Figure 4.20: Violin plots of the lift-off height of the Cabra flame for several co-flow temperatures, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF][START_REF] Wu | An experimental investigation of the turbulence structure of a lifted H-2/N-2 jet flame in a vitiated co-flow[END_REF][START_REF] Gordon | Further Characterisation of Lifted Hydrogen and Methane Flames issuing into a Vitiated Coflow[END_REF]. For enhanced readability, measurement uncertainties of the co flow temperature are only plotted for the experiment from Cabra et al.[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]. The error bars are the same for the other experiments.

  and 4.16 near ξ C D = 1.

Figure 4 . 21 :

 421 Figure 4.21: Comparison of the progress variable and hydroxyl mass fraction fields between the lowest flame simulated (left) and the nominal case (right)

  A 12 , ξ T , ξ C D ). Therefore, we can re-use the simulations used in Chapter 4: the results of 39 LES are used as training observations, and 20 more are used as validation observations. The averaged fields of hydroxyl mass fraction for the training samples are shown in Figure 5.1. The fields of every considered averaged LES variable for the training samples are shown in Annex F.

Figure 5 . 1 :

 51 Figure 5.1: Averaged Y OH fields of the training samples

Figure 5 . 2 .Figure 5 . 2 :

 5252 Figure 5.2: Variance explained by each mode of the PCA

Figure 5 . 3 :

 53 Figure 5.3: Mean and retained modes for Y OH

w 1 ,

 1 w 2 , w 3 are the scalar weights associated with the first three modes of the PCA. They depend on the position ξ in the uncertain space and are predicted in Section 5.3.

Figure 5 . 4 :

 54 Figure 5.4: Mean and retained modes for the gas temperature

Figure 5 . 5 :

 55 Figure 5.5: Y OH fields of the validation samples compressed by the PCA (right) compared to the LES (left)

Figure 5 . 6 :

 56 Figure 5.6: Summary plots of the surrogates used to predict each mode's coefficient

Figure 5 . 7 :Figure 5 . 8 :

 5758 Figure 5.7: Y OH fields of the validation samples predicted using the PCA modes and their coefficients predicted by their surrogates (right) compared with the LES (left)

Figure 5 . 9 :

 59 Figure 5.9: Summary plot of the lift-off height computed from the predicted fields, compared to the LES results

Figure 5 . 11 :

 511 Figure 5.11: Nominal, 50%, and 95% confidence interval of the O 2 mass fraction's axial profile in the LES, compared to the measurements from Cabra et al.[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure 5 . 13 :

 513 Figure 5.13: Nominal and 95% confidence interval of the OH mass fraction's radial profiles in the LES, compared to the measurements from Cabra et al.[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure 5 . 14 :

 514 Figure 5.14: Uncertain radial profile of Y H 2 O at X = 11D. Left: with the uniform distribution of C D . Right: With a Dirac of C D , at its nominal value.

Figure 5 .

 5 Figure 5.14 shows little differences between the uncertain profiles of Y H 2 O at X = 11D, aside from the location of the 95% confidence interval, although X = 11D is a location

Figure C. 1 :Figure C. 2 :Figure C. 3 :Figure C. 4 :Figure C. 5 :Figure C. 6 :Figure C. 7 :Figure C. 8 :Figure C. 9 :Figure C. 10 :Figure C. 11 :Figure C. 12 :Figure D. 1 :Figure D. 2 :Figure D. 3 :Figure D. 4 :Figure D. 5 :Figure D. 6 :Figure D. 7 :Figure D. 8 :Figure D. 9 :Figure D. 10 :Figure D. 11 :Figure D. 12 :Figure D. 13 :Figure D. 14 :Figure F. 1 :Figure F. 2 :Figure F. 3 :Figure F. 4 :Figure F. 5 :Figure F. 6 :

 1234567891011121234567891011121314123456 Figure C.1: Axial profile of the mixture fraction Z in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure G. 1 :Figure G. 2 :Figure G. 3 :

 123 Figure G.1: Mean and retained modes for Y H 2

Figure G. 4 :Figure G. 5 :

 45 Figure G.4: Mean and retained modes for Y OH

Figure G. 6 :Figure H. 1 :Figure H. 2 :Figure H. 4 :

 6124 Figure G.6: Mean and retained modes for Z

Figure I. 1 :

 1 Figure I.1: Nominal and 95% confidence interval of the H 2 mass fraction's axial profile in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 2 :

 2 Figure I.2: Nominal and 95% confidence interval of the H 2 O mass fraction's axial profile in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 3 :

 3 Figure I.3: Nominal and 95% confidence interval of the O 2 mass fraction's axial profile in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 4 :

 4 Figure I.4: Nominal and 95% confidence interval of the OH mass fraction's axial profile in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 5 :

 5 Figure I.5: Nominal and 95% confidence interval of the gas temperature's axial profile in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 7 :

 7 Figure I.7: Nominal and 95% confidence interval of the H 2 mass fraction's radial profiles in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 8 :

 8 Figure I.8: Nominal and 95% confidence interval of the O 2 mass fraction's radial profiles in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Figure I. 10 :

 10 Figure I.10: Nominal and 95% confidence interval of the OH mass fraction's radial profiles in the LES, compared to the measurements from[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]

Table 2 .

 2

	3: Major statistics of the auto-ignition delay time of the homogeneous reactor from
	Section 2.1.2
		Auto-ignition delay time τ (ms)
	mean median nominal 95% confidence interval
	45.4	15.3	14.4	[ 1.7 , 180.6 ]

Table 2 .

 2 4: Major statistics of the auto-ignition delay time of the 1D flame

		Auto-ignition delay time τ (ms)
	mean median nominal 95% confidence interval
	52.1	12.0	12.8	[ 2.6 , 340.2 ]

Table 3 .

 3 1: H 2 Cabra flame experimental conditions[START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] PDF simulation. This conclusion is shared by another experimental study byWu et al.

	Central jet		Co flow
	Q H 2 (slm) 25	Q H 2 (slm)	225
	Q N 2 (slm) 75	Q air (slm)	2100
	T jet (K)	305	T coflow (K)	1045
	V jet (m/s) 107	V coflow (m/s) 3.5
	Re jet	23600 Re coflow	18600
	D jet (mm) 4.57	d coflow (mm) 210
	X H 2	0.25	X H 2 O	0.1
	X N 2	0.75	X O 2	0.15
	t lips (mm) 0.89	X N 2	0.75
	Q: volumetric flow rate; X: mole fraction;
	Re: Reynolds number; D: diameter; t:thickness

Table 3 .

 3

		2: Discretization of the UFPV table
	Variable Discretization type	Num. of points Bounds
	z	Regular	400	[0, 1]
	S z C	Geometric (increasing steps) Double-geometric (increasing then decreasing steps)	20 153	[0, 1] [0, 1]
	χ z,st	Geometric (increasing steps)	10	[1.6, 555]

Table 4 .

 4 4: Properties of the surrogate model used for the Lift-off height in Span(ξ A 12 , ξ T )

	Surrogate Type	Regression Kriging
	Trend	Linear
	Kernel	Matérn 5/2
	Length-scales (θ) determination Maximum likelyhood approach
	Nugget	0.05

Table 4 . 5 :

 45 Properties of the best surrogate model found in Span(ξ A 12 , ξ T , ξ C D ) (Sec. 4.3), the response surfaces shown in Figures 4.9 and 4.15 are actually multiplied, not added, to obtain the aggregated surrogate of the flame lift-off. Its validation scores are the following :

	Surrogate Type	PC-Kriging
	Trend	Regression PCE, p = 3, q = 0.75
	Kernel	Matérn 5/2
	Length-scales (θ) determination Maximum likelyhood approach
	Nugget	0.05
		NRMSE PCK = 7.5%	(4.8)

  measured local values for Y OH , Y O 2 , Y H 2 , Y H 2 O , Z, and the temperature, a total of 33 million local surrogates should be constructed and 33 million PDFs evaluated.

  1 , . . . , Y OH,n nodes , Y O 2 ,1 , . . . , Y O 2 ,n nodes , (5.2) . . . Z 1 , . . . , Z n nodes Y is the vector of the mean of the training observations at each location for each field. Its shape is the same as Y prediction .

Table 5 .

 5 1: Error scores of the projection of the validation observations on the first three PCA modes, per field

	Y H 2	Y H 2 O	Y O 2	Y OH	T	Z
	NRMSE 2.27% 1.53% 1.21% 7.31% 1.16% 1.29%
	NMAE 5.92% 9.06% 5.19% 21.96% 6.49% 4.29%

Table 5 . 2 :

 52 Error metrics for the surrogate models

	Mode Surrogate Type PCE trend order NRMSE NMAE
	1	PC-Kriging	2	27.36% 49.10%
	2	PC-Kriging	2	32.32% 39.44%
	3	PC-Kriging	1	55.45% 67.10%

Table 5 .

 5 3: Error scores for the prediction of the validation observations, per field

	Y H 2	Y H 2 O	Y O 2	Y OH	T	Z
	NRMSE 3.98% 3.76% 3.15% 11.15% 2.49% 1.79%
	NMAE 7.97% 16.06% 14.37% 29.56% 11.69% 4.47%

Table A .

 A 1: H 2 /O 2 combustion mechanism. k = AT β exp (-E a /RT ), UF is the uncertainty factor of the pre-exponential constant A. Enhanced third-body efficiencies (relative to N 2 ): Ar=0.72, H 2 O=16.6, O 2 =1.0, H 2 =1.5, He=0.57 H 2 O 2 (+M)=OH+OH(+M) O=7.5, H 2 O 2 =7.7, O 2 =1.2, N 2 =1.5, He=0.65, H 2 =3.7 O+H

	Index Reaction	A	β	E a	UF
	H+H+M=H 2 +M	7.00E+17	-1	0	2
	Enhanced third-body efficiencies (relative to Ar):
	O 2 =0.0, N 2 =0.0, H=0.0, H 2 O=14.3
	H+H+H 2 =H 2 +H 2	1.00E+17 -0.6	0	2.5
	H+H+N 2 =O 2 +N 2	5.40E+18 -1.3	0	3.2
	H+H+H=H 2 +H	3.20E+15	0	0	3.2
	O+O+M=O 2 +M	1.00E+17	-1	0	2
	Enhanced third-body efficiencies (relative to Ar):
	O=28.8, 02=8.0, N 2 =2.0, H 2 O=5.0
	O+H+M=OH+M	6.75E+18	-1	0	3
	Enhanced third-body efficiencies:
	H 2 O=5.0
	H 2 O+M=H+OH+M	6.06E+27 -3.312 120770 2
	Enhanced third-body efficiencies (relative to Ar):
	H 2 O=0.0, H 2 =3.0, N 2 =2.0, O 2 =1.5, He=1.1
	H 2 O+H 2 O=H+OH+H 2 O	1.00E+26 -2.44 120160 2
	H+O 2 (+M)=HO 2 (+M)	4.66E+12 0.44	0	1.2
	Low pressure limit	1.23E+19 -1.2	0
	Fcent = 0.5
		2.00E+12	0.9	48750 1.5
	Low pressure limit	2.49E+24 -2.3	48750
	Fcent = 0.42
	Enhanced third-body efficiencies (relative to Ar):
	H 2 2 =OH+H	5.08E+04 2.67	6292	1.3
	H+O 2 =OH+O	1.04E+14	0	15286 1.1
	H 2 +OH=H 2 O+H	2.14E+08 1.52	3450	1.5
	OH+OH=H 2 O+O	2.67E+06 1.82	-1647 1.4
	HO 2 +O=OH+O 2	2.85E+10	1	-723.9	2
	H+HO 2 =OH+OH	7.08E+13	0	300	2
	H 2 O+O=H+HO 2	2.20E+08	2	61600	5
	H 2 +O 2 =H+HO 2	7.40E+05 2.43	53500 1.5
	HO 2 +OH=H 2 O+O 2	7.00E+12	0	-1093	2
	HO 2 +OH=H 2 O+O 2	4.50E+14	0	10930	2
	HO 2 +HO 2 =H 2 O 2 +O 2	1.03E+14	0	11040 2.5
	HO 2 +HO 2 =H 2 O 2 +O 2	1.94E+11	0	-1409 1.4
	H 2 O 2 +H=HO 2 +H 2	5.02E+06 2.07	4300	3
	H 2 O 2 +H=H 2 O+OH	2.03E+07 2.02	2620	3
	H 2 O 2 +O=HO 2 +OH	9.55E+06	2	3970	3
	H 2 O 2 +OH=HO 2 +H 2 O	1.74E+12	0	318	1.5
	H 2 O 2 +OH=HO 2 +H 2 O	7.59E+13	0	7269	1.5
	H+O 2 +H=H 2 +O 2	8.80E+22 -1.835	800	2
	H+O 2 +H=OH+OH	4.00E+22 -1.835	800	2
	H+O 2 +O=OH+O 2	7.35E+22 -1.835	800	2
	H+O 2 +OH=H 2 O+O 2	2.56E+22 -1.835	800	2

• The reduction of kinetic uncertainties in the last decade has led to a dramatic decrease in lift-off height uncertainty in the simulated H 2 Cabra flame.• The uncertainties on the turbulence-combustion interaction model are about as sensitive as the uncertainty on the co-flow temperature. Therefore, special care must be exercised when setting model parameters, in configurations similar to the H 2 Cabra flame, i.e. auto-igniting non-premixed lifted flames.

Remerciements

Takeaway points

• A set of methods has been selected for the forward uncertainty propagation in the context of very expensive observations.

• Most of these methods have been implemented from scratch in the in-house UQ code SURE, which was entirely developed during this thesis. A few others relied on on-the-shelf libraries. SURE is designed as a modular software which allows for combining UQ methods and pluging in any deterministic model. For the needs of this thesis, interface modules for AGATH (an in-house 0D and 1D flame simulation software) and AVBP [START_REF] Schonfeld | Steady and Unsteady Flow Simulations Using the Hybrid Flow Solver AVBP[END_REF] (reactive, compressible LES software) were developed.

• This chapter is written in the form of an accessible point of entry for future PhD students discovering uncertainty quantification.

NMAE PCK = 18.6% (4.9)

To improve the quality of the surrogate model, we tried to reduce again the uncertain dimension in order to reduce the distance between observations. The discovery of an active direction was attempted but to no avail. Then, we noticed that the problem was mostly additive, which means that interactions between ξ C D and the other uncertain variables were small enough to be neglected. This allows writing:

In this case, the 2D surrogate proposed in Section 4.3 was recycled and substracted from the observations in the 3D space. A 1D surrogate could then be constructed from the remaining data projected on the direction of ξ C D . The best surrogate model we found, in this case, was a regression Kriging with linear trend, presented in Table 4. [START_REF] Baudin | Open-TURNS: An industrial software for uncertainty quantification in simulation[END_REF] Using the additivity of the problem, both the 2D and 1D surrogates are finally aggregated to form a 3D surrogate in Span(ξ A 12 , ξ T , ξ C D ), as shown in Equation 4. [START_REF] Hamada | Experiments: Planning, Analysis, and Optimization[END_REF]. As a side note, because we submitted the observations to a logarithmic preconditioning therefore be exercised when analyzing simulations using C D values in the lower range.

Finally, the minor bulge in the response surface, which can be observed in Figures 4. [START_REF] Colin | Development of High-Order Taylor-Galerkin Schemes for LES[END_REF] Aside from the error metrics presented above, the aggregated surrogate can be visually verified against its training and validation observations using the summary plot presented in Figure 4.17, which compares the observations (on the x-axis) with the surrogate evaluations at the same coordinates in the uncertain space (on the y-axis). The surrogate is not perfect. In particular, it fails to accurately predict lift-off heights in the higher range, i.e. when the flame is lifted above 15 D. However, the general agreement of the surrogate with the training and validation observations is satisfactory.

Probability distribution of the lift-off height in the 3D space

The 3D surrogate constructed in Section 4.5.2 is sampled extensively using 2 19 evaluations at coordinates drawn from the 3D Sobol' sequence. The resulting PDF is presented in Figure 4.18.

The shape of the PDF in Span(ξ C D , ξ T , ξ 12 ) is very close to the one obtained in Span(ξ T , ξ 12 ), presented in Figure 4.11. However, the region of high probabilities around H/D = 9 is less peaked and the tail of the distribution is much heavier, due to the large lift-off heights permitted by the combination of low co-flow temperatures, slow-chemistry -which both slow down auto-ignition -and low dissipation coefficient of the mixture fraction variance -which enlarges the zone where auto-ignition does not even start due to a too high strain rate. of this range. This is also the cause of the bulge in the PDF around H/D = 20: lift-off heights that should have been in the tail of the distribution are wrongly predicted at the upper end of the range of the training range, which, in turn, leads to an overestimation of the probability at this location.

To improve the prediction of the lift-off height, and of the general flame shape for the rare events above H/D = 20, observations of such flames should be included in the training set of the PCA. However, an over-representation of rare events may tamper with the optimal variance reproduction property of the PCA. Due to a lack of time and computational budget to run more LES, we did not investigate this aspect any further. Furthermore, despite this issue, the proposed method is already very satisfying in the probable regions.

5.5.2

Comparison of the uncertain fields with the multi-scalar measurements from Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] The ability to predict time-averaged fields given a set of coordinates in the input uncertain space allows the determination of a PDF for every time-averaged quantity at every location in the original LES computational domain. Among other possible use cases, it allows more relevant comparisons with experimental measurements than simply comparing the deterministic simulation with a set of experimental readings. Figure 5.12: Nominal, 50%, and 95% confidence interval of the H 2 O mass fraction's radial profiles in the LES, compared to the measurements from Cabra et al. [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF] where the gas can be either ignited or not depending on the coordinates in the uncertain space. Similar conclusions can be reached by comparing other fields. This is consistent with Figure 4.18, in Chapter 4, which compared the PDF of the lift-off height of the Cabra flame in Span(ξ A 12 , ξ T , ξ C D ) and Span(ξ A 12 , ξ T ): The probability distributions of every variable of interest are more heavy-tailed in the 3D uncertain space than they are in the 2D space, but the probable regions are similarly populated.

Annex A

The Boiled-down Konnov (2019) mechanism Annex B

The UFPV table

0.000 0.005 0.010 0.015 0.020

0.000 0.005 0.010 0.015 0.020

0.000 0.005 0.010 0.015 0.020

0.000 0.005 0.010 0.015 0.020

0.000 0.005 0.010 0.015 0.020

0.000 0.005 0.010 0.015 0.020 Figure I.12: Nominal and 95% confidence interval of the mixture fraction's radial profiles in the LES, compared to the measurements from [START_REF] Cabra | Simultaneous laser raman-rayleigh-lif measurements and numerical modeling results of a lifted turbulent H2/N2 jet flame in a vitiated coflow[END_REF]