Mr Didier Lucor

Mr. Sébastien Mr Chris Lacor

Mr Stéphane Ducruix

Mr Stéphane Richard

Moreau

Mr Franck Nicoud

Nicole Boutet

Michèle Marie Labadens

availability, cooperation at all times, especially for the numerical data that was shared which helped to perform

Abstract

Thermoacoustic instabilities result from the interaction between acoustic pressure oscillations and flame heat release rate fluctuations. These combustion instabilities are of particular concern due to their frequent occurrence in modern, low emission gas turbine engines. Their major undesirable consequence is a reduced time of operation due to large amplitude oscillations of the flame position and structural vibrations within the combustor. Computational Fluid Dynamics (CFD) has now become a key approach to understand and predict these instabilities at industrial readiness level. Still, predicting this phenomenon remains difficult due to modelling and computational challenges; this is even more true when physical parameters of the modelling process are uncertain, which is always the case in practical situations. Introducing Uncertainty Quantification for thermoacoustics is the only way to study and control the stability of gas turbine combustors operated under realistic conditions; this is the objective of this work.

First, a laboratory-scale combustor (with only one injector and flame) as well as two industrial helicopter engines (with N injectors and flames) are investigated. Calculations based on a Helmholtz solver and quasi analytical low order tool provide suitable estimates of the frequency and modal structures for each geometry. The analysis suggests that the flame response to acoustic perturbations plays the predominant role in the dynamics of the combustor. Accounting for the uncertainties of the flame representation is thus identified as a key step towards a robust stability analysis.

Second, the notion of Risk Factor, that is to say the probability for a particular thermoacoustic mode to be unstable, is introduced in order to provide a more general description of the system than the classical binary (stable/unstable) classification. Monte Carlo and surrogate modelling approaches are then combined to perform an uncertainty quantification analysis of the laboratoryscale combustor with two uncertain parameters (amplitude and time delay of the flame response).

It is shown that the use of algebraic surrogate models reduces drastically the number of state computations, thus the computational load, while providing accurate estimates of the modal risk factor. To deal with the curse of dimensionality, a strategy to reduce the number of uncertain parameters is further introduced in order to properly handle the two industrial helicopter engines.

The active subspace algorithm used together with a change of variables allows identifying three dominant directions (instead of N initial uncertain parameters) which are sufficient to describe the dynamics of the industrial systems. Combined with appropriate surrogate models construction, this allows to conduct computationally efficient uncertainty quantification analysis of complex thermoacoustic systems.

Third, the perspective of using adjoint method for the sensitivity analysis of thermoacoustic systems represented by 3D Helmholtz solvers is examined. The results obtained for 2D and 3D test cases are promising and suggest to further explore the potential of this method on even more complex thermoacoustic problems.

Keywords: Thermoacoustic instabilities, Helmholtz equation, Computational fluid dynamics, Uncertainty Quantification, Risk Factor, Monte-Carlo, Surrogate modelling, Active Subspace, Adjoint method.

Résumé

Les instabilités thermo-acoustiques résultent de l'interaction entre les oscillations de pression acoustique et les fluctuations du taux de dégagement de chaleur de la flamme. Ces instabilités de combustion sont particulièrement préoccupantes en raison de leur fréquence dans les turbines à gaz modernes et à faible émission. Leurs principaux effets indésirables sont une réduction du temps de fonctionnement du moteur en raison des oscillations de grandes amplitudes ainsi que de fortes vibrations à l'intérieur de la chambre de combustion. La simulation numérique est maintenant devenue une approche clé pour comprendre et prédire ces instabilités dans la phase de conception industrielle. Cependant, la prédiction de ce phénomène reste difficile en raison de sa complexité; cela se confirme lorsque les paramètres physiques du processus de modélisation sont incertains, ce qui est pratiquement toujours le cas pour des systèmes réels. Introduire la quantification des incertitudes pour la thermo-acoustique est le seul moyen d'étudier et de contrôler la stabilité des chambres de combustion qui fonctionnent dans des conditions réalistes; c'est l'objectif de cette thèse.

Dans un premier temps, une chambre de combustion académique (avec un seul injecteur et une seule flamme) ainsi que deux chambres de moteurs d'hélicoptère (avec N injecteurs et des flammes) sont étudiés. Les calculs basés sur un solveur de Helmholtz et un outil quasianalytique de bas ordre fournissent des estimations appropriées de la fréquence et des structures modales pour chaque géométrie. L'analyse suggère que la réponse de la flamme aux perturbations acoustiques joue un rôle prédominant dans la dynamique de la chambre de combustion. Ainsi, la prise en compte des incertitudes liées à la représentation de la flamme apparaît comme une étape nécessaire vers une analyse robuste de la stabilité du système.

Dans un second temps, la notion de facteur de risque, c'est-à-dire la probabilité pour un mode thermo-acoustique d'être instable, est introduite afin de fournir une description plus générale du système que la classification classique et binaire (stable / instable). Les approches de modélisation de Monte Carlo et de modèle de substitution sont associées pour effectuer une analyse de quantification d'incertitudes de la chambre de combustion académique avec deux paramètres incertains (amplitude et temps de réponse de la flamme). On montre que l'utilisation de modèles de substitution algébriques réduit drastiquement le nombre de calculs initiales, donc la charge de calcul, tout en fournissant des estimations précises du facteur de risque modal. Pour traiter les problèmes multidimensionnel tels que les deux moteurs d'hélicoptère, une stratégie visant à réduire le nombre de paramètres incertains est introduite. La méthode «Active Subspace» combinée à une approche de changement de variables a permis d'identifier trois directions dominantes (au lieu des N paramètres incertains initiaux) qui suffisent à décrire la dynamique des deux systèmes industriels. Dès lors que ces paramètres dominants sont associés à des modèles de substitution appropriés, cela permet de réaliser efficacement une analyse de quantification des incertitudes de systèmes thermo-acoustiques complexes. with the initial burner length L i =0.125 [m]. In this case, the Crocco's value n=6.57

is fixed and the time delay τ is varying over a period T = 1

f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s.114
4.13 The approximate estimate of the burner parameter length L i for predicting the growth rate of the 1 st azimuthal mode of the system with N injectors. The Flame/Acoustic interactions are considered for analytical computation purpose.

In this case n = 6.57 and τ = 9.87 × 10 -4 s and the growth rate ω i is 8.81

5.4

The uniform and the β-PDF of an arbitrary random variable X with similar mean (µ) and standard deviation (σ), but with different ranges (R) . . . . . . . . . . . . n (where n is the nominal value of n) and the plot on the bottom represents 345 the PDF of the time delay for the dimensionless ratio τ τ (where τ is the nominal value of τ ). In both plots, 10,000 ATACAMAC computations were generated. . . . 6.2 Monte Carlo analysis performed with ATACAMAC solver using for 10,000 samples generated with a Uniform distribution. . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Uncertainty Quantification strategy applied to the real annular helicopter engine 350 with 2 × N uncertain parameters (2 uncertain parameters per flame (n and τ )).

Initially, the Active Subspace method is employed to reduce the system dimensionality to only 3 variables. Then, algebraic surrogate models for the complete and reduced probabilistic spaces are used to analyse the surface response of the system. Finally, the Risk Factor is computed using the low-order models and val-355 idated against the brute force Monte carlo Analysis with ATACAMAC on 10000 samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 Eigenvalues of the finite difference approximation to the growth rate gradient of the full annular system with 2 × N uncertainties. Finally, only the 3 first active variables are relevant when using the theoretical studies of [START_REF] Noiray | Investigation of azimuthal staging concepts in annular gas turbines[END_REF] and Bauerheim et al. (2014a). Convergence analysis with different samples are used to converge eigenvalues: M = {10, 20, 50, 100, 500, 1000} samples. The spectrum is associated to {Re(γ), Im(γ)}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Comparisons between the approximated fIm and the true f Im response surfaces using linear and quadratic surrogate models, different sample sizes and two different input space is presented: in the full-dimensional space and the reduced space (3 active variables determined from the Active Subspace method). . . . . . . . . . . . 370 6.7 Convergence of the the low dimensional active subspace model when it is fitted with M=50 samples (Q 50 3D ) and M=2000 samples (Q 2000 3D ). These models are replayed 100000 times to evaluate the Risk Factor variability when comparing to the Risk Factor obtained from the benchmark Monte Carlo database (RF M C ). An overall good agreement is found with a relative good error below 6% when fitting the LIST OF FIGURES 7.3 Comparison between eigenvalues obtained for the first acoustic mode of each geometry using the adjoint method (∂f AD Im /∂τ ) and finite differentiation (∂f F D Im /∂τ ).

385

A step size δτ = δϵ × τ with (δϵ = 1 × 10 -6 ) is used. The growth rate gradients are evaluated when τ varies over a period T= 1 f 0 . . . . . . . . . . . . . . . . . . . .

7.4

Comparison between eigenvalues obtained for the first acoustic mode of each geometry using the adjoint method (∂f AD Im /∂n) and finite differentiation (∂f F D Im /∂n). A step size δn = δϵ×n is used for which δϵ = 1×10 -6 . The global flame amplitude 390 n is varied as reported in Table . 7

List of Tables

Part I

General introduction

580

Chapter 1

The physics of combustion instabilities

History and phenomenology

The inherent features of oscillatory combustion process have been a long-standing concern for en-585 gineers. Research on combustion instabilities have been quite extensive during the recent period and much still so far a challenging topic in a range of engineering applications, see Fig. 1.1 (propulsion systems, rocket engines, domestic boilers, furnaces, rocket engines, gas turbine combustors etc.). This keen interest is encouraged, in particular, by the variety of physical phenomena involved in the combustion process such as thermodynamic properties of chemical reactions and fluid dynamics of the system.

The dynamics of combustion instabilities could be described as excited unsteady motions of the flame front that stem from the coupled interaction between resonant combustor acoustics (in terms of pressure and velocity) and flame heat release rate oscillations from the combustion process. These heat release fluctuations are generally delayed with respect to incident disturbances (noise, modulation of mixture fluctuations, convection of hydrodynamic processes etc.) and give rise to an unstable growth of pressure oscillations. thermoacoustic instabilities are generally observed in high performance and modern combustion chambers in which the flame confinement associated to turbulent flow oscillations lead to these significant heat oscillations coupled with noise. This was discussed by [START_REF] Candel | Flame interactions as a source of noise and combustion instabilities[END_REF], [START_REF] Schuller | Dynamics of and noise radiated by a perturbed impinging premixed jet flame[END_REF], [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], [START_REF] Ihme | On the generation of direct combustion noise in turbulent non-premixed flames[END_REF][START_REF] Ihme | On the generation of direct combustion noise in turbulent non-premixed flames[END_REF]O'Connor andLieuwen (2012a). Under perturbed operating conditions, flow oscillations would potentially make small disturbances grow exponentialy. When this happens, undesirable effects may occur such as the melting of engine materials, irregular high temperature changes, large amplitude pressure oscillations, flame flashback or large amplitude structural vibrations with well-defined frequencies close to the natural resonant modes of the combustor [START_REF] Lynch | Alrest high fidelity modeling program approach[END_REF], [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], [START_REF] Huang | Dynamics and stability of lean-premixed swirl-stabilized combustion[END_REF]).

Therefore, the operability of the engine is engaged because the flame/acoustic interaction may lead in extreme cases to the complete failure of the combustor itself. Devastating consequences of combustion instabilities are presented in Fig. 1.2. The literature on combustion instabilities is extensive but the works of [START_REF] Candel | Combustion instabilities coupled by pressure waves and their active control[END_REF], [START_REF] Dowling | Acoustic Analysis of Gas Turbine Combustors[END_REF], [START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF], [START_REF] Lieuwen | A mechanism of combustion instability in lean premixed gas turbine combustors[END_REF] and [START_REF] Lieuwen | Combustion Instabilities in Gas Turbine Engines. Operational Experience, Fundamental Mechanisms and Modeling[END_REF], Poinsot and Veynante (2011) may be cited among others. why some advanced studies on combustion instabilities were realised later on. In 1859, Rijke has highlighted vibratory combustion in a self-excited acoustic oscillator that consists of a cylindrical duct (opened at both ends), and a thermal energy source. Rijke investigations pointed out that whenever a thermal energy source is placed in the upper or lower half of a vertical tube, the response of acoustic oscillations is different. Indeed, at the upper half of the tube a dampening of the oscillations occurred while, when the thermal energy source was placed in the lower half part, self-excited thermo-acoustic oscillations were observed. By providing additional explanations on combustion instabilities, the Rijke tube turned out to be a good experimental support that allows analytical modelling of acoustic fluctuations in terms of sound pressure level measurements and acoustic modes assessment. Rijke explanations were an important landmark in the scientific study of combustion instabilities and it motivated another famous experimental study, by [START_REF] Mallard | Sur la vitesse de propagation de l'inflamation dans les mélanges explosifs[END_REF], on this topic.

According to the seminal work of [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF], instabilities are encouraged when heat release fluctuations are in phase with pressure oscillations. This theory is known under the famous Rayleigh criterion and constitutes the baseline interpretation of combustion instabilities:

«If heat be periodically communicated to, and abstracted from, a mass of air vibrating (for example) in a cylinder bounded by a piston, the effect produced will depend upon the phase of the vibration at which the transfer of heat takes place. If the heat is given to the air at the moment of greatest condensation, or be taken from it at the moment of greatest rarefaction, the vibration is encouraged. On the other hand, if heat be given at the moment of greatest rarefaction, or abstracted at the moment of greatest condensation, the vibration is discouraged.»

From a technical point of view, throughout the 1940s and 1950s, instabilities observed in solid-and liquid-propellant rockets, afterburners and ramjets generated many of reviews and articles on rocket instabilities [START_REF] Crocco | Theory of combustion instability in liquid propellant rocket motors[END_REF], [START_REF] Harrje | Liquid propellant rocket instability[END_REF]). Concomitantly, the study of instabilities became a central importance in industries which use combustion processes. In 1976, Culick significantly contributed to a quantitative prediction of combustion instabilities by establishing a mathematical formulation of the Rayleigh's criterion. His formulation relates the direct transfer of the thermal energy to the mechanical energy of acoustical motions. [START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems -an overview[END_REF], [START_REF] Culick | Some recent results for nonlinear acoustic in combustion chambers[END_REF] [START_REF] Zinn | A theoretical study of nonlinear combustion instability in liquid-propellant rocket engines[END_REF], [START_REF] Williams | Combustion Theory[END_REF], [START_REF] Raun | A review of Rijke tubes, Rijke burners and related devices[END_REF], [START_REF] Howe | Acoustics of Fluid-Structure Interaction[END_REF].

Because of the increasing and powerful computational resources [START_REF] Abramson | An Automatic Design Optimization Tool and its Application to Computational Fluid Dynamics[END_REF], [START_REF] Staffelbach | Highly parallel Large Eddy Simulations of multiburner configurations in industrial gas turbines[END_REF]), a common concern in the combustion community is the modelling of combustion instabilities. The scope is to characterize earlier the propensity of any combustion process to become unstable. Extensive experimental researches have been conducted to mimic the complex physics involved in the combustion process of real gas engines [START_REF] Poinsot | Analyse des instabilités de combustion de foyers turbulents prémélangés[END_REF], [START_REF] Palies | Dynamique et instabilités de combustion de flammes swirlées[END_REF], Palies et al. (2011a), [START_REF] Worth | Modal dynamics of self-excited azimuthal instabilities in an annular combustion chamber[END_REF], [START_REF] Meijia | Effets de la temperature de paroi sur la reponse de la flamme a des oscillations acoustiques[END_REF]). The ability to reproduce in laboratory-scale the combustion instabilities which appear in real gas turbine engines offers the opportunity to economically reduce industrial costs and offer a set of solutions to tackle them. However, the experimental reproduction of a complex system is not always feasible and the comparison of experimental data to real gas turbine engines results in the same operating conditions is not obvious. As aforementioned, the numerical study of combustion instabilities is a cumbersome task because the mechanisms leading to the excitation of acoustic oscillations are both various and dependent to the prevailing system complexity (Palies (2010), Palies et al. (2011a), [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF]). However, engineers are still progressing and even recently, numerical methods have proved their effectiveness to study combustion dynamics in complex industrial geometries [START_REF] Staffelbach | Large eddy simulation of self-excited azimuthal modes in annular combustors[END_REF], [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF], [START_REF] Hermeth | Mechanisms affecting the dynamic response of swirled flames in gas turbines[END_REF], [START_REF] Bourgouin | Self-sustained instabilities in an annular combustor coupled by azimuthal and longitudinal acoustic modes[END_REF]). Yet, as no universal method has been developed to determine combustion instabilities in the development cycle of gas turbine engines, it is crucial to minimize the computational cost and to instigate mechanisms that govern the complete instability process at technological readiness level.

Driving mechanisms of instabilities

Regarding the stringent environmental requirements, particularly in regards to Nitrous Oxides (N o x ) production, modern gas turbine engines for power generation have been optimized for low pollutants emissions. Hence, the rate of Nitrous Oxides (N o x ) production has been significantly reduced by operating the combustion process at low temperatures (about 1800 -2000 K) [START_REF] Bibliography Lefebvre | Lean premixed/prevaporized combustion[END_REF], [START_REF] Delabroy | Passive and active control of nox in industrial burners[END_REF], [START_REF] Cheng | Lean premixed burners[END_REF]). Specifically, the operat-CHAPTER 1. THE PHYSICS OF COMBUSTION INSTABILITIES ing mode consists in injecting a homogenous mixture of fuel/air inside the burner to operate in lean-premixed flame regime (Littlejohn et al. (2002), [START_REF] Ulhaq | Identification of the dynamics of technically premixed flames as multiple-input, single-output systems from les[END_REF]). Nevertheless, leanpremixed flames are very close to the flame extinction limits under lean conditions operating design. The flame speed being considerably reduced at lower equivalence ratio [START_REF] Lieuwen | The role of equivalence ratio oscillations in driving combustion instabilities in low nox gas turbines[END_REF], [START_REF] Sattelmayer | Influence of the combustor aerodynamics on combustion instabilities from equivalence ratio fluctuations[END_REF], [START_REF] Richardson | Effects of equivalence ratio variation on lean, stratified methane-air laminar counterflow flames[END_REF], [START_REF] Hermeth | LES evaluation of the effects of equivalence ratio fluctuations on the dynamic flame response in a real gas turbine combustion chamber[END_REF]), the flame would become sensitive to any perturbations, the system stability is altered thus prompting to combustion instabilities. The flame front dynamics is primarily impacted by upstream acoustic flow rate fluctuations, as well as equivalence ratio inhomogenities. However, flame/vortex interactions [START_REF] Poinsot | Vortex driven acoustically coupled combustion instabilities[END_REF], [START_REF] Mueller | Vorticity generation and attenuation as vortices convect through a premixed flame[END_REF], [START_REF] Bougrine | Fuel composition effects on flame stretch in turbulent premixed combustion: Numerical analysis of flame-vortex interaction and formulation of a new efficiency function[END_REF]), flame perturbation with the system boundaries [START_REF] Popp | A numerical study of laminar flame wall interaction with detailed chemistry: wall temperature effects[END_REF], [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF], Tay-Wo-Chong and Polifke ( 2013)), chemistry [START_REF] Quillatre | Large eddy simulation of turbulent premixed flames propagation in a small scale venting chamber: influence of chemistry and transport modelling[END_REF][START_REF] Popp | A numerical study of laminar flame wall interaction with detailed chemistry: wall temperature effects[END_REF], [START_REF] Selle | Largeeddy simulation of turbulent combustion for gas turbines with reduced chemistry[END_REF]) or unsteady strain rate [START_REF] Echekki | Unsteady strain rate and curvature effects in turbulent premixed methane-air flames[END_REF], [START_REF] Creta | Strain rate effects on the nonlinear development of hydrodynamically unstable flames[END_REF]) may lead to an important increase of pressure fluctuations as well as large unsteady heat fluxes. Knowingly, the modelling of the combustion process response to flow perturbations is a critical component to determine both the qualitative and quantitative dependence of combustor stability on geometrical parameters, fuel composition parameters and kinematic processes leading to the flame/acoustic interactions.

As mentioned by [START_REF] Lawn | A model for the thermoacoustic response of a premixed swirl burner, part i: Acoustic aspects[END_REF], further studies on the flame response mechanisms, even on a statistical point of view, are needed to understand the nonlinear combustion instability process. All these complexities show how Uncertainty Quantification analysis of the flame/acoustic coupling would be relevant in the field of combustion instabilities.

The onset of the self-sustained coupling between pressure oscillations and flame may be detailed as follows:

⋄ When the inherent incoming flow features are perturbed, this can result in inducing vortex shedding or fluctuations of the equivalence ratio of the fuel/oxidizer mixture. Subsequently, heat release fluctuations are generated as well as convective modes such as entropic waves.

⋄ Heat release rate oscillations will then create harmonic pressure waves that propagate within the combustion chamber and may reflect on the walls, inlet and outlet (nozzle exit at the downstream end of the combustor) of the cavity. Moreover, the flame/acoustic interaction can be interpreted as a transfer of energy: the system would become unstable when an excess of energy released from the flame during the 710 quasi-isobaric combustion process disturbs the energy balance of the acoustic system. In case of favourable phasing between heat release rate of the flame and acoustic pressure perturbations, the driving mechanisms of oscillations are amplified. A commonly used criterion for determining the stability of a combustion chamber is the Rayleigh criterion (Rayleigh (1878)), which reads :

CHAPTER 1. THE PHYSICS OF COMBUSTION INSTABILITIES ∫ T ∫ ∫ ∫ Ω p ′ q ′ dΩ dt > 0 (1.1)
where p ′ and q ′ represent the pressure and heat release fluctuations respectively, Ω is the flow domain. Depending on the phase of oscillation, the sign of the integral may vary. To establish the stability of the system at a given frequency, Eq. (1.1) is integrated over a period. To understand further the underlying physics of combustion instabilities, it is possible to extend the Rayleigh criterion to accommodate the system being studied. This point has been discussed in the article of [START_REF] Nicoud | Thermoacoustic instabilities: should the rayleigh criterion be extended to include entropy changes?[END_REF], in which, for example, the Rayleigh criterion has been extended

to account for entropy changes. Other studies of [START_REF] Motheau | A mixed acoustic-entropy combustion instability in a realistic gas turbine[END_REF], [START_REF] Motheau | Analysis and modeling of entropy modes in a realistic aeronautical gas turbine[END_REF] are mentioning the acoustic-entropy impact on combustion instability. Veynante ( 2011))

The study of instabilities may be also achieved by monitoring the time evolution of pressure oscillation in a combustion chamber. This is illustrated in In Regime III, due to saturation, the growth of pressure amplitude fluctuations drop-off and limit cycle oscillations appear in the combustion chamber. Between the linear and non-linear transition regimes, Regime II, an overshoot period is visible for which the limit cycle amplitude is lower than the amplitude of pulsation.

Large amplitude limit-cycle oscillations should be avoided to prevent combustor damage. At the limit cycle state, the growth rate of the acoustic pressure disturbances is equal to zero due to an increase of acoustic losses or to the time lag change between the flame responses to acoustic pressure perturbations. Therefore, characterizing properly the characteristic time scales in the overall combustion process is necessary. In this work, only the linear regime will be considered which corresponds to the Regime I. At this stage, acoustic is linear and the oscillation over the mean value of the pressure (p ′ /p) are small. Uncertainty Quantification analysis will be performed to characterize quantitatively the risk of the flame/acoustic coupling to destabilize the system by varying both the time lags between heat release fluctuations, pressure oscillations and the amplitude of flame response. To achieve this, suitable numerical tools will be used to identify firstly the key mechanisms leading to instabilities and last but not the least, their computational cost will be evaluated to perform affordable UQ analysis.

About suppression methods of combustion instabilities

The control of combustion instabilities relies mainly on suppressing the coupling phenomena between heat release perturbations and acoustic waves oscillations. However, this task is not easy when considering all the processes involved in the combustion dynamics. Besides, the control of instabilities is truly dependent on the system complexity because under a particular operating condition, several natural modes of the combustor may be excited simultaneously. This has been highlighted in the work of [START_REF] Gulati | Active control of unsteady combustion-induced oscillation[END_REF], [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], [START_REF] Schmid | Applications of the dynamic mode decomposition[END_REF]. It is then necessary to identify the role of each mode to better use an effective control approach.

Two methods to control combustion instabilities have been developed since the late 1940's.

⋄ Passive control techniques : in this case, acoustic dampers as Helmholtz resonators or acoustic liners may be used to master unstable modes of the combustion chamber. Furthermore, drastic changes on operating conditions may help to decrease the driving of oscillations:

by modifying the fuel delivery system or the mixture mass flow rate, the phasing between heat release fluctuations and acoustic pressure disturbances can be better controlled.

changing the system geometry (nozzle modifications, injection system, swirler design etc.) can also help to damp oscillatory phenomena.

Further detailed investigations on passive control techniques of combustion instabilities have been realised for example by [START_REF] Noiray | Passive control of combustion instabilities involving premixed flames anchored on perforated plates[END_REF], [START_REF] Evesque | Low-order acoustic modelling for annular combustors: Validation and inclusion of modal coupling[END_REF], [START_REF] Lieuwen | Analysis of acoustic wave interactions with turbulent premixed flames[END_REF], [START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF], [START_REF] Magri | A theoretical approach for passive control of thermoacoustic oscillations: Application to ducted flames[END_REF].

⋄ Active control techniques : here, the system is force in such a manner to alter the instability cycle by providing additional energy from an external source. By adding an extra source of energy, the system could be favourably perturbed so as to damp the oscillations.Many advances on active control techniques have been realised on a variety of combustor design. Among them the work of [START_REF] Mcmanus | A review of active control of combustion instabilities[END_REF], [START_REF] Poinsot | Suppression of combustion instabilities by active control[END_REF], [START_REF] Paschereit | Control of thermoacoustic instabilities in a premixed combustor by fuel modulation[END_REF], [START_REF] Candel | Combustion instabilities coupled by pressure waves and their active control[END_REF], [START_REF] Poinsot | Using cfd to understand active control of flow resonators[END_REF], Poinsot et al. (1988), Faivre and[START_REF] Faivre | Experimental and numerical investigations of jet active control for combustion applications[END_REF], [START_REF] Huang | Dynamics and stability of lean-premixed swirl-stabilized combustion[END_REF], [START_REF] Bauerheim | Experimental and numerical studies of active and passive control of combustion instabilities[END_REF], [START_REF] Meija | On the experimental determination of growth and damping rates for combustion instabilities[END_REF] may be cited.

As aforementioned, passive control techniques are very costly because they imply drastic changes in the development time of the engine, they are not suitable under low range of frequencies (typically few hundred Hertz), they might consider changing fuel delivery system or some other modifications on the system specificities. In this case, the offline testing needs to be done again to assess whether any changes in the control parameters are necessary. This is challenging in the context of industrial readiness control of combustion instabilities. Active control approaches are suitable under low range of frequencies, they are more practical and they have proved their effectiveness on different types of operating condition.

Tools to study combustion instabilities

Several approaches are available to model and simulate combustion instability mechanisms. Generally, the method chosen depends on the system complexity but also on the computational resources available. In this thesis, as the goal is to perform Uncertainty Quantification analysis of thermoacoustic modes, special attention needs to be paid to the system complexity. Indeed, the more the system is complex the more the number of uncertain parameters may increase. Also, the CPU time is a key element because UQ studies rely on performing many calculations at a time, which can rapidly become out of reach. Therefore, for each case a choice has to be made to perform affordable UQ analysis using the more adapted tool.

⋄ Analytical models: Significant efforts have been deployed in developing theoretical models to study combustion instabilities [START_REF] Williams | Combustion Theory[END_REF], [START_REF] Dowling | Acoustic Analysis of Gas Turbine Combustors[END_REF], [START_REF] Clavin | One-dimensional vibratory instability of planar flames propagating in tubes[END_REF], [START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF], Bauerheim et al. (2014a), [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF], [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF], [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF]). These analytical models are mostly adapted to simplified academic cases because many assumptions are considered to render the problem tractable.

⋄ Experiments: Experimental set-up have also been developed in order to study thermoacoustic instabilities. For example, advanced research has been done to analyse combustion dynamics in swirled stabilized combustors and to study the propagation of azimuthal and longitudinal waves in combustors [START_REF] Balachandran | Experimental investigation of the nonlinear response of turbulent premixed flames to imposed inlet velocity oscillations[END_REF] an academic annular configuration with swirled premixed flames was built to study several matters as the interaction between flames and the effect of mean swirl on the system stability as well as the nature of azimuthal modes [START_REF] Worth | Modal dynamics of self-excited azimuthal instabilities in an annular combustion chamber[END_REF], Bourgouin (014a), Bourgouin et al. (014b), [START_REF] Bourgouin | A new pattern of instability observed in an annular combustor: The slanted mode[END_REF]). Another experimental study has been also realised to study for example the effects of wall temperature on the flame response to acoustic oscillations [START_REF] Meijia | Effets de la temperature de paroi sur la reponse de la flamme a des oscillations acoustiques[END_REF]).

⋄ Large Eddy Simulation: By solving the filtered Navier-Stokes equations, Large Eddy

Simulations tools appear to be tremendously powerful to capture combustion instabilities dynamics in complex gas turbines combustors [START_REF] Staffelbach | Large eddy simulation of self-excited azimuthal modes in annular combustors[END_REF], [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF], [START_REF] Hermeth | LES evaluation of the effects of equivalence ratio fluctuations on the dynamic flame response in a real gas turbine combustion chamber[END_REF], [START_REF] Ghani | LES of longitudinal and transverse self-excited combustion instabilities in a bluff-body stabilized turbulent premixed flame[END_REF]). However, the modelling of combustion instabilities when using LES approaches depends on several operating conditions: the choice of boundary conditions, chemical models, wall temperatures, spray characteristics etc. It is necessary to identify which of these parameters have the most significant impact to accurately predict unstable modes of a given system. For example, when studying numerically thermoacoustic instabilities, is it necessary to:

account for the detailed geometry of the combustor, -use a very refined mesh,

use sophisticated chemistry model, -take into account heat losses.

Also, LES techniques are known to be CPU expensive because they require solving the 3D Navier-Stokes equations at high Reynolds number as well as taking into account several physical phenomena such as acoustics and combustion. These difficulties with LES techniques have been the forerunners of new scientific investigations on the study of combustion instabilities using low order models as acoustic network or Helmholtz solvers.

⋄ Low order modelling methods: They are based on linear acoustics and are ideal to provide phenomenological interpretations of the results provided by experiment or LES with affordable numerical resources and time. In this approach, the thermoacoustic system is represented as a network of acoustic elements inter-connected to each other [START_REF] Munjal | Acoustics of Ducts and Mufflers[END_REF], [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]). Each of these acoustic elements is connected by using mathematical transfer function matrices. Acoustic low order network tools have been successfully used to study acoustic modes in academic and complex industrial combustors [START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF], [START_REF] Stow | Modelling of circumferential modal coupling due to helmholtz resonators[END_REF], Bauerheim et al. (2014a[START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF], [START_REF] Mensah | Efficient computation of thermoacoustic modes in annular combustion chambers based on bloch-wave theory[END_REF]).

⋄ Acoustic solvers: 3D acoustic solvers such as Helmholtz solvers are extensively used to study thermoacoustic instabilities [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF], [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF], [START_REF] Benoit | Numerical assessment of thermo-acoustic instabilities in gas turbines[END_REF]). To do so, the set of Navier-Stokes equations for reactive flows are manipulated to obtain an inhomogeneous wave equation for acoustic pressure disturbances.

Therefore, the eigenfrequencies, the growth or the decay of the modes, the limit cycle amplitude of the oscillations of a given three-dimensional geometry can be calculated in the frequency domain. For small amplitude pressure disturbances p ′ (⃗ x, t) = p(⃗ x)e -iωt , the proper equation reads :

γ(⃗ x)p 0 ∇ • ( 1 ρ 0 (⃗ x) ∇p(⃗ x) ) + ω 2 p(⃗ x) = iω(γ(⃗ x) -1)q(⃗ x). (1.2)
where c 0 = γ(⃗ x)p 0 ρ 0 (⃗ x) is the mean speed of sound and ω the complex valued pulsation, ρ 0 the mean density and q(⃗ x) represents the unsteady heat release: q ′ (⃗ x, t) = q(⃗ x)e -iωt . The detailed development of this equation is given in Chapter 3. The flame response to acoustic perturbation at reference locations is modeled thanks to a n-τ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF].

This formulation may also be related to the Flame Transfer Function formalism. Besides the n-τ model, matrix transfer formulation [START_REF] Polifke | Determination of thermo-acoustic transfer matrices by experiment and computational fluid dynamics[END_REF], [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationnary computational fluid dynamics[END_REF]) may be used to account for the flame/acoustic coupling. When the flame response is modelled, Eq. (1.2) corresponds to a force Helmholtz equation which is solved in the frequency domain as a non-linear eigenvalue problem. This is achieved by using adapted discretization approaches with numerical algorithm [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]) or analytical tools. On top of providing the structure of all thermoacoustic modes of the combustor, the resolution of this equation provides the set of complex frequencies of the system. The real part of the complex pulsation ω r is related to the frequency of oscillation f r = ωr 2π while the imaginary part ω i represents the growth rate of the acoustic pressure disturbances. When ω i is negative, the mode is stable and when ω i is positive, the mode is unstable and needs to be controlled. This is illustrated in Fig. 1.5 showing a typical result of a thermoacoustic analysis, e.g. a set of modes, each with its own frequency (ω r ) and growth rate (ω i ).

Chapter 2

Uncertainty Quantification

Motivations and objectives

Noticeable efforts have been increasingly deployed to develop powerful computational resources in a capacity to inform decision-making. Consequently, important improvements have been made on devices designed in the past few decades, which spawned drastically the reduction of experimental costs. Computational simulation becomes now a routine and a crucial step necessary to reproduce the time evolution dynamics of engineering applications in a realistic point of view.

Besides reproducing the physical processes in engineering devices, it contributes to the validation of experimental observations and theoretical investigations. This large advancement of computational techniques has greatly improved the applicability of complex industrial systems in terms of modelling and simulation. Such techniques are generally based on mathematical models that are approximated under specific assumptions to represent the relevant physics of the complex system. Mathematical models take commonly the form of partial differential equations (PDEs) that incorporate miscellaneous effects related to geometrical scaling, initial and/or boundary conditions. Afterwards, these models are turned into operative computer codes for simulation purposes. Thus, the computational models performances and failures depend not only on the conceptual and mathematical modelling assumptions, but also on the numerical discretization of the mathematical model, implementation of the numerical algorithms, constitutive model inputs, CHAPTER 2. UNCERTAINTY QUANTIFICATION domain settings and tolerances, numerical approximations, convergence criteria.

In the modelling and numerical simulation of engineering devices, uncertainties are encountered because of the lack of knowledge of the physical processes and the difficulty to identify distinctively the numerous parameters that are governing the system dynamics [START_REF] Hoffman | Propagation of uncertainty in risk assessments: the need to distinguish between uncertainty due to lack of knowledge and uncertainty due to variability[END_REF]. Even the smallest change in the mathematical model may lead to huge changes on the scientific understanding of the system behaviour. Arguably, under these conditions, results computed by mathematical models may differ from reality or observations. Consequently, it is generally difficult to define a level of confidence on numerical simulations robustness [START_REF] Yu | Uncertainty quantification for chaotic computational fluid dynamics[END_REF], [START_REF] Lucas | Rigorous verification, validation, uncertainty quantification and certification through concentration-of-measure inequalities[END_REF], [START_REF] Riley | Quantification of model-form and predictive uncertainty for multi-physics simulation[END_REF]Grandhi (2011), Oberkampf (2005), [START_REF] Iaccarino | Quantification of uncertainty in flow simulations using probabilistic methods[END_REF].

A quantitative method for evaluating numerical simulation accuracy is therefore needed.

In this thesis, Uncertainty Quantification methodologies are applied in the context of thermoacoustic instabilities that originate from the two-way interaction between the flame dynamics and acoustic waves propagation in combustion chambers. Robust approaches, whether they are based on LES techniques or on pure acoustic theories, are rather accurate in predicting the growth rate of thermoacoustic modes developing in complex geometries. However, strategies to estimate the uncertainty of the underlying thermoacoustic flame model have not been investigated yet.

The interests are in the development and application of stochastic computational strategies and algorithms for the solution of several specific Uncertainty Quantification problems. Different methodologies are used to quantify uncertainties, from the traditional brute force Monte Carlo method to surrogate modelling techniques or even to reduced basis methods that are used to tackle the «curse of dimensionality» encountered in high-dimensional and complex applications. Before getting to the heart of the matter, a literature survey on Uncertainty Quantification techniques and a brief description on the state-of-the-art methodologies employed to tackle Uncertainty Quantification problems are discussed in this introductory chapter.

Literature survey and basic definitions

The field of Uncertainty Quantification is as old as the theory of probability and mathematical statistics. Its outstanding success is due to the combination of probability and statistics in the wide spread use of modelling, large-scale computations and experimental studies (Aposto-CHAPTER 2. UNCERTAINTY QUANTIFICATION lakis (1990), [START_REF] Helton | An exploration of alternative approaches to the representation of uncertainty in model predictions[END_REF], [START_REF] Roache | Quantification of uncertainty in computational fluid dynamics[END_REF], [START_REF] Mathelin | Stochastic approaches to uncertainty quantification in cfd simulations[END_REF], [START_REF] Chanstrami | Uncertainty quantification in simple linear and non-linear problems[END_REF]).

In computational fluid dynamics, the development of numerical simulation tools has further bolstered the use of Uncertainty Quantification in a wide range of disciplinary sciences such as aerodynamics [START_REF] Lin | Predicting shock dynamics in the presence of uncertainties[END_REF], [START_REF] Beran | Uncertainty quantification of limit-cycle oscillations[END_REF]), meteorology [START_REF] Rochoux | Towards predictive datadriven simulations of wildfire spread -part i: Reduced-cost ensemble kalman filter based on a polynomial chaos surrogate model for parameter estimation[END_REF]), structural dynamics [START_REF] Hasselman | A top-down approach to calibration, validation, uncertainty quantification and predictive accuracy assessment[END_REF]) among others.

The goal is to ease the quantification of input and response uncertainties in a computational framework to provide quantitative information of scientific phenomena. For example, let's consider a physical model whose expression is given by f(Y). In this model, Y is the vector containing the input parameters of the system, Y={Y Risk does not exist by itself. Risk is created when there is uncertainty. Therefore, accounting for quantitative risk analysis implies to know at first the kind of uncertainties that are involved in the computational simulations. Generally, uncertainties are divided in two groups, [START_REF] Hofer | An approximate epistemic uncertainty analysis approach in the presence of epistemic and aleatory uncertainties[END_REF], [START_REF] Oberkampf | Advanced simulation and computing workshop error estimation, uncertainty quantification, and reliability in numerical simulations[END_REF], [START_REF] Iaccarino | Quantification of uncertainty in flow simulations using probabilistic methods[END_REF], [START_REF] Eldred | Mixed aleatory-epistemic uncertainty quantification with collocation-based stochastic expansions and optimization-based interval estimation[END_REF]:

⋄ Aleatory uncertainty: Also called irreducible uncertainty, aleatory uncertainty is due to variability or randomness nature of the model input parameters. The latter are generated by intrinsic perturbations of a physical system or random measurement errors. Because of the random nature of the model parameters, different scenarii of the system behaviour must be taken into consideration in this case. This is the reason why aleatory uncertainty and the resulting risk are modeled with a Probability Distribution Function (uniform distribution, β-distribution, normal distribution etc.). Such a PDF describes all the possible values of the input parameters and how they would impact the output quantities of interest. As an example, aleatory uncertainties are related to the outcomes of tossing dice and drawing cards from a shuffled pack.

⋄ Epistemic uncertainties: In contrast to aleatory uncertainty, epistemic uncertainty is also called reducible uncertainty. This type of uncertainty concerns for instance the lack of knowledge about the physical system. Different causes can explain this:

-Incomplete or imprecise knowledge of the underlying processes of the system -Alternative point of view on the characteristics of the system etc.

This type of uncertainty is called reducible because further research or investigations would help to decrease or to overcome the lack of knowledge on the system. The modeling of epistemic uncertainties is generally achieved through margins analysis or evidence theories Helton et al.;[START_REF] Helton | A sampling-based computational strategy for the representation of epistemic uncertainty in model predictions with evidence theory[END_REF]2009), Swiler et al. (2009a), [START_REF] Swiler | Epistemic uncertainty in the calculation of margins[END_REF][START_REF] Diegert | Toward a more rigorous application of margins and uncertainties within the nuclear weapons life cycle -a sandia perspective[END_REF], [START_REF] Jakeman | Numerical approach for quantification of epistemic uncertainty[END_REF].

«How far is it possible to push research activities to get further information of the system behaviour?» : The answer of this question is a way of providing a brief distinction between aleatory and epistemic uncertainties towards risk management analysis. Once the type of uncertainties identified, efficient probabilistic approaches can be challenged to propagate uncertainties in the system and to derive meaningful uncertainty bounds of the model simulations. Indeed, not only is it important to quantify uncertainties but also one ought to account for decisive and sustained policies to calibrate and validate physical model for simulation-based predictions or design.

As mentioned above, uncertainties appear in mathematical models in various contexts. The specification of a well-posed mathematical model to represent the underlying phenomena of engineering applications is usually the starting point of any realistic analysis. Today's significant and relevant challenge for computational science and engineering is to make sure that these mathematical models are solved efficiently and accurately to provide the behaviour of the system.

Concomitantly, strategies for numerically solving the mathematical model on a computer imply significant approximations that would influence the range of validity of the subsequent model outputs. This means that uncertainty is an unavoidable aspect of modelling engineering application behaviours, whether the model is deterministic or stochastic:

⋄ Deterministic models: The output of the deterministic models is completely assessed by the exact values of the input parameters and the operating conditions initially stated in the problem. This is the case of Isaac Newton's dynamic laws for example.

⋄ Stochastic models: Stochastic models possess some intrinsic randomness input quantities sometimes due to the fact that the measurements are not sufficient to produce precise inputs.

Therefore, the range of validity of the outputs is large for the same set of parameter values and initial conditions; for example the Poisson model for describing wavelet expansions.

In this work, an uncertainty quantification analysis of n-τ model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF]) used to represent the Flame/Acoustic coupling (as mentioned above in Section 1.4) is conducted. Typically, quantifying uncertainties of the flame model in thermoacoustic system is crucial because small changes of the input parameters of this model are known to have non-negligible impacts on the stability of the system. Moreover, the flame parameters n and τ vary a lot from an experiment to another. The characteristics of the flame model are succinctly discussed in Section 2.4 and fully detailed in Section 3.

Even after a strategy for solving the set of the governing equations of the mathematical model is chosen, quantifying and characterizing the resulting output uncertainty is an important issue to anticipate the intrinsic variability and the lack of knowledge of the underlying phenomena occurring in the system. Thus, including Uncertainty Quantification in the entire physical/mathematical modelling process is fundamental to provide a probabilistic representation of the output uncertainties in numerical simulation as presented in Fig. 2.2. Under operability limits of the system (limit cycles in thermoacoustics, reaction to unusually high loads, temperatures, pressures, high Reynolds CHAPTER 2. UNCERTAINTY QUANTIFICATION number etc.) performing Uncertainty Quantification analysis is even more interesting. This step aims at propagating the sources of uncertainties in the system.

2.3

2 When dealing with high-dimensional and complex systems, the number of uncertain parameters may drastically increase thus making difficult the propagation of uncertainties through the simulation. The more the system dimensionality increases and the more the number of variable necessary to represent its behaviour grow exponentially. Therefore, several evaluations of the underlying model would be required to sample the uncertainty space thus leading to an intractable computation burden even on today's powerful computers. In some cases, reduced basis approximation methodologies could be used to bypass these issues of dimensionality by estimating the principal subspaces of input variations. However, the use of such methodologies it is not always intuitive and obvious.

3 Once the PDF of the main input uncertain parameters generated, the simulation of the computational model is performed for all the possible random values for the input parameters of the vector Y. The response surface of the output quantity of interest Z is then estimated. This is typically the Monte Carlo method, but other methods can be used to propagate uncertainties through the system.

Extensive studies in this aspect of Uncertainty Quantification approaches are more and more developed to reduce the computational effort and to address the challenges of probabilistic robust design and optimization in multidisciplinary CFD simulations. Such methods allow tackling numerically the propagation of uncertainties in space dynamics by either intrusive or nonintrusive techniques [START_REF] Reagan | Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection[END_REF], [START_REF] Sudret | Global sensitivity analysis using polynomial chaos expansions[END_REF], [START_REF] Beran | Uncertainty quantification of limit-cycle oscillations[END_REF], [START_REF] Acharjee | A non-intrusive stochastic galerkin approach for modeling uncertainty propagation in deformation processes[END_REF]):

1 Intrusive UQ approaches:

Intrusive Uncertainty Quantification methods require some changes in algebraic operators of the underlying model in the source code. This has to be done carefully to ensure a proper analysis of the system under uncertainty.

2 Non-intrusive UQ approaches:

Unlike the above intrusive approaches, they use a deterministic black-box (no modifications in the solver) for uncertainty propagation of input uncertainties of the model. These kinds of non-intrusive UQ methods interpolate samples in the range of the input distributions.

However, sampling methods based on non-intrusive techniques are rather difficult to use when the dimensionality of the system increases.

In this thesis both intrusive and non-intrusive methods are used.

In this section, the classical computational methods used to propagate uncertainties are briefly described:

⋄ Monte Carlo [START_REF] Bose | Uncertainty analysis of laminar aeroheating predictions for mars entries[END_REF], [START_REF] Reagan | Uncertainty quantification in reacting-flow simulations through non-intrusive spectral projection[END_REF].

⋄ Reduced basis approaches such as Polynomial Chaos [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF], Le Maître and Knio ( 2007), [START_REF] Marzouk | Dimensionality reduction and polynomial chaos acceleration of bayesian inference in inverse problems[END_REF], [START_REF] Raisee | Stochastic model reduction for polynomial chaos expansion using proper orthogonal decomposition[END_REF], Active Subspace methods Bauerheim et al.;[START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF]2014) ⋄ Reduced order modelling approaches: An exceedingly large number of scientific and engineering topics are confronted with the need of high computational resources to study complex, real world phenomena or to solve challenging design problems. Therefore, to overcome the roadblock of the simulation cost, the use of low-order modelling techniques is becoming increasingly popular. Different reduced order modelling techniques are described in this section and their advantages and drawbacks are discussed.

-Surrogate modelling techniques:

Surrogate models are used to generate an accurate approximation of a high-fidelity computational model while minimizing the computational cost. They are generally compact and cheap to evaluate, and they have proved their efficiency in a wide range of topics such as optimization, prototyping or sensitivity analysis. Consequently, in many fields there is great interest in tools and techniques that facilitate the construction of such regression models, while minimizing the computational cost and maximizing model accuracy. Building a good surrogate models is however not straightforward. For that purpose it is necessary to know a priori the physical behaviour of the system and to address the following questions:

1 How to couple the model with the reference simulation code ?

2 Which type of model should be appropriate to approximate the benchmark data (linear, quadratic, cubic etc.)?

3 How to run efficiently surrogate model simulations (locally or in parallel)?

4 Is it possible to estimate easily the model quality and to ensure a real estimation of the model outputs?

5 How to fit the surrogate model and how many samples do we need to collect to achieve this ?

The data collection aspect is worth emphasizing. Since data is computationally expensive to obtain and the optimal data distribution is not known initially, data points should be collected iteratively until covering reasonably the response surface of the high-fidelity model outputs. However, when the complexity of the system increases, the components of the surrogate models increase as well thus complicating the fitting process with reasonable number of samples. For these case, it is preferable to reduce the basis of the complex system at first before investigating on surrogate modelling techniques.

-The Polynomial Chaos technique:

Initially investigated by Norbert Wiener [START_REF] Wiener | The homogeneous chaos[END_REF]) before the advent of computers, the Polynomial Chaos method offers an efficient high-order accurate way of including non-linear effects in stochastic analysis. Several research activities, in a wide variety of topics, have been conducted using Polynomial Chaos technique. For example in CFD [START_REF] Lucor | Adaptive generalized polynomial chaos for nonlinear random oscillators[END_REF], [START_REF] Mathelin | Stochastic approaches to uncertainty quantification in cfd simulations[END_REF]), structural mechanics [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF], [START_REF] Ghanem | Spectral techniques for stochastic finite elements[END_REF]), nuclear engineering and design (Cooling et al. ( 2013)). The Polynomial Chaos technique has many attractive features which are potentially well suited for numerical computations and it is known to be more computationally efficient than the traditional stochastic Monte Carlo simulation. Among the attractive features of the Polynomial Chaos, two of them are very interesting:

1 Polynomial Chaos is a non-sampling method that is used to decompose a random function (or variable) into separate deterministic components. Therefore, the response surface of the model outcomes can be approximated by a sum of orthogonal polynomial series in the random uncertain parameters space.

2 The convergence of the Polynomial Chaos is much more efficient then Monte Carlo sampling method at least for simple geometries.

Following the theory of Polynomial Chaos, any stochastic quantity/equation can be approximated with a finite standard deviation using a truncated expansion. The solution of the stochastic equation can be represented as [START_REF] Wiener | The homogeneous chaos[END_REF]:

R(θ) = +∞ ∑ k=0 β k Ψ k (ξ(θ)) (2.1)
where β k represents the deterministic component e.g. the Polynomial Chaos coefficients of the stochastic equation R, Ψ k is the set of multidimensional polynomials, ξ(θ) is the vector containing the set of independent random variables with the given joint density

ρ(ξ 1 ) = ∑ ρ i (ξ i ).
The family Ψ k (ξ(θ)) satisfies the orthogonality relations:

⟨Ψ k , Ψ l ⟩ = 0 f or k ̸ = l, (2.2)
The property of orthogonality of the polynomial basis Ψ k is a very important characteristic in spectral analysis. It is mathematically expressed through the definition of the following inner product ⟨., .⟩:

⟨Ψ k , Ψ l ⟩ = ∫ Ψ k (ξ)Ψ l (ξ)ρ(ξ)dξ = δ kl ||Ψ k || 2 (2.3)
where δ kl is the Kronecker δ which is equal to 1 for j = k and equal to 0 otherwise and

||Ψ k || 2 = ⟨Ψ k , Ψ k ⟩.
For practical computation, the stochastic quantity R is approximated by a truncated expansion which depend on the number N of independent random variables of the stochastic equation R and the maximum degree of the Polynomials denoted p with respect to the following formula:

P = (N + p)! (N ! p!) (2.4)
When both the number of the polynomial order and the number of random parameters increase, the number of terms in the spectral expansion increases as well.

Now that the stochastic problem R has been replaced by a stochastic system for the Polynomial Chaos coefficients β k , intrusive or non-intrusive approaches can be used to solve the Polynomial Chaos system:

-The intrusive approach [START_REF] Acharjee | A non-intrusive stochastic galerkin approach for modeling uncertainty propagation in deformation processes[END_REF]Zabaras (2007), Tryoen et al. (2010)): This approach is known to be analytically cumbersome because it involves some algebraic manipulation of the underlying governing equations of the polynomial system. Therefore, additional implementation is needed to solve the novel set of equations derived.

The best-known intrusive method to solve Polynomial Chaos system is stochastic spectral Galerkin projection technique.

-The non-intrusive approach [START_REF] Raisee | A non-intrusive model reduction approach for polynomial chaos expansion using proper orthogonal decomposition[END_REF], Le Maître and Knio (2010), Zein et al. ( 2013)): There are two non-intrusive methods to construct the PCE approxi-mation: the projection method and the regression method. Unlike the previous intrusive approach, no modification of the system of equations is needed when using non-intrusive approach. Both of the projection and the regression method are black box methods that require a set of independent simulations for different values of the input parameters. As it was explained by [START_REF] Zein | An efficient samplingmethod for regression-based polynomial chaos expansion[END_REF], the regression method requires the definition of a design of experiments depending on the PCE polynomial function. When using the projection method fo example, the k th Polynomial Chaos coefficient is expressed by projecting the stochastic quantity R onto the polynomial basis in such a way that:

β k = ⟨R, Ψ k ⟩ ⟨Ψ k , Ψ k ⟩ (2.5)
Finally, Eq. (2.5) can be solved numerically with spectral projection and linear regression approach [START_REF] Eldred | Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification[END_REF]).

Uncertainty analysis from the computed Polynomial Chaos coefficients is therefore immediate as the expectation and the variance of the process are given respectively by Eq. (2.6) and Eq. (2.7).

E{R(θ)} = β 0 (2.6) V ar(R(θ)) = E [ (R(θ) -E [R(θ)]) 2 ] = +∞ ∑ k=1 β 2 k ||Ψ k || 2 (2.7)
A comparison between intrusive and non-intrusive Polynomial Chaos technique was investigated in the study of [START_REF] Onorato | Comparison of intrusive and nonintrusive polynomial chaos methods for cfd applications in aeronautics[END_REF] and some sensitivity analysis are performed using Polynomial Chaos technique in the work of [START_REF] Lucor | Sensitivity analysis of large-eddy simulations to subgrid-scale-model parametric uncertainty using polynomial chaos[END_REF] and [START_REF] Crestaux | Polynomial chaos expansion for sensitivity analysis[END_REF].

The cost of solving the Polynomial Chaos system grows at least proportionally to the number of terms in the truncated Polynomial Chaos expansion. Consequently, the method remains difficult to implement for high-dimensional systems and further investigations on this topic are still ongoing [START_REF] Raisee | Stochastic model reduction for polynomial chaos expansion using proper orthogonal decomposition[END_REF], [START_REF] Miranda | Adjoint-based robust optimization using polynomial chaos expansions[END_REF]).

-The Active subspace method:

Extensively discussed in the studies of [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF] and described in Chapter 6, the Active Subspace methodology is an emerging approach used to describe the strong variability of a model output (objective function) along the directions of the input parameters space. In this view, only the dominant one-dimensional subspace of the entire input parameter space is kept for future Uncertainty Quantification investigations. To identify the directions along with the variation of the model outputs is relevant, an eigenvalue decomposition of the gradients of the objective function is realised. Typically, uncentered covariance matrix C of the gradient vector of the model output is used. When considering a scalar function f of a column vector x, the covariance matrix C is expressed as the following:

C = E [ (∇ x f )(∇ x f ) T ] (2.8)
where E is the expectation operator and f the targeted scalar function e.g. the objective function. The elements of C are approximated with a sampling method (commonly a Monte Carlo), by randomly sampling gradient values in the parameter space. The approximated covariance matrix is therefore:

C = 1 M M ∑ i=1 (∇ x f i )(∇ x f i ) T (2.9)
where M is the number of samples, ∇ x f i = ∇ x f (x i ), x i follow a pre-defined distribution (uniform for example). Since this matrix is symmetric, positive, and semidefinite, it admits a real eigenvalue decomposition:

C = W ΛW T , Λ = diag(λ 1 , ..., λ m ), λ 1 ≥ ... ≥ λ m ≥ 0 (2.10)
where W is the eigenvector corresponding to the coefficients of a linear combination of input parameters (W T x) and are the eigenvalues which quantify the effect of the active variable W T x on the model output f(x): the larger λ i is, the more significant the active variable W T x is on the average output response. Consequently, the Active Subspace methodology dissociates the active to inactive subspaces to ease design optimization and surrogate modelling analysis. This method is generally compared to Principal Component Analysis (PCA), also known as Proper Orthogonal Decomposition (POD) but some differences remain between them:

1 PCA is typically used to either reduce the dimension of the output space, or the dimension of an input space that has initially conditioned by specific mathematical processes (pareto-front for instance [START_REF] Lukaczyk | Active subspaces for shape optimization[END_REF])).

2 Active subspace is different in that it reduces the number of input parameters based only on the model outputs and its corresponding gradients. No matrices conditioning is necessary as a first step.

⋄ Sensitivity based methods:

Sensitivity Analysis methodologies are used to quantify independent or correlated effects of input uncertainties and their subsequent impact on the model prediction. Typically, they help to address the following question:

Which of these input parameters have the most influence on the solution estimated from the model prediction?

To answer this question, sensitivity-based methods use the derivative of the model outcomes as a function of the model's input to quantify the ratio of output perturbations over the input perturbations.

The sensitivity derivative of an objective function f with respect to the random variable of y describing the sources of uncertainties is: ∂f ∂y .

The derivative of the objective function f can be assessed by numerical methods such as:

1 Finite difference implementation to calculate f (y 0 ) and f (y 0 + δy), where δ y stands for the perturbations on the input variables.

2 Adjoint based gradient calculation. Adjoint sensitivity analysis of incompressible flows was proposed by [START_REF] Hill | Panel discussion on computational combustion[END_REF] and developed further by [START_REF] Giannetti | Structural sensitivity of the first instability of the cylinder wake[END_REF] in order to reveal the region of the flow that causes a Von-Karman vortex street behind a cylinder. They used adjoint methods to calculate the effect that a small control cylinder has on the growth rate of oscillations, as a function of the control cylinder position downstream of the main cylinder. This control cylinder induces a force in the opposite direction to the velocity field. Gianetti and co-workers considered this feedback only on the perturbed fields but [START_REF] Marquet | Sensitivity analysis and passive control of cylinder flow[END_REF], extended this analysis to consider the cylinder effect on the base flow as well. Adjoint sensitivity analysis was also widely applied by [START_REF] Magri | Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach[END_REF], [START_REF] Magri | A theoretical approach for passive control of thermoacoustic oscillations: Application to ducted flames[END_REF], Magri and Juniper (2013a). They applied adjoint techniques to a time-delayed thermoacoustic system: a Rijke tube containing a hot wire. The idea was to calculate how the growth rate and frequency of small oscillations about a baseline state are affected either by a generic passive control element in the system (the structural sensitivity analysis) or by a generic change to its base state (the base-state sensitivity analysis).

Theoretically, adjoint techniques are described via two different approaches:

-Discrete Adjoint (DA): it operates on the numerically discretized system.

-Continuous Adjoint (CA): it operates on the continuous system as for partial differential equations.

The studies of [START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF] highlighted two new applications of adjoint methods in the study of thermo-acoustic instability. The first one relies on calculating gradients using the Active Subspace method previously presented. The second one relies on calculating the gradients in a non-linear thermo-acoustic Helmholtz solver. The latter task is an objective of this thesis.

Approximating the derivative of the function generally depends on the type of the solver being used. Generally finite difference methods are easier to handle with deterministic solvers because the implementation steps are rather straightforward. When dealing with 3D Finite Elements Methods and parallel solvers for example, the implementation of finite difference methods becomes more complex as the number of operations to achieve increases.

However, finite differences are known to produce inaccurate derivatives. On the contrary, adjoint techniques provide the exact derivative of the model outcomes. This is interesting when dealing with real time applications for instance.

About Uncertainty Quantification in the framework of thermoacoustics

For combustion engineers, a key challenge remains the development of accurate and predictive combustion response models to detect potential combustor instability. Indeed, effective modelling of the flame dynamics will certainly improve the understanding of processes such as nonlinear phenomena responsible for limit-cycle oscillations, the flame-acoustic coupling in industrial geometries, flame-vortices interactions and the interaction of flames with distributed reaction zones or well-stirred reactors. Due to the limited knowledge of all the aforementioned phenomena, introducing Uncertainty Quantification to analyse the probabilistic aspects of the simulation of combustion instabilities is interesting.

Flame models obtained experimentally or numerically are known to be highly dependent on the multiple input parameters whether geometrical or physical. One of the overriding concerns is the ability to address the sensitivity of thermoacoustic results with respect to the flame model input parameters, n and τ towards reliable predictions of unstable modes in gas turbine combustors:

Uncertainty Quantification will help in that sense. on n and τ propagate into uncertainties on the growth rate ω i and to determine the Risk Factor of the acoustic mode e.g. the probability for a mode to be unstable (ω i > 0):

Risk Factor(%) = 100

∫ ∞ 0 P DF (ω i )dω i (2.11)
where P DF (ω i ) stands for the probability density function of the growth rate of the acoustic disturbances. To fairly assess the Risk Factor, it is necessary to have a realistic statistical distribution of the input parameters n and τ , given by experimental data or early numerical results.

Aside from impedance boundary conditions and chamber design away from the flame, performing

Uncertainty Quantification analysis on the flame response parameters n and τ allows to account for uncertainties relevant to combustion chemistry, swirler design, wall heat transfer, inlet temperatures and spray characteristics. All these above mentioned uncertainties are the key elements that maintain the stability inside the combustor. To get the full statistics of the output quantity of interest, one critical issue is to define proper methodologies to propagate uncertainties in the system. Several techniques may be used to handle this task according to the number of input parameters involved. CHAPTER 2. UNCERTAINTY QUANTIFICATION In academic combustors, only one burner is generally present so that the shape and size of the uncertain regions depend only on a few uncertain parameters such as the inlet air temperature, the amplitude and phase of the flame response and the inlet/outlet boundary impedances.

The situation is more complex when dealing with industrial combustion chambers (as presented in Fig. 2.4). Such complex gas turbine engines contain a combustion chamber with an annular shape hosting several injectors as illustrated in Fig. 2.4. In such systems azimuthal thermoacoustic modes appear since the radial and longitudinal directions are shorter than the azimuthal one. Many studies on the effect of the nature of azimuthal modes in combustion chamber have been done with different tools [START_REF] Mensah | Efficient computation of thermoacoustic modes in annular combustion chambers based on bloch-wave theory[END_REF]). Moreover, the number of uncertain parameters may reach several tens since the gain n and time delay τ of each burner (and associated flame) are highly sensitive to manufacturing tolerances. The curse of dimensionality is thus becoming an issue when applying UQ to such systems. Moreover, the coupling between the combustion chamber, the burners and the upstream plenum is also rather complex as revealed by the recent experiment of [START_REF] Worth | Modal dynamics of self-excited azimuthal instabilities in an annular combustion chamber[END_REF], the numerical investigations of [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF] and [START_REF] Bourgouin | A new pattern of instability observed in an annular combustor: The slanted mode[END_REF]. Recent analytical descriptions of thermoacoustic instabilities in annular systems [START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF]), by taking into account burners heterogeneities (Bauerheim et al. (2014a), [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF]) open new perspectives regarding parametric studies and Uncertainty Quantification in these complex systems.

Objective and structure of the study

This thesis is a part of the European Project called UMRIDA (Uncertainty Management for Robust Industrial Design in Aeronautics), which started in October 2013. The objective of UMRIDA is to seek robust design optimization under uncertainties for industrial challenges.

This collaborative project aims at bridging the gap from current state-of-the-art at basic research to a technology readiness level where large numbers of simultaneous uncertainties can be treated in analysis and design. This thesis aims to bring new perspectives to quantify uncertainties in the thermoacoustic modelling of gas turbine combustors.

Objectives of the thesis

Various objectives are targeted in this work:

⋄ Predict combustion instabilities for academic and industrial combustion chambers to determine the frequency of oscillation, growth rate and structure of the eigenmodes.

⋄ Develop and introduce Uncertainty Quantification analysis in the framework of thermoacousticinstabilities to perform robust stability analysis of thermoacoustic systems. The use of Uncertainty Quantification aims at giving a consistent and industrially-realistic support by quantifying the confidence in the modelling of complex systems for risk assessment and decision making. Suitable algorithms are used to propagate uncertainties with respect to the flame model parameters and knowing that Large Eddy Simulation techniques are very CPU demanding, an Helmholtz solver and a quasi analytical tool are preferred for the studies.

Structure of the manuscript

The manuscript is structured in three parts that includes different chapters:

⋄ Part I: The current part is a general introduction on combustion instabilities and Uncertainty Quantification.

⋄ Part II: This part focuses on the study of thermoacoustic instabilities in combustors using low-order modelling techniques.

-Chapter 3 details the assumptions and the governing equations used to describe thermoacoustic instabilities in combustion chambers: from the Navier-stokes equations for a gas mixture to the linearized wave equation. The model used to represent the flame response to acoustic perturbations is also presented. The iterative procedure used to solve the discretized Helmholtz equation in a 3D Helmholtz solver is shown. It enables to provide eigenfrequencies and modal structures of the resonant modes of the system. Additionally, the mathematical framework and the basic concepts for using network modelling techniques to investigate thermoacoustic instabilities in industrial and annular combustors is presented.

-Chapter 4 aims at establishing the connectivity between LES and low-order modelling approaches to identify acoustic eigenmodes in large scale-geometries. The objective is to prepare the groundwork for the development and the application of computationally efficient Uncertainty Quantification approaches for complex industrial systems. to be faster, efficient and accurate in providing all thermoacoustic modes of the system. These tools provide a theoretical interpretation of the results given from Large Eddy Simulations and acoustic solvers. Moreover, such approach allows to ease the system modelling procedure because the interaction between combustion and acoustics can be essentially treated as a zero-dimensional process.

The literature confers numerous reviews and articles dedicated to the use of low-order analysis techniques for the study of thermoacoustic instabilities [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], [START_REF] Munjal | Acoustics of Ducts and Mufflers[END_REF]). Network of acoustic element was investigated by H.J. [START_REF] Merk | An analysis of unstable combustion of premixed gases[END_REF]) to characterize the unstable combustion process of premixed gases. Later, such methodology has been investigated by [START_REF] Bohn | An acoustical model to predict combustion driven oscillations[END_REF] who formalized a thermo-acoustic system into a set of network elements represented by specific transfer matrices. Other investigations on this topic CHAPTER 3. HELMHOLTZ SOLVERS AND NETWORK MODELS have been realised by [START_REF] Dowling | Nonlinear self-excited oscillations of a ducted flame[END_REF] and co-workers by taking care of non-linear effects, entropic waves, boundary conditions, mixture fraction oscillations and force oscillations due to flow instabilities has been also discussed. The implementation of low-order modelling techniques has mostly been realised for simple cases where a single burner is involved. Later on, such methodologies have been applied to annular combustion chambers by [START_REF] Keller | Thermoacoustic oscillations in combustion chambers of gas turbines[END_REF] and co-workers, Evesque

and [START_REF] Evesque | Low-order acoustic modelling for annular combustors: Validation and inclusion of modal coupling[END_REF] or even [START_REF] Kopitz | Thermoacoustic stability analysis of an annular combustion chamber with acoustic low order modeling and validation against experiment[END_REF] with a special care about the boundary conditions to impose in such complex configurations.

Modern gas turbine engines have a ring-shape structure and they are divided in different cavities that comprise a combustion chamber, an upstream air plenum and several injectors, typically from 10 to 25. These kinds of annular systems are widespread in helicopter and aircraft turbines because their design fits efficiently between the axial compressor and the turbine. In such complex systems, a constructive Flame/Acoustic coupling, occurring when heat release and acoustic pressure perturbations satisfy a phase difference relationship, as stated by the Rayleigh criterion (Rayleigh (1878)), favors the apparition of azimuthal acoustic waves. These azimuthal acoustic waves propagating inside the combustor are commonly observed for low fre-quency amplitudes and they represent a major issue for many industrial applications. An effective control of these modes is necessary to ensure the sustainability of modern combustion chambers and to supply the specific energy they require. Consequently, several research activities such as those of [START_REF] Lieuwen | Combustion Instabilities in Gas Turbine Engines. Operational Experience, Fundamental Mechanisms and Modeling[END_REF], [START_REF] Krueger | Prediction and measurement of thermoacoustic improvements in gas turbines with annular combustion systems[END_REF], [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], [START_REF] Leyko | Massively parallel les of azimuthal thermo-acoustic instabilities in annular gas turbines[END_REF] have been dedicated to the study of their structure and their complex physical mechanisms.

Until recently, only few experimental annular chambers have been built to study the physics of azimuthal modes [START_REF] Seume | Application of active combustion instability control to a heavy gas turbine[END_REF], [START_REF] Krebs | Thermoacoustic flame response of swirl flames[END_REF]). These experimental studies were cumbersome for a number of reasons including poor technological supplies to conceive realistic full annular combustors, limited optical access or even sustainable experimental costs. Applications were conducted on simplified and small-scale annular chambers thus leading to drastic modelling assumptions. As a result, that make difficult rigorous validations of experimental observations and theories on annular combustor engines behaviours. More recently, with enhancements of experimental means, the development of realistic laboratory-scale annular combustor has become more affordable and has shed some light on both the emergence and the nature of azimuthal thermoacoustic modes. Typically, they tend to develop as standing, turning (or spinning/mixed modes) or even rotating acoustic modes as it is detailed in Table . 3.1. Turning or spinning modes are characterized by pressure and velocity nodes traveling at the speed of sound whilst standing modes corresponds to fixed pressure nodes and wave modulations. These modes may be also represented as the combination of two waves A + and A -traveling in opposite directions. The ratio of the amplitude of the turning waves A + =A -determines the nature of the corresponding azimuthal mode. Rotating modes [START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF]) can be assimilated to standing modes for which the structure slowly rotates at the azimuthal convective speed. Although these types of modes are encountered in different experimental and numerical simulation studies, they are also observed in real gas turbine engine prototypes. Many non-linear and linear approaches [START_REF] Schuermans | Thermoacoustic modeling and control of multiburner combustion systems[END_REF], [START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF] Recent advances in computer software and hardware allow to combine state-of-the-art technologies and numerical methods to account for the physical processes involved in modern gas turbine combustors. Commonly, massively parallel 3D LES techniques are used to investigate and to control the dynamics of azimuthal thermoacoustic modes developing in annular systems.

Such techniques allow the study of other effects encountered in such complex combustors due to chemical aspects or even limit cycles [START_REF] Fureby | LES of a multi-burner annular gas turbine combustor[END_REF], Bourgouin et al.;[START_REF] Bourgouin | Self-sustained instabilities in an annular combustor coupled by azimuthal and longitudinal acoustic modes[END_REF]014b), Hermeth;[START_REF] Hermeth | Mechanisms affecting the dynamic response of swirled flames in gas turbines[END_REF][START_REF] Salas | Aspects numériques et physiques des instabilités de combustion dans les chambres de combustion annulaires[END_REF]). Because of the prohibitive computation time required by LES approaches, low-order modelling tools are preferred to study azimuthal modes [START_REF] Sensiau | Simulations numériques des instabilités thermoacoustiques dans les chambres de combustion aéronautiques[END_REF], [START_REF] Evesque | Spinning and azimuthally standing acoustic modes in annular combustors[END_REF], [START_REF] Bibliography Stow | A time-domain network model for non-linear thermoacoustic oscillations[END_REF], [START_REF] Morgans | Model-based control of combustion instabilities in annular combustors[END_REF]).

Going beyond computational time constraints of LES by using low order modelling tools allows to investigate fundamentally the pure acoustic of the system, to focus on other interesting mechanisms involved in annular configurations including the influence of transversal flame excitation [START_REF] Guirardo | Azimuthal instabilities in annular combustors: standing and spinning modes[END_REF]) or even the degree of interaction between the system cavities induced by flame response non-linearities [START_REF] Noiray | Investigation of azimuthal staging concepts in annular gas turbines[END_REF]). Moreover, Helmholtz solvers adapted to annular systems [START_REF] Benoit | Prédiction des instabilités thermo-acoustiques dans les turbines à gaz[END_REF], [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF], [START_REF] Sensiau | Simulations numériques des instabilités thermoacoustiques dans les chambres de combustion aéronautiques[END_REF]), are good candidates in predicting such annular combustor instabilities. However, the computation of such systems using Helmholtz solver may become expensive even if the solver is parallelized. Moreover, due to some difficulties in extracting phenomenological conclusions from Helmholtz solver computation, analytical network modelling techniques may be used to study physical processes involved in annular systems. Although providing theoretical interpretations of given solutions from Helmholtz solvers, network modelling techniques provide sustainable speed up of azimuthal mode computations. For Uncertainty Quantification purposes for which several runs could be required, the use of such techniques is very appealing.

In this chapter, the focus is on the study of azimuthal modes. To avoid expensive computation costs linked to LES techniques, the use of Helmholtz solvers and network modelling tools is preferred to investigate the stability and the control of azimuthal modes. Assuming harmonic time dependence, e -iωt and linear acoustics, mathematical/numerical models whose unknown is the (Fourier transformed) acoustic pressure p distribution over space can be derived. Two such models will be employed in this work:

1 A 3D Helmholtz solver called AVSP developed by CERFACS is used to account for all modes nature and complex geometry features [START_REF] Benoit | Prédiction des instabilités thermo-acoustiques dans les turbines à gaz[END_REF] This chapter will be organized as follows:

⋄ In section 3.1, the physical model used to represent thermoacoustic instabilities in combustors is presented. Initially, the derivation of the approximated linear wave equation for the small perturbations in reactive flows is performed. Then, the flame model, based on n-τ formalism, that is used to account for the coupling between acoustics and combustion is described. Once the Helmholtz equation is constructed, it is discretized on unstructured meshes, using a finite volume methodology. The latter leads to a complex nonlinear eigenvalue problem that is solved iteratively in the AVSP solver.

⋄ In section 3.2, the analytical theory used to study only azimuthal modes in annular systems is described. This analytical theory is based on a quasi-one-dimensional zero-Mach number formulation where many burners are connected to an upstream annular plenum and a downstream chamber. As for the 3D acoustic solver AVSP, the flame response is modeled using the n-τ formalism and is supposed to be compact. A methodology called Annular Network Reduction (ANR) is used to capture only azimuthal waves in the annular cavity network. The set of equations that results from this methodology allows to solve numerically a simple dispersion relation that furnishes a fair estimation of the frequency and the growth rate of all azimuthal modes of the combustors. This methodology is also useful to analyse other mechanisms as transverse forcing effects, symmetry breaking and mode nature.

Both tools used to study azimuthal modes appearing in annular combustors are complementary: the 3D Helmholtz solver AVSP provides qualitative interpretation of the behavior of the system while the analytical tool provides a theoretical interpretation of the results of AVSP solver.

Within a framework of Uncertainty Quantification analysis, the low order tool ATACAMAC has the advantage to be cheaper in CPU time than AVSP code besides furnishing quickly the risk associated to an azimuthal mode of the system to become unstable.

Thermoacoustic analysis using a Finite Volume Based Helmholtz Solver

The mechanisms of thermoacoustic instabilities is very complex due to the coupled interactions of acoustics waves and heat release fluctuations. Furthermore, the inherent non-linearities associated with the turbulent flow or chemical reactions can make the study of instabilities more complicated.

To analyze thermoacoustic instabilities, many simplifications are made to render the problem tractable:

weight and heat capacites.

2 The flame is modeled as a pure acoustic element.

3 Volume forces are neglected (the gravity for example).

4 Viscous effects are neglected.

As the validity of the assumptions used to study thermoacoustic instabilities are also case dependent, a well suited model is chosen to represent the coupling of heat release and acoustic wave propagations. In the framework of linear acoustics and under the assumptions cited above, Navier-Stokes equations [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]) can be manipulated to construct the wave equation for reactive flows that takes into account the interaction between the flame and the acoustic waves.

Mathematical formulation

The Euler equation for a gas mixture under the assumptions pre-cited in the above Section 3.1 reads [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Dρ Dt = -ρ∇.u, ρ Du Dt = -∇p, Ds Dt = rq p .
(3.1)

The system of Eq. (3.1) corresponds respectively to the equations of mass, momentum and entropy for a compressible inviscid flow (in absence of external forces). The parameters used in Eq. (3.1) are presented in The acoustic field is generally decomposed in terms of small amplitude perturbations that are superimposed on the mean flow field. When injecting this decomposition in the set of equations (3.1), and by keeping only first order terms, we get a set of linearized equations fitted by a specific term that accounts for the flame/acoustic interaction.

- r = C p -C v T Temperature [K] s Entropy [J/K]

The linear wave equation for reactive flows

Considering the simple case of large scale small amplitude fluctuations superimposed to a zero Mach number (u 0 ≈ 0) mean flow which depends only on space, the set of equations Eq. 3.1 can be decomposed in mean value (index 0) and low fluctuations (index 1). The zero Mach number assumption is valid as soon as the characteristic Mach number M = √ u 0 .u 0 /c 0 is small compared to the ratio between the thickness of the reaction zone and the typical acoustic wavelength λ 1535 [START_REF] Truffin | Comparison and extension of methods for acoustic identification of burners[END_REF], Poinsot and Veynante ( 2011)). In this case, ∇p 0 = 0 and q 0 = 0 and D Dt ≪ ∂ ∂t holds for any fluctuating quantity because, with u 0 ≈ 0, the non-linear convective terms are always of second order.

The instantaneous pressure, density, temperature, entropy, and velocity fields can then be CHAPTER 3. HELMHOLTZ SOLVERS AND NETWORK MODELS written as:

1540 p(⃗ x, t) = p 0 (⃗ x) + p 1 (⃗ x, t), ρ(⃗ x, t) = ρ 0 (⃗ x) + ρ 1 (⃗ x, t), s = s 0 (⃗ x) + s 1 (⃗ x, t), ⃗ u(⃗ x, t) = ⃗ u 0 (⃗ x) + ⃗ u 1 (⃗ x, t). (3.2)
Note that the quantities:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ p 1 (⃗ x,t) p 0 (⃗ x) , ρ 1 (⃗ x,t) ρ 0 (⃗ x) , s 1 (⃗ x,t) s 0 (⃗ x) , √ ⃗ u 1 (⃗ x, t) • ⃗ u 1 (⃗ x, t)/c 0 (⃗ x).
(3.3) are of order ϵ, where ϵ ≪ 1 and c 0 (⃗ x) = √ γp 0 (⃗ x)/ρ 0 (⃗ x) is the mean speed of sound. From the above, the set of linear equations for the fluctuating quantities ρ 1 (⃗ x, t), ⃗ u 1 (⃗ x, t) and p 1 (⃗ x, t), keeping only first order terms, reads :

∂ρ 1 (⃗ x, t) ∂t + ⃗ u 1 (⃗ x, t).∇ρ 0 (⃗ x) + ρ 0 (⃗ x)∇.⃗ u 1 (⃗ x, t) = 0, (3.4) ρ 0 (⃗ x) ∂⃗ u 1 (⃗ x, t) ∂t + ∇p 1 (⃗ x, t) = 0, (3.5) ∂s 1 (⃗ x, t) ∂t + ⃗ u 1 (⃗ x, t).∇s 0 (⃗ x) = rq 1 (⃗ x, t) p 0 (⃗ x) . (3.6)
Using the 2 nd Principle of thermodynamics, the entropy equation can be written as:

1545 Ds Dt = C v p Dp Dt - C p ρ Dρ Dt (3.7)
As the mean flow quantities are not time dependent, the mean entropy gradient reads:

∇s 0 = C v p 0 ∇p 0 - C ρ ρ 0 ∇ρ 0 (3.8)
As the flow is assumed to be at rest, the mean pressure gradient is equal to zero. Thus the entropy gradient, Eq. (3.8), becomes:

∇s 0 = - C ρ ρ 0 ∇ρ 0 (3.9)
When substracting Eq. (3.4) and Eq. (3.6), the following simplified system of equations is built:

1 γp 0 ∂p 1 ∂t + ∇.⃗ u 1 = 1 C v r γp 0 q 1 , (3.10) ∂⃗ u 1 ∂t + 1 ρ 0 ∇p 1 = 0. (3.11)
where q 1 stands for the fluctuating part of the heat release.

Taking the time derivative of Eq. (3.10) and adding the divergence of Eq. (3.11) allows finally to establish the linear wave equation for p 1 that describes the propagation of pressure fluctuations:

1 γ(⃗ x)p 0 ∂ 2 p 1 (⃗ x, t) ∂t 2 -∇ • ( 1 ρ 0 (⃗ x) ∇p 1 (⃗ x, t) ) = 1 C v γ(⃗ x) -1 γ(⃗ x)p 0 ∂q 1 (⃗ x, t) ∂t (3.12)
In Eq. (3.12), the left hand side term corresponds to a classic wave equation while the right hand side term takes into account the flame response to acoustic perturbations. However, the quantity ρ 0 (⃗ x) is not constant in space and it must be kept within the divergence because it accounts for temperature variations of the combustion process.

Therefore, the wave equation reads:

γ(⃗ x)p 0 (⃗ x)∇ • ( 1 ρ 0 (⃗ x) ∇p 1 (⃗ x, t) ) - ∂ 2 p 1 (⃗ x, t) ∂t 2 = -(γ(⃗ x) -1) ∂q 1 (⃗ x, t) ∂t (3.13)
It then proves useful to introduce harmonic variations in Eq. 3.13 in such a way that: (-iωrt) ] ,

p 1 (⃗ x, t) = e (ω i t) ℜ [ p(⃗ x)e (-iωrt) ] , (3.14) u 1 (⃗ x, t) = e (ω i t) ℜ [ û(⃗ x)e
(3.15)

q 1 (⃗ x, t) = e (ω i t) ℜ [ q(⃗ x)e (-iωrt) ] , (3.16)
where ω stands for the complex valued pulsation and is divided in two parts:

• ω r = ℜ(ω) = 2πf r , the frequency of oscillation (Hz),

• ω i = ℑ(ω) = 2πf i the growth rate of the acoustic pressure disturbances (s -1 ),

where

ω = ℜ(ω) + iℑ(ω) = ω r + iω i .
The combination of equations, Eq. 3.12, Eq. 3.14 and Eq. 3.16, yields the Helmholtz equation for the acoustic pressure disturbance which reads :

γ(⃗ x)p 0 (⃗ x)∇ • ( 1 ρ 0 (⃗ x) ∇p(⃗ x) ) + ω 2 p(⃗ x) = iω(γ(⃗ x) -1)q(⃗ x) (3.17)
In this equation the unknowns are p(⃗ x) the complex amplitude of the pressure disturbance, as well as the complex valued pulsation ω. Quantities ρ 0 and γ depend on the space x coordinates and must be provided as inputs. Modelling the right hand side of the equation Eq. (3.17) is the most difficult part when predicting thermoacoustic instabilities. In fact, this term is associated to the unsteady flame behavior and a well suited model must be used to express the unsteady heat release q(⃗ x).

Modelling of thermoacoustic instabilities using the Flame Transfer Function formulation

Several approaches have been proposed to predict resonant modes between acoustics and combustion [START_REF] Crighton | Modern methods in analytical acoustics[END_REF], [START_REF] Culick | Some recent results for nonlinear acoustic in combustion chambers[END_REF], [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationnary computational fluid dynamics[END_REF], [START_REF] Sattelmayer | Influence of the combustor aerodynamics on combustion instabilities from equivalence ratio fluctuations[END_REF], [START_REF] Selle | Compressible large-eddy simulation of turbulent combustion in complex geometry on unstructured meshes[END_REF]). Some other studies have been devoted to the description of the response of conical or V-shape premixed flames accounting for various phenomena such as stretching effects Generally, the flame response is characterized by its Flame Transfer Function which is defined as a linear relationship between incoming acoustic velocity fluctuations (generally located upstream of the flame front as it was discussed by [START_REF] Truffin | Comparison and extension of methods for acoustic identification of burners[END_REF] or [START_REF] Ducruix | Combustion dynamics and instabilities: elementary coupling and driving mechanisms[END_REF] and harmonic heat release rate perturbations. This idea was first introduced by [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors[END_REF] for compact flames, referred to as the n -τ formalism. The Flame Transfer Function is expressed as the ratio between the global heat released from the flame Q at time t to the time lagged acoustic velocity û measured in the cold gas region upstream of the flame front:

Q 1 (t) = ∫ V q 1 (t)dV = S ref γp 0 γ -1 × n × ⃗ u 1 (⃗ x ref , t -τ ). (3.18)
In Eq. (3.18), Q 1 (t) is the heat release integrated over the flow domain V , S ref is the cross section area of the burner mouth (see Fig. 3.2): In the frequency domain, the Flame Transfer Function becomes:

S ref = V f × δ f ,
Q1 = ∫ Ω q(⃗ x)dΩ = S ref γp 0 γ -1 × n × ⃗ û1 • ⃗ n ref e iωτ .
(3.19)

R = ∫ t ∫ Ω p 1 q 1 dΩ dt (3.20)
The classical Rayleigh criterion stipulates that Flame/Acoustics coupling induces the appearance of instabilities when R > 0 showing the importance of the parameter τ in the description and prediction of thermo-acoustic instabilities.

Using the global flame response modelling is convenient when the typical length of the flame region is small compared to the characteristic acoustic wavelength e.g. it is suitable for acoustically compact flames only. This condition is difficult to reach for industrial combustors but in experimental or analytical models using a global flame response (Eq. (3.18)) is more convenient.

Otherwise, it is possible to use the local flame response formulation [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF] to link the unsteady heat release emitted by the flame at time t to the acoustic velocity at an upstream reference point ⃗ x ref at an earlier time t -τ . In this case, heat release fluctuations are expressed by the following formula:

q 1 (⃗ x, t) q tot = n local (⃗ x) ⃗ u 1 [⃗ x ref , t -τ local (⃗ x)].⃗ n ref U bulk . (3.21)
where q tot stands for the total heat release and U bulk the bulk velocity. The parameter n(⃗ x)

has no dimension due to the scaling by q tot and U bulk . In the frequency domain heat release fluctuations are expressed as:

q1 (⃗ x) q tot = n local (⃗ x) û1 (⃗ x ref ) • ⃗ n ref U bulk e iωτ local (⃗ x) (3.22) (
∫ V f n local e iωτ local dV = U bulk q tot S ref γp 0 γ -1 ne iωτ . (3.23)
where V f is the flame volume. Therefore the connection between the local flame formulation to the 1D flame formulation of Crocco is done following the formula Eq. (3.24):

n local = γp 0 (γ -1) U bulk q tot S ref V f
× n and τ local = τ.

(3.24)

The three-dimensional finite volume based acoustic solver AVSP

To solve Eq. 3.17, it is necessary to provide at first ρ 0 (⃗ x), γ(⃗ x) and the fields of the Flame Transfer Function parameters n and τ . In this work, these data are extracted in two different ways:

⋄ from an experimental combustion chamber [START_REF] Palies | Dynamique et instabilités de combustion de flammes swirlées[END_REF]).

⋄ from LES computations [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF]).

The boundary conditions hereafter described can be used to solve Eq. (3.17) with AVSP:

⋄ Dirichlet type boundary condition:

p = 0. (3.25)
This corresponds to fully reflecting boundary conditions at the outlets.

⋄ Homogeneous Neumann type boundary condition:

∇p • ⃗ n = 0, (3.26)
where ⃗ n is the wall's normal vector. This boundary condition corresponds to fully rigid walls or reflecting inlets.

⋄ Robin type boundary condition:

∇p • ⃗ n -i ω c 0 (⃗ x)Z(ω) p = 0, (3.27)
where Z = p ρ 0 c 0 û•⃗ n is the local reduced complex impedance (generally extracted from LES computations) and c 0 the mean sound speed.

Once the sound speed and the flame parameters fields are provided, the Helmholtz equation is discretized using a finite volume formulation on unstructured tetrahedral meshes thus leading to a nonlinear complex eigenvalue problem. Therefore, Eq. (3.17) is turned into the following matrix form [START_REF] Sensiau | Simulations numériques des instabilités thermoacoustiques dans les chambres de combustion aéronautiques[END_REF], Salas (2013)):

A p + B(ω)p + ω 2 p = F (ω)p, (3.28)
where A is the matrix containing the discretization of the operator ∇• is the structure of the corresponding acoustic mode.

The complex nonlinear eigenvalue problem Eq. 3.28 is then solved in a 3D acoustic solver called AVSP developed at CERFACS. AVSP is based on a finite volume methodology and it is used to fully discretize all the geometrical features of the combustion chamber. It solves, in the frequency domain, the discretized formulation of the Helmholtz equation Eq. 3.17 by assuming harmonic variations at frequency f = ω 2π for the velocity (Eq. (3.15)), the pressure (Eq. (3.14)) and the local heat release fluctuations (Eq. (3.16)).

In the AVSP solver, Eq. 3.28 can be either solved in a steady flame regime or an active flame regime:

1 Steady flame regime: In this case, the unsteady flame response is neglected; the right hand side term of Eq. (3.17) is set to zero so as F p in Eq. (3.28). Consequently, the problem is drastically simplified into an eigenvalue problem depending only on the complex valued pulsation ω. From a physical point of view, steady flame computations are performed to get an idea of the natural acoustic modes in the combustion chamber. Under the assumption that the unsteady flame response acts as a small perturbation of the modes without combustion, a linear expansion technique can be developed to assess the imaginary part of ω and hence the stability of the perturbed modes [START_REF] Mcmanus | A review of active control of combustion instabilities[END_REF], [START_REF] Sensiau | A comparison of solvers for quadratic eigenvalue problems from combustion[END_REF]).

CHAPTER 3. HELMHOLTZ SOLVERS AND NETWORK MODELS

2 Active flame regime: In this case, the unsteady flame response is not neglected and may lead to significant changes of the frequencies inside the combustor. Therefore, an iterative 1675 process based on a fixed point strategy [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF], [START_REF] Sensiau | Simulations numériques des instabilités thermoacoustiques dans les chambres de combustion aéronautiques[END_REF], [START_REF] Salas | Aspects numériques et physiques des instabilités de combustion dans les chambres de combustion annulaires[END_REF]) is used to solve iteratively, the non-linear eigenvalue problem of Eq. (3.28). This iterative procedure is used to solve the following discretized eigenvalue problem:

A p + B(ω + k )p + ω +2 k p = F (ω + k )p (3.30)
where ω + k is the output solution of the problem. This algorithm is sketched in Fig. 3.4 and it can be summarized by the following relation:

1680 ω k+1 = αω + k + (1 -α)ω k , (3.31)
The set of parameters in Eq. (3.31) are detailed in Table . 3.3. ω k → Input of the computation at the 0 th iteration (k=0) provided from the resolution of Eq. (3.28) in the passive flame regime:

ω 0 = ω k ω + k → Output solution of the eigenproblem α →
The relaxation coefficient which is used to smooth the iteration process in case of convergence problems.

Table 3.3: Definition of the parameters of Eq. (3.31) that represents the fixed point algorithm [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]).

As studied by Miguel-Brebion ( 2017), the relax parameter can be fixed or imposed dynamically to optimize the convergence process.

The fixed point algorithm is repeated until the successive solutions of the sequence of linear eigenproblems converge to the sought nonlinear eigenvalue ω. This arises when |ω k -ω + k | < ϵ, where ϵ is the prescribed tolerance. Overall, the convergence of the iterative process depends on the complexity of the system being studied.

Analytical description of thermoacoustic instabilities in annular combustors with network modelling techniques

Network models allow the study of annular configurations as a network of interconnected acoustic elements (chamber, plenum, flame tube, nozzle for example) communicating by means of jump conditions [START_REF] Schuermans | Thermoacoustic modeling and control of multiburner combustion systems[END_REF]) or scattering matrices. The coupling relations for the unknowns across an element are combined into the transfer or scattering matrix of the element.

The transfer matrix coefficients of all network elements are combined to form the complete matrix of the network that can be solved by hand or numerically [START_REF] Polifke | Determination of thermo-acoustic transfer matrices by experiment and computational fluid dynamics[END_REF] When comparing such analytical results to those given by the full 3D Helmholtz solver AVSP, a very good agreement is found in terms of frequencies and growth rate of acoustic modes of the system. Such a methodology opens the path to predict and control azimuthal modes in annular acoustic systems using a fully analytical approach. However the BC type configuration does not fully reflect realistic and modern annular combustors that are linked not only to an annular chamber but also an upstream plenum (see Fig. In this thesis as the main focus is about performing Uncertainty Quantification of thermoacoustic instabilities developing in realistic combustion chamber affordably, the analytical approach of Bauerheim et al. (2014a), [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF] for PBC type of configuration will be mostly used.

Theoretical description

The basic aspects of the analytical model of Bauerheim et al. (2014a), [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF] are briefly described in this Section. More details on the theoretical developments are available in Bauerheim et al. (2014a), [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF], [START_REF] Salas | Aspects numériques et physiques des instabilités de combustion dans les chambres de combustion annulaires[END_REF]. This model is based on an Annular Network Reduction (ANR) methodology that allows to simplify the system complexity by solving an analytical dispersion relation which is implicit and non-linear in the frequency domain. Therefore, this equation may be solved either analytically (under additional assumptions) or numerically and its solutions provide the complex angular frequency ω = ω r + ω i . When the imaginary part of the angular frequency ω i is positive (ω i > 0), the mode is unstable and conversely when the imaginary part of the angular frequency ω i is negative (ω i < 0) the mode is in a stable regime. The ANR methodology allows to recast the system cavities into independent where the value of c 0 depends on the location (c 0 u in the burners and plenum, but c 0 b in the chamber, Fig. 3.7). Thus, using Eq. (3.32), the azimuthal propagation in the i th sector of the annular plenum and chamber can be combined to form a propagation matrix R i (ω) such that:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s i+1 ) = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ e jku2Lp/N 0 0 0 0 e -jku2Lp/N 0 0 0 0 e jk b 2Lc/N 0 0 0 0 e -jk b 2Lc/N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s i ) = [R i ] ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s i ) (3.33)
where w p stands for the acoustic wave propagating in the plenum, w c is the acoustic wave propagating in the chamber. In Eq. 3. Each of the N sectors of annular plenum is linked to the annular chamber through a burner.

Therefore, the interaction between the i th burner and the annular chamber is characterized model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors[END_REF], [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF]) for the unsteady heat release Q ′ produced by the flame (Q ′ = n i e jωτ i u ′ , where n i and τ i are the gain and the time-delay for the i th Flame Transfer Function), an interaction matrix [T i ] is deduced by Bauerheim et al. (2014a). It relates acoustic quantities before the i th junction (coordinate s - i ) to the ones after the junction (s + i ):

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s + i ) = [P ] -1 ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 Γ i,1 1 Γ i,2 0 0 0 1 0 Γ i,3 0 Γ i,4 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ [P ] ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s - i ) = [T i ] ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (s - i ) (3.34)
where [P ] is the matrix relating the Riemann invariants w ± to the acoustic pressure and velocity, and Γ i,k=1..4 are the coupling parameters derived by Bauerheim et al. (2014a):

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Γ i,1 = -S i 2Sp cotan(k u L i ) Γ i,2 = S i 2Sp 1 sin(kuL i ) Γ i,3 = S i 2Sc ρ 0 b c 0 b ρ 0 u c 0 u 1+ne jωτ sin(kuL i ) Γ i,4 = -S i 2Sc ρ 0 b c 0 b ρ 0 u c 0 u (1 + ne jωτ )cotan(k u L i ) (3.35)
where L i is the i-th burner length and S i its cross Section. These coupling parameters have been deduced by assuming that the flames are located exactly at the burner/chamber junction. This location plays a crucial role for plenum modes. These coupling parameters are also obtained in longitudinal configurations [START_REF] Schuller | Acoustic decoupling of longitudinal modes in generic combustion systems[END_REF]) and characterize how cavities are coupled and interact (Fig. 3.8). Decoupling can be achieved using a large section change at the burner junction, but it can be also affected by the flame itself (e.g., by n i and τ i ). Note that if Γ i,1 = Γ i,2 = 0 for all junctions i = 1..N , then the annular plenum is disconnected from the rest of the system. Using the propagation R i (ω)and interaction matrices [T i ] to connect the annular sectors, the annular periodicity leads to the equation governing the acoustic modes behavior in the annular plenum and chamber:

( N ∏ i=1 [R i ][T i ] ) ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ w + p w - p w + c w - c ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
(3.36) Equation (3.36) has non-trivial solutions if and only if the determinant is null, which yields the dispersion relation to be solved: det ⋄ Non-symmetric annular combustors: in this case all sectors and flames are different.

( N ∏ i=1 [R i ][T i ] -I d ) = 0 (3.
The coupling parameters Γ i may differ from a burner to another. Contrary to axisymmetric annular combustors, an implicit analytical dispersion relation for the pulsation ω should be derived as performed by Bauerheim et al. (2014a).

Despite this apparent simplicity, annular configurations containing a chamber and a plenum can exhibit complex lock-in and veering phenomena, for which the active flames are a key ingredient.

1 Under the null coupling assumption or fully decoupled case;

In this case all coupling parameters are zero, Γ 1 = Γ 2 = Γ 3 = Γ 4 = 0. Consequently, the plenum and the chamber are fully decoupled from the burners and flames. As a result, eigenfrequencies are f 0 p = pc 0 u /2L p (pure azimuthal decoupled mode in the plenum) or f 0 c = pc 0 b /2L c (pure azimuthal decoupled mode in the chamber). Since the fresh mixture and hot gases have different temperatures, and the half-perimeter of the plenum and chamber are different, eigenmodes in the plenum and chamber are typically distinct.

f c = pc 0 b 2L c - c 0 N Γ 0 4 4πL c and f p = pc 0 u 2L p - c 0 b N Γ 0 1 4πL p (3.39)
where Γ 0 1 (respectively Γ 0 4 ) is the value of the coupling parameter Γ 1 (respectively Γ 4 ) at the frequency f = f 0 p (respectively f = f 0 c ): these modes are called weakly coupled.

3 Under strong coupling assumption;

The two annular cavities can couple and oscillate at the same frequency, even if f 0 p and f 0 c do not match: the burners and flames tune one of the two cavities so that they can both resonate. In this case, the acoustic mode cannot be identified strictly to belong either to the annular plenum or the annular chamber because the whole combustor is resonating.

The mathematical framework described previously will be applied on two perfectly axisymmetric annular combustion chambers typical of helicopter engines. The objective is then to use suitable probabilistic Uncertainty Quantification approaches to investigate uncertainties related to the Flame Transfer Function in high dimensional systems using the analytical network modelling tool ATACAMAC. It will contribute to determine the Risk Factor of the predominant azimuthal mode of the system namely its probability to be unstable. To minimize the computational cost, only the 3D reactive LES of a single sector is performed to achieve these tasks.

2 Helmholtz computations with AVSP solver: Once the local mean flame fields and sound speed have been extracted, they are used as inputs for AVSP to solve Eq. (3.17).

3 Network modelling with ATACAMAC tool: Solutions of Helmholtz computations are taken as reference to fit the quasi-analytical tool ATACAMAC to push further the thermoacoustic analysis of the system by bringing phenomenological interpretations of the combustor dynamics. However, ATACAMAC requires at first a geometrical fitting of the full-scale combustor. Direct geometrical adjustments of the combustor limit the predictive character of such analytical tool and this is the reason why it is recommended to fit them to 3D results obtained with AVSP by accounting for the whole complexity of the combustion chamber. A good calibration of the network model ATACAMAC with respect to LES and Helmholtz solutions will provide substantial speedups for thermoacoustic calculations and an appealing perspective for Uncertainty Quantification analysis. In other words, ATACA-MAC can be seen as a surrogate model for LES or Helmholtz solvers which then allows to perform UQ studies. The procedure described in Fig. 4.1 is applied on two multi-burner combustion chambers typical of industrial helicopter engines. Safran Helicopter Engines provided the two combustors within the European project UMRIDA.

Description of the 1 st annular combustor of interest with N burners

The industrial system investigated in this section corresponds to a reverse full annular helicopter combustion chamber composed of N circumferentially arranged and identical burners.

The schematic view of the single sector used in the Large Eddy Simulation is shown in Fig. 4.2a

and the full annular engine used for Helmholtz solver computations is presented in Fig. 4.2b.

Each sector of the annular system features an upstream casing where the airflow coming from 1860 the compressor is injected and a downstream combustion chamber where the combustion process takes place. In the primary zone of the combustion chamber, fuel is injected through the swirler and the cooling of burnt gases as well as the thermal protection of the combustion walls are ensured by multi-perforated plates and dilution holes [START_REF] Mendez | Large-eddy simulation of a bi-periodic turbulent flow with effusion[END_REF], [START_REF] Lahbib | Aerothermal modeling of mutiperforated plate cooling for Large Eddy Simulations[END_REF]). The combustion chamber ends with a choked nozzle that is used to release burnt gases and to conserve the sonic state of the stator.

The LES of the single sector and the full annular configuration have been initially performed in the work of [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF], [START_REF] Wolf | Azimuthal instabilities in annular combustion chambers[END_REF] with the LES solver AVBP developed at CERFACS. AVBP is a hybrid (structured/unstructured) and compressible solver that includes chemical aspects and variable heat capacities used to solve the Navier-Stokes equations for reactive flows. It relies on centered-spatial schemes and explicit time-advancement that allow proper control of the numerical-dissipations/filter to accurately resolve all relevant multi-scale of complex industrial systems and acoustic effects [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF]). On top of studying the dynamics of the flow inside the engine, the single sector pulsated LES has been performed in order to extract the input parameters γ(⃗ x), the mean density ρ 0 (⃗ x) and the mean sound speed c 0 (⃗ x) and the TURBINE COMBUSTION CHAMBERS Flame Transfer Function parameters n and τ to account for the Flame/Acoustic interactions.

Afterwards, these fields are injected in the Helmholtz solver AVSP to perform thermoacoustic calculations in the 360 degrees configuration. The objective is to determine azimuthal thermoacoustic modes that are prone to develop in such annular system. Such a study allows to construct the stability map for the combustion chamber and to analyse deeply the mode structure.

In this work, the focus is mainly on the first azimuthal mode of the system, initially identified by Wolf et al. (2012a) and [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF], as an unstable standing mode that slowly rotates at convective velocity controlled by the mean swirl velocity. He shown that a reduction of the FTF delay combined with modification in the chemistry would overcome the unstable effects. However it should be interesting to ensure these conclusions by providing at least a quantitative estimation on the risk of this first azimuthal thermoacoustic mode to become unstable. To do so, its Risk

Factor will be computed by following an Uncertainty Quantification methodology adapted to multi-burner systems. Performing such Uncertainty Quantification studies is a highly challenging undertaking in terms of input models uncertainties and computation resources. The idea in this section is to establish a methodology that allows to considerably speed-up thermoacoustic mode computations using analytical modelling techniques in view of the non-negligible computational time required by Helmholtz solvers (hours of computation using 64 cores) and the prohibitive CPU time required by LES techniques (3,000,000 CPU hours to compute only 30ms physical time on 4,096 cores).

Having only a limited information on the progress of the thermoacoustic simulations performed by Wolf and co-workers, the global thermoacoustic analysis of the system is completely re-done in this work. Two types of simulations were conducted in this study to classify all thermoacoustic modes of the combustor: passive flame and active flame computations. The structure of the system indicates longitudinal and azimuthal waves propagating inside the different cavities.

However, the Helmholtz solver AVSP, which is used for the thermoacoustic analysis, does not allow to identify clearly in which zone acoustic modes belong or to provide the coupling degree between all acoustic cavities of the combustor. Therefore, to push further the acoustic analysis of the system, the analytical model of Bauerheim et al. (2014a) and [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF] described in Section 3.2.1 is used.

Acoustic mode computations of the annular system with N injectors with Helmholtz solver

Thermoacoustic mode computations of the system are realised in this study using a mesh composed of 69019 nodes and 336135 cells for the single sector system and 1010370 nodes and 5 042025 cells for the full annular configuration. As shown in Table 4.1, the mesh size for the Helmholtz computations is drastically reduced compared to that used for Large Eddy Simulations. When performing LES computations, the quality and type of mesh used for the discretization of the Navier-Stokes equations in the computational volume play a crucial role on both accuracy (in term of solutions) and CPU cost. Moreover, to enable a good resolution of the flame front, a well resolved LES mesh is mandatory to sufficiently capture the flame changes dynamic as discussed in early works of [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] and [START_REF] Selle | Joint use of compressible large-eddy simulation and Helmholtz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF]. Such a grid resolution is TURBINE COMBUSTION CHAMBERS not mandatory to capture pure acoustic eigenmodes of the combustor with Helmholtz solver and spectral analysis methodologies would help to define the number of cells required in this case.

To perform thermoacoustic calculations with AVSP solver, it is necessary to provide at first γ(x), the mean density ρ 0 (x), the mean sound speed field c 0 (x) and the local fields of the flame parameters n local (x) and τ local (x). A constant adiabatic coefficient γ and identical sectors and flames (the system is considered to be axisymmetric) are considered for the thermoacoustic analysis. These inputs come from the time-averaged reactive compressible Large Eddy Simulations of a single sector for the operating conditions presented in Table 4. Longitudinal modes are found and dual frequencies correspond to degenerate azimuthal modes that are typical to industrial combustors [START_REF] Lieuwen | Combustion Instabilities in Gas Turbine Engines. Operational Experience, Fundamental Mechanisms and Modeling[END_REF]). The growth rate of each thermoacoustic modes of the combustor is null (ω i = 0.0[s -1 ]) because there is no flame response suggests a coupling activity exists between the cavities of the combustion chamber. This means that the acoustic pressure developing inside the casing and the combustion chamber are both linked to the axial distance and the radial coordinate. To go further, it is interesting to seek how the system would evolve when accounting for the flame effects. a way to proceed allows to ease the exploitation of thermoacoustic solutions and to assess any potential changes in the system response when the flame parameters are perturbed. The global interaction index n and time delay τ injected in AVSP are: n=1486.43[J/m] and τ = 9.87×10 -4 s.

The corresponding global Crocco's values where also determined : n = 6.57 and τ = 9.87×10 -4 s.

Computing the full annular system requires a proper definition of a reference upstream position

x ref , in each sector, to relate the local unsteady heat release to the complete acoustic field. Generally, this point is located at few millimeters upstream the burner mouth in the cold gas area [START_REF] Truffin | Comparison and extension of methods for acoustic identification of burners[END_REF]).

Under the above operating conditions, numerical simulations were conducted by imposing a homogeneous Neumann boundary condition on all walls, inlets and outlets. The first nine eigen-frequencies computed in active flame regime are listed in ⋄ When ℑ(ω) is below 0 the mode is stable and when ℑ(ω) is above 0 the mode is unstable.

Eigenmodes of the system shift from the stable to unstable regime for a value of τ equal to τ = τ 0 = 8.8367 × 10 -4 s approximately equal to a half of the period T = 1 Therefore, when the rotational symmetry of the system is not conserved, two azimuthal counterrotating eigen-pairs appear as it is the case for this annular system. This explains why very close azimuthal modes are computed and remain different in terms of structure.

An observation of the entire pressure field of the system shown that the acoustic activity is present between the chamber and the upper front of the plenum. At this step, clearly stating on the nature of azimuthal mode and being able to quantify rigorously the coupling phenomena of each part remains difficult with the Helmholtz solver. At this point, analytical modelling techniques similar to the one described in section 3.2 are more adapted to push further the thermoacoustic analysis in terms of phenomenological interpretations of azimuthal thermoacoustic modes. Therefore, the 1D analytical tool ATACAMAC is used to deal with azimuthal modes of the combustor. For Uncertainty Quantification purpose, this tool will be also used due to its affordable computational time in determining azimuthal modes. However, providing a good fitting of the industrial 3D geometry to ATACAMAC tool is the first stage of the study and the principal concern of the next Section 4.4.

Acoustic mode computations of the annular system with N injectors using network modelling tool

The whole annular combustion chamber of interest has been studied in the early work of Wolf the plenum and the burners on one hand and between the burners and the chamber on the other hand are provided.

Such methodology is here applied for the first time to typical real industrial combustion chambers of full annular helicopter engines. At first, the objective is to provide a good fitting of the full-scale gas turbine combustion chamber into a 1D thermoacoustic network representation.

To achieve this, a study of acoustic propagations through the complex 3D geometry is first conducted using the full 3D Helmholtz solver AVSP. This was done in Section 4.3. A good fitting of the industrial system is found when eigenmodes and acoustic pressure perturbations estimated from the analytical model and the Helmholtz solver are in good agreement. This explains why the full scale complex system should be modelled as a network of acoustic interconnected elements based on the 3D results obtained with AVSP and the functional operating conditions of the combustor. As ATACAMAC is a simple 1D acoustic network model, the geometrical fitting process may require to be optimized to represent the acoustics of the industrial geometry. An ill-posed setting of the combustor parameters would certainly bias the description of the target mode frequency and growth rate thus impacting the correct representation of the system stability when modifying the flame time delay τ .

The 1 st azimuthal thermoacoustic mode of the combustor in active flame regime, which appears to be the predominant mode, is targeted in this study, for both passive and active flame computations (see Table 4.3 and Table 4.4). A two-step process has been followed to ensure an appropriate fitting of the real industrial combustor and hence a good predictive representation of the system eigenmodes: The 1 st azimuthal mode of the chamber and the plenum are presented in The results suggest that the acoustic activity of the whole combustor is located in the chamber cavity. The above results are then used to calibrate the network tool ATACAMAC.

To achieve this task, the following formula used to compute the k th azimuthal mode of a simple annular cavity is used:

f kb 0 = pc b 0 2πR c and f ku 0 = pc u 0 2πR p (4.1)
where f 0 stands for the acoustic mode computed in steady flame regime, k stands for the mode number of the cavity, R c corresponds to the radius of the chamber cavity and R p stands for the radius of the plenum. The equation. Table 4.8: Parameters used for numerical applications of the annular system with N injectors. L c represents the half perimeter of the chamber and L p is the half perimeter of the plenum. L i stands for the initial burner length used for acoustic computation. S i is the section of the injector, S c stands for the section of the chamber and S p represents the section of the plenum.

Once the parameters needed to fit the network tool ATACAMAC are assessed, an eigenvalue analysis is performed to predict the stability characteristics and pulsating amplitudes of the Moreover, the stability map of the first thermoacoustic mode of the system with N injectors is investigated both with the 3D Helmholtz solver AVSP and the 1D analytical tool ATA-

a period T = 1 f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s.
Under the operating conditions stated in Table . 4.8, the network-modelling tool is not able to represent appropriately the behaviour of the physical system in active flame regime. The sign of the growth rate is not well predicted by the analytical model and the eigenfrequency is underestimated (see Fig. 4.12). This shows how the modelling process of the network model fitting is highly correlated to the geometrical parameters estimation.

Generally, simple corrections on the burner Length L i and its section S i need to be incorporated to capture 3D effects. Commonly, these two parameters are not easy to extract from the real CAD and subsequently they do not coincide with the absolute values of the industrial combustor burner. The 3D effects near the burner/chamber junctions can be accounted for [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF], [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF]) using a standard length correction in the low-frequency range for a flanged tube [START_REF] Silva | The behavior of subgrid-scale models near the turbulent/nonturbulent interface in jets[END_REF]) which is applied at the downstream burner's end (∆L i ≈ 0.4 √ 4S i π). Improper selection of these parameters would certainly bias the description of the targeted mode frequency and growth rate thus impact-TURBINE COMBUSTION CHAMBERS ing the correct representation of the system stability when modifying the flame time delay τ . Therefore, as the burner section S i has been successfully extracted from the industrial system geometry by [START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF], only the burner length L i is investigated to match with 3D Helmholtz calculations.

In this work, the range of growth rate obtained when varying the burner length is displayed in Fig. 4.13. To reach the growth rate of the first acoustic mode of interest, the burner length is estimated as

L * i = 0.231[m]
. Therefore, a posterior analysis of the growth rate disturbances accounting for the new burner length L * i is conducted (L * i corresponds to the corrected length used to fit ATACAMAC results to the 3D Helmholtz solver results).

The first azimuthal mode computed with ATACAMAC is shown in Table .

L * i = 0.231[m]
. The global interaction index n is fixed, n=1486.43[J/m] (the Crocco's value is n = 6.57 for ATACAMAC computations), and τ is varying over a period T = 1

f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s. TURBINE COMBUSTION CHAMBERS
The thermoacoustic analysis of the first azimuthal mode of the combustor is pushed further by investigating the possibility of strong coupling activities between the plenum and the combustion chamber. Beyond evaluating the maximum likelihood estimation of interactions between downstream annular chamber and upstream annular plenum, the goal is to capture the steep bifurcation of modes. This corresponds to the strongly coupled regime discussed in Section 3.2.

For that, the stability of the first azimuthal mode is constructed by varying the interaction index n and the time delay τ over a period T = 1

f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s. Knowing that the interaction index n is T b 0 T u 0 -1 = 2
.0 in the low-frequency limit, the Crocco's interaction index n is varied from n=2 to n=14. These values are taken identical for all N sectors. The corresponding stability map is shown in Fig. 4.15. Figure 4.15 shows that no major changes of frequencies in the annular plenum and the annular chamber are observed for the ranges of n and τ considered.

This indicates that the two cavities behave independently, at least to first order. The coupling parameters Γ (see Eq. (3.35)) were also evaluated for each values of the interaction index n=2 to n=14 and the time delay τ over a period T = 1

f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s.
Results are presented in Table . 4.12 which shows that the coupling parameters in the annular combustion chamber and the annular plenum are very small. Results show also that a coupling phenomenon does exist between the cavities of the combustor but most of the acoustic activity is located in the combustion chamber. In the 19 burner configuration studied by [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF], the strongly coupled regime was reached and the coupling parameters were significantly larger: about 10 orders of magnitude when comparing to Γ i presented in Table . 4.12.

As it was explained by [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF], the length and the cross section area of the burner play a predominant role on the coupling parameter (see Eq. (3.35)). Typically, Γ i goes to infinity as the burner length L i tends to zero. The burner length of the configuration studied by [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF] was much larger than in the present study: L i = 0.6m in his case against L * i = 0.23m for the system being investigated in this chapter. Additionally, the cross section area of the burner is very small in the present study when comparing to the 19 burners problem: S i = 9.9 × 10 -5 m 2 for the former and S i = 1 × 10 -2 m 2 in the latter. Consequently, Γ i of Eq. (3.35) computed in the present study are very small because of the burner parameters and thus a bifurcation of eigenfrequencies (strongly coupled regime) is never observed. The thermoa-

n |Γ i,1 | × 10 -3 |Γ i,2 | × 10 -3 |Γ i,3 | × 10 -3 |Γ i,4 | × 10 -3 CHAMBER 0 -1.

Part III

Uncertainty Quantification methods for the study of thermoacoustic 2145 instabilities in combustors STABILIZED COMBUSTOR EXPERIMENT to determine the PDF of the growth rate and the Risk Factor of the 1 st thermoacoustic mode of the system.

Monte Carlo analysis generally require a large number of model evaluations thus increasing potentially the computational burden even when combined with parallel numerical simulation tools. Therefore, for substantial computational savings, a reduced approach for Uncertainty Quantification analysis is adopted to deal with thermoacoustic systems.

The procedure is hereinafter detailed:

1 Surrogate modelling techniques are developed and introduced based on the two input uncertain flame parameters n and τ . Such surrogate modelling methods are widely used in

Computational Fluid Dynamics and have proved their efficiency at optimizing computationally expensive problems [START_REF] Rochoux | Towards predictive datadriven simulations of wildfire spread -part i: Reduced-cost ensemble kalman filter based on a polynomial chaos surrogate model for parameter estimation[END_REF]).

2 The optimal surrogate models coefficients are then determined with just a few Helmholtzbased thermoacoustic simulations arbitrarily selected from the Monte Carlo database. This task is achieved with a least mean squares methodology.

3 Once well fitted, a Monte Carlo analysis with surrogate models can replace time consuming AVSP computations to speedup by orders of magnitude the modal Risk Factor assessment.

The Uncertainty Quantification analysis is applied to a single injector, swirled stabilized combustor experiment. This system developed and built at EM2C laboratory was devoted to the study of the non linear behaviour of swirled flame dynamics accounting for changes of the acoustic environment. Section 5.2 presents the experimental set-up as well as the early experimental and numerical stability analysis conducted by [START_REF] Palies | Dynamique et instabilités de combustion de flammes swirlées[END_REF] and [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF].

The Uncertainty Quantification analysis methods are then presented:

The first UQ analysis is conducted in Section.5.3 by using a standard Monte Carlo method that is described in Section.5. This experiment is handy and practical because it was thought and conceived in such a way 2205 that both the upstream manifold (l 1 ) and the combustion chamber (l 3 ) may take respectively three and four different lengths. Hence, this simple system leads to twelve possible geometries as summarized in Table . 5.1. Thus, from twelve possible geometries, the system offers the advantage to investigate finally 24 different operating conditions. Also, acoustic losses of the system were measured during the experimental phase. From a practical point of view, measuring acoustic dissipations of a system is difficult and a global experimental strategy has not been defined to capture them. Therefore, to evaluate the acoustic damping of the experimental system, an acoustic wave has been sent through the combustion chamber to measure the response of the flame for a range of frequencies around resonance. These losses are expressed for both types of flames with an uncertainty of ∆α = ±10s -1 : α A = 82s -1 for flame A and α B = 125s -1 for flame B. The numerical acoustic modelling of the swirled combustor and its associated linear stability analysis has been realised

Cases studied

l 3 =100 l 3 =150 l 3 =200 l 3 =400 Expe./Simu. l 1 =96.0 C01 C02 C03 C04 Expe./Simu. l 1 =160.0 C05 C06 C07 C08 Expe./Simu. l 1 =224.0 C09 C10 C11 C12
by [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] by considering very small acoustic velocity perturbations for the flames.

The study was conducted with the 3D Helmholtz solver AVSP Nicoud et al. (2007). However, no intrinsic dissipation is accounted for in the Helmholtz equation and thus in the numerical simulation of the combustor. In this case, the numerical stability analysis is performed by taking 125 CHAPTER 5. UNCERTAINTY QUANTIFICATION OF A SWIRLED STABILIZED COMBUSTOR EXPERIMENT into account acoustic losses measured experimentally for both flames. Hence, the system is considered to be stable when the growth rate ω i is smaller than the damping rate α and similarly, when the computed growth rate is larger than the damping rate, the system is considered to be unstable. Moreover, accounting for the error ∆α, leads to the subsequent classification:

⋄ Stable S : ω i < α -∆α ⋄ Unstable U : ω i > α + ∆α ⋄ Marginal S/U : α -∆α < ω i < α + ∆α
Experimentally, a mode is denoted S/U when a low amplitude frequency of oscillation is detected, S if no fluctuation appears and U if a large amplitude limit cycle is observed. Numerical A good agreement is found in most of the cases when comparing the numerical and experimental stability analysis. Only three partial disagreements are observed because the experiment predicts marginal stability (S/U) while the computation gives an instability or conversely. Such a methodology which consists in classifying thermoacoustic modes in a stable or unsta-2245 ble regime does not deliver quantitative information about the risk of a mode to be unstable.

Performing an Uncertainty Quantification analysis would help to account for risk in quantitative analysis, thus offering a continuous classification of the thermoacoustic modes of the combustor.

The objective of the study is to focus mainly on the partial disagreements of Table . 5.3 and to compute the Risk Factor of the first longitudinal acoustic mode (its probability to become unsta-

f U X = 1 ||x max -x min || for x min ≤ x ≤ x max (5.1)
Therefore, the mean µ U X and the variance v U X are:

µ U X = x min + x max 2 and v U X = 1 12 (R U µ U X ) 2 (5.2)
where R U represents the normalized range xmax-x min µ U X of the uniform distribution : here

R U = 10%.
⋄ The β-distribution : The β-distribution is characterized by its density function:

f ζ Y = B(α, ζ) -1 y α-1 (1 -y) ζ-1 for 0 ≤ y ≤ 1 (5.3) where B(α, ζ)= Γ(α)Γ(ζ) Γ (
α = µ ζ Y ( µ ζ Y (1 -µ ζ Y ) v ζ Y -1 ) (5.4) 129 CHAPTER 5. UNCERTAINTY QUANTIFICATION OF A SWIRLED STABILIZED COMBUSTOR EXPERIMENT and ζ = 1 -µ ζ Y ( µ ζ Y (1 -µ ζ Y ) v ζ Y -1 ) (5.5)
To close the problem, the reduced variable Y in [0, 1] is related to the desired random variable X in [x min , x max ]:

X = µ ζ X (1 + R ζ [2Y -1]) (5.6)
Taking the mean and variance of the previous equation leads to the following relations between characteristics of X and Y:

µ ζ Y = 1/2 and v ζ Y = ν ζ X 4R 2 ζ (µ ζ X ) 2 (5.7)
Consequently, the mean value of Y is fixed and its variance can be deduced by imposing that the Beta and uniform PDFs have the same characteristics, e.g.

µ ζ X = µ U X and ν ζ X = ν U X .
However, the range of the β-PDF appears in (µ ζ X )(Eq. ( 5 In Fig. 5.7a, each point corresponds to a Helmholtz simulation in the complex domain. The horizontal solid lines denotes the acoustic losses α: 115 s -1 < α B < 135 s -1 . The stable or unstable regions are evaluated using the difference ω i -α: 2305 1. ω i -115 s -1 < 0 corresponds to a stable system (S).

2. ω i -135 s -1 > 0 corresponds to a unstable system (U).

3. 115 s -1 < ω i < 135 s -1 corresponds to a situation where the system is marginal (neither stable nor unstable) (S/U).

The M samples are then classified as follows: stable regime (S), unstable regime (U) and marginal 2310 regime (S/U). In Fig. 5.7b, the PDF of the growth rate (ω i ) is presented and shows that most of the thermoacoustic modes found by the Helmholtz solver are in the stable regime. This leads to a Risk Factor close to 24 % thus meaning that the acoustic mode has 24 % of chance to be unstable. Following a similar methodology as for the uniform distribution, 4000 runs have been performed using the Helmholtz solver considering this time a β-distribution for the input parameters n and τ . Results are presented in Fig. 5.8. The Risk Factor obtained from the β-distribution is close to the one obtained by the uniform distribution: 24 % for the uniform distribution against 22 % for the β-distribution. This shows that UQ results are weakly affected by the distributions chosen for the input parameters n and τ for the study of such academic cases which suggests that assessing the Risk Factor of a mode without a clear knowledge of the uncertainties on the input data is relevant. Moreover, the Risk Factor being 22 -24 %, this simple UQ analysis shows that the computation is actually consistent with the experimental data. Indeed, accounted for a realistic 10 % uncertainty in the flame response, this Risk Factor value means that the mode of interest is computationally found stable in approx. 76 -78 % of the cases (recall that the mode of Case 07-Flame B was observed stable in the experiment, see Table . 5.1). In the rest of the study only the uniform distribution is kept for the UQ analysis.

For this simple system, one Helmholtz simulation took approximately 24 minutes on 16 cores. STABILIZED COMBUSTOR EXPERIMENT This run time seems to be not prohibitive but may quickly become so when performing a Monte Carlo analysis when the complexity of the system increases, which typically is a computationally intensive undertaking. In light of this, investigating suitable surrogate modelling methodologies would help to reduce this computational cost.

Surrogate modelling techniques

In this section, an Uncertainty Quantification strategy based on reduced-order models approach is proposed and described in Fig. 5.9. Reduced-order models are developed and introduced to determine the growth rate variation of the system. Such surrogate models are tailored to tackle uncertainties related to the Flame Transfer Function parameters n and τ . A linear regression method is then used to determine the models coefficients on the basis of the 4000 Helmholtz simulations previously generated in Section.5.3.1. Further analysis are conducted to evaluate the statistical efficiency of these models as well as their level of accuracy in approximating the Risk Factor of the first longitudinal mode of the system.

Linear regression

Because Eq. (1.2) is an eigenvalue problem which is nonlinear in ω i , the response surface

ω i = ω i (n,
ω n-τ i = ζ 0 + ζ 1 n + ζ 2 τ (5.8)
2 LM F T F : based on the Flame Transfer Function evaluated at ω = ω 0 , where ω 0 corresponds to the mode without flame coupling (corresponding to n=0). The Flame Transfer Functions incorporate here physical non-linearities into the model:

ω FTF i = ζ 0 + ζ 1 ℜ(ne jω 0 τ ) + ζ 2 ℑ(ne jω 0 τ ) (5.9)
3 QM F T F : is a quadratic model based on the Flame Transfer Function also evaluated at ω = ω 0 . Here, the physical non-linearities are taken into account into the model.

ω Q FTF i = ζ 0 + ζ 1 ℜ(ne jω 0 τ ) + ζ 2 ℑ(ne jω 0 τ ) + ζ 3 ℜ(ne jω 0 τ ) 2 (5.10) +ζ 4 ℑ(ne jω 0 τ ) 2 + ζ 5 (ℜ(ne jω 0 τ ) × ℑ(ne jω 0 τ )) (5.11)
The models LM n-τ , LM F T F and QM F T F can be written in linear algebra notation as follows:

ω i = Xζ + ϵ (5.12)
where Xζ is the matrix-vector product, ζ=

[ζ 0 , ζ 1 , ζ 2 , ζ 3 , ζ 4 , ζ 5 ]
T corresponds to the regression coefficients of the model. These coefficients represent the mean change in the response variable for one unit of change in the predictor variable. ω i is considered to be a N × 1 dimensional vector containing the growth rate ω i determined from N Helmholtz computations, X is the matrix containing:

⋄ 1, n and τ when using LM n-τ , ⋄ 1, ℜ(ne jω 0 τ ), ℑ(ne jω 0 τ ) for the linear model LM F T F , ⋄ 1, ℜ(ne jω 0 τ ), ℑ(ne jω 0 τ ), ℜ(ne jω 0 τ ) 2 , ℑ(ne jω 0 τ ) 2 and (ℜ(ne jω 0 τ )×ℑ(ne jω 0 τ )) for the quadratic model QM F T F , and this for each sample and ϵ the N×1 vector of residuals: STABILIZED COMBUSTOR EXPERIMENT

ω i = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ω i 1 ω i 2 . . . ω i N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , X = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 n 1 τ 1 1 n 2 τ 2 . . . . . . . . . 1 n N τ N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ζ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ζ 0 ζ 1 ζ 2 ζ 3 ζ 4 ζ 5 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ and ϵ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ϵ 1 ϵ 2 . . . ϵ N ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦
A least squares methodology is used to determine the coefficients ζ of the three models which minimize the error ϵ:

ζ = ( X t X ) -1 X t ω i (5.13)
where ζ corresponds to the estimated parameters from the least squares, (X t X) -1 is called the "information matrix" and X t corresponds to the transpose of the X matrix. The predicted values ωi for the mean of ω i of the three models are then determined as follows:

ωi = X ζ = X ( X t X ) -1 X t ω i (5.14)
The idea is now to use the surrogate models formulated above to approximate the results found in section .5.3.1. Such a validation process is achieved through the following steps:

1 The ζ-coefficients of each model are found using the full set of 4000 Helmholtz simulations of the Monte Carlo database. These coefficients are computed using Eq. (5.13) and displayed in Table . 5.4. STABILIZED COMBUSTOR EXPERIMENT

ζ-coefficients LM n-τ LM F T F QM F T F ζ 0 -0.0312×10 3 -4.5014 -7.5811 ζ 1 0.0 -0.0160 -0.0142 ζ 2 4.9897×10 3 -0.0152 -0.0264 ζ 3 - - 4.8176×10 -6 ζ 4 - - -1.8057×10 -6 ζ 5 - - -1.0596×10 -5
Table 5.4: ζ-coefficients computed for surrogate models LM n-τ , LM F T F and QM F T F using the 4000 samples of the Monte Carlo database.

2 Once the ζ-coefficients computed, the Pearson Correlation Coefficient is computed to provide an index of the degree of correlation between the surrogate models outcome and the reference Monte Carlo database.

R = E[(ω i -E(ω i ))(ω model i -E(ω model i ))] σ ω i σ ω model i (5.15)
In Eq. (5.15), E is the expectation, ω i corresponds to the reference growth rate, ωi is the growth rate issued from linear least squares fitting and σ corresponds to the standard deviation from the reference growth rate and the estimated growth rate from linear least The regression analysis shown that LM F T F (Eq. (5.9)) and the quadratic model QM F T F (Eq. (5.11)), are able to reproduce respectively 98% and almost 100% of the growth rate variation whereas the model LM n-τ reproduced 95% correlation of the growth rate variations.

3 Therefore, the algebraic surrogate models LM F T F and QM F T F should be rather accurate to mimic the actual response surface of the system and to estimate, with a minimum error, the Risk Factor of the mode. To assert this, a Monte Carlo analysis is applied to the surrogate models LM n-τ , LM F T F and QM F T F to construct the PDF of the growth rate and to estimate the modal Risk Factor. Fig 5 .11 shows the PDF of the growth rate determined from surrogate models and Table 5.6 shows the corresponding Risk Factor estimated.

Surrogate model Risk Factor in

% LM n-τ 21 LM F T F 23 QM F T F 24
Table 5.6: Risk Factor and computation time estimated from from surrogate models. The whole set of Helmholtz simulations (4000) were used.

The surrogate models evaluations are here almost instantaneous and provide good trends of the growth rate distribution and a good estimation of the modal Risk Factor of interest. Among all surrogate models, the model QM F T F appears to be more accurate in predicting the Risk Factor of the mode when comparing to the reference Risk Factor obtained with Helmholtz solver (≈ 24%).

So far, whole sets of Helmholtz simulations (4000) obtained with the 3D AVSP solver have 2400 been used to tune the surrogate models. For the sake of Uncertainty Quantification analysis, it would be interesting to seek how to tune these models to estimate accurately, with just a few STABILIZED COMBUSTOR EXPERIMENT Helmholtz simulations, the Risk Factor of the first acoustic mode of the system (relying on much less then 4,000 Helmholtz simulations to fit the surrogate models).

Risk Factor estimation with reduced-order models

In this section, a reduced Uncertainty Quantification strategy which combines few Helmholtz simulations and surrogate modelling is employed. This UQ strategy has distinctive features:

⋄ To avoid CPU-intensive Helmholtz simulations with the 3D and parallel solver AVSP, the surrogate models are tuned using only a limited number of Helmholtz simulations. A Monte Carlo analysis is then achieved with the surrogate models to get an estimation for both the PDF of the growth rate ω i and the modal Risk Factor of interest.

⋄ However, the subset of Helmholtz simulations required to fairly estimate the modal Risk

Factor with the surrogate models needs to be determined. To do so, several evaluations of the surrogate models are realised based on different subsets of randomly selected Helmholtz computations from the full Monte Carlo database. As a consequence, the mean Risk Factor and its standard deviation are evaluated for each subset of Helmholtz simulations used.

This allows to get an insight on the variability of the Risk Factor for each size of Helmholtz samples. Moreover, the confidence intervals for the mean Risk Factors are computed which in addition provides a deduction on the number of Helmholtz simulation required to approximate the modal Risk Factor with the surrogate models.

⋄ Finally, the impact of each uncertain parameter (n and τ ) on the growth rate variations is discussed after deriving the surrogate models.

The surrogate models developed in section 5.3.2.1 are used to ease the construction of the growth rate distribution. Only a small dataset of Helmholtz-based thermoacoustic simulations to provide an unbiased estimate of the modal Risk Factor. The large number of runs required for accurate predictions is necessarily not compatible with costly computational tools based for example on finite/volume element models or complex industrial systems, and this even when high-performance computing platforms are at hand.

Quantifying the impact and accuracy of such models is necessary to producing defensible claims in the context of reliable Risk Factor approximation. One approach is to choose a subset of Helmholtz simulations to determine the ζ-coefficients of the surrogate models using the least mean squares fitting method described in section 5.3.2.1. Once models fitted, the Risk Factor is evaluated for a Monte Carlo analysis based on the models. This process should be then repeated by increasing gradually the number of Helmholtz samples to the model fitting procedure until adequate convergence of the Risk Factor is reached (comparable to the reference Risk Factor obtained from the Monte Carlo analysis with AVSP ≈ 24%). However, each of the Helmholtz samples are added without replacement otherwise this would biased information in the Risk Factor approximation for each subset. Maximizing the number of Helmholtz samples will provide better coverage in the growth rate design space and should provide locally the level of accuracy of the surrogate. For completeness, monitoring the Risk Factor estimated for each subset of Helmholtz simulation is interesting to determine the error between the surrogate model and the deterministic model evaluation (with AVSP). As efficient computational surrogate models are used in this work, the computer cost is not a stumbling block to perform several surrogate model evaluations. This provides the standard deviation of the Risk Factor for each subset of Helmholtz simulations and an indication of their corresponding confidence intervals. Finally, this will provide the minimum number of Helmholtz simulations required to get a fair estimation of the modal Risk Factor with reduced-order models.

From a practical point of view, the Uncertainty Quantification analysis goes through the steps presented in Fig. 5.12 and hereafter detailed:

1

Step 1: The work achieved in section 5.3.2.1 proved that LM F T F and QM F T F are better correlated with the full Monte Carlo database than LM n-τ . Thus only surrogate models LM F T F and QM F T F are kept for UQ analysis purpose in the rest of the study. For each model, the goal is to determine their regression coefficients at reduced cost thus relying only on a few samples of Helmholtz simulations instead of the 4,000 initially performed in section 5.3.1. Therefore, for surface fitting of each surrogate model, a subset of 3, 5, 6, 10, 20, 40 and 100 Helmholtz simulations are randomly collected (sampling without and QM F T F using a different number of Helmholtz simulations from the full MC database.

4

Step 4: Moreover, the mean Risk Factors and associated standard deviation are investigated for each size of Helmholtz samples used (from 3 to 100 samples). Results are summed up in Table 5.7 and Fig. 5.15 describes the evolution of the standard deviation when using LM F T F (black) and QM F T F (red). For both surrogate models, the standard deviations exhibit a significant drop for lower subset of Helmholtz samples (from 3 to 10 Helmholtz samples). Then, the variation of the standard deviations becomes very weak until being almost imperceptible as shown in Fig. 5.15. This suggests that only a few tens of Helmholtz simulations is enough to converge towards a good estimate of the modal Risk

Factor when using such surrogate models. Another way to ensure these observations is to provide a prediction confidence interval (CI) with the surrogate models to evaluate the confidence for the mean Risk Factors obtained with the different size of the Helmholtz samples.

These confidence intervals are computed by the following formula: The UQ strategy followed in this work shows that combining surrogate models with a limited number of Helmholtz simulations allows to capture, to a satisfactory degree, the Risk Factor of the mode with a good predictive confidence interval. The use of such surrogate modelling techniques allows to overcome the impediment of time consuming by orders of magnitude.

CI = µ RF ± z * σ √ n ( 5

Investigation of the other cases

This section aims at investigating the other partial disagreements of Table 5.3: the configuration 11 Flame A and the geometry 11 Flame B. Instead of performing an expensive Monte Carlo analysis with the 3D Helmholtz solver AVSP, Uncertainty Quantification studies are pursued based on reduced-order models developed and introduced for the previous geometry 07 Flame B. For the latter case, the standard deviation decreases as the number of Helmholtz samples increases. Moreover, the decrease in the average relative error of the standard deviation is not large when the number of Helmholtz samples varies from 10 up to 100 and there is not a significant improvement in the reliability of the modal Risk Factor when larger sample are used. Based on these observations, only a hundred of Helmholtz simulations are sampled from a uniform distribution using AVSP solver. STABILIZED COMBUSTOR EXPERIMENT Initially, the overall hundred computations are used to fit the ζ-coefficients of the surrogate models LM F T F and QM F T F and to approximate the modal Risk Factor of the system. Then, as for the geometry 07 Flame B, a sensitivity analysis on the Risk Factor is investigated through different tuning of the ζ-coefficients of the surrogate models.

Test case 2: The configuration 11-Flame A

The Uncertainty Quantification analysis of the 1 st acoustic mode of the geometry 11 Flame A is now investigated. The objective is to seek the probability of the mode to be unstable (f 0 = ω 0 /2π = 104 Hz) namely its Risk Factor. For this operating point, the experimental stability analysis predicted a marginal regime while a stationary state has been concluded numerically.

The objectives are to investigate if:

⋄ reduced-order models provide good fits to the entire data set made of 100 Helmholtz samples ⋄ small relative errors on the Risk Factor estimation are found when the sampling size is drastically reduced to 10 Helmholtz runs. For this, 5 subsets composed of 10 Helmholtz runs each are constructed based on the entire data set. Then, for each scenario, 100 Monte Carlo model evaluations are performed to determine if a reduced sampling size of 10 is enough to obtain reliable estimates of the variability in the growth rate and hence in the modal Risk Factor of the system.

The statistical analysis is carried out using only the models LM F T F and QM F T F which shown better results in the previous case. The range of uncertainty used are similar to those of the geometry 07 Flame B: ∆n n = ∆τ τ = ±10%. To propagate uncertainties, a uniform distribution is used to generate random perturbations of the flame parameters n and τ . Based on the findings of the case 07 Flame B, the choice of the PDF has not an important impact as much on the Risk Factor estimation. However, since the realistic growth rate distribution of the mode is unknown, the accuracy of the growth rate estimates would be determined by how well the surrogate models fit the Helmholtz database.

At first, 100 Helmholtz simulations are performed using the Helmholtz solver AVSP. The overall Helmholtz runs performed are then used to tune the surrogate models LM F T F and QM F T F sis is then performed using the surrogate models to get the PDF of the growth rate and hence an estimation of the Risk Factor of the first thermoacoustic mode of the configuration. The ζcoefficients, defined by Eq. (5.13) and calculated for both surrogate models LM F T F and QM F T F , are presented in Table . 5.8. The least mean squares fitting as well as the Pearson's correlation coefficients computed using Eq. 5.15 are shown in Fig. 5.17 and merged in Table . 5.9. The

ζ-coefficients LM F T F QM F T F ζ 0 5.6 2.4 ζ 1 -3.6 ×10 -3 -4.7 ×10 -3 ζ 2 -3.5×10 -3 -6.6 ×10 -3 ζ 3 -4.9 ×10 -7 ζ 4 -3.9 ×10 -7 ζ 5 -1.4 ×10 -8
Table 5.8: ζ-coefficients determined for surrogate models LM F T F and QM F T F based on the 100 samples computed with AVSP code for the geometry 11 of Flame A. Results show that there is a risk of 96%, within ±1% depending on the surrogate model used,

for the 1 st acoustic mode to become unstable under these operating conditions.

To further investigate the effect of the Helmholtz sample size, a sensitivity analysis of the Risk Factor predicted with the surrogate models is conducted using a set of 10 Helmholtz simulations 2550 STABILIZED COMBUSTOR EXPERIMENT the accuracy of the surrogate models LM F T F and QM F T F in predicting the modal Risk Factor of the system. This complements and comes to reinforce the results of the statistical analysis conducted for the configuration 07 of the flame B.

For this configuration 11 Flame A, the experimental results could not provide a clear evidence of the mode regime. The Uncertainty Quantification study helped to refine the thermoacoustic analysis by confirming that this operating point is most probably unstable. The reason why a strong instability was not detected experimentally remains unclear.

Test case 3: The configuration 11-Flame B

For this last case, the stability analysis with AVSP predicted a stable regime while a marginal regime was found from the experimental stability analysis. The Uncertainty Quantification analysis that combines reduced-order modelling techniques and few Helmholtz samples is used once again.

As shown in the previous sections, choosing only a few tens of Helmholtz simulations is enough to get an accurate estimate of the modal Risk Factor. Besides, the UQ analysis conducted for the geometry 11 Flame A highlighted that 15 Helmholtz samples are enough to tune the surrogate models LM F T F and QM F T F . On the basis of this, only 15 Helmholtz computations were performed for this operating point and used to determine the ζ-coefficients of both surrogate models. Results are presented in Table . 5.12. The least mean squares fitting obtained by using the 15 Helmholtz samples is displayed in Fig 5 .19 and the Pearson's correlation coefficients computed using Eq. 5.15 are presented in Table . 5.13.

The results show a good correlation between the surrogate models and the Helmholtz samples computed from AVSP. The Risk Factor computed when using the surrogate models LM F T F and To ensure the results obtained when using 15 Helmholtz samples, 100 and 1000 additional runs were performed to get an estimation of the modal Risk Factor. The Pearson's correlation coefficient computed when using 100 and 1000 samples to fit the surrogate models LM F T F and 2585 QM F T F are presented in The Risk Factors estimated when using either 100 or 1000 Helmholtz runs is similar to those obtained when using only 15 Helmholtz samples. This means that 15 Helmholtz runs are enough to fit both surrogate models and to reproduce the growth rate variations of the system.

QM F T F are summarized in
For this configuration 11 flame B, assuming uncertainties on the Flame Transfer Function 2590 STABILIZED COMBUSTOR EXPERIMENT using a Monte Carlo approach based on several Helmholtz simulations of a single experimental operating point but with random perturbations on the Flame Transfer Function parameters.

Then, a two-step UQ strategy was used to deal with thermoacoustics in such a system: (i) First, three surrogate models were tuned from a moderate number of Helmholtz solutions (ii) Then, these algebraic models were used to perform a Monte Carlo analysis affordably and to approximate the Risk Factor of the mode. The study proves that analytical surrogate models can be used to predict the Risk Factors within good predictive confidence intervals.

The modal Risk Factor assessed for each geometry is hereinafter summarized:

⋄ The configuration 07 Flame B: For this case, the experiment predicted a stable regime while the numerical stability analysis predicted a marginal regime. When accounting for uncertainties on the flame model parameters, the Risk Factor associated to the first acoustic mode of the geometry is approximated to 24%, meaning that the mode has 24% of chance to be unstable when accounting for a 10% uncertainty on the flame model input parameters. In other words, the partial disagreement between the experimental and the numerical stability analysis can be partially explained by uncertainties on the flame model parameters.

⋄ The configuration 11 Flame A: For this geometry, the experiment predicted a marginal regime while the numerical stability analysis predicted an unstable regime. A 10% uncertainty on the flame model parameters lead to 99% of probability or the mode to stay unstable. This means that the mode is found unstable numerically, even if the flame parameters are quite uncertain. Thus, the partial disagreement between the experimental and the numerical stability analysis can hardly be explained by the limited knowledge of the flame response and other explanations must be sought.

⋄ The configuration 11 Flame B: For this case, the experiment predicted a marginal regime while the numerical stability analysis predicted a stable one. When accounting for uncertainties on the flame model parameters, there is no probability for the acoustic mode to be unstable. In other words, the stability of the mode is not altered when accounting for a 10% uncertainty on the flame model parameters. As for the configuration 11 flame A, the partial disagreement found between the numerical and the experimental stability analysis STABILIZED COMBUSTOR EXPERIMENT could not be explained by uncertainties on the flame model input parameters.

In the work of [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF], the stability analysis of the combustor was investigated by accounting for the amplitude of the velocity perturbation by using the Flame Describing Function formulation. Typically, the Flame Describing Function formulation is used to describe the nonlinear flame response to harmonic velocity perturbations over a range of forcing frequencies.

Therefore, this method allows to predict the amplitude and frequency of limit cycle oscillations in non-linear feedback systems. The Flame Describing Function is defined as:

F (|û|, ω) = n(|û|, ω)e iωτ = Q(|û|, ω)/Q tot û/U bulk (5.17)
where |U bulk | stands here for the amplitude of acoustic perturbations (see Section 3).

Typically, the work of [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] was achieved in two steps:

⋄ The numerical stability analysis of the system was performed by considering only the smallest value of the acoustic perturbations |û| ūA and |û| ūB for the two flames A and B.

⋄ Then, the frequencies and the growth rate variations of the modes were investigated as a function of the amplitude of the acoustic oscillations.

The results of [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] show that when increasing the amplitude of acoustic velocity perturbations, the growth rate of the acoustic modes decreases before reaching a limit cycle when the growth rate equals the damping rate of the system. It means that the flame function parameters n and τ are not the only sources of uncertainties that control the stability of the system. Indeed, small variations of the amplitude of the acoustic velocity perturbations |û| may also modify the growth rates. Therefore, the idea would be to investigate the uncertainties related to the amplitude of the velocity perturbations in the Flame Describing Function model.

This would help to complete the UQ analysis by measuring for example the effects of these acoustic perturbations on the modal growth rates at least for the two partial disagreements of ⋄ Collocation formulations are slightly less computationally expensive than Polynomial chaos methods as discussed by [START_REF] Dwight | Efficient uncertainty quantification using gradient-enhanced kriging[END_REF].

⋄ Dimension reduction approaches through gradient-based global sensitivity analysis are proposed to reduce the number of parameters in the system and to ease scalability to highdimensional problems. Active subspace method [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF]) is one of these approaches.

This chapter intends to highlight the potential of dimension reduction by exploiting active subspaces to quantify uncertainty. These approaches are applied to the realistic annular helicopter engine studied in Section 4 that features N circumferentially arranged and identical burners. Each burner is described by two uncertain input parameters used to represent the flame response n and τ . Therefore, we are facing the famous «curse of dimensionality» as no less than thirty independent uncertain parameters are involved in this case. The Uncertainty Quantification analysis is performed using the 1D Analytical tool ATACAMAC detailed in Section 3.2. This tool has been retained because it encompasses the essential features of azimuthal modes developing in complex annular combustors. Furthermore, it does not require heavy computational resources since only an algebraic model is evaluated to provide azimuthal eigenmodes (about few minutes of computation against hours with 3D Helmholtz solver and days with LES techniques). This allows extensive and quick comparison of different Uncertainty Quantification strategies: (i) the brute force Monte Carlo method and (ii) the Active subspace technique combined with surrogate modelling approaches are used for the study. Moreover, this tool has been successfully employed ACTIVE SUBSPACE METHOD recently to develop a novel Uncertainty Quantification approach combining Active Subspace and Adjoint towards the study of symmetry breaking effects of azimuthal modes in annular combustors [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF], [START_REF] Magri | Stability analysis of thermoacoustic nonlinear eigenproblems in annular combustors. (part 2: Uncertainty quantication)[END_REF]).

To work around the dimensionality issue towards Uncertainty Quantification analysis, the following tasks are performed:

1 At first, the brute force Monte Carlo is applied on the full parameter space. To achieve this task, the least biased uniform distribution is employed to generate random perturbations of the flame input parameters n and τ . Uncertainties are then propagated through the system to determine the PDF of the growth rate ω i and to approximate the Risk Factor of the first azimuthal mode e.g. its probability to become unstable.

2 The Active Subspace method is then used to capture and exploit the relevant subspaces of the system along which the growth rate variations are important. To do so, an eigenvalue decomposition of the gradients of the growth rate must be performed. Numerically, finite difference techniques are then used to approximate the derivatives of the growth rate and thus the active subspace of the system. Hence, the system dimensionality is drastically reduced.

3 Linear and quadratic surrogate models are built, based on the active variables discovered from the Active Subspace method. Such models proved satisfactory in cheaply and accurately estimating the Risk Factor of a mode as discussed in Section. 5 and [START_REF] Ndiaye | Uncertainty quantification of thermoacoustic instabilities in a swirled stabilized combustor[END_REF]. Such surrogate being inexpensive to evaluate, exhaustive sampling is realised to determine the PDF growth rate and subsequently the modal Risk Factor of the system.

These are then compared to the results obtained with the brute force Monte Carlo method performed in the first task.

Analysis with Monte Carlo method

One established solution and widely used method for risk management under uncertainties is Monte Carlo. Therefore, taking advantage of the affordable computation with ATACAMAC, n (where n is the nominal value of n) and the plot on the bottom represents the PDF of the time delay for the dimensionless ratio τ τ (where τ is the nominal value of τ ). In both plots, 10,000 ATACAMAC computations were generated.

The Active Subspace approach

In this section, the definition of the Active Subspace is reviewed from [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF].

Recently, this method has been applied by [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF] to explore symmetry breaking effects in a simplified annular combustor.

Problem formulation

Active subspace method is an emerging approach that gives insight into the relevant directions in the input parameter space; the relative change in each component of the input space along these directions generate the largest change of the output quantities of interest. This method relies upon the fact that the response tends to vary more prominently in a few dominant directions.

The latter are defined by linear combinations of the original model's inputs.

Consider a differentiable and square-integrable function f Im ∈ R in such a way that:

f Im = f Im (x).
(6.1)

In the present case, f Im is the objective function describing the growth rate response of the system for which the inputs are x = {n i , τ i } i=1...D (D is the number of dimensions in the problem).

Denote the gradient f Im by ∇f Im ∈ R L with partial derivatives ∂f Im ∂x i . Evaluation of ∇f Im might be achieved in different ways e.g. finite differences, adjoint method or automatic differentiation (typically, an active subspace for f Im will be a linear subspace for which f Im change a lot more on average along direction in the active subspace than along those in the complementary inactive ACTIVE SUBSPACE METHOD subspace). By considering that all the partial derivatives of f Im are square-integrable, an average derivative functional expressed as the matrix C ∈ R L×L , also called the uncentered covariance matrix, can be defined by:

C = E [ (∇f Im (x))(∇f Im (x)) T ] , (6.2)
where E is the expectation operator.

As the matrix C is symmetric, positive semi-definite it admits the following real eigenvalue decomposition: (6.3) where W ∈ R L×L is an orthogonal matrix whose columns w 1 , ..., w L are the eigenvectors of C .

C = W ΛW T ,
Consequently, W T i (x) are the reduced coordinates e.g. the active variables. Λ ∈ R L×L is a diagonal matrix with diagonal entries diag(λ 1 , ..., λ L ), λ 1 ≥ ... ≥ λ i ≥ 0, that include eigenvalues of the matrix C .

The eigenvalue λ i that relates the effect of the active variables W T i (x) on the growth rate response f Im , is in fact the mean-squared value of the directional derivative of f Im in the direction w i :

λ i = w T i C w i = w T i E [ (∇f Im )(∇f Im ) T ] w i = E [ (∇f Im • w i ) 2 ] .
(6.4)

The partitioning of the eigenvalues in Eq. 6.4 can be used to define a new coordinate system: the more λ i is important and the more significant the active variable W T i x is on the average output response. Therefore, the strongest active variables can be isolated.

As explained by [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF] and [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF], when only a few linear combinations of the input parameters are relevant (a few eigenvalues are much larger than any others) the system dimensionality can be reduced to just a few. For this reason, exploring such low-dimensional subspace is extremely valuable for Uncertainty Quantification analysis and this is the interest of the study.

The Uncertainty Quantification strategy applied to the realistic annular combustor with 2×N uncertain parameters (2 uncertain parameters per flame (n and τ )) is sketched in Fig. 6.3:

1 Active Subspace method is used to reduce the system dimensionality to only 3 dimensions. Risk Factor is computed. Risk Factors approximated with surrogate models are compared against the Risk Factor estimated from the brute force Monte Carlo analysis (Section. 6.3). Finally, the Risk Factor is computed using the low-order models and validated against the brute force Monte carlo Analysis with ATACAMAC on 10000 samples.

Identification of Active Subspaces

The numerical approximation of the Active Subspace can be realised using the Monte Carlo method [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF]. Therefore, ∇f Im = ∇ k x f Im for the k th sample must be 2805 computed using the following Monte Carlo approximation to the covariance matrix C :

C = E [ (∇ x f Im )(∇ x f Im ) T ] ≈ 1 M M ((∇ x f Im )(∇ x f Im ) T ) = W Λ W T , (6.5) ACTIVE SUBSPACE METHOD
of magnitudes of the eigenvalues components and it shows that M=50 samples are enough to converge Λ correctly; for smaller samples, the eigenvalues are scrummed and difficult to identify.

The following conclusions can be drawn from the spectrum analysis:

⋄ The first eigenvalue is a good metric for evaluating the global sensitivity of the combustor to the input uncertainties x = {n i , τ i } i=1...D .

⋄ At first, it appears that the Uncertainty Quantification problem can be reduced to a 5dimensional problem ( instead of the 2 × N initial dimensions).

Moreover, Bauerheim and co-workers [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF]) investigated how to reduce the dimension of the problem with the Active subspace method when the eigenvalues are difficult to determine, when the physical behaviour of the system become complex and bifurcation of modes occur in the combustor. This is typically the case when eigenmodes of the combustor are strongly coupled. In this case, instead of increasing the number of gradient evaluations, an alternative is to perform a change of variables to ease the physical interpretation of active variables W T i x and to improve the eigenvalue decomposition of the matrix C . To achieve this, the theoretical studies of [START_REF] Noiray | Investigation of azimuthal staging concepts in annular gas turbines[END_REF] and Bauerheim et al. (2014a) for annular system without plenum, i.e. Γ i,1 = Γ i,2 = 0, can be used. These theories stipulate that the complex frequencies of the mode of order p for weakly coupled modes are:

f ± c = pc 0 b 2L c - c 0 b 4πL c (Σ 0 ± S 0 ) , (6.6)
where L c is the chamber length and c 0 b the sound speed in hot gases (see in Section. 3.2) In Eq. 6.6, Σ 0 is the «coupling strength» defined as:

Σ 0 = N ∑ i=1 Γ 0 i (6.7)
This parameter is the sum of all the coupling parameters of the system, and is independent of the pattern used to distribute the burner uncertainties along the annular chamber. It corresponds to a symmetric effect.

The parameter ±S 0 is the «splitting strength» which distinguishes the two azimuthal mode frequencies f + c and f - c . A convenient form of this parameter is obtained by using the spatial CHAPTER 6. UNCERTAINTY QUANTIFICATION USING THE ACTIVE SUBSPACE METHOD Fourier transform of the coupling parameter distribution γ:

S 0 = √ γ(2p)γ(-2p) where γ(k) = N ∑ i=1 Γ 0 i e -j2kπi/N (6.8)
Note that the «coupling strength» can be also be recast in this form, i.e., Σ 0 = γ(0). It shows that only few specific patterns can affect the azimuthal mode stability. They correspond to the 0 th and the ±2p th Fourier coefficients γ of the coupling parameter or heat release distribution [START_REF] Noiray | Investigation of azimuthal staging concepts in annular gas turbines[END_REF]). Unlike the coupling strength Σ 0 , the splitting parameter S 0 can be changed by modifying the pattern of the burner types along the annular chamber. Such a modification can be intended as when controlling devices are introduced, or unintended when turbulence or uncertainties affect randomly the flame response to acoustics. In a UQ perspective, the explicit solution of Eq. (6.8) allows the CPU cost to be drastically reduced since only patterns associated with γ(0) and γ(±2p) can be retained [START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF]). Recently, [START_REF] Ghirardo | Stability criteria for standing and spinning waves in annular combustors[END_REF] also shown that non-linearities of the flame response itself can produce a splitting effect [START_REF] Ghirardo | Stability criteria for standing and spinning waves in annular combustors[END_REF][START_REF] Bauerheim | Experimental and numerical studies of active and passive control of combustion instabilities[END_REF]2016)). The azimuthal mean flow induced by swirlers or modern effusive plates can also promotes such a splitting Bauerheim et al. (2014a).

The above theoretical asserts are used in this work to incorporate phenomenological interpretation of the active variables through the Fourier Transform of the Flame Transfer Function such as:

{n i , τ i } → {Re(γ), Im(γ)} . (6.9)

Eigenvalues spectrum determined using Eq. (6.9) and the corresponding gradient matrix (again computed by finite differences) are presented in Fig. 6.5.

The results show that the eigenvalues convergence is quicker when using the Fourier transform formalism and the spectrum Λ is accurately predicted when using only M = 20 samples. Moreover, it is observed that the system reduces from 5D to only a 3D parameter space in this case thus meaning that only the 3 first active variables are relevant and lead to the strong perturbations of the growth rate in the combustor.

Exploiting Active Subspaces to Quantify Uncertainty

In the above section, a technique for discovering the possible dependence of the growth rate response to a lower-dimensional active subspace was addressed. This lower-dimensional subspace is based upon a small subset of the original design full-space dimension. The procedure enables ACTIVE SUBSPACE METHOD to reduce significantly the dimension of the problem from full-dimensional space space space to a 3D active space involving physical quantities associated to the Fourier transform of the Flame Transfer Function. The inactive variables of the system having been chased down, the objective is now to take advantage of the low-dimensional active subspace discovered. Thus, physics-based reduced order models are proposed to get insight of the growth rate variations when accounting for uncertainties on the flame response parameters n and τ .

Training surrogate models can be difficult for complex problems because of the amount of evaluation-time needed to provide a good fit. Typically, the number of simulations required depends mostly on the characteristics of the surrogate (i.e. the polynomial order) and the dimensionality of the input parameter space. Fortunately in this work, the total simulation time needed to provide eigenmodes of the system is well affordable (few minutes of computation with ATACAMAC solver) and subsequently surrogate models of different complexities are investigated. Although focusing on the model's response along active directions, a «whole» polynomial representation of the problem over the full-dimensional space space space is constructed and evaluated.

Four types of surrogate models are studied:

fIm = ζ 0 + D ∑ j=1 α j W j    Linear (L) + D ∑ j=1 D ∑ k=1 β j,k W j W k    Quadratic (Q) (6.10)
⋄ Linear models:

-L F D : The first linear model is constructed in the full probabilistic space.

-L 3D : The second linear model is spanned along the reduced subspace with the 3 active variables discovered with Active Subspace method.

⋄ Quadratic models:

-Q F D : The first quadratic model is constructed in the full probabilistic space.

-Q 3D : The second quadratic model is built on the reduced subspace with the 3 active variables discovered with Active Subspace method.

A summary of the different surrogate models investigated is presented in Fig. 6.2. ACTIVE SUBSPACE METHOD 2 The resulting approximated surrogate models are then evaluated randomly from 100 to 100000 times on a Monte Carlo dataset. Convergence tests prove that 25000 simulations are enough to reach a converged estimation of the modal Risk Factor with surrogate models.

3 Finally, Risk Factors computed with surrogate models are compared to the one obtained from the brute force Monte Carlo in Section. 6.3.

The fitting procedure

This section explores further the fitting procedure of the surrogate models by evaluating:

⋄ the total number of coefficients required by the four surrogates.

⋄ the total number of evaluation points needed for approximating correctly the modal Risk

Factor with surrogates.

The surrogate forms of interest are linear and quadratic. therefore, the total number of coefficients needed in a D dimensional space is Φ L (D) = (D + 1) for linear models and Φ Q (D) = (D+1)(D+2) 2 for quadratic ones. The number of coefficients required for each model is summed up in Accounting for the complexity of each surrogate models, it is important to investigate the number of model evaluations needed to fit linear and quadratic models in the full and the reduced ACTIVE SUBSPACE METHOD basis. To do so, a least mean squares method is applied with different samples sizes, M = {20, 50, 100, 500, 1000, 2000, 3000} and the Pearson's Correlation Coefficients are computed to provide an index of the degree of correlation between surrogate models fIm and the true response surface of the system f Im . These are computed using the formula of Eq. 5.15 and results are merged in Table . 6.4. In Fig. 6.6, a comparison between the approximated fIm and the true f Im response surfaces using linear and quadratic surrogate models, different sample sizes and two different input space is presented (in the full space and the reduced space from the Active Subspace method).

Model Type M=20 M=50 M=100 M=500 M=1000 M=2000 M=3000 The following observations can be made from Table . 6.4 and Fig. 6.6:

⋄ Linear models: Less than a hundred samples are not enough to approximate the growth rate variations when using the linear surrogate model L FD . However, the growth rate starts to be adequately approximated when tuning the model L 3D with only 50 samples (80%). Above a thousand samples, the predictions are enhanced but a lack of accuracy in the growth rate approximation is particularly noteworthy (between 82% and 84% with the model L 3D ).

⋄ Quadratic models: The least mean square regressions fail when the surrogate model Q FD is fitted with less than a hundred samples. Above a thousand samples, the growth rate variations are correctly approximated (92% with 1000 samples and 95% with 2000 samples). Meanwhile, when using the model Q 3D , these variations are quite well captured 2940 with only 50 samples (92%) and even better with a thousand samples (95% with 1,000 samples and 97% with both 2,000 and 3,000 samples).

Because of the non-linearities induced by the Flame Transfer Function, linear surrogate models cannot fully capture the response surface of the system within a relative error bound. For better accuracy, it is necessary to increase the complexity of the models by using the quadratic surrogate 2945 ACTIVE SUBSPACE METHOD models even if this implies tuning more coefficients. Obviously, it is expensive in high dimensions as 465 additional coefficients need to be tuned in the full-dimensional space but it is extremely beneficial in the reduced active subspace as there are only 6 additional coefficients to tune. Even better, when the quadratic model is spanned along the active directions, an accurate response surface is obtained when evaluating the model with only 50 samples as illustrated in Fig. 6.6. It highlights how reduced basis methods such as active subspace can lead to efficient Uncertainty Quantification strategies for high dimensional thermoacoustic problems.

Risk Factor estimation

Throughout this section, the Risk Factor of the first azimuthal mode of the combustor is investigated. To achieve this task, the following steps are followed:

1 At first, quadratic models Q 3D , that provides better correlation to the real response surface of the system with M={50, 2000} samples, are used. In the rest of the study Q 50 3D and Q 2000 3D will stand respectively for the reduced quadratic model fitted with M=50 samples and M=2000 samples. Factor estimated values are merged in Table . 6.5.

An overall good agreement is found when comparing the Risk Factor assessed from ATACA-MAC and surrogate models in Table . 6.5. Particularly, the low dimensional models are rather accurate in mimicking the actual response surface of the system. As expected, quadratic models provide better estimations of the Risk Factor than linear models. Globally, the Uncertainty Quan-2980 tification strategy adopted, which consists in combining a reduced basis technique and surrogate modelling approach, can be used to provide an accurate estimation of the modal Risk Factor in high dimensional thermoacoustic problems. ACTIVE SUBSPACE METHOD rate response, 10,000 Helmholtz simulations of the random inputs parameters n and τ were collected using the 1D Analytical Tool ATACAMAC. These random perturbations are generated using a uniform distribution. Therefore, the probability density function of the growth rate ω i is constructed and hence the modal Risk Factor of the system is approximated.

2 Then, the Active Subspace method is proposed as an interesting alternative towards the quantification of uncertainties in high dimensional problems. This technique, is based on the definition of a reduced basis able to catch most of the variation information of the system by exploiting the gradient of the growth rate with respect to the input parameters.

This gradient information is provided using Finite Difference discretization technique. The system dimensionality is reduced to only 3 variables.

3 Finally, linear and quadratic surrogate models are built over the full and the reduced spaces to approximate response surfaces of the problem. To appreciate the reliability and the accuracy of these models in predicting the Risk Factor of the system, a validation against the benchmark brute force Monte Carlo analysis is performed. The Risk Factor is accurately estimated when fitting a quadratic surrogate model based on only 3 active variables with only 50 ATACAMAC simulations (with a statistical error less than 6%).

This UQ method can be applied to other configurations and tools such as the 3D Helmholtz solvers AVSP instead of the ATACAMAC tool. Therefore, to avoid heavy gradient computation by finite difference method, gradient information can be obtained by using perturbative approaches such as Adjoint Sensitivity analysis procedure [START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF]). This is discussed in further details in Chapter 7. FOR THERMOACOUSTIC INSTABILITIES dient computations are realised by Finite Difference method and the corresponding results are confronted to the gradients obtained from the continuous adjoint approach.

Finally, concluding remarks and perspectives on the study are given in Section 7.6.

The adjoint method: Motivations

Computational Fluid Dynamics tools represent core elements in the design and development process of complex engineering devices. However, these techniques are expensive and time consuming specially for large-scale applications. Consequently, the direct calculation of uncertainties is unfeasible because the design under uncertainty may require the equivalent of many CFD computations. Therefore, the challenge is to approximate only the important physical phenomena of the system in a meaningful but tight CPU cost way. One method for overcoming the CPU limitation of high-fidelity computational models is to use surrogate based methods as discussed in Chapter 5. Nevertheless, surrogate models may not be able to faithfully represent some of the relevant features present in thermoacoustic systems. For example, the ATACAMAC model used in Chapter 4 can neither represent modes with a longitudinal component nor the effect of multi-perforated liners in complex geometries. Another challenge when dealing with realistic combustors is the presence of many swirler and associated flames, each of them being modeled by at least two uncertain parameters. In terms of UQ, this brings the curse of dimensionality into play. In order to break the curse, dimension reduction strategies, such as the Active Subspace methodology presented in Chapter 6, can be used to incorporate gradient information into reduced-order models thus extending their applicability for Uncertainty Quantification analysis.

Yet, gradients can be computed in a variety of ways. Traditional methods consist in using finite difference method that are relatively straightforward to implement, but at the expense of accuracy and far outweigh computational time to evaluate the model's output derivatives [START_REF] Bibliography Martins | The connection between the complex-step derivative approximation and algorithmic dierentiation[END_REF]). Such a way to compute the gradients was not an issue in Chapter 6 because simple 1D analytical network tool and algebraic models were employed. However, gradient computation by finite difference is a major bottleneck when dealing with more complex and parallel CFD solver such as LES or 3D acoustic code such as AVSP. FOR THERMOACOUSTIC INSTABILITIES

The use of adjoint methods was initially triggered in the late 1950's particularly in the framework of optimal control theory [START_REF] Lions | Laboratory investigation of an ultralow nox premixed combustion concept for industrial boilers[END_REF]). In the framework of fluid dynamics, gradient computations by adjoint-based methods were initially investigated by [START_REF] Pironneau | On optimum design in fluid mechanics[END_REF] who derived a continuous adjoint formulation of the Navier-Stokes equations. Numerous other studies were also conducted to perform sensitivity analysis towards aerodynamic design optimization.

Among them, one can cite the work of [START_REF] Jameson | Aerodynamic design via control theory[END_REF], [START_REF] Jameson | Optimum aerodynamic design using cfd and control theory[END_REF] who applied the adjoint Euler equations to transonic two-dimensional airfoils and Navier-Stokes equations to optimize a three-dimensional and aeronautical wing. Extensive studies in the same context are provided in [START_REF] Newman | Overview of sensitivity analysis and shape optimization for complex aerodynamic configurations[END_REF] and [START_REF] Giles | An introduction to the adjoint approach to design[END_REF].

Adjoint CFD solvers are still gaining in maturity in several scientific studies for the development of high-fidelity gradient-based optimization algorithms. Typically, they allow to get a broad insight on the variability of the system when all the model's input parameters are perturbed. There are two types of adjoint methods:

1 The continuous approach for which the adjoint equations are derived from the governing computational model and then subsequently discretized.

2 The discrete adjoint method for which the adjoint equation are directly derived from the discretized governing computational model. Discrete adjoint formulation, that are built on top of the discretized direct equation, should match exactly to the direct solutions. They would potentially be more suitable and accurate in the case of gradient estimations. Recall that the AVSP solver is an iterative, matrix-free solver because in the case of realistic problems, the matrix arising from the discretization of the Helmholtz equation may be very large (O(10 6 )) and storage becomes very undesirable for memory reasons. Therefore, developing a discrete adjoint algorithm in the AVSP solver would not be easy as the matrixvector products should be stored iteratively for gradient calculations purpose. As it would be hardly manageable to handle a discrete adjoint formulation because it would invert the operations in the differentiated code in a counterintuitive way, the continuous adjoint formulation is preferred in this study.

As discussed in Chapter 1, thermoacoustic oscillations occur due to feedback between heat FOR THERMOACOUSTIC INSTABILITIES release rate fluctuations and acoustic pressure fluctuations in confined spaces. These oscillations may lead to excessive vibrations, higher heat transfer to the walls and mechanical failures. The use of adjoint methods for gradient computation and sensitivity analysis of thermoacoustics allows to evaluate how all acoustic modes of the system would be potentially affected by any changes with respect to model's parameters. This is interesting for meaningful validation of computational models and prediction uncertainties. Recent studies of [START_REF] Magri | Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach[END_REF], [START_REF] Magri | A theoretical approach for passive control of thermoacoustic oscillations: Application to ducted flames[END_REF], Magri and Juniper (2013a)

Continuous adjoint approach the Helmholtz equation for thermoacoustic instabilities

In this section, we are interested in the continuous adjoint formulation for the Helmholtz equation Eq. 3.17 detailed in Section 3. At first, a brief explanation on the formulation of the problem is given in Section 7.3.1. For more mathematical details and functional analysis, refer to the work of Juniper and co-workers [START_REF] Juniper | The structural sensitivity of open shear flows calculated with a local stability analysis[END_REF], [START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF], [START_REF] Magri | Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach[END_REF], [START_REF] Magri | A theoretical approach for passive control of thermoacoustic oscillations: Application to ducted flames[END_REF], Magri and Juniper (2013a)).

Formulation of the problem

The direct Helmholtz problem Eq.(3.17) can be expressed as:

(V {ω, q ω }) p = 0 (7.1) FOR THERMOACOUSTIC INSTABILITIES where V is the matrix acting on the eigenfunction p; ω is one of the associated complex eigenvalue; q ω is the vector containing the parameters of the problem (geometrical parameters, n -τ parameters, speed of sound, ...).

The adjoint [START_REF] Griffiths | Consistent Quantum Theory[END_REF]) of the compact linear operator V , denoted V † , is the conjugate transpose of the operator V also called Hermitian adjoint to V . Similarly, the adjoint eigenfunction p † is the conjugate transpose of the operator p also called Hermitian adjoint to p. In an orthonormal basis, the adjoint eigenfunction p † and adjoint operator V † are obtained from that of p and V by complex conjugation and transposition with respect to the Hermitian inner product:

⟨ p † , (V {ω, q ω })p ⟩ = ⟨ (V {ω, q ω }) † p † , p⟩ , (7.2)
where ⟨, ⟩ is the inner product defined as:

⟨f, g⟩ = ∫ Ω f * g dΩ, ( 7.3)
for any functions f and g defined in the flow domain Ω. f * denotes the complex conjugate of f so that:

⟨f, g⟩ * = ⟨g, f ⟩ . (7.4)

In other words, the adjoint operator is defined through the following formula:

∫ Ω ( p † * (V {ω, q ω })p ) dΩ = ∫ Ω (( V {ω, q ω }) † p † ) * p) dΩ (7.5)
Finally, to find the adjoint operator relevant to the continuous formulation, integrations by parts of Eq. (7.5) need to be performed. As it will be made clear, the operators V and V † differ mainly because of the contribution of the flame and boundary conditions.

In the following, more focus is put upon the derivation of the adjoint Helmholtz equation with respect to the boundary conditions implemented in AVSP and detailed in Section 3. Later, sensitivity derivatives are screened to see how the coupling between the direct and the adjoint equations is achieved.

Derivation of adjoint Helmholtz equations

This section describes the derivation of the adjoint Helmholtz equation, its implementation and validation within the 3D Helmholtz solver AVSP. To achieve this task, the inner product of the Helmholtz equation and adjoint pressure is first formed:

⟨ p † (⃗ x), ∇. ( 1 ρ 0 (⃗ x) ∇p(⃗ x) ) + ω 2 γ(⃗ x)p 0 p(⃗ x) -iω γ(⃗ x) -1 γ(⃗ x)p 0 n(⃗ x)e iωτ (⃗ x) ∇p(⃗ x ref ).⃗ n ref ⟩ = 0 (7.6)
Which is also equivalent to:

3140 ⟨ p † (⃗ x), ∇.
( (

1 ρ 0 (⃗ x) ∇p(⃗ x) ) + ω 2 γ(⃗ x)p 0 p(⃗ x) ⟩    Term I = ⟨ p † (⃗ x), iω γ(⃗ x)-1 γ(⃗ x)p 0 n(⃗ x)e iωτ (⃗ x) ∇p(⃗ x ref ).⃗ n ref ⟩    Term II (
1 ρ 0 (⃗ x) ∇p(⃗ x) ) + ω 2 γ(⃗ x)p 0 p(⃗ x) ⟩    Term I :
When using the inner product definition of Eq. 7.3, Term I becomes:

∫ Ω p † * (⃗ x)∇ • ( 1 ρ 0 (⃗ x) ∇p(⃗ x) ) dΩ    A + ∫ Ω p † * (⃗ x) ω 2 γ(⃗ x)p 0 p(⃗ x) dΩ    B . (7.8)
When integrating by parts the first term of Eq. (7.8) labelled A , the following volume and surface integrals appear:

∫ Ω p † * (⃗ x)∇ • ( 1 ρ 0 (⃗ x) ∇p(⃗ x) ) dΩ    A = ∫ ∂Ω p † * (⃗ x) 1 ρ 0 (⃗ x) ∇p(⃗ x) • ⃗ n dS - ∫ Ω 1 ρ 0 ∇p(⃗ x) • ∇p † * (⃗ x) dΩ, (7.9)
Further integrating the second term of Eq. (7.9) leads to:

∫ Ω p † * (⃗ x)∇ • ( 1 ρ 0 ∇p(⃗ x) ) dΩ    A = ∫ ∂Ω p † * (⃗ x) 1 ρ 0 (⃗ x) ∇p(⃗ x) • ⃗ n dS - ∫ ∂Ω p(⃗ x) 1 ρ 0 (⃗ x) ∇p † * (⃗ x) • ⃗ n dS + ∫ Ω ∇ • ( 1 ρ 0 (⃗ x) ∇p † * (⃗ x) ) p(⃗ x) dΩ.
(7.10) By substituting Eq. ( 7.10) into Eq. (7.8), Term I is thus transformed into:

Term I = ∫ ∂Ω p † * (⃗ x) 1 ρ 0 (⃗ x) ∇p(⃗ x) • ⃗ n -p(⃗ x) 1 ρ 0 ∇p † * (⃗ x) • ⃗ n dS + ∫ Ω ∇ • ( 1 ρ 0 (⃗ x) ∇p † * (⃗ x) ) p(⃗ x) + p † * (⃗ x) ω 2 γ(⃗ x)p 0 p(⃗ x) dΩ.
(7.11)

By taking the complex conjugate of Eq. (7.11), the following expression is obtained:

Term I * = ∫ ∂Ω p † (⃗ x) 1 ρ 0 (⃗ x) ∇p * (⃗ x) • ⃗ n -p * (⃗ x) 1 ρ 0 (⃗ x) ∇p † (⃗ x) • ⃗ n dS + ∫ Ω ∇ • ( 1 ρ 0 (⃗ x) ∇p † (⃗ x) ) p * (⃗ x) + p † (⃗ x) ω * 2 γ(⃗ x)p 0 p * (⃗ x) dΩ.
(7.12)

The surface integral term in Eq. (7.12) automatically vanishes as soon as any combination of Neumann (∇p • ⃗ n = 0) and Dirichlet (p = 0) boundary condition is used for the direct 3155 Helmholtz problem.

When a complex impedance boundary is used:

Z = p(⃗ x) ρ 0 (⃗ x)c 0 û(⃗ x) • ⃗ n = iω p(⃗ x) c 0 (⃗ x)∇p(⃗ x) • ⃗ n , (7.13)
a proper boundary condition must be chosen for the adjoint problem in order to cancel the surface integral term of Eq. (7.12). This is typically the case when:

p † (⃗ x) ∇p † (⃗ x) • ⃗ n = p * (⃗ x) ∇p * (⃗ x) • ⃗ n = Z * c 0 -iω * . ( 7

.14) FOR THERMOACOUSTIC INSTABILITIES

In the case where Eq. (7.14) is selected as a boundary condition, Term I is thus such that:

Term I * = ⟨ p, ∇ • ( 1 ρ 0 (⃗ x) ∇p * (⃗ x) ) + ω * 2 γ(⃗ x)p 0 p * (⃗ x) ⟩ (7.15)
Note that Term I is self adjoint since the operator acting on p † is simply:

∇ • ( 1 ρ 0 (⃗ x) ∇p † (⃗ x) ) + ω * 2 γ(⃗ x)p 0 p † (⃗ x) . (7.16)
Due to the self-adjoint nature of the state equations, the adjoint equations have the same differential operators and the adjoint pulsation ω † is the complex conjugate of the direct pulsation ω * (ω † = ω * ). It constitutes a very important statement which both eases the derivation of adjoint equations and the validation of adjoint algorithms in the AVSP solver.

From Eq. (7.14), it also means that the proper boundary impedance for the adjoint problem is -Z * when Z is used for the direct problem.

2 Adjoint formulation for When incorporating Eq. (7.17) in Term II, one obtains:

⟨ p † (⃗ x), iω γ(⃗ x)-1 γ(⃗ x)p 0 n(⃗ x)e iωτ (⃗ x) ∇p(⃗ x ref ).⃗ n ref ⟩   
∇p(⃗ x ref ).⃗ n ref ≡ ∫ Ω ∇p(⃗ x) 1 σ √ π e ( - (⃗ x-⃗ x ref ) 2 σ 2 )    f G (⃗ x-⃗ x ref ) •⃗ n ref dΩ, ( 7
⟨ p † (⃗ x), iω γ(⃗ x) -1 γ(⃗ x)p 0 n(⃗ x)e iωτ (⃗ x) ∫ Ω f G (⃗ x -⃗ x ref )∇p(⃗ x) • ⃗ n ref dΩ ⟩ . (7.18)
For sake of simplicity, the term iω γ-1 γp 0 n(⃗ x)e iωτ (⃗ x) is noted F(ω) in the rest of the study.

Inverting the two integrals in Eq. ( 7.18) leads to:

3185 ∫ Ω ⟨ p † (⃗ x), F(ω) ⟩ f G (⃗ x -⃗ x ref )∇p(⃗ x) • ⃗ n ref dΩ. ( 7
∫ ∂Ω ⟨ p † (⃗ x), F(ω) ⟩ p(⃗ x)f G (⃗ x -⃗ x ref )⃗ n ref • ⃗ n dS    =0 - ∫ Ω ⟨ p † (⃗ x), F(ω) ⟩ p(⃗ x) ∇f G (⃗ x-⃗ x ref )⃗ n ref dΩ.
(7.20)

The surface integral term of Eq. (7.20) is zero as soon as the flame region does not reach the boundary ∂Ω, which is the case in practice. Therefore only the volume integral term of Eq. (7.20) remains.

3190

By using the inner product relation of Eq. ( 7.3), Eq. (7.20) reads:

Term II = - ∫ Ω ⟨ p † (⃗ x), F(ω) ⟩ p(⃗ x) ∇f G (⃗ x -⃗ x ref ) • ⃗ n ref dΩ = ⟨⟨ p † (⃗ x), F(ω) ⟩ * ∇f G (⃗ x -⃗ x ref ).⃗ n ref , p(⃗ x) ⟩ . (7.21)
Term II is thus such that:

Term II * = - ⟨ p(⃗ x), ⟨ p † , F(ω)) ⟩ * ∇f G (⃗ x -⃗ x ref ) • ⃗ n ref ⟩ (7.22)
When gathering Eq. 7.15 and Eq. 7.22, it follows that: without the flame effects for a chosen eigenpair (ω 0 , p0 ).

⟨ p(⃗ x), ∇. ( 1 ρ 0 (⃗ x) ∇p † (⃗ x) ) + ω * 2 γ(⃗ x)p 0 p † (⃗ x) ⟩ - ⟨ p(⃗ x), ⟨ p † (⃗ x), F(ω) ⟩ * ∇f G (⃗ x -⃗ x ref ).⃗ n ref ⟩ = 0. ( 7
-Find the adjoint eigenmode p † 0 by solving the discretized adjoint Helmholtz equation without flame coupling (Eq. (7.25)) for a chosen eigenpair (ω * 0 , p † 0 ).

2-Active Flame resolution:

-Set ω = ω 0 , ω * = ω * 0 and k = 1 to initiate the fixed point iteration algorithm (see Section 3.1.4) for the direct and the adjoint problems.

-Solve both Eq. (3.28) and Eq. (7.25) using the fixed point method, the k th iteration consisting in solving the following eigenproblem in ω k and ω * k defined as: (7.29) where θ = ω 2 . Following the approach of [START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF], when the matrix A is perturbed by δA, in which ||δA|| ≈ ϵ ≈ o(1), the shift in the converged eigenvalue ω K is given by: (7.30) where K is the number of fixed point iterations to reach the convergence and δA K = ∂A ∂ρ 0 = δA 0 denotes the perturbations on the mean density δρ 0 .

A p = θp,
δθ K = - ⟨ p † K , δA K pK ⟩ ⟨ p † K , pK ⟩ ≡ ⟨⟨p K , δA K ⟩⟩ ,
When accounting for both impedance boundary conditions and the flame/acoustic coupling, the discretized and unperturbed direct Helmholtz equation reads: where δT is the perturbation on the discretized matrix T containing the time delay τ in its diagonal. When using the fixed point iterative procedure, Eq. (7.34) becomes:

δL k = δA 0 + δB(θ k ) -(δN)Φ(θ k-1 )G -NΦ(θ k-1 )(δG) -iΦ(θ k-1 )NG [ 1 2 θ -1/2 k-1 Tδθ k-1 + θ 1/2 k-1 δT ] (7.35)
Finally, when substituting Eq. (7.35) in Eq. (7.32) and considering that K the number of the fixed point iterations to reach convergence, the shift in the converged eigenvalue θ K , is:

δθ K = ⟨⟨p K , δA 0 ⟩⟩ + ⟨⟨p K , δB K ⟩⟩ -⟨⟨p K , δNΦ(θ K-1 )G⟩⟩ -⟨⟨p K , NΦ(θ K-1 )(δG)⟩⟩ -iθ 1/2 K-1 ⟨⟨p K , δTΦ(θ K-1 )NG⟩⟩ -i 1 2 θ -1/2
K-1 ⟨⟨p K , Φ(θ K-1 )NTG⟩⟩ δθ K-1 (7.36) FOR THERMOACOUSTIC INSTABILITIES For convenience, the last term is denoted ξ K = i 2 θ -1/2 K-1 ⟨⟨p K , Φ(θ K-1 )NGT⟩⟩ and therefore the eigenvalue shift is: δθ K = ⟨⟨p K , δA 0 ⟩⟩ + ⟨⟨p K , δB K ⟩⟩ -⟨⟨p K , δNΦ(θ K-1 )G⟩⟩ -iθ 1/2 K-1 ⟨⟨p K , δTΦ(θ K-1 )NG⟩⟩ -⟨⟨p K , NΦ(θ K-1 )(δG)⟩⟩ + ξ K δθ K-1 (7.37) Eq. (7.37) is repeated until the right hand side contains δθ K-K , which is known to be zero.

Each step of the gradient iteration process implies one forward solution of the direct equation and one backward solution of the adjoint equation. Therefore, both eigenvalues and eigenvectors from the direct and adjoint equations must be stored at each iteration step of the point fixed algorithm.

The estimated initial conditions are then updated using the computed gradient direction. This process is not expensive since the gradient computations are completely independent of the number of input variables. The next step consists in comparing the gradients estimated by adjoint method with gradients calculated from a forward finite difference calculations.

Gradients calculations by finite difference method:

In order to measure the accuracy of the gradients computed by the adjoint method, a first order finite difference approximation of the growth rate function f Im is used:

∂f Im ∂x i = f Im (x i + δϵ i ) -f Im (x i ) δϵ i + O(δϵ i ), (1 ≤ i ≤ m) (7.38)
where δϵ i is the input parameter step perturbation, x i is the set of input parameter of the system and m is the number of input parameter. As it was mentioned in Section 7.2, the function f Im needs to be calculated once at point x i and further m times at f Im (x i + δϵ i ) for 1 ≤ i ≤ m. This results in m + 1 evaluations of the growth rate function f Im . Consequently, the computational effort for the gradient approximation using finite differentiating method is proportional to the number of input parameters. A sketch of the procedure to compute the gradients by finite difference approach is presented in Fig. and δn = δϵ × n [J/m] is applied. Note that the amplitude of the perturbation is δϵ = 1 × 10 -6 .

At first, only a perturbation on the flame time delay is applied and the gradients computed for the first acoustic mode of the systems are gathered in 7.9: Comparison between the gradients computed by adjoint method (∂f AD Im /∂τ ) and finite difference approximation (∂f F D Im /∂τ ). Only a perturbation on the time delay τ is taken into account with a step size δτ = δϵ × τ [s -1 ]. The amplitude of the perturbation is δϵ = 1 × 10 -6 .

The results prove satisfactory as the adjoint gradients are estimated within a reasonable error of 1%. Moreover, the gradients are well estimated when the time delay τ is varied over a period

Part IV

General conclusions

Since the 90's, there have been increasingly stringent regulations on pollutants emitted out of gas turbines. These have led engine manufacturers to operate combustors with lean premixed fuel and air thus allowing to control the temperature during the combustion process and hence the concentration of emissions. However, the major drawback in the use of lean premixed combustion is the emergence of thermoacoustic instabilities in gas turbine combustors. These instabilities occur because of the coupling between heat release rate and acoustic oscillations. They are frequently encountered in both aircraft and land-based power generation engines. The understanding and the control of this coupling phenomenon is key to the reliable and robust operation of gas turbine engines. Accounting for the uncertainties in the input parameters in any models for thermoacoustics is also required in order to reach a roubust prediction of the related instabilities.

In this thesis, we have provided a procedure to represent, characterize, and analyse the uncertainties for thermoacoustics to investigate and control the stability of gas turbine combustors.

Typically, we have developed and analysed computational strategies and algorithms based on both classical Uncertainty Quantification methods and model order reduction techniques, in order to improve the reliability of simulation-based analysis of gas turbine combustors. To convey a comprehensive understanding of the work achieved, generic conclusions and perspectives of further research and application possibilities are drawn in the following.

⋄ One objective of this thesis was to prepare the groundwork for an efficient development and implementation of Uncertainty Quantification methods (see in Chapter 4). This is necessary in order to solve high-dimensional thermoacoustic problems within an affordable computation time which remains an important requirement when performing Uncertainty Quantification analysis. A step-by-step methodology that bind Large Eddy Simulation Techniques, a Helmholtz solver and a quasi 1D analytical tool have been established to provide an estimate of the frequency and modal structures of two industrial helicopter engines (with N injectors and flames). The methodology is based on a model-fitting procedure that allows to represent easily the industrial geometry as a network of inter-connected acoustic elements by using the forward LES and Helmholtz solver solutions. This procedure proved satisfactory in predicting the stability characteristics and pulsating amplitudes of the industrial systems. Besides, thermoacoustic modes of the system were assessed with affordable computational effort without sacrificing the numerical accuracy.

⋄ The Uncertainty Quantification analysis of a mono-injector combustor with two uncertain parameters have been first investigated in Chapter 5. The thermoacoustic analysis of the mono-injector have been realised experimentally and numerically (with a 3D Helmholtz solver) in different settings. The comparison of the experimental and the numerical stability analysis appeared to be in good agreements except for three operating points that were expected to be more sensitive to the flame response to acoustic perturbations. To unravel the stability analysis of the systems, a continuous description of thermoacoustic modes has been adopted. This description is based on the definition of the modal Risk Factor that corresponds to the probability for a mode to be unstable given the uncertainties on the input parameters. To predict the modal Risk Factor of the geometries, a hybrid algorithm based on the «brute-force» Monte Carlo method and surrogate modelling techniques have been investigated. In particular, to reduce the computational cost in Monte Carlo Sampling that requires full solves of the underlying model, only a few Helmholtz simulations are used to fit the surrogate models. A Monte Carlo has been then applied on these surrogate models to provide an accurate estimate of the modal Risk Factor for each operating points.

A comparison between the Risk factor estimated by the Monte Carlo of the underlying model and the approximate Risk Factor obtained from the surrogate models show a good agreement. Although gaining further benefit on approximating the Risk Factor of the mode at low cost, the global error analysis has been conducted for more evaluation of the failure probability when using such algebraic surrogate models. The results have certified the efficiency and accuracy of the surrogate models in determining the Risk Factor of the system within a reasonable error, with remarkable applications in solving uncertainty quantification problems for thermoacoustics.

⋄ A large-scale and high-dimensional Uncertainty Quantification analysis have been conducted for two helicopter engines with N injectors and flames. To avoid heavy computational burden of the full system with LES techniques and Helmholtz solvers, the step-by-step methodology developed in Chapter 6 is harnessed. Thus, the Monte Carlo method is straightforwardly applied to provide an accurate estimate of the modal Risk Factor of the geometries.

To accelerate the Uncertainty Quantification analysis, a reduced basis method called «Active Subspace» is employed to reduce the full-dimensional space subspace to just a few. This technique detects the directions of the strongest variability using evaluations of the gradient and subsequently exploits these directions to construct a response surface on a low-dimensional subspace. In this work, the gradients were computed using Finite Differences approximations thus allowing to identify only 3 dominant directions (instead of the initial 2 × N directions), which are enough to describe the dynamics of the industrial systems. A posteriori analysis that combines the three dominant active variables and surrogate modelling techniques achieve a good computational performance in estimating the modal Risk Factor of the industrial systems. The latter is compared against the benchmark Risk Factor estimated from the «brute-force» Monte Carlo method and a good agreement was found. Besides, the global error analysis of the surrogate models was proved satisfactory thus highlighting the potential of the Active Subspace method to handle high dimensional Uncertainty Quantification problems.

⋄ In this work, another gradient-based method, namely the adjoint approach, has been investigated to deal with thermoacoustic problems when using a 3D Helmholtz solver (see Chapter 7). Adjoint methods are known to be computationally economical in providing accurate gradient estimations independently of the number of uncertain parameters of the system. In this work , the continuous adjoint Helmholtz equation has been developed and implemented (with respect to in-and outflow boundary conditions) in a Helmholtz solver for sensitivity analysis of the growth rate of the acoustic pressure disturbances. The treatment of high-dimensional and large scale thermoacoustic problems with adjoint method have not been addressed in this work. Its applicability requires the uptake and further robust developments to better handle the parallel processing of the 3D adjoint solver. Therefore, the implementation of the adjoint equation has been validated on different two-and threedimensional design problems. The growth rate gradients were evaluated with respect to the flame response parameters. The accuracy of the gradients evaluated by adjoint method was then validated against a first order Finite Difference approximation. Good agreements were found and it appears that less computational effort is required to evaluate the gradient by adjoint technique when perturbing the flame parameters. Moreover, the numerical convergence of the continuous adjoint equation is quicker for all the cases comparing to the direct equation resolution. In light of the results obtained, Uncertainty Quantification analysis using adjoint method is encouraging and albeit promising to handle more complex thermoacoustic systems with 3D Helmholtz solvers. code AVBP described in Chapter 4. Although avoiding performing expansive tests based on pressure and heat release records, performing Large Eddy Simulations provide interesting insight on the dynamics of turbulent flames and their interactions with the acoustic waves of the combustor.

The Large Eddy Simulations of the configuration was performed at Safran Helicopter Engines using the operating conditions displayed in Table A.1. To reduce uncertainties on boundary 3490 conditions the chamber casing is also simulated. The computational domain starts after the inlet diffuser and ends between the high-pressure stator and rotor. In this section, the flow is choked, allowing an accurate acoustic representation of the outlet. The results of the active flame computations are presented in Table . A.6. performed. These parameters are not reported in th manuscript for confidentiality reasons. As it was discussed in Chapter 4, ATACAMAC relies on a simplified description of the combustor geometry and therefore adjustment of some geometrical parameters has first to be performed to fit 3D results from AVSP. Such an adjustment is based on the objective of reproducing both the real and imaginary part of the targeted eigenmode for a number of imposed time delays, which is a key parameter for flame instabilities prediction. In practice, this is mostly done by slightly varying the burner length, since the burner (or injector) geometry is complex and its acoustic length is not easy to extract from a CAD. Typically, this adjustment is done based on the standard length correction in the low-frequency limit for a flanged tube ∆L i ≈ 0.4 √ 4S i π [START_REF] Silva | The behavior of subgrid-scale models near the turbulent/nonturbulent interface in jets[END_REF], [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF]). Consequently, the parametric analysis of the burner length suggests that a correction L * i = 9.45 × 10 -3 m should represent correctly the azimuthal mode of interest and match with the one computed with AVSP. The comparison is made in Table . A.7 and good agreements are found.

Finalement3

  , on examine la perspective d'utiliser la méthode adjointe pour analyser la sensibilité des systèmes thermo-acoustiques représentés par des solveurs 3D de Helmholtz. Les résultats obtenus sur des cas tests 2D et 3D sont prometteurs et suggèrent d'explorer davantage le potentiel de cette méthode dans le cas de problèmes thermo-acoustiques encore plus complexes. Mots clés: Instabilités thermoacoustiques, equation d'Helmholtz, Simulation numérique, Quantification d'incertitudes, Facteur de Risque, Monte-Carlo, Modèle de substitution, and phenomenology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Driving mechanisms of instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 1.3 About suppression methods of combustion instabilities . . . . . . . . . . . . . . . . 1.4 Tools to study combustion instabilities . . . . . . . . . . . . . . . . . . . . . . . . . 2 Uncertainty Quantification 2.1 Motivations and objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Literature survey and basic definitions . . . . . . . . . . . . . . . . . . . . . . . . . 135 2.3 State-of-the-art methodologies for Uncertainty Quantification analysis in CFD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 About Uncertainty Quantification in the framework of thermoacoustics . . . . . . . 2.5 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Structure of the manuscript . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140 II Low-order analysis tools for thermoacoustic instabilities in combustion CONTENTS chambers Helmholtz solvers and Network models 3.1 Thermoacoustic analysis using a Finite Volume Based Helmholtz Solver . . . . . . 3.1.1 Mathematical formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 3.1.2 The linear wave equation for reactive flows . . . . . . . . . . . . . . . . . . 3.1.3 Modelling of thermoacoustic instabilities using the Flame Transfer Function formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 The three-dimensional finite volume based acoustic solver AVSP . . . . . . 3.2 Analytical description of thermoacoustic instabilities in annular combustors with 150 network modelling techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Theoretical description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Thermoacoustic analysis of annular gas turbine combustion chambers 4.1 Towards the network modelling of industrial annular combustion chambers . . . . 4.2 Description of the 1 st annular combustor of interest with N burners . . . . . . . . .

  mode computations of the annular system with N injectors with Helmholtz solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Steady flame calculation of the full annular combustor with N injectors using the 3D Helmholtz solver AVSP . . . . . . . . . . . . . . . . . . . . . . 4.3.2 Active flame calculations of the full annular combustor with N injectors CONTENTS A.2.4 Active Flame computation of the system with N injectors using the 3D Helmholtz solver AVSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2.5 Acoustic computations using the quasi-analytical tool ATACAMAC . . . . 210 A.3 Uncertainty Quantification study . . . . . . . . . . . . . . . . . . . . . . . . . . . . side and after the instability on the right hand side). Picture b represents a damaged liquid-rocket engine after combustion instabilities. . . . . . . . . . . . . . 1.3 Feedback mechanism of combustion instabilities, inspired from Noiray et al. (2008). 1.4 Monitoring of pressure oscillations over the time in a combustion chamber (from Poinsot and Veynante (2011)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220 1.5 Typical result of the study of acoustic modes of a combustion chamber. Acoustic modes are considered to be stable when ω i < 0 (modes 2 , 3 , 5 and 7 in the bottom area in blue) and unstable when ω i > 0 (modes 1 , 4 , and 6 on the top area in red). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Uncertainty Quantification analysis: example of the PDF of model outcomes. The 225 risk associated to each part of the PDF is estimated (in %) to account for potential model deficiency or system failure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Conceptual view of the physical modelling process: from empirical observations to fine statistic analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 2.3 Uncertainty quantification analysis of thermoacoustic modes in a combustion cham-230 ber. Each mode belongs to an admissible region of the frequency plane with an associated Risk Factor to be unstable. . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Uncertainty Quantification using different set of thermoacoustic tools: cost evaluation with analytical tool, Helmholtz solver or LES techniques. . . . . . . . . . . . 3.1 Annular combustion chamber (right hand side picture, from combustor from Safran 235 Helicopter Engines and left hand side picture from CFM International). . . . . . . 3.2 Sketch of Crocco's flame model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Representation of the Flame/Acoustic coupling within a combustion chamber. The vector ⃗ u 1 represents the incoming force acoustic perturbation generated through the injector inlet, ⃗ x ref corresponds to the reference position where the velocity 240 fluctuations are measured, Q1 is the global heat release fluctuation integrated over the flame volume and q1 is the local heat release fluctuation per unit flame volume. 3.4 Representation of the fixed point algorithm implemented in AVSP solver. . . . . . 3.5 BC configuration to study azimuthal modes in annular combustor. . . . . . . . . . 3.6 PBC configuration to study azimuthal modes in annular combustor. . . . . . . . . 245 3.7 Representation of an annular combustion chamber connecting burners to an annular plenum. Because of the flame, the annular plenum and burners contain a fresh mixture characterized by a density ρ 0 u and sound speed c 0 u , whereas hot products with ρ 0 b and c 0 b are located in the combustion chamber. . . . . . . . . . . . . . . . . 3.8 H-junction: connections of each N sectors of the plenum to the combustion cham-250 ber through the i th burner. The analytical derivation by Bauerheim et al. (2014a) leads to four coupling parameters Γ i=1..4 . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Procedure for the fitting of real industrial annular combustion chamber and Uncertainty Quantification analysis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . (provided by Safran Helicopter Engines). . . . . . . . . . . . . . . . . . . . . . . . . 4.3 3D unstructured meshes for Helmholtz computation for the system with N injectors: the single sector on the left hand side and the full annular system on the right hand side. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Sound speed field c 0 (⃗ x) extracted from a LES time-average solution and used for 260 Helmholtz computations with AVSP solver. . . . . . . . . . . . . . . . . . . . . . . 4.5 Acoustic pressure field of the first azimuthal mode of the full annular helicopter combustion chamber with N injectors found from passive flame computation with AVSP solver.: f=612.0 Hz. The FTF parameters n and τ are set to 0. . . . . . 4.6 Frequencies and growth rates of acoustic modes with active flame (squares) and 265 modes with passive flame (diamonds). . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Map of stability for the first thermoacoustic mode of the system with N injectors in active flame regime with AVSP solver. The global value of the interaction index n is fixed to n=1486.43[J/m]. The time delay τ is varying over a period 64 × 10 -3 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 270 4.8 Structure of the first azimuthal mode of the full annular helicopter combustion chamber with N injectors found from active flame computation with AVSP solver. The global value of the interaction index n is n=1486.43[J/m] and the time delay τ is τ = 9.87 × 10 -4 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 Structure of the second azimuthal mode of the full annular helicopter combustion 275 chamber with N injectors found from active flame computation with AVSP solver. The global value of the interaction index n is n=1486.43[J/m] and the time delay τ is τ = 9.87 × 10 -4 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Sketch of the downstream chamber computed with AVSP solver to determine the first acoustic mode of the system with N injectors. . . . . . . . . . . . . . . . . . . LIST OF FIGURES 4.11 Sketch of the upstream plenum computed with AVSP solver to determine the first acoustic mode of the system with N injectors. The sector is duplicated N times to obtain the geometry full annular plenum. . . . . . . . . . . . . . . . . . . . . . . . 109 4.12 Map of stability of the first thermoacoustic mode of the combustor with N injectors: ATACAMAC computation (losanges) vs Helmholtz solver computation (squares)

  [s -1 ] . 116 4.14 Stability map of the first thermoacoustic mode of the combustor with N injectors: ATACAMAC computation (losanges) vs Helmholtz solver computation (squares) using the corrected length L * i = 0.231[m]. The global interaction index n is fixed, n=1486.43[J/m] (the Crocco's value is n = 6.57 for ATACAMAC computations), and τ is varying over a period T = 1 f 0 1 = 2Lc pc 0 ≈ 1.64 × 10 -3 s. . . . . . . . . . . . . 116 4.15 Stability map of the full annular helicopter combustor with N=N injectors when varying the interaction index n from 2.0 to 14 and the time delay τ over a period 0 ≈ 1.64 × 10 -3 s. WCC corresponds to the weakly coupled case modes chamber regime, WCP the weakly coupled case modes plenum area and SC represent the strongly coupled modes area. The latter case is never observed. . . . 118 5.1 The swirled combustor experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 5.2 Trace of the flame chemiluminescence in the symmetry plane of the burner. From Palies (2010) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 LIST OF FIGURES 5.3 Linearized stability prediction.The gray bounds indicate the marginally stable region defined by ∆α = ±10[s -1 ]. Empty symbols indicate agreement with experimental results while filled symbols represent partial agreement. From Silva et al. (2013). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  Monte Carlo with the 3D Helmholtz solver AVSP: sampling method workflow. 5.6 Convergence study of the risk factor ( in %); Top: Uniform distribution, Bottom: β-distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 (a) Monte Carlo results using M = 4000 Helmholtz-based thermoacoustic samples and a uniform distribution. (b) Histogram and Kernel density estimations of the growth rate. The Risk Factor is evaluated to 24%. . . . . . . . . . . . . . . . . . . 5.8 (a) Monte Carlo results using M = 4000 Helmholtz-based thermoacoustic samples and a β-distribution. (b) Histogram and Kernel density estimations of the growth rate. The Risk Factor is evaluated to 27.4%. . . . . . . . . . . . . . . . . . . . . . . 5.9 force Monte Carlo with the reduced-order model evaluation: sampling method workflow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 The least mean squares fitting of the geometry 11 Flame B. . . . . . . . . . . . . . 5.11 Histogram of the growth rate constructed with surrogate models. . . . . . . . . . . 5.12 Workflow for estimating the variability of the modal Risk Factor for a given size of Helmholtz samples randomly selected from the reference Monte Carlo AVSP database. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 Risk Factor estimated from a Monte Carlo analysis using the linear model LM F T F . . . . 5.14 Risk Factor estimated from a Monte Carlo analysis using the linear model QM F T F . . . . 5.15 Evolution of the standard deviation of the mean Risk Factor when using LM F T F (black) and QM F T F (red) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 5.16 On the left hand side: Evolution of the confidence interval of the mean Risk Factor when using LM F T F . On the right hand side: Evolution of the standard deviation of the mean Risk Factor when using LM F T F (black) and QM F T F (red). . . . . . . 5.17 The least mean squares fitting of the geometry 11 Flame A. . . . . . . . . . . . . . 5.18 (a) Uncertainty region for the first acoustic mode for a uniform PDF with 10% 335 uncertainty on the flame amplitude n and the flame time delay τ . (b) Histogram of the growth rate of acoustic disturbance for 100 Helmholtz samples computed using a Uniform PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.19 The least mean squares fitting of the geometry 11 Flame B using 15 Helmholtz samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 340 5.20 Histogram of the growth rate of acoustic disturbance for 15 Helmholtz samples computed using a Uniform PDF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Representation of the uniform distribution followed by the flame parameters: the plot on the top represents the PDF of the flame amplitude for the dimensionless ratio n

  375 model with 50 simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.1 Growth rate computed for each geometry in Table. 7.1 when accounting for a Gaussian formulation in the AVSP solver. The dotted black line represents the growth rate computed using the Dirac formulation with one reference point. When the standard deviation σ decreases, the growth rates are similar to those found with the Dirac formulation as expected. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Procedure to compute gradients by finite difference approach for m number of input parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  .11. . . . . . . . . . . . . . . . . . . . . . . . . . . A.1 Sketch of the full annular helicopter engine equipped with N injectors (provided by Safran Helicopter Engines). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 View of the temperature field (a) and the velocity and heat release fluctuations evolution over time from LES computations of Safran Helicopter Engines. . . . . . 395 A.3 Sound speed field c 0 (⃗ x) extracted from a LES time-average solution and used for Helmholtz computations of the helicopter engines using AVSP solver. . . . . . . . . A.4 3D unstructured meshes for LES and Helmholtz computation of the second helicopter engine investigated. View of the single sector. . . . . . . . . . . . . . . . . . A.5 Structure of the first longitudinal mode of the full annular helicopter combustion 400 chamber with N injectors found from passive flame computation: ω r = 495.55 Hz. A.6 Structure of the first azimuthal mode of the full annular helicopter combustion chamber with N injectors found from passive flame computation: ω r = 683.2 Hz. . A.7 Map of stability for the first thermoacoustic mode of the system with N injectors in active flame regime with AVSP solver. The global value of the interaction 405 index n is fixed to n=7612.0 J/m. The time delay τ is varying over a period 6 × 10 -3 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 LIST OF FIGURES A.8 Stability map of the first thermoacoustic mode of the combustor with N injectors: ATACAMAC computation (losanges) vs Helmholtz solver computation (squares) using the corrected length L * i = 9.45 × 10 -3 m. The Crocco's interaction index n 410 is fixed, n=3.92, and τ is varying over a period T = 1 1.6 × 10 -3 s. . . . A.9 Stability map of the full annular helicopter combustor with N injectors when varying the Crocco's interaction index n from 3.92 to 12 and the time delay τ 1.6 × 10 -3 s. Azimuthal plenum modes begin to change direction at n=5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 415 A.10 Monte Carlo analysis performed with ATACAMAC solver for the system with N injectors and flames. 10,000 samples were generated with a Uniform distribution. The Risk Factor is approximately 51%. . . . . . . . . . . . . . . . . . . . . . . . . . A.11 Eigenvalues of the finite difference approximation to the growth rate gradient of the full annular system with 2 × N uncertainties. The system is reduced to only a 3D space using 500 ATACAMAC calculations involving physical quantities associated with the Fourier transform of the Flame Transfer Function. . . . . . . . . . . . . . A.12 Comparisons between the approximated and the true response surfaces using the quadratic surrogate models Q 50 3D and Q 1000 3D . Different sample sizes are used for the fitting procedure: 50 and 1000 samples. . . . . . . . . . . . . . . . . . . . . . . . .

3. 1

 1 Azimuthal modes classification. From Wolf et al. (2012b). . . . . . . . . . . . . . . 3.2 Parameters in the mass conservation and momentum equations for a compressible viscous fluid, in absence of external forces (Eq. (3.1)). . . . . . . . . . . . . . . . . 3.3 Definition of the parameters of Eq. (3.31) that represents the fixed point algorithm 430[START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]). As studied by[START_REF] Miguel-Brebion | Joint numerical and experimental study of thermo-acoustic instabilities[END_REF], the relax parameter can be fixed or imposed dynamically to optimize the convergence process. . . . . . 4.1 Computational domains and grids used for LES and Helmholtz simulations. . . . . 4.2 Operating conditions for the LES and Helmholtz computations . . . . . . . . . . . 4.3 Frequency and decay rate of the first 9 eigenfrequencies of the 3D annular combus-435 tor with N injectors in passive flame regime. Computations realised with AVSP solver. All azimuthal modes are degenerate. . . . . . . . . . . . . . . . . . . . . . . 4.4 Frequency and decay rate of the first 9 eigenfrequencies of the 3D annular combustor with N injectors in active flame computation with AVSP solver. The global values n=1486.43[J/m] and τ = 9.87 × 10 -4 s where used to account for the flame 440 effects for AVSP computations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 First azimuthal frequencies computed with AVSP solver when the chamber and the annular plenum of the combustor are treated independently. . . . . . . . . . . LIST OF TABLES 4.6 Comparisons of the 1 st azimuthal mode computed in the full annular combustor (Section 4.3.1) and the one computed only in the annular chamber cavity. Computations are realised in steady flame regime with AVSP solver. . . . . . . . . . . . 110 4.7 Mean sound speed and radius used to determine analytically the first acoustic mode of the upstream plenum and the downstream chamber in passive flame regime.111 4.8 Parameters used for numerical applications of the annular system with N injectors.L c represents the half perimeter of the chamber and L p is the half perimeter of the plenum. L i stands for the initial burner length used for acoustic computation. S i is the section of the injector, S c stands for the section of the chamber and S p represents the section of the plenum. . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.9 First azimuthal chamber mode determined analytically (Eq. 4.1), with AVSP Helmholtz solver and the network model tool ATACAMAC. . . . . . . . . . . . . . 112 4.10 Eigenfrequency and growth rate of the first azimuthal mode of the system with N injectors: comparison between AVSP and ATACAMAC prediction for the Crocco's values n=6.57 and τ = 9.87 × 10 -4 s. . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.11 Eigenfrequency and growth rate of the first azimuthal mode of the system with N injectors: comparison between AVSP and ATACAMAC prediction. In this case the global interaction index is n=1486.43[J/m] and τ = 9.87 × 10 -4 s for the AVSP calculation and the Crocco's parameters n = 6.57 and τ = 9.87 × 10 -4 was used for ATACAMAC computations. The corrected length L * i = 0.231[m] was employed to determine the acoustic modes with ATACAMAC tool. . . . . . . . . . 115 4.12 Coupling parameters when increasing the interaction index from n=2 to n=14 but for a constant value of the time delay τ = 9.87 × 10 -4 [s -1 ]. . . . . . 119 5.1 Twelve different configurations explored: l 1 indicates the upstream manifold length and l 3 corresponds to the combustion chamber length. Dimensions are given in millimeters. From Silva et al. (2013). . . . . . . . . . . . . . . . . . . . . . . . . . . 124 LIST OF TABLES 5.2 Operating conditions used and eigenmodes computed using the 3D Helmholtz 470 solver AVSP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Linear stability analysis of flame A and flame B. Comparison between experimental and numerical results. (S)) Stable, (S/U)) Marginally stable/unstable, (U)) Unstable. The geometrical configurations C01 to C12 are defined in Table. 5.1. The three operating point with partial disagreement are highlighted. . . . . . . . . 475 5.4 ζ-coefficients computed for surrogate models LM n-τ , LM F T F and QM F T F using the 4000 samples of the Monte Carlo database. . . . . . . . . . . . . . . . . . . . . 5.5 Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. The sample size with the surrogate models is 4,000 samples. 5.6 Risk Factor and computation time estimated from from surrogate models. The 480 whole set of Helmholtz simulations (4000) were used. . . . . . . . . . . . . . . . . . 5.7 Risk Factors and their associated standard deviations computed by the Monte Carlo and surrogate models LM F T F and QM F T F using a different number of Helmholtz simulations from the full MC database. . . . . . . . . . . . . . . . . . . 5.8 ζ-coefficients determined for surrogate models LM F T F and QM F T F based on the 485 100 samples computed with AVSP code for the geometry 11 of Flame A. . . . . . . 5.9 Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Risk Factor estimated from surrogate models for the geometry 11 flame A. . . . . . 5.11 Risk Factors and their associated standard deviations computed by the Monte 490 Carlo surrogate models evaluations using LM F T F and QM F T F . 5 subsets of 10 Helmholtz samples each, randomly extracted from the full Helmholtz runs database, were used for the Risk Factor estimation. . . . . . . . . . . . . . . . . . . 5.14 Risk Factor estimated from surrogate models for the geometry 11 flame B using 15 Helmholtz samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 495 5.12 ζ-coefficients determined for surrogate models LM F T F and QM F T F based on the 15 samples computed with AVSP code for the geometry 11 of Flame B. . . . . . . 155 5.13 Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. 15 Helmholtz samples were used to provide these coefficients.155 5.15 Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. 100 and 1000 Helmholtz samples were used for the calculations.156 5.16 Risk Factor estimated from surrogate models for the geometry 11 flame B when using 100 and 1000 Helmholtz samples. . . . . . . . . . . . . . . . . . . . . . . . . 157 6.1 Mean Flame Transfer Function parameters considered in this study. . . . . . . . . 163 6.2 Summary of the surrogate models investigated to approximate the response surface of the annular combustor with N injectors. . . . . . . . . . . . . . . . . . . . . . . . 173 6.3 Summary of the number of coefficients for each surrogate model in the full-dimensional space and the 3D low-dimensional active subspace. . . . . . . . . . . . . . . . . . . 174 6.4 Pearson's correlation coefficients computed for surrogate models L FD , L 3D , Q FD and Q 3D using M = {20, 50, 100, 500, 1000, 2000} samples. The subscript «-» denotes the number of samples for which the Pearson's correlation coefficients cannot be computed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 6.5 Risk Factor estimated with the different surrogate models. These are compared to the Risk Factor determined from the benchmark Monte Carlo database (RF=84%). 179 7.1 Geometries investigated for the validation of the adjoint Helmholtz equation in the 3D Helmholtz solver AVSP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 7.2 Operating conditions of each of the geometries in Table.7.1 that are used to validate the implementation of the adjoint Helmholtz equation in the AVSP solver: l is the length of the geometry, h denotes the height and R is the radius of the Cylinder. The global interaction index is denoted n and τ stand for the flame time delay of the Flame Transfer Function. . . . . . . . . . . . . . . . . . . . . . . . . . 192 LIST OF TABLES 7.3 Standard deviations used to compute of the First acoustic modes of each of the geometry in Table.7.1 using the Gaussian formulation. . . . . . . . . . . . . . . . . 7.4 Pulsations and growth rates computed for the Dirac and the Gaussian formulation of the pressure gradient ∇p(⃗ x ref ).⃗ n ref in the AVSP solver. . . . . . . . . . . . . . . 525 7.5 Boundary conditions used to validate the implementation of the continuous adjoint Helmholtz equation in the AVSP solver: D denotes a Dirichlet boundary condition, N an Homogeneous Neumann and -Z * a complex impedance boundary condition. 7.6 Eigenmodes computed when solving the direct Helmholtz equation and the continuous adjoint Helmholtz equation in the AVSP solver. Homogeneous Neumann, 530 Dirichlet and complex impedance boundary conditions are used for the computations. Results proved satisfactory as the direct and adjoint eigenvalues should be complex conjugates of each other. The round off error is much less than 1% for the eigenvalues estimated. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.7 Growth rate derivatives ∂f Im /∂n and ∂f Im /∂τ computed for all the geometries 535 by finite difference when the amplitude of the perturbation δϵ is varied from from 1 × 10 -12 to 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.8 Definitions of the of the input parameter step perturbation used to compute the growth rate gradients by adjoint and Finite Difference approximation. . . . . . . . 7.9 Comparison between the gradients computed by adjoint method (∂f AD Im /∂τ ) and 540 finite difference approximation (∂f F D Im /∂τ ). Only a perturbation on the time delay τ is taken into account with a step size δτ = δϵ × τ [s -1 ]. The amplitude of the perturbation is δϵ = 1 × 10 -6 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.10 Comparison between the gradients computed by adjoint method (∂f AD Im /∂n) and finite difference approximation (∂f F D Im /∂n). Only a perturbation on the global 545 flame amplitude n is taken into account with a step size δn = δϵ×n. The amplitude of the perturbation is δϵ = 1 × 10 -6 . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :

 11 Figure 1.1: Examples of power systems where combustion instabilities can take place.

Figure 1 . 2 :

 12 Figure 1.2: Drawbacks of combustion instabilities. Picture a shows an injector system damaged after the instability in the combustor (before the instability on the left hand side and after the instability on the right hand side). Picture b represents a damaged liquid-rocket engine after combustion instabilities.

CHAPTER 1 .

 1 THE PHYSICS OF COMBUSTION INSTABILITIES ⋄ The reflection and the propagation of acoustic waves could perturb the flow back upward to the flame where the combustion process is taking place. Therefore, the flow may be 705 perturbed again and heat release fluctuations are re-generated. Generally, the feedback mechanism of combustion instability is described as closed instability loop, as displayed on Fig. 1.3.

Figure 1 . 3 :

 13 Figure 1.3: Feedback mechanism of combustion instabilities, inspired from Noiray et al. (2008).

Figure 1 . 4 :

 14 Figure 1.4: Monitoring of pressure oscillations over the time in a combustion chamber (from Poinsot and

  Fig. 1.4 where an instability is triggered at t=0. Initially in Regime I, linear oscillations of acoustic pressure appear (e.g. triggered CHAPTER 1. THE PHYSICS OF COMBUSTION INSTABILITIES by low-level turbulent fluctuations). Under favourable operating conditions, the amplitude of oscillations grows exponentially until reaching saturation. At this point, the combustion source terms overcome acoustic losses.

  , Palies (2010), Palies et al. (2010), Palies et al. (2011b), Palies et al. (2011c), Schuller et al. (2012)). Recently,

Figure 1 . 5 :

 15 Figure 1.5: Typical result of the study of acoustic modes of a combustion chamber. Acoustic modes are considered to be stable when ω i < 0 (modes 2 , 3 , 5 and 7 in the bottom area in blue) and unstable when ω i > 0 (modes 1 , 4 , and 6 on the top area in red).

Figure 2 . 1 :

 21 Figure 2.1: Uncertainty Quantification analysis: example of the PDF of model outcomes. The risk associated to each part of the PDF is estimated (in %) to account for potential model deficiency or system failure.

Figure 2 . 2 :

 22 Figure 2.2: Conceptual view of the physical modelling process: from empirical observations to fine statistic analysis.

⋄

  Sensitivity based approaches as Adjoint-based gradients techniques[START_REF] Putko | Approach for input uncertainty propagation and robust design in cfd using sensitivity derivatives[END_REF],[START_REF] Magri | Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach[END_REF],Magri and Juniper (2013a),[START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF].More information on these methods is given in the next chapters of the manuscript.⋄ The Monte Carlo method: Brute force Monte Carlo methodology is a widely used method for uncertainty analysis in multi-disciplinary applications. It is used to quantify the uncertainty on model outputs resulting from uncertainties on the model input parameters or input experimental data. Monte Carlo methods imply random sampling from the distributions of the uncertain inputs and the model is evaluated successively until a desired statistically significant distribution of outputs is obtained. Monte Carlo is conceptually simple and straightforward in term of implementation but requires a large number of model evaluations, e.g. large number of simulations, of the computational model to generate the output response surface of the system. It is inappropriate to full-scale complex applications because this would require a non-negligible parallel high performance computing. To overcome the issue, alternative methods such as Reduced Basis approaches (the Proper Orthogonal Decomposition, the Polynomial Chaos Expansion or the Active Subspace method etc.) can be used to decrease at first the dimensionality of the system. Numerous investigations have been conducted to reduce the number of Monte Carlo simulation runs effectively (Latin Hypercube sampling for example). In spite of the improved efficiency of the Monte Carlo methods, a well-established convergence criterion to complete the computations at a desired level of accuracy is still missing. Investigating on reduced-order techniques would help to determine the maximum number of simulations required to get an accurate estimate of the output quantities.

Figure 2 . 3 :

 23 Figure 2.3: Uncertainty quantification analysis of thermoacoustic modes in a combustion chamber. Each mode belongs to an admissible region of the frequency plane with an associated Risk Factor to be unstable.

Figure 2 . 4 :

 24 Figure 2.4: Uncertainty Quantification using different set of thermoacoustic tools: cost evaluation with analytical tool, Helmholtz solver or LES techniques.

⋄

  Part III: In this part, various Uncertainty Quantification methods are applied on a laboratory scale combustor (with only one injector and flame) as well as two industrial helicopter engines ( with N injectors and flames). The thermoacoustic analysis of the systems are conducted with an Helmholtz solver and a network modelling tool to determine eigenmodes of the geometries. The results suggests that the flame response plays an important role on the stability of the system and thus Uncertainty Quantification analysis on the flame model parameters would help to get more insight on the system behavior. CHAPTER 2. UNCERTAINTY QUANTIFICATION -Chapter 5 presents the Uncertainty Quantification study performed on the academic combustor to determine the probability for the first acoustic mode of the combustor to be unstable. The thermoacoustic analysis of the system is conducted with an Helmholtz solver and Monte Carlo methods and surrogate modelling techniques are combined for Uncertainty Quantification analysis purposes. Although reducing drastically the number of state computations, it is shown that algebraic surrogate models are efficient in providing accurate estimate of the modal risk factor.-Chapter 6 tackle the Uncertainty Quantification of the annular helicopter engines with N injectors and flames (The second helicopter engine is treated in Appendix A). A quasi 1D analytical tool is used for both the thermoacoustic and the Uncertainty Quantification of the problem. At first, the dimensionality of the system is reduced using the Active Subspace methodology (from 2×N uncertainties (2 uncertain parameters per flame (n and τ )) to only 3). Then, the Uncertainty Quantification study is conducted with appropriate surrogate models that are based only on the active variables assessed from the Active Subspace approach. The results proved satisfactory when comparing to a forward Monte Carlo analysis.-Chapter 7 focuses on the application of adjoint method for thermoacoustic problems. A derivation of the adjoint Helmholtz equation using a continuous adjoint approach is presented. The implementation aspects on a 3D Helmholtz solver and the validation on twoand three-dimensional test cases are shown. The results obtained are promising and open the perspective of further exploring the potential of adjoint method for the Uncertainty Quantification of thermoacoustic problems. ⋄ Part IV: This part proposed further discussions and the future perspectives of this work. the multi-physics phenomena involved in combustion instabilities is very challenging. Generally, the methodology used to study the system behavior is highly dependent on the combustor design complexity. Low order tools and theories on simplified geometries (Sensiau (2008), Bauerheim et al. (2014a), Bauerheim et al. (2014b), Bauerheim et al. (2014b), Mensah and Moeck (2015), Parmentier et al. (2012), Salas (2013)) have been spread out and turned out

Figure 3 . 1 :

 31 Figure 3.1: Annular combustion chamber (right hand side picture, from combustor from Safran Helicopter Engines and left hand side picture from CFM International).

  ,[START_REF] Benoit | Numerical assessment of thermo-acoustic instabilities in gas turbines[END_REF],[START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF],[START_REF] Sensiau | Simulations numériques des instabilités thermoacoustiques dans les chambres de combustion aéronautiques[END_REF]).2 A low order tool called ATACAMAC developed by CERFACS based on geometry simplifications is used to capture only azimuthal modes in annular configurations(Bauerheim et al. (2014a),[START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF],[START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF]). The outcomes of ATA-CAMAC solver are then used to extract phenomenological analysis of the results of AVSP code.

(

  [START_REF] Wang | Les prediction of wall-pressure fluctuations and noise of a low-speed airfoil[END_REF],[START_REF] Shin | Flame wrinkle destruction processes in harmonically forced, laminar premixed flames[END_REF]), effects of the type of velocity perturbation impinging the flames(Schuller et al. (2002a),[START_REF] Schuller | A unified model for the prediction of laminar flame transfer functions: comparisons between conical and v-flames dynamics[END_REF]), non-linearities effects(Schuller et al. (2002a),[START_REF] Preetham | Dynamics of laminar premixed flames forced by harmonic velocity disturbances[END_REF]), multiple flame effects(Duchaine and Poinsot (2011),[START_REF] Kornilov | Experimental assessment of the acoustic response of laminar premixed bunsen flames[END_REF]).

  where V f is the flame volume and δ f stands for the flame thickness. The vector ⃗ u 1 denotes the velocity vector of the main flow which feeds the flame. The global parameter n, also called the interaction index, measures the amplitude of the flame response to acoustic perturbations and the global parameter τ corresponds to the phase time lag between acoustic perturbation (at an upstream reference point ⃗ x ref ) and the flame response.

Figure 3 . 2 :

 32 Figure 3.2: Sketch of Crocco's flame model.

  a) Global formulation of the FTF. (b) Local formulation of the FTF.

Figure 3 . 3 :

 33 Figure 3.3: Representation of the Flame/Acoustic coupling within a combustion chamber. The vector ⃗ u 1 represents the incoming force acoustic perturbation generated through the injector inlet, ⃗ x ref corresponds to the reference position where the velocity fluctuations are measured, Q1 is the global heat release fluctuation integrated over the flame volume and q1 is the local heat release fluctuation per unit flame volume.

  containing the impedances when using a Robin type boundary condition (this term is null when setting either Neumann or Dirichlet boundary conditions). The matrix F includes the discretization of the right hand side term of the Helmholtz equation representing the flame/acoustic coupling in such a way that: F p = (NΦG)p, (3.29) where N is the matrix containing the flame amplitude n(⃗ x) at each grid point, Φ contains the exponential e iωτ (⃗ x) and the matrix G includes the gradient of the pressure measured at the reference point and along the reference direction ⃗ n ref : ∇p(⃗ x ref ) • ⃗ n ref . Therefore, the system features discrete non-linear eigenpair (ω, p(⃗ x)) for which ω represents an eigenfrequency and p(⃗ x)

Figure 3 . 4 :

 34 Figure 3.4: Representation of the fixed point algorithm implemented in AVSP solver.

,

  [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationnary computational fluid dynamics[END_REF]). Recently, a methodology to incorporate the effect of non-purely acoustic mechanisms into Helmholtz solvers has been developed by[START_REF] Ni | Accounting for acoustic damping in a helmholtz solver accounting for acoustic damping in a helmholtz solver[END_REF] with transfer matrices measured from experiments and large-eddy simulation. The use of network models allows to investigate different processes that are related to the coupling between acoustic cavities, input uncertainties or even symmetry breaking effects. It offers the opportunity to capture the leading mechanisms affecting the modes nature and to get 85 CHAPTER 3. HELMHOLTZ SOLVERS AND NETWORK MODELS an insight to control them at the early design stage. Parmentier et al. (2012) developed a 1D Analytical Tool used to Analyze and Control Azimuthal Modes in Annular Chambers. This tool is based on the linearized acoustic equations with a steady and uniform azimuthal mean flow. This technique is efficient in representing analytically azimuthal eigenmodes in a BC type configuration (Burner + Chamber configuration) connected by several injectors (see Fig. 3.5).

Figure 3 . 5 :

 35 Figure 3.5: BC configuration to study azimuthal modes in annular combustor.

  3.6) that delivers the air. Further studies on PBC type configuration (Plenum + Burner + Chamber configuration, Fig. 3.6), were performed and proved effective solutions in mimicking industrial annular combustors behaviour

Figure 3 . 6 :

 36 Figure 3.6: PBC configuration to study azimuthal modes in annular combustor.

Figure 3 . 7 :

 37 Figure 3.7: Representation of an annular combustion chamber connecting burners to an annular plenum. Because of the flame, the annular plenum and burners contain a fresh mixture characterized by a density ρ 0 u and sound speed c 0 u , whereas hot products with ρ 0 b and c 0 b are located in the combustion chamber.

  33, N corresponds to the number of sectors, the perimeter of the annular combustion chamber and the annular casing are respectively noted 2L c = 2πR c and 2L p = 2πR p . The wave numbers in the cold and hot gases reads k u = ω/c 0 u and k b = ω/c 0 b each.

Figure 3 . 8 :

 38 Figure 3.8: H-junction: connections of each N sectors of the plenum to the combustion chamber through the i th burner. The analytical derivation by Bauerheim et al. (2014a) leads to four coupling parameters Γ i=1..4 .

  by an H-junction (O'Connor and Lieuwen (2012b), O'Connor and Lieuwen (2012a), Blimbaum et al. (2012)) as shown in Fig. 3.8. Consequently, the pressure p ′ and the velocity u ′ in the chamber are related to those in the plenum. Based on jump conditions (Dowling (1995), Bauerheim et al. (2014b)), the acoustic propagation in the burner described by Eq. (3.32), and a n -τ

  37) where I d is the 4-by-4 identity matrix. This dispersion relation (3.37) is non-linear in ω. Numerical solvers can efficiently solve Eq (3.37) (Newton Raphson algorithm, say), but explicit expressions are still useful to understand key mechanisms controlling combustion instabilities. The ANR methodology differs according to respective symmetrical aspects of the combustor: ⋄ Axisymmetric annular combustors: in this case, all sectors and flames are identical. In the analytical model, all matrices [R i ] and [T i ] are similar (the subscript i can be ommited) thus leading to the following explicit dispersion relation: det({[R][T ]} N -I d ) = 0. This equation can be recast asN ∏ p=1 det([R][T ] -e j2pπ/N I d ) = 0 ⇔ det([R][T ] -e j2pπ/N I d ) = 0 for p = 1..N (3.38)This simplification highlights that in axisymmetric configurations, each sector has the same acoustic behavior: the stability of the system can be deduced by considering only one sector (matrix [R][T ]) which necessarily acts as a pure phase-lag 2pπ/N , where p corresponds physically to the azimuthal order.

2

  The coupling factors are not null but satisfy |Γ k=1..4 | ≪ 1; In this case, solutions are close to the fully decoupled case. Consequently, they can be searched as f c = f 0 c + δf and f p = f 0 p + δf . A Taylor expansion of the dispersion relation yields the solutions in the case where the two annular cavities are not naturally coupled, viz. namely when f 0 p and f 0 c are not multiple of each other:

  experimentation and computational modelling of unsteady reactive flow open new opportunities to get insight about the physical phenomena relevant to engineering applications. Even though there are still numerous open theoretical questionings related to numerical approaches for thermoacoustic instabilities, the computation cost related to numerical tools remains one of the major roadblocks. This chapter is preparing the groundwork for the development of Uncertainty Quantification methods for large-scale systems within a reasonable numerical timeframe. The overall process is sketched in Fig. 4.1 and consists in establishing the connectivity between Large Eddy simulation techniques and low-order modelling approaches described in Chapter 3 with the aim to provide a way to identify pure acoustic eigenmodes in complex geometries: 1 LES solver: At first, Large Eddy simulations are performed to retrieve the mean sound speed of the system and the local fields of the Flame Transfer Function parameters n local (x) TURBINE COMBUSTION CHAMBERS and τ local (x).

Figure 4 . 1 :

 41 Figure 4.1: Procedure for the fitting of real industrial annular combustion chamber and Uncertainty Quantification analysis.

Figure 4 . 2 :

 42 Figure 4.2: Schematic representation of the full annular helicopter engine fed by N injectors (provided by Safran Helicopter Engines).

Figure 4 . 3 :

 43 Figure 4.3: 3D unstructured meshes for Helmholtz computation for the system with N injectors: the single sector on the left hand side and the full annular system on the right hand side.

Figure 4 . 4 :

 44 Figure 4.4: Sound speed field c 0 (⃗ x) extracted from a LES time-average solution and used for Helmholtz computations with AVSP solver.

  Acoustic pressure evolution over the azimuthal angle.

Figure 4 . 5 :

 45 Figure 4.5: Acoustic pressure field of the first azimuthal mode of the full annular helicopter combustion chamber with N injectors found from passive flame computation with AVSP solver.: f=612.0 Hz. The FTF parameters n and τ are set to 0.

Table 4. 4 :

 4 Frequency and decay rate of the first 9 eigenfrequencies of the 3D annular combustor with N injectors in active flame computation with AVSP solver. The global values n=1486.43[J/m] and τ = 9.87 × 10 -4 s where used to account for the flame effects for AVSP computations.

Figure 4 . 6 :

 46 Figure 4.6: Frequencies and growth rates of acoustic modes with active flame (squares) and modes with passive flame (diamonds).

Figure 4 . 7 :

 47 Figure 4.7: Map of stability for the first thermoacoustic mode of the system with N injectors in active flame regime with AVSP solver. The global value of the interaction index n is fixed to n=1486.43[J/m]. The time delay τ is varying over a period T = 1 f 0 1

  × 10 -3 s. As shown in the stability chart, accounting here for the Flame/Acoustic coupling has destabilizing effects on the first azimuthal mode of interest. In active flame regime, the two first eigenmodes computed exhibit almost the same frequencies TURBINE COMBUSTION CHAMBERS and growth rate of the acoustic pressure perturbations. Their structure are hereinafter investigated and shown in Fig.4.8 and Fig.4.9.(a) Modulus of the acoustic pressure. (b) Phase of the acoustic pressure.(c) Reconstruction of the pressure field.

  Acoustic pressure evolution over the azimuthal angle.

Figure 4 . 8 :

 48 Figure 4.8: Structure of the first azimuthal mode of the full annular helicopter combustion chamber with N injectors found from active flame computation with AVSP solver. The global value of the interaction index n is n=1486.43[J/m] and the time delay τ is τ = 9.87 × 10 -4 s.

  Acoustic pressure evolution over the azimuthal angle.

Figure 4 . 9 :

 49 Figure 4.9: Structure of the second azimuthal mode of the full annular helicopter combustion chamber with N injectors found from active flame computation with AVSP solver. The global value of the interaction index n is n=1486.43[J/m] and the time delay τ is τ = 9.87 × 10 -4 s.

  et al. (2012b) based on the analytical method of[START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF]. The latter method was devoted to the study of 1D acoustic waves propagating in annular combustion chambers connected to several burners. In the formulation of[START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF], the network-based TURBINE COMBUSTION CHAMBERS model does not account for an upstream annular plenum, and thus does not truly represent the design of real gas turbine combustors.In the present work, the network model to study thermoacoustic oscillations in real industrial combustors introduced byBauerheim et al. (2014a),[START_REF] Bauerheim | An analytical model for azimuthal thermoacoustic modes in an annular chamber fed by an annular plenum[END_REF] and named ATA-CAMAC is used. As detailed in Section 3.2, this methodology allows to reduce the size of the full scale acoustic problem to a simple 4-by-4 matrix containing all information of the resonant modes combustor. Therefore, explicit dispersion relations for Plenum + Burner + Chamber configurations are obtained and exact forms of the coupling parameters for azimuthal modes between

  the plenum cavities independently with the 3D Helmholtz solver AVSP. The meshes used for thermoacoustic computations of the downstream annular plenum and the upstream annular chamber are displayed in see Fig 4.10 and Fig 4.11.(a) Annular chamber. (b) Single sector of the chamber.

Figure 4 . 10 :

 410 Figure 4.10: Sketch of the downstream chamber computed with AVSP solver to determine the first acoustic mode of the system with N injectors.

( a )

 a Single sector of the plenum.

Figure 4 . 11 :

 411 Figure 4.11: Sketch of the upstream plenum computed with AVSP solver to determine the first acoustic mode of the system with N injectors. The sector is duplicated N times to obtain the geometry full annular plenum.

  (4.1) is used to determine the radius of the chamber and the plenum of the combustor, see Fig. 3.7). Fields of the mean sound speed in the plenum c u 0 and the chamber c b 0 are directly extracted from AVSP computations. Results are presented in

2080CAMAC.

  This stability chart is built by varying the time delay τ = τ = 9.87 × 10 -4 s over a period of the first azimuthal mode T = 1 Growth rate vs Time delay.

Figure 4 . 12 :

 412 Figure 4.12: Map of stability of the first thermoacoustic mode of the combustor with N injectors: ATACAMAC computation (losanges) vs Helmholtz solver computation (squares) with the initial burner length L i =0.125[m]. In this case, the Crocco's value n=6.57 is fixed and the time delay τ is varying over

  4.11. A good agreement is found when comparing to the first azimuthal mode computed with AVSP code. Eigenfrequency and growth rate of the first azimuthal mode of the system with N injectors: comparison between AVSP and ATACAMAC prediction. In this case the global interaction index is n=1486.43[J/m] and τ = 9.87 × 10 -4 s for the AVSP calculation and the Crocco's parameters n = 6.57 and τ = 9.87 × 10 -4 was used for ATACAMAC computations. The corrected length L * i = 0.231[m] was employed to determine the acoustic modes with ATACAMAC tool.Moreover, the stability map of the system has been studied by fixing the value of the interaction index n=6.57 and by varying the time delay τ over a period a period T = 1 1.64 × 10 -3 s. Results are presented in Fig. A.8 and good trends of the growth rate variations are predicted by the analytical model ATACAMAC.

Figure 4 . 13 :

 413 Figure 4.13: The approximate estimate of the burner parameter length L i for predicting the growth rate of the 1 st azimuthal mode of the system with N injectors. The Flame/Acoustic interactions are considered for analytical computation purpose. In this case n = 6.57 and τ = 9.87 × 10 -4 s and the growth rate ω i is 8.81[s -1 ]

Figure 4 . 14 :

 414 Figure 4.14: Stability map of the first thermoacoustic mode of the combustor with N injectors: AT-ACAMAC computation (losanges) vs Helmholtz solver computation (squares) using the corrected length

Figure 4 . 15 :

 415 Figure 4.15: Stability map of the full annular helicopter combustor with N=N injectors when varying the interaction index n from 2.0 to 14 and the time delay τ over a period T = 1 f 0 1 = 2Lcpc 0 ≈ 1.64 × 10 -3 s. WCC corresponds to the weakly coupled case modes chamber regime, WCP the weakly coupled case modes plenum area and SC represent the strongly coupled modes area. The latter case is never observed.

Figure 5 . 1 :

 51 Figure 5.1: The swirled combustor experiment.

l 3

 3 corresponds to the combustion chamber length. Dimensions are given in millimeters. From[START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF]. STABILIZED COMBUSTOR EXPERIMENT To measure the flame response, a loudspeaker is placed at the back end of the system. Moreover, two experimental conditions corresponding to two different air flow rates were experimentally tested corresponding to flames A and B, with larger power in the latter ( QA = 1.94kW ) than in the former ( QB = 3.03kW ). These two operating points have the same equivalence ratio equal to 0.7 but with different bulk flow velocities in the injector equal to ūb = 2.67ms -1 for the flame A and ūb = 4.13ms -1 for the flame B.

Figure 5 . 2 :

 52 Figure 5.2: Trace of the flame chemiluminescence in the symmetry plane of the burner. From Palies (2010)

  computations of[START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] have been redone for Uncertainty Quantification purpose. The operating conditions used and the numerical results are presented in Table.5.2. The stability map of all thermoacoustic modes of the geometries studied is presented in Fig.5.3 and the global comparative study between the experimental and numerical stability results of[START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] is summed up in Table.5.3.

Figure 5 . 3 :

 53 Figure 5.3: Linearized stability prediction.The gray bounds indicate the marginally stable region defined by ∆α = ±10[s -1 ]. Empty symbols indicate agreement with experimental results while filled symbols represent partial agreement. From Silva et al. (2013).

  α+ζ) denotes the Beta function, Γ(.) is the Gamma function, and α and 2275 ζ are two free parameters. Note that f ζ Y is only defined for a reduced random variable Y on [0, 1]. The parameters α and ζ which characterize the β-PDF are deduced from the desired mean µ ζ Y and variance ν ζ Y of this reduced variable Y:

  .7). If this range is chosen equal to the range of the previous uniform PDF (e.g. R ζ = R U = 10%) then the ζ-distribution degenerates to the previous uniform PDF. Consequently, the range R ζ is an additional free parameter. For this study, this range is fixed to R ζ = 30% leading to the characteristic values α = ζ = 2.87.[⋄]

Figure 5 . 5 :

 55 Figure 5.5: force Monte Carlo with the 3D Helmholtz solver AVSP: sampling method workflow.

Figure 5 . 6 :

 56 Figure 5.6: Convergence study of the risk factor ( in %); Top: Uniform distribution, Bottom: βdistribution

Figure 5

 5 Figure 5.7: (a) Monte Carlo results using M = 4000 Helmholtz-based thermoacoustic samples and a uniform distribution. (b) Histogram and Kernel density estimations of the growth rate. The Risk Factor is evaluated to 24%.

Figure 5 . 8 :

 58 Figure 5.8: (a) Monte Carlo results using M = 4000 Helmholtz-based thermoacoustic samples and a β-distribution. (b) Histogram and Kernel density estimations of the growth rate. The Risk Factor is evaluated to 27.4%.

Figure 5 . 9 :

 59 Figure 5.9: force Monte Carlo with the reduced-order model evaluation: sampling method workflow.

  τ ) is implicit and non-linear. To speed up the Uncertainty Quantification analysis, it is worth investigating if this response surface designed from the full Monte Carlo database in section .5.3.1 can be estimated by explicit surrogate models. Linear and quadratic models based on the uncertainties on the Flame Transfer Function parameters n and τ are investigated: 134 CHAPTER 5. UNCERTAINTY QUANTIFICATION OF A SWIRLED STABILIZED COMBUSTOR EXPERIMENT 1 LM n-τ : a linear model based on the parameters n and τ of the Flame Transfer function :

Figure 5 . 10 :

 510 Figure 5.10: The least mean squares fitting of the geometry 11 Flame B.

Figure 5 . 11 :

 511 Figure 5.11: Histogram of the growth rate constructed with surrogate models.

Figure 5 . 12 :23Figure 5 . 13 :Figure 5 . 14 :

 512513514 Figure 5.12: Workflow for estimating the variability of the modal Risk Factor for a given size of Helmholtz samples randomly selected from the reference Monte Carlo AVSP database.

Figure 5 . 15 :Figure 5 . 16 :

 515516 Figure 5.15: Evolution of the standard deviation of the mean Risk Factor when using LM F T F (black) and QM F T F (red)

Figure 5 . 17 :CHAPTER

 517 Figure 5.17: The least mean squares fitting of the geometry 11 Flame A.

  Fig. 5.18 and the Risk Factor estimated are presented in Table. 5.10.

Figure 5 . 19 :

 519 Figure 5.19: The least mean squares fitting of the geometry 11 Flame B using 15 Helmholtz samples.

Figure 5 . 20 :

 520 Figure 5.20: Histogram of the growth rate of acoustic disturbance for 15 Helmholtz samples computed using a Uniform PDF.

  the configuration 11 of the flame A and 11 of the flame B. ACTIVE SUBSPACE METHOD methodologies have been focused on forward uncertainty propagation, including Monte Carlo methods, adaptive sparse and Generalized Polynomial Chaos for Galerkin and collocation formulations or even Active Subspace methods. However all these techniques become hardly implementable on high fidelity CFD solvers for very large scale systems: ⋄ As discussed in Section. 2.3, Polynomial Chaos Expansion models are expensive to derive unless the number of terms in the expansion is moderate, which requires a relatively small number of uncertain variables and a low degree of expansion.

Figure 6 . 1 :

 61 Figure 6.1: Representation of the uniform distribution followed by the flame parameters: the plot on the top represents the PDF of the flame amplitude for the dimensionless ratio n n (where n is the nominal value of n) and the plot on the bottom represents the PDF of the time delay for the dimensionless ratio τ

  2760

Figure 6 . 2 :

 62 Figure 6.2: Monte Carlo analysis performed with ATACAMAC solver using for 10,000 samples generated with a Uniform distribution.

Figure 6 . 3 :

 63 Figure 6.3: Uncertainty Quantification strategy applied to the real annular helicopter engine with 2 × Nuncertain parameters (2 uncertain parameters per flame (n and τ )). Initially, the Active Subspace method is employed to reduce the system dimensionality to only 3 variables. Then, algebraic surrogate models for the complete and reduced probabilistic spaces are used to analyse the surface response of the system.

Figure 6 . 5 :

 65 Figure 6.5: Eigenvalues of the finite difference approximation to the growth rate gradient of the full annular system with 2 × N uncertainties. Finally, only the 3 first active variables are relevant when using the theoretical studies of Noiray et al. (2011) and Bauerheim et al. (2014a). Convergence analysis with different samples are used to converge eigenvalues: M = {10, 20, 50, 100, 500, 1000} samples. The spectrum is associated to {Re(γ), Im(γ)}.

Figure 6 . 6 :

 66 Figure 6.6: Comparisons between the approximated fIm and the true f Im response surfaces using linear and quadratic surrogate models, different sample sizes and two different input space is presented: in the full-dimensional space and the reduced space (3 active variables determined from the Active Subspace method).

2Figure 6 . 7 :

 67 Figure 6.7: Convergence of the the low dimensional active subspace model when it is fitted with M=50 samples (Q 50 3D ) and M=2000 samples (Q 2000 3D). These models are replayed 100000 times to evaluate the Risk Factor variability when comparing to the Risk Factor obtained from the benchmark Monte Carlo database (RF M C ). An overall good agreement is found with a relative good error below 6% when fitting the model with 50 simulations.

1

  7.7) Term I and Term II of Eq. (7.7) are investigated by taking into account the following boundary conditions implemented in the Helmholtz solver AVSP: ⋄ Dirichlet boundary condition (see Eq. (3.25)). ⋄ Homogeneous Neumann boundary condition (see Eq. (3.26)). ⋄ Robin boundary condition (see Eq. (3.27)). 3145 Adjoint formulation for ⟨ p † (⃗ x), ∇.

  [START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF], the right hand side term of Eq. (7.7), labelled Term II, needs to be derived carefully to avoid extreme sensitivity at the reference point, where the acoustic velocity is measured. To make the adjoint problem well posed, the Dirac distribution δ(⃗ x -⃗ x ref ) which is used to generate ∇p(⃗ x ref ).⃗ n ref from the pressure field gradient is regularized as a Gaussian distribution noted f G (⃗ x -⃗ x ref ). Indeed, the eigenvalue is extremely sensitive to the velocity eigenfunction at the reference point thus affecting the numerical resolution of the adjoint problem. The above heat release model is therefore approximated as:

  Geometries investigated for the validation of the adjoint Helmholtz equation in the 3D Helmholtz solver AVSP.

Figure 7 .Figure 7 . 1 :

 771 Figure 7.1 presents the evolution of the growth rates when decreasing the standard deviation

  7.5. Only the first acoustic eigenmodes of each geometry will be targeted in this work. FOR THERMOACOUSTIC INSTABILITIES -Find the direct eigenmode p0 by solving the discretized direct Helmholtz equation (Eq. (3.28))

A

  p + B(ω k-1 )p + ω 2 k p = F (ω k-1 )p (7.27) A † p † + B † (ω * k-1 )p † + ω * 2 k p † = F † (ω * k-1 )p † (7.28) -Iterate on k until |ω k -ω k-1 | < tol and |ω * k -ω * k-1 | < tol,where tol is the tolerance desired.The gradient calculations are realised in a post processing step by starting from the discretized and unperturbed direct Helmholtz equation without impedances and flame effects:

  (A + B -(NΦG)) p = θp. (7.31) FOR THERMOACOUSTIC INSTABILITIES where N is a diagonal matrix containing the flame amplitude n(⃗ x) at each grid point, Φ includes the exponential e iωT in which T is the diagonal matrix containing the time delay τ (⃗ x) and the matrix G contains the gradient of the pressure measured at the reference point and along the reference direction ⃗ n ref : f G (⃗ x -⃗ x ref ). When defining L = A + B(θ) -NΦ(θ)G, a perturbation on the matrix L by δL, in which ||δL|| ≈ ϵ ≈ o(1), leads to the following eigenvalue drift: δL = δA 0 + δB(θ) -[(δN)Φ(θ)G + N(δΦ(θ))G + NΦ(θ)(δG)] , (7.33) where δA 0 , δB, δN, δΦ and δG are respectively the perturbation of the discretized matrices A, B, N, Φ and G.Further developing Eq. (7.33) leads to: δL = δA 0 + δB -

  Growth rate derivatives ∂f Im /∂n and ∂f Im /∂τ computed for all the geometries by finite difference when the amplitude of the perturbation δϵ is varied from from 1 × 10 -12 to 1Definitions of the of the input parameter step perturbation used to compute the growth rate gradients by adjoint and Finite Difference approximation. method for a posteriori comparison. As for the previous analysis with finite difference method, the global flame amplitude n and the time delay τ are perturbed. An increment of δτ = δϵ×τ [s -1 ]

Figure 7 . 3 :

 73 Figure 7.3: Comparison between eigenvalues obtained for the first acoustic mode of each geometry using the adjoint method (∂f AD Im /∂τ ) and finite differentiation (∂f F D Im /∂τ ). A step size δτ = δϵ × τ with (δϵ = 1 × 10 -6 ) is used. The growth rate gradients are evaluated when τ varies over a period T= 1 f0 .

Figure 7 . 4 :

 74 Figure 7.4: Comparison between eigenvalues obtained for the first acoustic mode of each geometry using the adjoint method (∂f AD Im /∂n) and finite differentiation (∂f F D Im /∂n). A step size δn = δϵ × n is used for which δϵ = 1 × 10 -6 . The global flame amplitude n is varied as reported in Table.7.11.

  7.11. Very good agreements are found for all the cases when the global flame amplitude n and the 3335 time delay τ are independently varied. The results shows that the growth rate gradients are more sensitive to the perturbations on the time delay τ than to the global flame amplitude n. Further investigations are then conducted by simultaneously increasing the value of the global 205 CHAPTER 7. ON THE APPLICATION OF THE ADJOINT METHOD FOR THERMOACOUSTIC INSTABILITIES flame amplitude n while varying the time delay τ over a period T= 1 f 0 . The results are presented in Fig. 7.12.

  Operating conditions for the LES computation of the annular system with N injectors. annular plenum. However the stability of these modes remains unclear and this is the reason why active flame computations are conducted to get insight on the system behaviour. and the second (683.2 Hz) predominant mode of the combustor observed in the passive flame computations are targeted. To achieve this, the local fields of the flame parameters n(⃗ x) and τ (⃗ x) are first extracted from the single sector pulsated LES. The values used in the acoustic calculations with AVSP are gathered in Table.AValues for the flame interaction n and the time delay τ used to compute eigenmodes of the annular system with N injectors in active flame regime.

  ℜ(ω) Hz ℑ(ω)[s -1 ] ℜ(ω) Hz ℑ(ω)[rad/s] Frequency and decay rate of the first and the second acoustic modes computed in the active flame regime for the second helicopter engine investigated. The global values n=7612.0 J/m and τ = 1.46 × 10 -3 s where used to account for the flame effects for the AVSP computations. The first longitudinal mode at 495.0 Hz identified in the passive flame computations is found unstable thus confirming the previous LES observations. Active flame computations also show an unstable acoustic activity of the predominant azimuthal mode of the combustor. To further investigate the stability of this azimuthal mode, additional thermoacoustic calculations of the full-scale geometry are realised with AVSP by varying the flame time delay τ over a period 3555 6 × 10 -3 s. The objective is to evaluate the variation of the azimuthal mode frequency and growth rate. The results are presented in Fig. A.7.

  Eigenfrequency and growth rate of the first azimuthal mode of the system with N injectors: comparison between AVSP and ATACAMAC prediction. In this case the global flame amplitude n=7612.0[J/m] (the Crocco's interaction index being n = 3.92) and τ = 1.47 × 10 -4 s. The corrected length L * i = 9.45 × 10 -3 m was used to determine the acoustic modes with ATACAMAC tool.Moreover, the stability analysis of the system has been conducted using the ATACAMAC tool by varying the time delay τ over a period T = 1 6 × 10 -3 s. The results are then compared to the forward stability analysis conducted with the AVSP solver in Fig. A.7 and the growth variations are well represented. Growth rate vs Time delay.

Figure A. 8 :

 8 Figure A.8: Stability map of the first thermoacoustic mode of the combustor with N injectors: ATA-CAMAC computation (losanges) vs Helmholtz solver computation (squares) using the corrected length L * i = 9.45 × 10 -3 m. The Crocco's interaction index n is fixed, n=3.92, and τ is varying over a period T = 1 f 0 1 = 2Lc pc 0 ≈ 1.6 × 10 -3 s.

  1 , Y 2 , ..., Y k }. The model response denoted Z

	is computed using the input data of the vector Y in such a way that Z=f(Y). The Uncertainty
	Quantification analysis of the model f (Y) starts by generating random perturbation of the in-
	put quantities using a well suited Probability Density Function. Then, to construct uncertainty
	bounds for the model response Z, a sampling method is used to propagate input uncertainties
	through the model (for example Monte Carlo). That is to say, instead of looking for a single
	result by running the physical model f(Y) only once, Uncertainty Quantification explores the
	range of findings provided by running the same model multiple times, each time with different
	set of values of its corresponding key input parameters distributions. This leads to a probabilistic
	representation of the output Z thus providing the other alternative and plausible scenarii of the
	phenomena represented by f(Y). The statistical representation of Z is then interpreted to account
	for risk in quantitative analysis as it is shown in Fig. 2.1.

State-of-the-art methodologies for Uncertainty Quantifica- tion analysis in CFD simulations

	Uncertainty Quantification increases the reliability and robustness of high-fidelity CFD simula-
	tion of industrial systems by accounting for variability in operating conditions. Common input
	factors of these variability are transient forcing functions, boundary conditions, stated assump-
	tions, chemical kinetics aspects, parametric uncertainties (simplification of the geometry and/or
	limitation of the domain studied, leading edge, blade shapes, roughness, etc.), no-modelled phys-
	ical processes or forms of the physical models (e.g. turbulence modelled as an extra diffusivity),
	turbulence modelling uncertainties, etc. Uncertain inputs may also be theoretically constant or
	follow known relationships but may have some inherent uncertainty. These factors may vary
	in large, tractable but unknown ways and this is even more cumbersome to handle when deal-
	ing with realistic applications. Consequently, to quantitatively measure the effects of the above
	model uncertainties in CFD simulations, the use of efficient computational methods for Uncer-
	tainty Quantification analysis is required.
	Let's recall the mathematical model f(Y) defined earlier in Section 2.2. Denoting Y={Y 1 , Y 2 , ..., Y k }
	the vector containing the uncertain input parameters of the system and Z the output response of
	the model f(Y), the Uncertainty Quantification analysis is realised as follows:
	1 The joint Probability Density Function of the vector Y is quantified by using a discretized
	random process to generate random perturbations of the input parameters Y 1 , Y 2 , ..., Y k .

Table 3

 3

.1: Azimuthal modes classification. From

[START_REF] Wolf | Acoustic and large eddy simulation studies of azimuthal modes in annular combustion chambers[END_REF]

.

  Table. 3.2.

	Quantity	Definition	Units
	ρ	Density	[kg/m 3 ]
	u	Velocity vector	[m/s]
	p	Pressure	[Pa]
	q	Volumetric heat release [W/m 3 ]
	r	Perfect gas constant:

Table 3 .

 3

2: Parameters in the mass conservation and momentum equations for a compressible viscous fluid, in absence of external forces (Eq. (3.1)).

Domain Number of nodes Number of tetrahedral cells

	Single Sector LES	518 649	2 819 176
	Full Annular LES	7 694 265	42 287 640
	Single Sector Helmholtz	69 019	336 135
	Full Annular Helmholtz	1 010 370	5 042 025

Table 4 .

 4

1: Computational domains and grids used for LES and Helmholtz simulations.

  Table 4.3.

	1935
		Steady Flame
	Mode Number ℜ(ω) Hz ℑ(ω)[s -1 ]	Mode description
	1.	612.0	0.0	1 st Azimuthal mode
	2.	612.0	0.0	2 nd Azimuthal mode
	3.	849.8	0.0	1 st Longitudinal mode
	4.	1147.3	0.0	3 rd Azimuthal mode
	5.	1147.3	0.0	4 th Azimuthal mode
	6.	1312.2	0.0	5 th Azimuthal mode
	7.	1312.2	0.0	6 th Azimuthal mode
	8.	1597.5	0.0	7 th Azimuthal mode
	9.	1597.5	0.0	8 th Azimuthal mode
	Table 4.3: Frequency and decay rate of the first 9 eigenfrequencies of the 3D annular combustor with
	N injectors in passive flame regime. Computations realised with AVSP solver. All azimuthal modes are
	degenerate.

Table .

 .

	1.	622.2	8.8	612.0	0.0	1 st Azimuthal mode
	2.	623.3	7.3	612.0	0.0	2 nd Azimuthal mode
	3.	848.0	3.5	849.8	0.0	1 st Longitudinal mode
	4.	1137.1	-11.5	1147.3	0.0	3 rd Azimuthal mode
	5.	1139.7	-10.7	1147.3	0.0	4 th Azimuthal mode
	6.	1313.3	-2.2	1312.2	0.0	5 th Azimuthal mode
	7.	1313.3	-2.0	1312.2	0.0	6 th Azimuthal mode
	8.	1598.8	3.6	1597.5	0.0	7 th Azimuthal mode
	9.	1599.2	5.5	1597.5	0.0	8 th Azimuthal mode
						4.4.

Table .

 .

	CHAPTER 4. THERMOACOUSTIC ANALYSIS OF ANNULAR GAS
	TURBINE COMBUSTION CHAMBERS
	First Chamber mode:f 1C 0	First Plenum mode:f 1P 0
	AVSP frequency Hz	614.17	413.45
	Table 4.5: First azimuthal frequencies computed with AVSP solver when the chamber and the annular
	plenum of the combustor are treated independently.
	Moreover, the 1 st azimuthal acoustic mode assessed in the chamber cavity is very close to
	the one determined in the steady flame computation in Section 4.3.1 for the full annular
	combustor, see Table. 4.6.
	1 st azimuthal mode computed with AVSP
	Full annular combustor	612.0
	Annular chamber	614.17
			4.5.
	109

Table 4.6: Comparisons of the 1 st azimuthal mode computed in the full annular combustor (Section 4.3.1) and the one computed only in the annular chamber cavity. Computations are realised in steady flame regime with AVSP solver.

  Table.4.7.

	Chamber
	Half perimeter	L i /L c	0.20
	Section	S i /S c	2.06 ×10 -3
	Plenum
	Half perimeter	L i /L p	0.20
	Section	S i /S p	1.50 ×10 -2
	Fresh gases
	Mean Density	ρ u 0	4.92 [kg/m 3 ]
	Mean Sound Speed	c u 0	480.43 [m/s]
	Hot gases
	Mean Density	ρ b 0	2.08 [kg/m 3 ]
	Mean Sound Speed	c b 0	706.72 [m/s]
	Flame parameters
	Crocco's interaction index	n	6.57
	Time delay	τ	varying
	Flame Thickness	δ f	1 ×10 -3

  As shown in Table.4.9, when using the parameters in Table.4.8, the targeted azimuthal chamber mode of the combustor is very well estimated. The next step consists in further investigating both the system behaviour when taking into account the flame effects and the coupling between the chamber and the plenum cavities with ATACAMAC.

	TURBINE COMBUSTION CHAMBERS
	2075
	2 Step 2: Computations with ATACAMAC are performed using the operating conditions of
	Tab 4.8 in active flame regime. Results are presented in Table. 4.10 and compared against
	the first azimuthal mode computed with AVSP in Section 4.3.2.
	3D Helmholtz solver result (AVSP) 1D Model Result (ATACAMAC)
	622.24+8.81i		567.98-12.85i
	Analytical result Hz AVSP Hz ATACAMAC Hz
	1 st chamber mode f 1C 0	614.17	614.17	614.17
	1 st plenum mode f 1P 0	413.15	413.15	413.15
	Table 4.9: First azimuthal chamber mode determined analytically (Eq. 4.1), with AVSP Helmholtz solver and the
	network model tool ATACAMAC.

2070

industrial combustion chamber. Results are presented in Table

.

4.9. Table 4.10: Eigenfrequency and growth rate of the first azimuthal mode of the system with N injectors: comparison between AVSP and ATACAMAC prediction for the Crocco's values n=6.57 and τ = 9.87 × 10 -4 s.

Table 4 .12:

 4 Coupling parameters when increasing the interaction index from n=2 to n=14 but for a constant value of the time delay τ = 9.87 × 10 -4 [s -1 ].

		92	7.15	6.13	-1.65
	2	-1.92	7.15	-4.79	1.28
	4	-1.92	7.15	-1.57 × 10 1	4.23
	6.57	-1.92	7.15	-2.26	7.12
	8	-1.91	7.15	-3.76	1.01
	10	-1.91	7.15	-4.85 × 10 -2	1.30
	12	-1.91	7.15	-5.94 × 10 1	1.60 × 10 1
	14	-1.91	7.15i	-7.04 × 10 1	1.89 × 10 1
			PLENUM
	0	2.26	7.25	6.21	1.93
	2	2.26	7.25	-4.86	-1.52
	4	2.26	7.25	-1.59	-4.98
	6.57	2.26	7.25	-3.01 × 10 1	-9.42
	8	2.27	7.25	-3.80 × 10 1	-1.19 ×
	10	2.27	7.25	-4.91 × 10 1	-1.53
	12	2.27	7.25	-6.02 × 10 1	-1.88 ×
	14	2.27	7.25	-7.13 × 10 1	-2.22 ×

Table 5 . 1 :

 51 Twelve different configurations explored: l 1 indicates the upstream manifold length and

Table 5 . 2 :

 52 Operating conditions used and eigenmodes computed using the 3D Helmholtz solver AVSP.

	126

Table 5 .3: Linear

 5

	Case	Flame A				Flame B
		C01	C02 C03 C04	C01	C02 C03 C04
	Experiment	S	S	S	U	S	S	S-U	U
	Simulation	S	S	S	U	S	S	S-U	U
		C05	C06 C07 C08	C05	C06 C07 C08
	Experiment	S	S	S-U	U	S	S	S	UU
	Simulation	S	S	S-U	U	S	S	S-U	U
		C09	C10 C11 C12	C09	C10 C11 C12
	Experiment	S	S	S-U	U	S	S	S-U	U
	Simulation	S	S	U	U	S	S	S	U

stability analysis of flame A and flame B. Comparison between experimental and numerical results. (S)) Stable, (S/U)) Marginally stable/unstable, (U)) Unstable. The geometrical configurations C01 to C12 are defined in Table. 5.1. The three operating point with partial disagreement are highlighted.

Table 5 . 5 :

 55 Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. The sample size with the surrogate models is 4,000 samples.

Table 5 .

 5 Table 5.14. The Risk Factors computed are null in this case. 14: Risk Factor estimated from surrogate models for the geometry 11 flame B using 15 Helmholtz samples. STABILIZED COMBUSTOR EXPERIMENT

	Surrogate model Risk Factor in %
	LM F T F	0
	QM F T F	0

Table 5 .

 5 12: ζ-coefficients determined for surrogate models LM F T F and QM F T F based on the 15 samples computed with AVSP code for the geometry 11 of Flame B.

Table 5 .

 5 13: Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. 15 Helmholtz samples were used to provide these coefficients.

Table 5 .

 5 Table. 5.15 and the Risk factors computed are shown in Table. 5.16. 15: Correlation coefficients of the surrogate models and the full Monte Carlo database computed from AVSP. 100 and 1000 Helmholtz samples were used for the calculations.

			Correlation
	Models N = 100 Samples N = 1000 Samples
	LM F T F	96.89%	97.39%
	QM F T F	97.13%	98.56%

of coefficients Linear models L FD 31 L 3D 4 Quadratic models Q FD 496 Q 3D 10Table 6 . 3 :

 63 Table. 6.3. Summary of the number of coefficients for each surrogate model in the full-dimensional space and the 3D low-dimensional active subspace.

	Model Type	Number

Table 6 .

 6

	Linear models
	L FD	-	0.63	0.78	0.78	0.81	0.81	0.81
	L 3D	-	0.80	0.80	0.82	0.82	0.83	0.84
	Quadratic models
	Q FD	-	-	-	-	0.92	0.95	0.95
	Q 3D	0.83	0.92	0.92	0.92	0.95	0.97	0.97

4: Pearson's correlation coefficients computed for surrogate models L FD , L 3D , Q FD and Q 3D using M = {20, 50, 100, 500, 1000, 2000} samples. The subscript «-» denotes the number of samples for which the Pearson's correlation coefficients cannot be computed.

  have proved how adjoint sensitivity analysis can be efficiently applied to an electrical heated Rijke tube by taking into account the effect of the mean-flow temperature jump in the acoustics. Later,[START_REF] Juniper | Sensitivity analysis of thermoacoustic eigenproblems with adjoint methods[END_REF] presented two different methods for Uncertainty Quantification of thermoacoustic instabilities for nonlinear Helmholtz eigenvalue problems. The methods allow to compute gradients a thousand times faster than finite difference methods. Based on this, the present study is initiated to enhance and complement the Uncertainty Quantification analysis performed in Chapter 6. The objective is to speed up the gradients computations using adjoint methods when the AVSP solver is used to model the thermoacoustics instead of the 1D analytical network tool ATACAMAC used in Chapter 4.

  .17) where σ is the standard deviation of the Kernel and ⃗ x ref stands for the nominal coordinates of the reference point. Note that when σ goes to zero, f G tends to the Dirac distribution and the integral in Eq. (7.17) is exactly ∇p(⃗ x ref ).⃗ n ref . Otherwise, with a finite value of FOR THERMOACOUSTIC INSTABILITIES σ, it is a regularized version of this quantity, more suitable for further developments and numerical implementation.

Table 7 . 2 :

 72 Operating conditions of each of the geometries in Table.7.1 that are used to validate the implementation of the adjoint Helmholtz equation in the AVSP solver: l is the length of the geometry, h denotes the height and R is the radius of the Cylinder. The global interaction index is denoted n and τ stand for the flame time delay of the Flame Transfer Function.As a first step, the implementation of the Gaussian formulation to measure the pressuregradient f G (⃗ x-⃗ x ref ) inthe AVSP solver is investigated for each geometry. The standard deviation 3210 of the Gaussian function used to compute eigenmodes is presented in Table. 7.3. FOR THERMOACOUSTIC INSTABILITIES

	Geometry	σ[m]
	Simple 2D Tube	1.3 × 10 -2
	2D mono-injector 1.0 × 10 -2
	3D Cylinder	1.4 × 10 -1

Table 7 . 3 :

 73 Standard deviations used to compute of the First acoustic modes of each of the geometry in Table.7.1 using the Gaussian formulation.The first acoustic modes computed for each geometry using the Gaussian formulation are summarized inTable. 7.4; (i) ωr and ωi stands for the growth rate obtained by the Dirac formulation, (ii) ω r and ω i are those obtained with the Gaussian one.

		Dirac Formulation Gaussian Formulation
	Geometry	ωr [Hz]	ωi [s -1 ]	ω r [Hz]	ω i [s -1 ]
	Simple 2D Tube	342.2	+0.6	355.2	+3.4
	2D Mono-injector 2802.4	+4.4	2802.4	+5.3
	3D Cylinder	2632.9	-0.1	2633.0	-1.9

Table 7 . 4 :

 74 Pulsations and growth rates computed for the Dirac and the Gaussian formulation of the pressure gradient ∇p(⃗ x ref ).⃗ n ref in the AVSP solver.

Table

  Table. 7.9:

		Adjoint	Finite difference	Error
		Real part Imag. part Real part Imag. part Relative
	Simple 2D Tube	2997.58	1129.23	2970.10	1140.18	0.9%
	2D mono-injector	1.42	0.83	1.34	0.84	5%
	3D Cylinder	1.59	-1.39	1.61	-1.31	3%

3325

Perturbation on

τ : δτ = δϵ × τ [s -1 ]

(a) The single sector for LES computation. (b) Full annular system computed with Helmholtz solver.

Step 1: At first, no flame effects are considered (Steady flame computation). Having access to mesh generation data of the industrial combustor, the annular chamber and the annular plenum cavities are decoupled. This allows to compute acoustic modes in the chamber and

(a) Response surface of the first acoustic mode. (b) Histogram of the growth rate of acoustic disturbance.

(a) The surrogate model LMn-τ (b) The surrogate model LMF T F (c) The surrogate model QMF T F

(a) Response surface of the growth rate. (b) Histogram of the growth rate.

Algebraic surrogate models are built in the full dimension space and over the low-dimensional subspace.

(a) M = {10, 20, 50} samples (b) M = {100, 500, 1000} samples

of the Gaussian function.

Acknowledgements

Campassens and Chantal Nasri. The work presented here was supported by the European Project called UMRIDA (Uncertainty Management for Robust Industrial Design in Aeronautics).

List of Figures

Then, the study is performed for different operating conditions in Section 5.4.1 and Section 5.4.2.

Discussions and conclusions are given in section 5.5.

Experimental set-up description

The laboratory-scale experiment used in this study corresponds to a single swirled stabilized 127 CHAPTER 5. UNCERTAINTY QUANTIFICATION OF A SWIRLED STABILIZED COMBUSTOR EXPERIMENT ble) for each operating condition: C11 for the flame A, C07 and C11 for the flame B. By doing this, it is expected to explain the disagreement found between the experimental and the numerical stability analysis by the lack of knowledge on the flame input parameters n and τ . These parameters have generally an important impact on the stability prediction of thermoacoustic systems.

Uncertainty Quantification inquiries will begin with the case 07 Flame B then the geometry 11

Flame A and finally the geometry 11 Flame B.

Test case 1: Configuration 07-Flame B

Monte Carlo analysis with 3D Helmholtz solver

At first, the range of uncertainty for the flame parameters n and τ is investigated by collecting quantitative data from two independent experimentalists groups at EM2C (Paris) and IMFT (Toulouse). From these datasets, a 10% uncertainties on both n and τ parameters was selected:

This range of uncertainty is applied to the following nominal experimental value of the global value of the interaction index n = 1079J/m and τ = 4.73ms. Also, the type of distribution followed by the FTF parameters is not known and it is necessary to make sure that the shape of the PDF has only a limited impact on the computed Risk Factor value. Consequently, two typical distributions, namely a Uniform Distribution and a β-distribution (Fig. parameters n and τ does not impact the stationary state of the fundamental acoustic mode.

Therefore, the partial disagreement found between the numerical and the experimental stability analysis is not related to the present Flame Transfer Function model. Extrapolating the range of uncertainty kept for n and τ (a 10% uncertainty for each) would certainly perturbed the modal growth rate but this should not be consistent with the range of uncertainty observed by experimentalists.

Conclusions and discussions

Surrogate modelling techniques have been designed in this study for Uncertainty Quantification analysis. This approach has been applied in the context of thermoacoustic analysis of a single swirled combustor experiment. All eigenmodes of the combustor have been assessed by means of a parallel Helmholtz solver. The Flame Transfer Function measured experimentally has been used as a flame model to feed the Helmholtz solver. The frequency of oscillation as well as the growth rate of the first thermoacoustic mode were computed for 24 different operating points and the stability analysis of the system has been performed by [START_REF] Silva | Combining a Helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF]. Numerical predictions are coherent with the experimental observations of the combustor, except in 3 cases (out of 24) where the agreement is only partial. Introducing Uncertainty Quantification allows a more accurate mode classification than the usual binary one (stable or unstable), and thus a more reliable comparison between experimental observations and numerical predictions. As a consequence, a continuous classification of the thermoacoustic modes is adopted based on the probability of a mode to be unstable given the uncertainties on the flame response, also called Risk Factor. At first the Risk Factor associated to the first acoustic mode of the combustor was assessed Chapter 6

Uncertainty Quantification using the Active Subspace method

Introduction

This chapter investigates the effects of uncertainties on the thermoacoustics of annular combustor with several swirlers and flames. The Active Subspace method mentioned in Section 2.3 is combined with efficient surrogate techniques to determine the statistical output of the growth rate of the acoustic pressure disturbances and thus the modal Risk Factor of the system. An overview of the UQ strategy in this work is presented in Section 6.2. The brute force Monte Carlo used to get insight on the response of the system is detailed in Section 6.3. The dimension reduction realised by mean of the Active Subspace method is discussed in Section 6.4 and Section 6.5 details the surrogate methods constructed to provide the modal Risk Factor at low cost. Finally, discussion and perspectives on the Uncertainty Quantification strategy developed in the study are discussed in Section 6.6 and the possibility to settle it on the 3D Helmholtz solver AVSP is broached.

Overview

Various computational methods have been proposed and developed during the last few decades to solve high dimensional Uncertainty Quantification problems. The majority of the theories and ACTIVE SUBSPACE METHOD the study is initiated by generating an ensemble of random perturbations of the Flame Transfer Function parameters n and τ . These are drawn using a uniform probability distribution and the bounds considered to parameterise the latter distribution are set to σn n = στ τ = 10%, where n = 6.57J/m and τ = 9.87 × 10 -4 s are the nominal values respectively for the interaction index n and the time delay τ (see Fig. 6.1). Furthermore, all injectors and flames are considered to be statistically identical and the operating conditions are similar to those reported in Table . 4.8.

A preliminary convergence diagnostics is performed (e.g. mean and standard deviation) to

Identical Flames 6.57 9.84 × 10 -4 ensure uniformly distributed statistical input parameters and thus a well-established convergence of the Monte Carlo database. This task is achieved by using an increasing refinement of the probabilistic space discretization. Results are shown in Fig. 6.1. The convergence analysis shows that performing 10,000 deterministic calculations with ATACAMAC is enough to provide the PDF of the growth rate and subsequently a sufficient converged estimate for the modal Risk

Factor of the combustor.

Monte Carlo results are presented in Fig. 6.2 and the Risk Factor computed for the first azimuthal mode of the combustor is 84%.

The brute force Monte Carlo approach can be used without difficulty when the system is represented by ATACAMAC. In cases where a more complete description like a 3D Helmholtz solver must be used (to account for example for modes which are non fully azimuthal), the Monte Carlo approach would not be feasible. Hence, the purpose is to take advantage from the analytical tool ATACAMAC to investigate an efficient UQ strategy that will be applicable prospectively to more complex solvers to approximate the response surface of the system. That is why the Active Subspace method is examined as an alternative solution to determine as a first step the subspace of inputs that most strongly affect the growth rate response, and to reduce the dimension of the input space. ACTIVE SUBSPACE METHOD where M stands for the number of the gradient evaluations. ATACAMAC provides the growth rate ω i by finite differences. In the case considered, there are 2 × N uncertain parameters (2 uncertain parameters per flame (n and τ )) characterizing the growth rate response of the full annular combustion chamber. A Finite Difference approximation of the gradients is realised For linear surrogate models, the number of basis functions increases linearly with the number of input parameters. However, for quadratic models, the number of basis functions (monomials with a degree of at most 2) evolves quadratically with the number of parameters. Besides, surrogate models are referred to an approximate model fitting sample data meaning that a sufficient number of simulations is required to approximate accurately the statistics of the model's output e.g. Eq. (6.10). Moreover, to sample the full-dimensional space, the number of points should be increased as the number of model's coefficients increases. The use of such high dimensional surrogate models becomes quickly unmanageable even when using the top-notch high-fidelity CFD solver (based on LES techniques for example) and consequently, building up a surrogate model by iteratively fitting along the active subspace is highly desirable. But, by reducing the input space dimensionality, a slight penalty in the accuracy of the surrogate model is accounted in exchange for the opportunity to tackle the high dimensional problem. Illustrating the potential of dimension reduction towards Uncertainty Quantification analysis is the main interest in this work. To this end, the following strategy is adopted:

1 The surrogate models reported in Table .6.2 are used to approximate the response surface of the system and hence to compute the modal Risk Factor. These are fitted using a least mean squares method (see chapter.

Discussions and perspectives

Dealing with complex industrial system, like a full annular combustion chamber, implies the need for the development of proper simulation tools for safety analysis and contribute to rational design policies. Several coupled physical mechanisms are involved when modelling such complex systems and thus a large number of uncertain parameters are implied. Therefore, the question of the reliability of these simulations must be addressed. Consequently, innovative Uncertainty Quantification methodologies must be used to tackle the «curse of dimensionality» which makes the technique often infeasible when increasing the size of the problem.

Uncertainty Quantification strategy has been applied to the thermoacoustic stability of a realistic full annular helicopter engine to determine its Risk Factor, defined as the probability of the first azimuthal chamber mode to be unstable. The system contains N burners and flames in a weakly coupled regime as it was discussed in chapter. Finally, the continuous adjoint Helmholtz equation is:

(7.24)

Implementation of the continuous adjoint Helmholtz equation in the AVSP solver

In this section, the implementation of the continuous adjoint Helmholtz equation in the 3D solver AVSP is investigated. Such an adjoint capability makes the calculations of the growth rate sensitivities accessible when the input parameters of a system are perturbed.

The key changes necessary to implement the continuous adjoint Helmholtz equation in the AVSP solver consist of:

1 Introducing the Gaussian formulation (see, Eq.(7.17)) to measure the pressure gradient at the reference location.

2 Constructing only the second term of Eq.(7.21) to make the adjoint problem well posed, the first term of Eq.(7.7) being self-adjoint.

To validate the implementation of the continuous adjoint equation in the AVSP solver, different geometries are used. Each of these configuration is presented in Table . 7.1 and the operating conditions used for the AVSP calculations are shown in Table . 7.2. The continuous adjoint eigenvalues are compared against the direct eigenvalues in Table . 7.6.

Boundary condition

Geometry

Additionally, information on the relative error, ||ω-ω † * || ||ω|| , between both direct and adjoint eigenvalues is shown.

Direct resolution Adjoint resolution

Error

Gradient estimations by adjoint method in the 3D Helmholtz solver AVSP

In this section, the continuous adjoint method is used to compute the gradients of the growth rate ∇f Im with respect to the flame input parameters n and τ . The accuracy of the approach is first assessed by comparison with finite difference estimates. The computational costs required to compute the gradients with both methods are then compared.

-Gradients calculations by adjoint method:

As for the direct Helmholtz equation Eq. (3.17), Eq. (7.24) is discretized using finite volume method thus leading to the following matrix formulation:

In absence of complex valued boundary condition and heat release, Eq. (7.25) reduces to: (7.26) thus leading to a linear eigenproblem in p † easy to solve in AVSP solver. When accounting for the flame effects or non trivial boundary condition, Eq. (7.25) is solved with the same fixed point iterative algorithm described in Section 3.1.4 to determine the discrete non-linear adjoint eigenpair (ω * , p † ).

To evaluate the growth rate gradients, both the direct and adjoint eigenmodes must be first provided by solving the discretized direct and adjoint Helmholtz equations (Eq. (3.28) and

Eq. (7.25)). Typically, the following iterative algorithm is used: For all the cases, a plateau appears where:

1-Passive Flame resolution: FOR THERMOACOUSTIC INSTABILITIES

⋄ ∂f Im /∂n is independent on δϵ in the range {1 × 10 -7 , 1 × 10 1 } ⋄ ∂f Im /∂τ is independent on δϵ in the range {1 × 10 -7 , 1 × 10 -4 }.

For smaller increments, the difference between f Im (τ ) and f Im (τ + δτ ) or either f Im (n) and f Im (n + δn) is very small and sensitive to numerical errors so that the derivative estimate is not robust. For larger increment, the thermoacoustic system does not behaves linearly on both the ranges [n, δn] and [τ, δτ ] and the finite difference approximation is not accurate. From Fig. 7.7, δϵ = 1 × 10 -6 can be used to provide accurate and robust estimates of ∂f Im /∂n and ∂f Im /∂τ .

Comparisons between gradients by adjoint and finite difference methods:

The computation of the gradients of the first acoustic mode for each of the geometry in Good agreements are also found between the gradients estimated by adjoint method and those computed by finite difference approach when the global flame amplitude n is perturbed.

The gradients are also computed when increasing the flame amplitude n. The results are presented in Fig. 7.4 and the ranges of variation for the flame amplitude n are reported in

Gradient by Adjoint

Gradient by Finite Difference

Simple 2D Tube 2D Mono Injector 3D Cylinder

Table 7.12: Growth rate derivatives ∇f Im computed by adjoint and Finite difference method. Both the global flame amplitude n and the time delay τ are varied: the global flame amplitude n is increased while the time delay τ varies over a period T= 1 f0 .

CHAPTER 7. ON THE APPLICATION OF THE ADJOINT METHOD FOR THERMOACOUSTIC INSTABILITIES

When both the flame parameters are perturbed, gradient calculations made by finite difference and adjoint method match. The computation cost when computing the gradient by adjoint is less demanding than the finite difference approximations because no additional solutions of the Direct equation are required. In this work the study has been focused on the flame parameters because they are know to have a non-negligible impact on the stability of thermoacoustic systems.

However, the same analysis can be conducted by varying more parameters such as the complex impedance Z, the mean sound speed c 0 (⃗ x), the mean density ρ 0 (⃗ x), the heat capacity ratio γ(⃗ x), the mean pressure P 0 or even the geometrical parameters of the systems. This shows that the adjoint method provides an efficient framework to evaluate accurately the gradients and would be suitable to account for more uncertain parameters on complex geometries. Although being an alternative to compute the gradients, adjoint would contribute to further optimize complex gas turbine combustors.

Concluding remarks and perspectives

Continuous adjoint equations have been derived and implemented in the three dimensional Helmholtz solver AVSP. This adjoint method was developed to allow for the calculation of the thermoacoustic eigenmode gradients by solving only a second set of equations, the so-called adjoint equations. Combining the results from the solution of the adjoint equation and the direct one allows to compute the gradients with respect to the input parameter of the system.

The treatment of high-dimensional and large scale thermoacoustic problems with adjoint method have not been realised in this work. Its applicability requires the uptake and further robust developments to better handle the parallel processing of the 3D adjoint solver. Therefore, the algorithm have been validated on two-and three-dimensional test cases. A complimentary finite difference method have been constructed and used as a benchmark to validate the accuracy of the gradients computed by the adjoint method. Overall, a good agreement is found.

Several conclusions can be made from the study:

⋄ To measure the pressure gradient using a Gaussian formulation, the standard deviation of the Gaussian function must be selected wisely according to the minimum mesh size of the

CHAPTER 7. ON THE APPLICATION OF THE ADJOINT METHOD FOR THERMOACOUSTIC INSTABILITIES

geometry.

⋄ The step size perturbation of the input parameter need to be carefully selected otherwise potential numerical errors would appear. This would be impacting for the gradient estimations by both finite difference and adjoint method.

In this work, it was observed that the continuous-adjoint equation requires generally less resolution and usually converges more quickly than the direct equation. Therefore, considering the gradient computations by adjoint method would be far more interesting to tackle high dimensional problems.

Appendices

Appendix A

The second annular helicopter engine

A.1 Description of the annular geometry

The second industrial configuration targeted in this study is a full annular helicopter combustion 3475 chamber designed by Safran Helicopter Engines. The combustion chamber is made up with a downstream annular combustor and an upstream annular casing that are connected to N injectors.

Instead of the first annular case investigated in Chapter 6, this system features less injectors and flames. Each burner is composed of swirler in whom fuel is injected to efficiently mix kerosene with air prior to combustion. A sketch of the helicopter combustion chamber is presented in

A.2 Thermoacoustic analysis of the second full annular combustor with N injectors and flames

A.2.1 Large Eddy Simulation of the second annular helicopter engine

The Large Eddy simulation of the annular helicopter engine has been conducted using the LES

The air flowing in the casing feeds the combustion chamber through the swirler, cooling films and dilution holes, all of those being explicitly meshed and resolved. Multi-perforated walls used to cool the liners are taken into account by a homogeneous boundary condition. Such a condition is not suited to account for acoustic damping at the combustor wall, resulting in a zero dissipation of acoustic waves at the combustor liner, thus often leading to an overestimation of the acoustic activity in the combustion chamber.

The analysis of LES results has revealed strong acoustic oscillations at a frequency close to 500 Hz. At this frequency, the pressure fluctuations grow in amplitude and lead to acoustic In order to get a better understanding of the system behaviour, a similar study to that of the N burners configuration realised in Chapter 4 is conducted:

⋄ At first, the pulsated single sector LES calculations are used to extract the input parameters c 0 (⃗ x), γ(⃗ x), ρ 0 (⃗ x) as well as the flame parameters fields n(⃗ x) and τ (⃗ x).

⋄ These inputs are then used to perform pure acoustic calculations using the AVSP solver.

Single sector and full annular computations are performed to determine both the structure and the growth rate of the thermoacoustic modes developing inside the combustor. The objective is to identify the unstable mode observed in the single sector LES computations (500 Hz) and to deal with unstable azimuthal modes that would potentially expand inside the configuration.

⋄ The 3D results obtained with AVSP are then fitted to the quasi-analytical tool ATACAMAC to get insight of the coupling phenomena and the nature of the unstable azimuthal mode developing in the system. This allows conducting computationally efficient Uncertainty Quantification analysis to determine the Risk Factor of the predominant azimuthal mode of the combustor.

A.2.2 Acoustic computations using the Helmholtz solver AVSP

As for the N-burner configuration, the AVSP calculations are performed in the steady and the In both computations, an acoustic mode at 495.5 Hz is observed and its structure is presented in Fig Uncertainty Quantification to determine the probability of this mode to stay unstable. For that purpose, it is however important to perform the thermoacoustic computations within a reasonable computational timeframe because the computational cost with the AVSP solver is about 60000 CPU hours per simulation on 120 processors. Therefore, the same procedure used to investigate the stability of the N injectors configuration (see Chapter 4) is reused. Typically, the 3D AVSP results are used to fit the analytical network modelling tool ATACAMAC to focus on the coupling between the system cavities and to perform Uncertainty Quantification analysis at low cost.

A.2.5 Acoustic computations using the quasi-analytical tool ATACAMAC

The functional operating conditions used for numerical applications of the system with N injectors were provided by Safran Helicopter Engines and extracted from the forward AVSP calculations

A.3 Uncertainty Quantification study

In this section, the Uncertainty Quantification of the full annular combustor with N injectors and flames is investigated. To achieve this task, the analytical tool ATACAMAC is used to determine the probability of the first azimuthal mode, reported in Table . A.7, to be unstable (namely its 3600 Risk Factor its determined). The Uncertainty Quantification strategy conducted in this work is similar to the one employed for the first annular system studied in Chapter 6.

The approach is the following:

1 Taking advantage of the computational efficiency of the ATACAMAC tool, a Brute Force Monte Carlo is first performed on the 2×N uncertain input parameters (2 uncertain parameter n and τ per flame) to get insight of the growth rate surface response of the combustor.

Uncertainties are propagated through the system by using a uniform distribution to generate random perturbation of the flame input parameters n and τ . The same uncertainty range as those of the UQ analysis of the N burners configuration are kept: 2 The objective is now to use the Active Subspace method described in Chapter 6 to perform the UQ analysis of the combustor using only the relevant low-dimensional subspaces of the full annular geometry. The Active Subspace method [START_REF] Constantine | Active subspace methods in theory and pratice: applications to kriging surfaces[END_REF]) is used to reduce the dimension of the full parameter space to just a few. To find active variables of the system, the method requires gradient evaluation to detect which directions in the parameter space lead to strong variations of the growth rate. Other directions leading to a flat response surface are not useful for describing the combustor stability and hence can be disregarded. As discussed in Chapter 6 and [START_REF] Bauerheim | On a novel uncertainty quantification approach combining active subspace and adjoint: application to symmetry breaking of azimuthal modes[END_REF], using the physical quantities associated with the Fourier transform of the FTF is better to ease the physical interpretation of active variables as well as improve the accuracy of the gradient calculations. Based on these asserts, the eigenvalues spectrum is drawn in It suggests the existence of a 3D active subspace, the first one being a constant corresponding to an equi-weighted linear combination, i.e., associated with the mean Flame Transfer Function over the N burners.

3 The idea is now to use the 3D active variables to fit algebraic surrogate models to determine the Risk Factor of the first azimuthal mode of the combustor with much less than 10000 ATACAMAC simulations. Linear and quadratic models are used to investigate the response surface of the system (see Table .A.8):

⋄ Linear models:

-L F D : The first linear model is constructed in the full dimension and initial probabilistic space.

-L 3D : The second linear model is spanned along the reduced subspace with the 3 active variables discovered with Active Subspace method.

⋄ Quadratic models:

-Q F D : The first quadratic model is constructed in the full dimension and initial probabilistic space.

-Q 3D : The second quadratic model is built on the reduced subspace with the 3 active variables discovered with Active Subspace method. In the second industrial helicopter engine investigated, the use of M=50 and M=1000 AT-ACAMAC calculations was enough to provide accurate estimate of the modal Risk Factor.

Model

Based on this, these sample sizes are used to fit the algebraic surrogate models of the Nburner configuration. To provide the correlation between the low-order models and the true response surface of the system, the least means squares method is used to determine the Pearson's correlation coefficients with Eq. (5.15). The results are gathered in Table . A.9. The results prove that only 50 ATACAMAC simulations are enough to provide a correct approximation of the response surface of the combustor. However, the quadratic models Q F D and Q 3D are better correlated with the true response surface. As discussed in the above study of the N burner geometry in Section 6.5.1 a better accuracy of the system To ensure an accurate Risk Factor estimation, 25000 replays of the quadratic surrogate models Q 50 3D and Q 1000 3D are realised. The resulting Risk Factor approximated is then compared to the one determined from the Brute force Monte Carlo study in Table . A.10. and % with 54.01 when using Q 50 3D . Particularly, when using only the 3D active variables of the combustor, it is shown that the quadratic surrogate model Q 50 3D provides reasonable approximation of the Risk Factor within an error below 5%. The use of Active Subspace method proves again satisfactory in reducing the system dimensionality and by providing 3670 accurate estimate of the modal Risk Factor of the combustion chamber.