Emmanuel Motheau

Thierry Poinsot

Pr Simone

Pr Aimee Morgans

Dr Michel Cazalens -Snecma

Accounting for mean flow effects in a zero-Mach number thermo-acoustic solver: Application to entropy induced combustion instabilities

Keywords: Combustion instabilities, Large-Eddy Simulation, indirect combustion noise, mixed-modes, low-order thermoacoustic models

Virtually all combustion chambers are subject to instabilities. Consequently there is a need to better understand them so as to control them. A possibility is to simulate the reactive flow within a combustor with the Large-Eddy Simulation (LES) method. However LES results come at a tremendous computational cost. Another route is to reduce the complexity of the problem to a simple thermoacoustic Helmholtz wave equation, which can be solved in the frequency domain as an eigenvalue problem. The coupling between the flame and the acoustics is then taken into account via proper models. The main drawback of this latter methodology is that it relies on the zero-Mach number assumption. Hence all phenomena inherent to mean flow effects are neglected. The present thesis aims to provide a novel strategy to introduce back some mean flow effects within the zero-Mach number framework. In a first part, the proper way to impose high-speed elements such as a turbine is investigated. The second part focuses on the coupling between acoustics and temperature heterogeneities that are naturally generated at the flame and convected downstream by the flow. Such phenomenon is important because it is responsible for indirect combustion noise that may drive a thermoacoustic instability. A Delayed Entropy Coupled Boundary Condition (DECBC) is then derived in order to account for this latter mechanism in the framework of a Helmholtz solver where the baseline flow is assumed at rest. In the last part, a realistic aero-engine combustor that features a mixed acoustic/entropy instability is studied. The methodology developed in the present thesis is tested and compared to LES computations. It is shown that computations with the Helmholtz solver can reproduce a complex combustion instability, and that this latter methodology is a potential tool to design new combustors so as to predict and avoid combustion instabilities.

Résumé en Français

Pratiquement toutes les chambres de combustion présentent des instabilités. Par conséquent, il est nécessaire de mieux les comprendre afin de les contrôler. Une possibilité est de simuler l'écoulement réactif à l'intérieur d'une chambre de combustion grâce à la Simulation aux Grandes Echelles (SGE). Cependant la SGE est très coûteuse en terme de capacité de calcul. Une autre possibilité est de réduire la complexité du problème à une simple équation d'onde thermoacoustique (équation dite de Helmholtz), qui peut être résolue fréquentiellement comme un problème aux valeurs propres. Le couplage entre l'acoustique et la flamme est alors prise en compte au travers des modèles appropriés. Le principal problème de cette méthode est qu'elle repose sur l'hypothèse d'un nombre de Mach nul. Tous les phenomènes lieś à l'écoulement moyen sont donc négligés. La présente thèse propose une nouvelle stratégie pour prendre en compte certains effets de l'écoulement dans un contexte à Mach nul. Dans une première partie, la manière la plus judicieuse d'imposer un élement présentant un écoulement très rapide est étudiée. La seconde partie se focalise sur le couplage entre l'acoustique et les hétérogenéites de temperature qui sont générées par la flamme et naturellement convectées par l'écoulement moyen. Ce phénomène est important car il est responsable du bruit indirect de combustion qui peut conduire à une instabilité thermoacoustique. Un nouveau type de condition limite (DECBC) est proposé afin de prendre en compte ce méchanisme dans un contexte de résolution de l'équation de Helmholtz à Mach nul. Dans la dernière partie, une chambre de combustion aéronautique présentant une instabilité mixte acoustique/entropique est etudiée. Le bénéfice des méthodes développées dans la présente thèse est testé et comparé à des calculs avec la SGE. Il est montré que les calculs avec un solver de Helmholtz peut reproduire une instabilité de combustion complexe, et que cet outil s'avère avoir le potentiel prédire les instabilités afin de concevoir de nouvelles chambres de combustion.

Mots clefs : Instabilités de combustion, Simulation Grandes Echelles, bruit indirect de combustion, modes mixtes, modélisation thermoacoustique basordre.

Jamais aucun effort administratif ou scolaire ne remplacera les miracles du hasard auquel on doit les grands hommes.

Honoré de Balzac

Etrange parcours que celui conduisant à la réalisation de cette thèse. Il y a environ une dizaine d'années, en 2005, j'escaladais les voûtes de l'église Saint-Louis-en-l' Île afin d'y glisser des micros pour réaliser la prise de son du tout nouvel orgue Aubertin avec Zig-Zag Territoires. C'est durant cette période qu'Alban Moraud me fit découvrir -et je l'en remercie -ses cours d'acoustique qu'il suivait au CNAM, sous le regard dubitatif de Franck Jaffrès qui préférait s'en remettre à ses oreilles (et il a bien raison !) plutôt qu'aux fonctions de Green. Quelques mois plus tard, je réalisais mon premier stage en laboratoire universitaire, dans le domaine de l'acoustique musicale. Joies de la synthèse sonore par modèles physiques ! L'année suivante, changement de décor, j'étais à DCNS pour faire du rayonnement acoustique de sous-marins, et l'année suivante c'était en toute logique que je découvrais la combustion et la thermo-acoustique avec Thomas Lederlin à Turbomeca, me permettant ainsi de clore mon cursus d'ingénieur en acoustique. J'en profite pour remercier Pr. Alexandre Garcia, titulaire de la Chaire d'Acoustique au CNAM et Manuel Melon, maître de conférence, pour la transmission de leur savoir de l'acoustique sous toutes ses formes et leur aide précieuse, surtout pour les tracasseries administratives ! C'est ainsi que de fil en aiguille on en arrive à cette thèse CIFRE sous financement Snecma. Je remercie tout d'abord mes deux directeurs de thèse Thierry Poinsot et Franck Nicoud pour m'avoir confié le sujet (le code AVSP, quoi), et surtout pour leur aide et conseils tout au long de ces trois années. Je ne pensais pas au tout début que cette simple idée de modéliser la convection d'entropie dans un contexte à Mach nul nous conduirait à autant de développements. Merci aussi de m'avoir emmené à Stanford pour le Summer Program, la piscine quotidienne c'était super. Plaisanteries mises à part, ce fut une expérience très forte et enrichissante. A ce sujet, je remercie Florent Duchaine pour nous avoir trouvé cette magnifique demeure à deux pas du campus ! A ma femme Ling-Yin, à mon fils Léandre.

x Nomenclature

General introduction

Les hommes se souviennent mal qu'on leur a défendu, étant enfants, d'ouvrir le ventre des pantins... (c'est déjà un crime de lèse-mystère): ils continuent à vouloir fourrer leur esthétique nez là où il n'a que faire. S'ils ne crèvent plus de pantins, ils expliquent, démontent et, froidement, tuent le mystère: c'est plus commode et alors on peut causer.

Monsieur Croche (Claude Debussy)

As recalled by [START_REF] Candel | Advances in combustion and propulsion applications[END_REF] in a review of advances in combustion for propulsion applications, despite the fact that man mastered fire half a million years ago, the science of combustion is much more recent. In his own words concerning the evolution of aviation: One may wonder how we got from there to here in just 100 years, from the 30 m Wright brothers flight to the transportation of an astounding 4000 billion revenue passenger-km (RPK) per year. However, the theories are far from delivering all the secrets of combustion and engineers did not wait for Science to develop machines for the incredible achievements cited above.

Combustion is used first because the chemical reactions that take place in a combustor generate a very large energy. As recalled by [START_REF] Sawyer | Science based policy for adressing energy and environmental problems[END_REF] combustion provides more than 90 % of the energy in the world, not only for transportation but also for industrial, agrarian and, of course, domestic applications. Although alternative sources are emerging, such as wind or solar energies, their efficiency is limited and the production of energy will remain combustion in the next decades. However this source of energy is facing several problems:

• burning any kind of fuels generates air pollution that directly affects health,

• although this is a matter subject to controversy, the climate is changing and combustion does not seem unconcerned to this phenomenon,

• there is no debate that fossil resources are not infinite and that they will be depleted one day.

While new forms of fuels or energy sources are currently under study for the future, existing systems must be improved so as to face the challenges cited above. Concerning the aeronautical domain, since 2001 the Advisory Council for Aeronautics Research in Europe (ACARE) has set several objectives to 2050 through the Flightpath 2050 project. To cite a few:

• reduce fuel consumption and CO2 emissions by 75 % per passenger kilometer,

• reduce nitrogen oxides (NOx) emissions by 90 %,

• reduce perceived noise by 65 %.

Most of the aeronautical propulsion relies on gas turbine engines, because they provide a large power compared to their size and weight and can be embedded on planes. Moreover they are reliable machines and operate on a wide variety of fuels. Note that gas turbines are not only devoted to aircraft engines but also cover a broad range of applications including electricity generation or water and sewage pumping stations. The power output also varies depending on the machine. For example the SIEMENS largest ground-based gas turbine for electricity generation can provide up to 350 MW and has a weight of approx. 440 t. For an order of magnitude, in the aeronautical context the A400M turboprop TP400 has a maximum power of approx. 8 MW for a weight of less than 2000 kg.

The main idea of a gas turbine is to convert the motion of a gas into a mechanical motion. For the story, the premises of this kind of machine go back to the Ancient Greek civilisation around 150 BC, with the Hero's aeolipile (see Fig. 1, left). The application of this invention is not really known, maybe it was for education purpose. A much more practical use of the gas turbine principle was provided by Leonardo Da Vinci with his Chimney Jack (see Fig. 1, right). The hot air rising from a fire under a roasting pork passes through an axial turbine in an exhaust duct, which convert the energy into a motion that turns the roasting spit. It is a perpetual motion that is self-sustaining as long as there is combustion. Apart from these anecdotes, the modern industrial gas turbine goes back to the beginning of the 20th century. It is mainly composed of three major components: a compressor, a combustion chamber and a turbine. The air passes through the compressor so as to get a high pressure and enters in the combustor where it is mixed with a fuel. The burned gas are ejected downstream through an exhaust. To increase the efficiency, a turbine recovers the energy so as to transfer it through the compressor which is on the same axis.

In order to achieve the objectives of reduction of pollution and consumption, manufacturers are facing a great challenge because gas turbines have achieved a high level of complexity and efficiency optimisation, and now the smallest change can destabilise the overall system. One of the new concepts for combustion technology recently developed is lean combustion (see [START_REF] Huang | Dynamics and stability of lean-premixed swirlstabilized combustion[END_REF][START_REF] Lieuwen | Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling[END_REF]. Fuel and air are premixed upstream of the combustor to avoid the flame to burn at the stoichiometry. The flame takes place in a zone where the air is in excess, which has the effect to reduce the flame temperature and consequently also the thermal NOx formation. Due to a higher combustion efficiency a reduction of formed CO is also achieved. Unfortunately, these lean-burn conditions often result in the occurrence of instabilities, since combustion tends to occur near the lean blow-out limit [START_REF] Venkataraman | Mechanism of combustion instability in a lean premixed dump combustor[END_REF] which is more susceptible to unstable combustion dynamics.

Combustion Instabilities

Instabilities are resulting from a driving process which generates a perturbation in the reacting flow that couples back with a resonant interaction. As simply explained by [START_REF] Huang | Dynamics and stability of lean-premixed swirlstabilized combustion[END_REF], two fundamental causes are responsible for the occurrence of instabilities in combustion systems:

• combustion chambers are almost entirely closed and the damping internal processes are too weak to attenuate unsteady motions;

• only a small fraction of the energy released by chemical reactions (less than 0.001%) is sufficient to drive unsteady motions having peak amplitudes of the order of the mean chamber pressure.

Basically this mechanism is cyclic and involves time lags. Depending on the phase shift between a driving process and the time response of the flow, a constructive coupling can lead the overall system to accumulate energy that may lead to strong, sustained fluctuations. Virtually all kinds of combustion engines such as gas-turbines, aeronautical turbines or rocket engines experience instabilities at some stage of their design or life-cycle.

Combustion instabilities have many undesirable effects such as large-amplitude structural vibrations, flame flashback or blowoff, or an abnormally high temperature of the wall of the combustor. The consequences of combustion instabilities range from loss of performance or premature fatigue of materials to a possible mechanical failure or destruction of the system. Figure 2 presents pictures of a combustor before and after being damaged by an instability. Combustion instabilities are a very broad field and an important review of the driving mechanisms may be found in the literature [START_REF] Kotake | On combustion noise related to chemical reactions[END_REF][START_REF] Candel | Combustion dynamics and control: progress and challenges[END_REF] Following [START_REF] Barrère | Comparison of combustion instabilities found in various types of combustion chambers[END_REF], combustion instabilities may be distinguished between three classes:

• Intrinsic instabilities which are inherent to motions induced by the flame itself on its own field, like Darrieus-Landau or thermo-diffusive instabilities. They are observed in the absence of any external influences.

• Chamber instabilities which are essentially related to the acoustic modes of the chamber that interact with the flame. Flame/acoustic coupling is a process that has been observed for a long time, either visually by noting that the flames in a fish-tail gas burner enjoy the music of Beethoven (see [START_REF] Leconte | On the influence of musical sounds on the flame of a jet of coal-gas[END_REF] or either by hearing the sound produced by a flame in a tube at a particular location (see [START_REF] Rijke | Notice of a new method of causing a vibration of the air contained in a tube open at both ends[END_REF]. This latter device was first studied by [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF]. An extensive review of the socalled Rijke tube can be found in the literature (see [START_REF] Raun | A review of rijke tubes, rijke burners and related devices[END_REF]. Basically the idea is that the flame is a source of acoustic-pressure fluctuations [START_REF] Kotake | On combustion noise related to chemical reactions[END_REF], and when placed in a location where fluctuations of pressure may occur such as near a pressure anti-node, the acoustic eigenmodes of the system can be excited and oscillate as long as the flame exists. In order to determine if a system is stable or not, Rayleigh (1878) has developed a well-known Rayleigh criteria. Formally it is often written as

V p ′ q ′ dV > 0 ( 1
)
where V is the volume of the combustor while p ′ and q ′ are the fluctuations of pressure and heat release rate, respectively. Basically this criterion states that the combustor is unstable when the pressure and heat release rate fluctuations over the combustor volume are in phase. This criterion is somewhat simplistic and the stability of a thermoacoustic system should be rather studied with a more rigorous theoretical framework based on acoustic or disturbance energy balance (see §1.2 for more details). Moreoever, one should notice that not only flame/acoustic coupling can exist in a chamber. Vortex roll-up due to shear-flow instability near the wall of the chamber may also be a source of flame perturbation [START_REF] Renard | Dynamics of flame/vortex interactions[END_REF]. This latter mechanism is most often encountered in rockets and ramjets, and most common combustion instabilities encountered in gas turbines rely on a thermo-acoustic path.

• System instabilities involve various components of the configuration. In practical systems, the combustion chamber is not alone but coupled to a casing by complex injection systems like a swirler, or by multiperforated cooling plates. A compressor and a turbine are also present on both sides of the chamber. All of these components have a potential interaction with the flame and the overall stability may be affected. For example, temperature fluctuations generated by the unsteady flame can be convected by the mean flow to the downstream nozzle. It is well known [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] that temperature fluctuations impacting a nozzle generate pressure fluctuations. This driving process may couple back with the flame. This latter mechanism is one of the main subjects of the present thesis and will be studied in detail in the following chapters.

It is worth noting that manufacturers want to avoid combustion instabilities and need to understand them so as to control them. Experiments on practical industrial configurations are possible, but they are costly and do not permit to access to the fundamental processes of the combustion. Moreover, the experiments allow especially to observe from a macroscopic point of view the existence of problems rather than to predict them. The prediction of combustion instabilities at an early design stage is a challenging task, especially because all mechanisms as well as their interactions in a practical context are not fully understood.

The control of combustion instabilities a posteriori is possible and rely on the addition of passive acoustic dampers such as liners or Helmholtz resonator so as to enhance the rate of absorption of acoustic energy (see [START_REF] Dupère | The use of helmholtz resonators in a practical combustor[END_REF][START_REF] Noiray | Passive control of combustion instabilities involving premixed flames anchored on perforated plates[END_REF][START_REF] Zhao | Tuned passive control of combustion instabilities using multiple helmholtz resonators[END_REF]. It is also possible to act on the flame dynamics by modification of the manifold that injects reactants inside the combustor. The combustion region can becomes less susceptible to perturbations. Active control techniques proved to be efficient on academic combustors (see [START_REF] Bloxsidge | Active control of reheat buzz[END_REF][START_REF] Mcmanus | A review of active control of combustion instabilities[END_REF][START_REF] Annaswamy | Active control of combustion instability: theory and practice[END_REF][START_REF] Candel | Combustion dynamics and control: progress and challenges[END_REF][START_REF] Dowling | Feedback control of combustion oscillations[END_REF][START_REF] Morgans | Model-based control of combustion instabilities[END_REF]. However the use of such methods are difficult because of certification constraints. Indeed, active control methods are based on measurements techniques that detect the instability so as to act on different parts of the system, such as the modulation of the fuel mass-flow rate or the modification of acoustic boundary conditions by a loudspeaker. The addition of a such a closed feedback loop in the overall system is critical because the risk of failure due to model-based approximations may enhance an instability [START_REF] Annaswamy | Impact of linear coupling on the design of active controllers for the thermoacoustic instability[END_REF].

Consequently, there is a need for a better understanding of combustion instabilities at the fundamental level, but also the development of tools that can predict them at an early design stage. The next section is devoted to a review of the different methods that can help to understand and predict combustion instabilities.

Prediction of combustion instabilities

Before describing the different methods that allow to predict combustion instabilities, one should introduce first the typical behaviour of an instability in a combustor. The time evolution of the pressure signal recorded in a combustor when a combustion instability is triggered at t = 0 is depicted in Fig. 3. The overall evolution of the instability can be decomposed into three zones:

• a linear zone where the oscillations appear and follow an exponential growth;

• the amplitudes do not grow indefinitely and reach a saturation characterised by a limit-cycle, the non-linear effects are here essential;

• the transition between the linear and non-linear zone quite often exhibit an overshoot period.

Note that in certain thermoacoustic systems, hysteresis and bifurcation effects may be encountered [START_REF] Mariappan | Modelling nonlinear thermoacoustic instability in an electrically heated rijke tube[END_REF]. Depending on the amplitude and the frequency of a perturbation, the response of the system may not follow a linear evolution (Balachandran et al., 2005; Palies et al., 2011). Moreover, a linearly stable combustor, i.e. that does not pulse spontaneously, may be triggered by ambient noise or high amplitude perturbations and reach an unstable behaviour [START_REF] Balasubramanian | Thermoacoustic instability in a rijke tube: Non-normality and nonlinearity[END_REF][START_REF] Juniper | Triggering in the horizontal rijke tube: non-normality, transient growth and bypass transition[END_REF].

In principle, a direct numerical resolution of the unsteady reacting Navier-Stokes equations allows one to predict combustion instabilities. Due to the augmentation of computational resources and advances in higher order numerical methods so as to get extremely precise schemes, the Direct Numerical Simulation (DNS) of the equations is possible in three dimensions but still limited to laboratory scale flames and to low Reynolds flows. DNS has been used extensively to investigate fundamental turbulence-chemistry interactions in combustion topics so as to understand and validates theoretical models [START_REF] Trouvé | The evolution equation for the flame surface density[END_REF][START_REF] Poinsot | Application of direct numerical simulation to premixed turbulent combustion[END_REF][START_REF] Vervisch | Direct numerical simulation of non premixed turbulent flames[END_REF][START_REF] Echekki | Direct numerical simulation of autoignition in nonhomogeneous hydrogen-air mixtures[END_REF][START_REF] Hawkes | Direct numerical simulation of hydrogen-enriched lean premixed methane-air flames[END_REF]). At the opposite end of the spectrum, the Reynolds-Averaged-Navier-Stokes (RANS) methods only compute averaged values without solving any turbulent quantities. The central difficulty of the RANS method is that it mostly relies on models for turbulence and flow-flame interactions (Poinsot & Veynante, 2011). Although the accuracy of such an approach when dealing with reacting turbulent flows is highly questionable, this method has found great interest in engineering due to the coarse mesh involved and small restitution times for results. Halfway between DNS and RANS methods, Large Eddy Simulation (LES) constitutes an interesting alternative for numerical computations of reacting flows. The mean idea of this latter method is to simulate the large scale eddies in the flow which contain most of the energy, while the small scales are represented by subgrid models. Hence, the filtering operation is directly related to the grid mesh resolution. One difficulty raised by the LES is that the spatial filtering of the equations leads to unresolved terms such as the Reynolds stresses or the heat release and reaction rates. These terms need to be modelled so as to resolve the flame front on a LES grid. However, as thermo-acoustic instability involves large scale structures which are calculated explicitly, the LES provides a relevant description of turbulence/combustion/acoustic interactions that occur in the combustion chamber. This means that the LES covers all the zones described in Fig. 3 and is virtually able to reproduce all non-linear effects. Several studies using Large Eddy Simulations have shown the potential of this method, even for realistic engine geometries and operating points [START_REF] Pitsch | Large eddy simulation of turbulent combustion[END_REF] However LES results still come at a tremendous computational cost, impractical for the study of design or operating conditions variations and faster tools are required. A natural approach is to reduce the complexity of the problem by using a linearised formulation of the equations. Basically the main idea is to decompose the flow into two components: a mean part and a fluctuating part. As the conservation equation for mass, momentum and energy are satisfied for the mean flow field, one may deduce a set of equations for the perturbations containing only first order terms. This approach focuses on solving the perturbations equations and the mean flow field is not a part of the solution and must be determined beforehand. Moreover, the interaction of the flame with the acoustic field needs to be described by an appropriate model. Basically it relies on the observation of the response of the flame to a perturbation. One may then deduce a so-called flame transfer function (FTF) that relates the unsteady heat release to acoustic quantities at reference locations (Schuller et al., 2003b;[START_REF] Truffin | Comparison and extension of methods for acoustic identification of burners[END_REF]. The flame is then often modelled as a purely acoustic element thanks to a nτ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF], n being a gain while τ is a time delay (which may also be replaced by a phase). Note that this model relies on a linear description of the flame, the parameters n and τ being frequency dependent. More details about the nτ active flame model are available in Appendix B. The FTF can be obtained either experimentally [START_REF] Durox | Combustion dynamics of inverted conical flames[END_REF][START_REF] Kornilov | Experimental assessment of the acoustic response of laminar premixed bunsen flames[END_REF], numerically (Polifke et al., 2001b;[START_REF] Giauque | Validation of a flame transfer function reconstruction method for complex turbulent configurations[END_REF] or by semi-analytical models such as the resolution of the G-equation (Schuller et al., 2003a). Note that the nτ model is rather a simple linear representation of the flame and has been extended to a weakly non-linear characterisation with the concept of flame describing functions (FDF) by [START_REF] Noiray | A unified framework for nonlinear combustion instability analysis based on the flame describing function[END_REF]. The main idea of FDF is that parameters n and τ depend both on the amplitude and the frequency of the perturbations. The FDF framework has been extensively used to investigate the non-linear dynamic of flames, for example by [START_REF] Durox | Experimental analysis of nonlinear flame transfer functions for different flame geometries[END_REF] for several different flame geometries and by [START_REF] Kim | Effects of nonuniform reactant stoichiometry on thermoacoustic instability in a lean-premixed gas turbine combustor[END_REF] who focus on the effects of nonuniform reactant stratification.

As viscous effects can be neglected in the description of the acoustic problem, it is possible to solve either in the frequency or temporal domain a set of 3D Linearised Euler Equations (LEE) derived from the Navier-Stokes equations. Five equations are needed so as to represent the reactive flow where the local unsteady heat release appears as a forcing term, which is responsible for combustion noise and thermoacoustic instabilities. The LEE approach gives an accurate representation of the physical thermoacoustic phenomena [START_REF] Nicoud | About the zero mach number assumption in the calculation of thermoacoustic instabilitie[END_REF][START_REF] Yu | Effects of mean flow, entropy waves, and boundary conditions on longitudinal combustion instability[END_REF] at a reasonable computational cost. However, solving LEE requires numerical schemes that are subject to instability and specific numerical techniques must be used, such as the addition of artificial viscosity or filtering operations [START_REF] Agarwal | Calculation of sound propagation in nonuniform flows: suppression of instability waves[END_REF][START_REF] Rao | Use of finite element methods in frequency domain aeroacoustics[END_REF].

A natural approach is to characterise the solutions of the problem by its stable/unstable modes in the frequency domain. Low-order methods focus on the acoustic problem by modelling the configuration as a network of unidimensional ducts where the wave propagation is solved analytically [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF] or by signal processing techniques (Polifke et al., 2001b). Changes in the geometry or mean flow properties, such as temperature, density or sound velocity, are taken into account by jump conditions at the interface of each segment. The acoustic quantities are related to the amplitudes of the forward and backward waves, also known as Riemann invariants, which are determined to satisfy all jump relations and boundary conditions. The set of equations can be solved as an eigenvalue problem, giving a discrete set of complex frequencies which are the roots of a dispersion relation in the complex plane. The real part gives information about the frequency of the oscillation, while the imaginary part represents the amplitude of the damping or amplification of the instability. Despite the fact that the implementation and the use is rather straightforward even for 3D configurations [START_REF] Parmentier | A simple analytical model to study and control azimuthal instabilities in annular combustion chambers[END_REF], making it an adaptive method for the industry, the strong approximations that must be made to the geometry make this one-dimensional network approach hardly predictive in complex industrial combustors.

Another way to simplify the LEE approach while keeping an accurate description of the 3D geometry is to assume that the mean flow is at rest, often referred to as the zero-Mach number approximation. Under this assumption, an approximate linear wave equation for the amplitude p(x) of the pressure perturbations p ′ (x, t) = p(x) exp(-jωt) in reacting flows may be derived from the Navier-Stokes equations (see Poinsot & Veynante, 2011) and reads:

ρ (x) c 2 (x) ∇ ⋅ 1 ρ(x) ∇p(x) + ω 2 p(x) = jω (γ(x) -1) q(x) (2)
where q(x) is the amplitude of the unsteady heat release q ′ (x, t) = q(x) exp(-jωt), ρ is the density, c the speed of sound and ω is the complex angular frequency. In order to close the problem, the flame is modelled as a purely acoustic element thanks to a nτ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF]; Eq. ( 2) is referred as the inhomogeneous Helmholtz equation and then corresponds to a non-linear eigenvalue problem which can be solved by using appropriate discretisation techniques and numerical algorithms [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]. The main advantage of this method is its small computational burden while keeping a realistic geometry. The main drawback is that only acoustic perturbations are accounted for, since vorticity and entropy waves do not propagate when the mean flow is at rest. The effects of convection on the acoustic waves are also neglected. Nevertheless, this method offers great flexibility for physical modelling and systematic variation of input parameters. For example, the flame transfer function that represents flame/acoustic interaction or the acoustic absorption by multi-perforated plates have been successfully implemented into a three-dimensional Helmholtz solver [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF][START_REF] Gullaud | Effect of multiperforated plates on the acoustic modes in combustors[END_REF].

Combined with LES, this approach proved useful to gain understanding on the structure and nature of the instabilities observed in industrial burners [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]Wolf et al., 2012b). One should notice that frequency domain methods are devoted to assess the growth/decay rate of an unstable mode, meaning that they are restricted to the linear zone of the thermoacoustic system (see Fig. 3). However, a recent study conducted by [START_REF] Silva | Combining a helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF] has shown that the non-linear behaviour of an academical combustor and the limit-cycle amplitudes of the instability can be reproduced with the help of the FDF framework [START_REF] Noiray | A unified framework for nonlinear combustion instability analysis based on the flame describing function[END_REF].

About the zero-Mach number assumption

When seeking the thermoacoustic fluctuations in a combustor, the most natural strategy consists of computing the whole engine (cf. Fig. 4.a). In this case, one should use the LES/LEE approach because the Mach number, M , is usually not small in the compressor and the turbine. An alternative solution is to use the Helmholtz approach Eq. (2) in the combustion chamber, where M << 1 while accounting for the upstream and downstream elements through appropriate boundary conditions (cf. Fig. 4.b), namely a complexvalued impedance noted Z. These impedances can be deduced from transfer functions describing the response of acoustic elements to acoustic or entropic perturbations, either under the compact hypothesis [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] or by solving the LEE numerically [START_REF] Lamarque | Boundary conditions for acoustic eigenmodes computation in gas turbine combustion chambers[END_REF] or analytically (Duran & Moreau, 2013). At least two difficulties arise when using a Helmholtz solver (thus assuming M = 0) for the combustion chamber and complex impedances as boundary conditions:

• The impedance must be consistent with the zero-Mach number assumption although it represents an acoustic element where the mean flow is not at rest. Moreover, this impedance must be imposed at a location in the combustion chamber where the Mach number is small.

• As shown by many authors (see [START_REF] Morfey | Amplification of aerodynamic noise by convected flow inhomogeneities[END_REF][START_REF] Ffowcs-Williams | The generation of sound by density inhomogeneities in low mach number nozzle flows[END_REF][START_REF] Howe | Contributions to the theory of aerodynamic sound, with application to excess jet noise and the theory of the flute[END_REF][START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF], entropy fluctuations accelerated in a region where the mean flow is not homogeneous, which is the case when the combustor opens onto a high pressure distributor, generate acoustic waves. Acoustic waves transmitted through the distributor produce indirect combustion noise, while acoustic waves travelling back to the flame may trigger low-frequency resonant modes called rumble (Polifke et al., 2001a;[START_REF] Bloxsidge | Active control of reheat buzz[END_REF]. This important mechanism for combustion instabilities is not accounted for in a Helmholtz solver, since entropy fluctuations do not propagate if the mean flow is at rest.

It has been demonstrated by [START_REF] Nicoud | About the zero mach number assumption in the calculation of thermoacoustic instabilitie[END_REF] and [START_REF] Wieczorek | Numerical study of mach number effects on combustion instability[END_REF] that the zero-Mach number assumption can lead to significant errors for the prediction of both the frequency and the growth rate of combustion instabilities, suggesting that improvements are required for the two items cited above. Since accounting for the non zero-Mach number terms in the equations leads to a drastic increase in the complexity of the problem [START_REF] Rao | Use of finite element methods in frequency domain aeroacoustics[END_REF], the aim of this thesis is to propose models that account for these effects while keeping the simplicity and flexibility of a Helmholtz solver.

About indirect combustion noise and entropy/acoustic mixed modes

A classical mechanism for combustion instability is a constructive coupling between acoustic waves and the unsteady combustion that arises when pressure and heat release fluctuations are in phase (see [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF][START_REF] Lieuwen | Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms, and Modeling[END_REF][START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF]. Another mechanism that may also support self-sustained instabilities relies on the acoustic perturbations induced by entropy spots (temperature and/or mixture heterogeneities) being generated in the flame region and evacuated through the downstream nozzle. This latter mechanism is called indirect combustion noise when its contribution to the whole noise outside the combustion chamber is studied [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] to ease the technological issues. This means that the combustor chamber is not connected to a choked nozzle, eliminating any contribution of the indirect noise effect ab initio.

Knowledge about mixed acoustic-entropy combustion instabilities goes back to the history of rocket engines and ramjet in the 1950s. Earlier works focused on the downstream nozzle that ends the flame tube. [START_REF] Tsien | The transfer functions of rocket nozzles[END_REF] formulated a transfer function at the entrance of the nozzle via the derivation of an impedance, modelling the damping or increasing contribution of the device. The problem was reformulated later by [START_REF] Zinn | Longitudinal mode acoustic losses in short nozzles[END_REF] and [START_REF] Morfey | Amplification of aerodynamic noise by convected flow inhomogeneities[END_REF] in terms of acoustic energy balance. Multiple works dedicated to the characterisation of the nozzle impedance emerged in the 70's. [START_REF] Bohn | Response of a subsonic nozzle to acoustic and entropy disturbances[END_REF] for the subsonic case, as well as [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] for both subsonic and supersonic cases, derived the reflection/transmission coefficients for impinging acoustic and entropic waves under the zero frequency limit (compact nozzle assumption). This approach was later extended to the non-compact cases by [START_REF] Moase | The forced response of choked nozzles and supersonic diffusers[END_REF] and [START_REF] Goh | Phase prediction of the response of choked nozzles to entropy and acoustic disturbances[END_REF], and generalised for both subsonic and choked flows with and without shock wave by Duran & Moreau (2013).

It must be noted that the influence of the nozzle on the combustion zone located upstream is most of the times not accounted for. For example, experiments [START_REF] Bell | Experimental and theoretical determination of the admittances of a family of nozzles subjected to axial instabilities[END_REF][START_REF] Zinn | Experimental determination of three-dimensional liquid rocket nozzle admittances[END_REF] or numerical studies [START_REF] Bloy | The pressure waves produced by the convection of temperature disturbances in high subsonic nozzle flows[END_REF] focused on the acoustic response of nozzles for several flow parameters (influence of the Mach number or amplitudes of temperature spot) and its sensitivity to the geometry details (sharpness, area contraction effects). Purely mathematical models (Ffowcs-Williams & Howe, 1975) addressed the problem from the noise radiation point of view via the derivation of Green functions. Concerning the interaction mechanism between the combustion and the indirect noise generated at nozzles, to the authors knowledge the first attempt to develop an analytical theory goes back to the works of [START_REF] Abouseif | Ramjet rumble: The lowfrequency instability mechanism in coaxial dump combustors[END_REF], closely followed by [START_REF] Keller | Thermally induced low-frequency oscillations[END_REF]. They found that the time lag characteristic of temperature spots detached from the flame and convected to the nozzle could add a new family of eigenmodes affecting the overall stability of the system. However, as recalled by [START_REF] Yu | Low-frequency pressure oscillations in a model ramjet combustor[END_REF], there was a lack of experimental evidence on the contribution of a mixed, convective mode to the flame instability. Despite the fact that this latter work focused on the interaction between vortex and nozzle, the results can be extrapolated to the interaction between entropy and acoustic modes and constitute a first step to outline that the instability mechanism relies on the combination of acoustic and convective time lags. Contrary to the indirect combustion noise, the vortex shedding and its interaction with the acoustics is an important part of the stability of solid rocket engines and many studies are devoted to this subject [START_REF] Vuillot | Vortex-shedding phenomena in solid rocket motors[END_REF][START_REF] Chedevergne | Direct numerical simulation and biglobal stability investigations of the gaseous motion in solid rocket motors[END_REF]. The experimental study that enlightened the role of the indirect noise on combustion instability was performed by [START_REF] Macquisten | Low-frequency combustion oscillations in a model afterburner[END_REF], who observed that the flame burning in a simple tube could be strongly affected by the addition of a nozzle. However it was just an observation which was not followed by a modelling effort. The works cited above were motivated by the study of rocket engines and ramjets. The subject has taken a renewed gain of interest with the development of lean premixed combustor for gas turbines, especially in the aeronautical application where these engines cover a wide range of operating points and are more prone to combustion instabilities due to the technological choices. The complexity of such systems has naturally led to the development of low-order methods to predict thermo-acoustic instabilities. More comprehensive works about low-order modelling of combustion instability [START_REF] Keller | Thermoacoustic oscillations in combustion chambers of gas turbines[END_REF] (2006), one important issue is that the influence of indirect combustion noise is not systematic and a low-frequency instability may appear even without the presence of a nozzle, and vice-versa. Some elements of response were provided by Polifke et al. (2001a) who showed that entropy and acoustic waves may have a constructive or destructive phase dependency, as well as [START_REF] Sattelmayer | Influence of the combustor aerodynamics on combustion instabilities from equivalence ratio fluctuations[END_REF] who argued that entropy spots are submitted to strong dilution (mixing) due to the highly turbulent nature of the flow in practical combustors. Recent numerical study by [START_REF] Morgans | The dissipation and shear dispersion of entropy waves in combustor thermoacoustics[END_REF] suggests however that the amount of dissipation of entropy spots convected by a turbulent flow may remain moderate. The importance of the entropy wave dissipation in mixed modes introduced a new mechanism to consider in models which did not have to be considered for purely acoustic unstable modes which are only marginally affected by wave dissipation. The effects of dissipation and dispersion on unstable mixed modes have been investigated in a recent study by Goh & Morgans (2013). They confirmed that entropy waves can destabilise a stable combustor, but can also have a stabilising effect on an unstable configuration. Moreover, some configurations can exhibit a mode switch between unstable modes. However, despite the fundamental interest of these studies, they still rely on 1D low-order methods and simple academic configurations.

Only a few numerical studies have been performed to reproduce a mixed entropyacoustic combustion instability in an air-breathing combustor. [START_REF] Zhu | Self-excited oscillations in combustors with spray atomizers[END_REF] and [START_REF] Yao | Combustion rumble prediction with integrated computational-fluid-dynamics/low-order-model methods[END_REF] considered an academic combustor and applied the Reynolds Averaged Navier-Stokes (RANS) methods although the accuracy of this approach when dealing with reacting turbulent flows is highly questionable. Unsteady RANS computations have been performed by [START_REF] Hochgreb | Forced and self-excited instabilities from lean premixed, liquid-fuelled aeroengine injectors at high pressures and temperatures[END_REF] in an academic combustor. These authors found that the mechanism of convection of entropy spots to the downstream nozzle has a keyrole on the establishment of a self-sustained thermoacoustic oscillation, whose computed frequency match the experiments. Large-Eddy Simulation (LES) is now recognised as appropriate (although CPU demanding) to deal with flame dynamics related issues, with or without coupling with acoustics (see [START_REF] Pitsch | Large eddy simulation of turbulent combustion[END_REF][START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF]Wolf et al., 2012b;[START_REF] Gicquel | Large eddy simulations of gaseous flames in gas turbine combustion chambers[END_REF]. LES was used to investigate the noise generated by the impact of acoustic and entropy waves on nozzles or turbine stages (see [START_REF] Bodony | Scattering of an entropy disturbance into sound by a symmetric thin body[END_REF]Duran et al., 2013). However to the authors knowledge no LES study deals specifically with a mixed entropy-acoustic instability.

Thesis objectives and organisation

This thesis is devoted to the analysis and the modelling of mean flow convective effects on combustion instabilities. It is expected that Large-Eddy Simulations are able to reproduce such effects in realistic configurations, but the restitution time for results is very long, especially when dealing with low-frequency instability that requires recording a significant time signal. The strategy depicted in Fig. 4.b is then a good candidate to get similar results in a very short computational time, because the joint use of a Helmholtz solver combined with analytical models does not require a huge computational effort.

The Helmholtz equation (2) relies on the so-called zero-Mach number assumption stating that the mean velocity is very small compared to the speed of sound. [START_REF] Nicoud | About the zero mach number assumption in the calculation of thermoacoustic instabilitie[END_REF] suggest that the domain of validity of the zero-mean-flow assumption might be rather small. Nevertheless, this somewhat restrictive assumption is necessary to derive a wave equation for the thermoacoustic perturbations superimposed to a non isentropic baseline flow, the alternative being to use the complete set of Linearised Euler Equations. Unfortunately, this would make the computational effort needed to compute the thermoacoustic modes significantly larger (five coupled equations being solved for) than what is required when dealing with Eq. (2). Being able to (partly) account for the non zero-Mach number effects without relying on the LEEs is therefore desirable.

Chapter 1 focuses on the purely aeroacoustic effects induced by the presence of a mean flow in a configuration. It is shown that a mean flow has a damping effect on the acoustic modes. The theoretical derivation of a conservation energy balance for the acoustic perturbation shows that the energy absorption originating from the mechanical work during the convection process mostly appears through the energy flux formulation. When using a Helmholtz solver, it is shown that the zero-Mach number assumption introduces an energy flux mismatch at the boundary conditions. A theoretical study based on an analytical asymptotic expansion of the energy flux shows that it is possible to reduce the error introduced by the zero-Mach number assumption by imposing the right energy flux through boundaries via a proper formulation of the impedances. Comparisons between a LEE and a Helmholtz solver show that the damping rates of the eigenmodes is well estimated for different Mach number when using the proper impedance formulation.

Eq. (2) does not support entropy waves. Thus the acoustic generation due to the entropy spot being accelerated in the nozzle/turbine located downstream of the combustion chamber is not accounted for. As explained above the production of sound by acceleration of entropy fluctuations is a key phenomenon when dealing with combustion noise and neglecting this acoustic source when studying thermoacoustic instabilities is questionable. Chapter 2 is devoted to the development of a methodology to introduce the coupling between entropy and acoustic waves through boundary conditions. The methodology is decomposed into two parts: first, as the entropy variable does not exist in the computational domain, it must be modelled by a so-called transfer function. Then, the entropy variable is related to the acoustic variable via the imposition of a proper boundary condition, which has been named Delayed Entropy Coupled Boundary Condition (DECBC).

In chapter 3, a realistic gas turbine combustor provided by SAFRAN is studied. During experiments, this combustion chamber exhibited a strong low-frequency instability for some operating points and was stable otherwise. Large-Eddy Simulations are performed on two operating points, one which exhibits an instability and another one which is stable. LES is able to reproduce the behaviour of the combustor and results in terms of frequency and pressure amplitudes are close to the experiments. In order to understand the mechanism that underlies the combustion instability, a Dynamic Mode Decomposition is performed and shows that the instability is driven by the convection of temperature pockets from the flame to the exit nozzle. Hence, this configuration is a good candidate to apply the DECBC method within a zero-Mach number framework. The Helmholtz computation coupled with the DECBC approach are able to reproduce the instability in a short computational time. This demonstrates that this instability is indeed driven by the noise generated by strong entropy fluctuations being convected through the nozzle. To the authors knowledge, this is the first time that such instability is identified in an aeroengine configuration. Finally, a sensitivity study with the Helmholtz solver coupled with DECBC is performed to optimise the shapes of the downstream nozzle. It is shown that for certain geometry of the nozzle, the growth rate of the unstable mode can be reduced. LES computations with an extended geometry of the nozzle confirm that the combustion instability can be significantly attenuated with the proposed methodology.

A schematic representation of the organisation of the thesis is depicted in Fig. 5. Note that some parts of the geometry are hidden for confidentiality reasons. In order to assess combustion instabilities in a gas turbine engine, the strategy developed in the present thesis relies on the use of a Helmholtz solver so as to compute the thermoacoustic problem inside the combustor. Such a method solves the acoustic pressure equation as an eigenvalue problem and provides both the frequency and growth rate of the associated instability. However the zero-Mach number assumption used to derive the Helmholtz equation is not applicable to the high-speed elements surrounding the combustor, as for example the compressor or the turbine. As stated in the introduction, the chosen strategy is to model such elements under the form of an impedance that will be applied as a boundary condition to the Helmholtz domain (cf. The purpose of this chapter is to assess the proper way to formulate an impedance at the interface of the zero/non-zero Mach number domains. Before any mathematical development, the properties of acoustics in a non-quiescent flow are briefly reviewed in §1.1 by taking some key examples from the aeroacoustic literature. Equations and impedance formulations are proposed in §1.2, while §1.3 is devoted to the modelling of a simple application case. Finally results are presented and discussed in §1.4.

About acoustics in a non-quiescent flow

Acoustics is a discipline that is found in practically all scientific and engineering applications as well as in the everyday life of humans. As recalled by [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF] acoustics is the science of sound, where the term sound is not only about the human sensation of hearing, but is referred to any mechanical wave motion. Basically he defines the notion of a wave as an oscillatory disturbance that propagates from a source. This somewhat rudimentary statement goes back to the ancient Greek civilisation and the legend says that philosophers such as Pythagoras (550 B.C.) were teaching mathematics with the help of musical acoustic theories. The physical laws of acoustics were first established in the 17th century by Galileo and Marin Mersenne, this latter one also referred as the father of acoustics (see [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF][START_REF] Bohn | Environmental effects on the speed of sound[END_REF]. The advances in mathematics and calculus techniques during the 19th century allowed the emergence of strong theories of acoustics, especially with the monumental works of Helmholtz (1877) and [START_REF] Rayleigh | The Theory of Sound . Mac Millan[END_REF]. The industrial revolution and the increasing technological developments have allowed acoustics to spread in all areas and be an actor of major revolutionary advances, such as sound recording and transmissions with the telephone or the radio, detection of submarines, medical care or, for our concern, knowledge in combustion instabilities.

Most practical applications deal with the basic aspects of acoustics where the waves propagate in a medium at rest. Several reference textbooks (see [START_REF] Morse | Theoretical acoustics[END_REF][START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF] are dedicated to the underlying theories but the propagation of acoustic waves in a moving media is succinctly treated. When the effects of a flow become non negligible as to interact with the acoustics, physical theories and applications move to the so-called branch of aeroacoustics (see [START_REF] Howe | Theory of vortex sound[END_REF][START_REF] Rienstra | An introduction to acoustics[END_REF]. A close look at the literature shows that the acoustics is split between theories that take into account the mean flow and others that neglect it. This is particularly of our interest because the present study concerns the interface where a zero-Mach number domain coincide with a finite Mach number element.

The separation between communities that take into account the mean flow or not is particularly put in evidence by [START_REF] Ingard | Effect of flow on the acoustic resonances of an open-ended duct[END_REF], who explain in a simple manner the effects of a mean flow on the acoustics of a simple tube with an open end. To the authors knowledge, this is one of the first theoretical and experimental work that investigate the mean flow effects in term of damping of eigenmodes as well as the impact of the Mach number on the acoustic reflection coefficient (ratio between a reflected and an impacting acoustic wave) at the open end. This latter point is crucial, because as already mentioned by [START_REF] Meyer | Physical and applied acoustics: an introduction[END_REF], the acoustic reflection coefficient takes values superior to unity when the Mach number increase. This is particularly surprising since it would mean that an acoustic wave traveling in the tube to the open end output comes back amplified. The 1.1 About acoustics in a non-quiescent flow authors suggest that the vorticity generated at the edge of the end of the tube could be converted into noise that enters back into the tube. It is also argued that the reflection coefficient should be taken in the sense of energy rather than pressure. The energy reflection coefficient is then formulated by the multiplication of the pressure reflection coefficient by a term (1 -M ) (1 + M ) where M is the Mach number and should not exceed unity. This is partially true but it will be demonstrated later either experimentally by [START_REF] Peters | Damping and reflection coefficient measurements for an open pipe at low mach and low helmholtz numbers[END_REF] or theoretically by [START_REF] Munt | Acoustic transmission properties of a jet pipe with subsonic jet flow: I. the cold jet reflection coefficient[END_REF] that depending on the geometry of the end of the tube, the energy reflection coefficient can also take values superior to unity under certain flow conditions.

These observations are controversial because a confusion is made between the sound generated by the flow, by turbulence for example, and the fact that the mean flow itself has an impact on the wave propagation. The correction term (1 -M ) (1 + M ) refers to this latter point. This is also why some authors [START_REF] Ingard | Effect of flow on the acoustic resonances of an open-ended duct[END_REF][START_REF] Davies | Reflection coefficients for an unflanged pipe with flow[END_REF][START_REF] Davies | Plane wave reflection at flow intakes[END_REF] apply empirical correction terms to their formulations of the reflection coefficient to capture the sound generation coupled with the flow and to fit experiments. In the works of [START_REF] Munt | Acoustic transmission properties of a jet pipe with subsonic jet flow: I. the cold jet reflection coefficient[END_REF] the sound generated by the jet outside the tube is explicitly a source term in the formulation of the reflection coefficient and is responsible for the peak response superior to unity, despite the correction term (1 -M ) (1 + M ) that models the damping flow effects inside the tube.

On the other hand, the introduction of flow effects in the acoustics has also been widely debated during the establishment of a theory of the aerodynamic sound. As reviewed by [START_REF] Howe | Contributions to the theory of aerodynamic sound, with application to excess jet noise and the theory of the flute[END_REF], the first major contribution goes back to the works of [START_REF] Lighthill | On sound generated aerodynamically: I. General theory[END_REF][START_REF] Lighthill | On sound generated aerodynamically: II. Turbulence as a source of sound[END_REF] which derived a stationary, scalar inhomogeneous wave operator where the sound generation is shifted in a source term. The main drawback is that the mean flow in not taken into account in the wave operator, the fluid being assumed at rest and extended to the entire space, and the interactions between acoustics and the flow are encapsulated in a distribution of equivalent sources terms which must be known a priori. The mean flow has been introduced in the wave operator by [START_REF] Phillips | On the generation of sound by supersonic turbulent shear layers[END_REF] which derived an inhomogeneous convected wave equation for the pressure where the moving flow velocity appears in the time derivative. Unfortunately, as pointed out by [START_REF] Lilley | On the noise from jets[END_REF], this latter approach is based on the hypothesis that the convected wave operator correctly describes the propagation of small disturbances through a medium of variable mean properties, which is not usually correct especially in shear-flow problems. Indeed, theses cases contain sources terms that are linear in the acoustic perturbation and should be included in the propagation operator.

The derivation of a generalised acoustic field theory is a controversial subject, especially when dealing with reactive flow and is not the matter of the present thesis. The reader is invited to investigate the literature (see [START_REF] Powell | Some aspects of aeroacoustics: From rayleigh until today[END_REF][START_REF] Tam | Jet noise: Since[END_REF][START_REF] Goldstein | A generalized acoustic analogy[END_REF][START_REF] Rienstra | An introduction to acoustics[END_REF]). An important point that emerges from the theories derived from the Lighthill's analogy is that they are convenient when dealing on the sound produced by the flow, but fails to considers the counter part of a the flow perturbed by the acoustic field. When such a feedback is significant (for example when dealing with homentropic low Mach number flow), the concept of vortex sound introduced by Powell (1964) is most appropriate. Based upon this statement, [START_REF] Howe | Contributions to the theory of aerodynamic sound, with application to excess jet noise and the theory of the flute[END_REF] proposes to define the acoustic field as an unsteady scalar potential flow. From this point the pressure variable does not satisfy the inhomogeneous convected wave equation and [START_REF] Howe | Contributions to the theory of aerodynamic sound, with application to excess jet noise and the theory of the flute[END_REF] proposes to use rather the fluctuating total stagnation enthalpy as state variable. Under that form, the vorticity is directly responsible for the generation of sound and appears explicitly in the source terms of the wave equation.

The use of the fluctuating total enthalpy as a state variable to derive a generalised acoustic field has been mathematically demonstrated by [START_REF] Doak | Fluctuating total enthalpy as the basic generalized acoustic field. Theoret[END_REF]. The choice of such a state variable appears in a systematic manner in aeroacoustic studies. To cite a few, Hofmans et al. What emerges from this brief overview is that the choice of the fluctuating pressure variable, which is natural when dealing with acoustics in a zero-Mach number domain such as a Helmholtz solver, is questionable when computing the impedance of high-speed elements such as a turbine or a compressor. As explained in this overview, the subject is controversial and the use of another set of state variables is a good candidate to solve this question. However, there is a need to develop a theoretical background that supports this statement. In the following section, the problem of connecting non-zero and zero Mach number elements is then addressed with the help of an energy flux analysis which suggests a method to introduce non zero-Mach effects at nozzles terminating combustion chamber in a zero Mach number solver.

Energy considerations

An important property of boundary conditions in an acoustic solver is to ensure that the energy flux that crosses the frontiers is correctly predicted. An interesting starting point is to go back to the definition of the principle of the energy conservation, which may be written into the general form:

∂E ∂t + ∇ ⋅ F = P (1.1)
where E represents the energy density evolution in a medium, F the energy flux through boundaries and P a source term for the rate of energy production or dissipation per unit volume. The form of Eq. (1.1) is convenient, because it becomes a conservation law when integrated over an arbitrary fixed volume V and the term ∇ ⋅ F can be re-expressed as a surface integral by means of the Gauss' theorem:

V ∇ ⋅ F dV = S F ⋅ n dS (1.2)
where n is the outward unit vector that defines the normal direction of the surface S. As pointed out by [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF], such equation is an energy corollary, which means that it is just an intellectual concept that reformulates the fundamental acoustic equations. This corollary adds nothing to the physics but facilitates the description of the problem by a phenomenological modelling. In other words, Eq. (1.1) represents a temporal balance between a source or a sink of energy inside the domain and the incoming or leaving energy flux through the boundaries. Such a balance is then a useful tool to derive a stability criterion of the system, especially in combustion to predict if the coupling between a flame and the acoustics could lead to an uncontrolled accumulation of energy inside a combustor.

Forms of energy corollary

However the main drawback of this energy corollary is that it depends on what is used for E and that multiple possibilities exist for E [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF][START_REF] Nicoud | Thermoacoustic instabilities: should the rayleigh criterion be extended to include entropy changes ?[END_REF].

Basically the main idea is to start from the conservation equations of mass, momentum and energy. The reader is invited to consult the classical textbooks for more details on the origin of these equations (see [START_REF] Batchelor | An introduction to fluid dynamics[END_REF]Landau & Lifchitz, 1987). In the absence of body forces, it gives the following set of equations:

∂ρ ∂t + ∂(ρu i ) ∂x i = 0 (1.3) ∂u j ∂t + u i ∂u j ∂x i + 1 ρ ∂p ∂x j = 1 ρ ∂τ i,j ∂x i (1.4) ∂s ∂t + u i ∂s ∂x i = Q ρ (1.5)
where ρ is the density, u i,j are the components of the velocity vector u, p is the thermodynamic pressure and s is the specific entropy. Q is an energy source term related to the heat addition to the volume and τ i,j denotes the components of the viscous stress tensor.

As the acoustic discipline focuses on perturbations in a flow, an intense discussion on the appropriate expression of the energy corollary has shaken the community during the last decades, especially because the different proposed expressions do not hold for all applications. Two major corollaries emerged:

• The principle of acoustic energy has been developed first by [START_REF] Cantrell | Interaction between sound and flow in acoustic cavities: Mass, momentum, and energy considerations[END_REF] for homentropic, irrotational flows in a cavity, and generalised later by [START_REF] Morfey | Acoustic energy in non-uniform flows[END_REF] for non-isentropic, heat-conducting viscous flows. The idea is to linearise the set of conservation laws governing the fluid motion described at Eqs. (1.3) to (1.5). Each quantity φ(x, t) is decomposed into a mean and a fluctuating part φ(x, t) = φ(x) + φ ′ (x, t) where the hypothesis made is that fluctuations are small compared to the mean values (i.e. φ ′ << φ) and that mean values are not time dependant. The linearisation at the first order of Eqs. (1.3) to (1.5) then becomes:

∂(ρ + ρ ′ ) ∂t + (u + u ′ )∇(ρ + ρ ′ ) + (ρ + ρ ′ )∇ ⋅ (u + u ′ ) = 0 (1.6) (ρ + ρ ′ ) ∂(u + u ′ ) ∂t + (ρ + ρ ′ )(u + u ′ ) ⋅ ∇(u + u ′ ) + ∇(p + p ′ ) = Ψ ′ (1.7) ∂(s + s ′ ) ∂t + (u + u ′ )∇(s + s ′ ) = r q + q ′ p + p ′ (1.8)
where Ψ is defined such that Ψ j = (1 ρ)(∂τ i,j ∂x i ). The energy source term Q in Eq. (1.5) has been reformulated as Q = q T , where T the absolute temperature, q the rate of heat release per unit of volume and r the universal gas constant that appears after application of the state equation p = ρrT . Of course the mean flow field is supposed to obey to the conservation law, leading to the following set of equations:

Dρ Dt = -ρ∇ ⋅ u (1.9) ρ Du Dt = -∇p + Ψ (1.10) Ds Dt = r q p (1.11) where D(⋅) Dt is the particular derivative ∂(⋅) ∂t + u ⋅ ∇(⋅).
It is then possible to combine Eqs. (1.6) to (1.8) so as to reformulate the following second order acoustic energy corollary:

∂E a ∂t + ∇ ⋅ F a = P a (1.12) with E a = 1 2 p ′ 2 ρ c 2 + 1 2 ρu ′ 2 + u ⋅ u ′ p ′ c 2 (1.13) F a = p ′ u ′ + u ρ c 2 p ′ 2 + u 2 c 2 p ′ u ′ + ρu ⋅ u ′ 2 (1.14)
One must note that Morfey (1971) splits the disturbance velocity field into irrotational and solenoidal components which are defined as the acoustic (u) and the unsteady vortical (w) motions, respectively. At boundaries, the perturbations in the energy flux F a only involve the fluctuating pressure and acoustic velocity variables. Any resulting entropy or vortical disturbances are shifted to the source term P a , not described here for sake of simplicity, but extensively detailed in the study of [START_REF] Bloxsidge | Active control of reheat buzz[END_REF]. However, this definition is problematic, especially in reacting flows, because velocity disturbances may be associated with convective disturbances such as entropy heterogeneities. One should note that the Rayleigh criterion term p ′ q ′ appears in P a . Hence, from Eq. (1.12), the stability of a thermoacoustic system is driven by the balance between energy sources/sink terms contained in P a and energy flux F a at boundaries.

• However, as pointed out by [START_REF] Candel | Acoustic conservation principles and an application to plane and modal propagation in nozzles and diffusers[END_REF], this principle may be not completely satisfactory, especially when dealing with geometrical acoustic in nozzles or reacting configurations that involve thermal related energy terms. This leads to the other principle of disturbances energy, derived by [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF] and [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF] for zeromean flow and generalised by [START_REF] Myers | Transport of energy by disturbances in arbitrary steady flows[END_REF] for arbitrary disturbances in steady flows.

The main idea of this principle is that the equation for the energy associated with disturbances is not derived by handling algebraic combination of the linearised mass and momentum equations, but directly by expanding and linearising the following general energy equation for the fluid:

∂(ρJ -p) ∂t + ∇ ⋅ (Jm) -m ⋅ Ψ -T Q = 0 (1.15)
which introduces the total stagnation enthalpy J and mass flux m, defined as:

J = h + u 2 2 and m = ρu (1.16)
where h = e + p ρ is the specific enthalpy and e the internal energy.

After quite heavy algebra it then comes the following expression of the disturbance energy density E d and energy flux F d :

E d = 1 2 p ′ 2 ρ c 2 + 1 2 ρu ′ 2 + u ⋅ u ′ p ′ c 2 + ρ T s ′ 2 2C p (1.17)
F d = m ′ (J ′ -T s ′ ) -m ′ τ ρ ′ + m T ′ s ′ + T ′ Q T ′ (1.18)
The disturbances energy principle considers that not only acoustic perturbations are part of the disturbances energy, but that a proper energy also involves correlations terms between entropy and acoustic fluctuations that appear into the energy density and energy flux.

Currently, the establishment of an universal formulation for the energy of fluctuations has not yet reached a consensus. It is especially true in the combustion community because this raises the question of the choice of a proper criterion to study the stability of thermoacoustic systems. [START_REF] Nicoud | Thermoacoustic instabilities: should the rayleigh criterion be extended to include entropy changes ?[END_REF] argued that the so-called Rayleigh criterion, derived from an energy balance within the zero-Mach number assumption and based on the product of fluctuations of pressure and rate of heat release, should be wrong and must be extended to includes fluctuations of entropy. This idea was also demonstrated by [START_REF] Brear | Disturbance energy transport and sound production in gaseous combustion[END_REF] who extended the works of [START_REF] Myers | Transport of energy by disturbances in arbitrary steady flows[END_REF] to reacting flows by including species and chemical non-equilibrium, and computed the disturbance energy balance with the help of Direct Numerical Simulations of a premixed flame. In addition, whatever the formulation of the Rayleigh criterion, it must be noted that it represents an incomplete description of the overall stability of the thermoacoustic system and must be completed by a close look at the energy flux leaving or entering the domain.

Recall that the purpose of this chapter is to assess the proper way to formulate an acoustic impedance at the interface of zero/non-zero Mach number domains and that the present study is restricted to homentropic irrotational flows. As shown by [START_REF] Myers | Transport of energy by disturbances in arbitrary steady flows[END_REF], in such as particular case, the disturbances energy principle reduces to the acoustic energy principle of [START_REF] Cantrell | Interaction between sound and flow in acoustic cavities: Mass, momentum, and energy considerations[END_REF] because all terms related to entropy or vorticity are dropped to zero. In other words, Eq. (1.18) reduces to Eq. (1.14). As Eq. (1.1) holds at each location and time, it is however convenient to consider integrated versions of this equation. Since we restrict ourselves to the linear approximation, the instantaneous value of the acoustic energy is not as important as its evolution after one period of oscillation T . Hence the energy flux through a boundary with the presence of a mean flow takes the following form:

F M ≠0 = T 0 S 0 p ′ u ′ + u ρ c 2 p ′ 2 + u 2 c 2 p ′ u ′ + ρu ⋅ u ′ 2 ⋅ n dS dt (1.19)
Introducing first J ′ and m ′ the linearised isentropic fluctuating parts of the total enthalpy (J = C p T + u 2 /2) and the mass flow rate (m = ρu), respectively:

J ′ = p ′ ρ + u ⋅ u ′ (1.20) m ′ = ρ u ′ + u p ′ c 2 (1.21)
Eq. (1.19) may then be recast to the following expression:

F M ≠0 = T 0 S J ′ m ′ ⋅ n dS dt (1.22)
Considering acoustics in the zero-Mach number limit (u = 0) leads to the following expression, specific to the Helmholtz equation:

F M =0 = T 0 S p ′ u ′ ⋅ n dS dt (1.23)

Impedances and flux mismatch

Following the approach described in Fig. 4.b, only the combustion chamber is solved for using a Helmholtz solver while finite-Mach-number effects in the compressor, turbine, nozzle, etc. are accounted for only through boundary conditions by solving the LEE either analytically (Duran & Moreau, 2013) or numerically [START_REF] Lamarque | Boundary conditions for acoustic eigenmodes computation in gas turbine combustion chambers[END_REF]. Equation (1.14) suggests that the natural variables for a Helmholtz solver are p ′ and u ′ while Eq. (1.22) indicates that J ′ and m ′ are more attune for LEE and provide the right formulation for acoustic losses, including non zero-Mach effects. Note that these two sets of independent variables are strictly equivalent and one may write the LEE of Helmholtz equations with any set. An impedance, Z may also be defined for each set of variables.

Energy considerations

Note that impedances being defined in the frequency domain, any fluctuating quantity g ′ (x, t) must first be recast as g ′ (x, t) = R ĝ(x)e -jωt , where ĝ is the complex amplitude, and ω the angular frequency.

Impedances formulation

An energy flux or an impedance represent physical quantities across a boundary and involve the velocity which is defined along a direction. Introducing n a unit vector normal to the boundary, the reduced impedance Z may be defined along the n direction for each set of variables. Using p and û one has:

Z p,u = p ρcû ⋅ n (1.24)
while using Ĵ and m one may define:

Z J,m = ρ c Ĵ m ⋅ n (1.
Z J,m = ρ p ρ + u ⋅ û c ρû ⋅ n + u ⋅ n c 2 p (1.26) = p + ρ(u ⋅ û) ρc(û ⋅ n) + u ⋅ n c p (1.27)
Let's divide all terms by ρc(û ⋅ n). One then has:

Z J,m = p ρcû ⋅ n + u ⋅ û cû ⋅ n 1 + u ⋅ n c p ρcû ⋅ n (1.28)
Using Eq. (1.24) leads to the following formulation of the impedance Z J,m as a function of Z p,u :

Z J,m = Z p,u + u ⋅ û cû ⋅ n 1 + u ⋅ n c Z p,u (1.29)
Care must be taken with the term (u ⋅ û) (c û ⋅ n). Considering first a 3D space and introducing n the unit vector normal to the boundary and t 1 ,t 2 two unitary vectors such that (n, t 1 , t 2 ) forms an orthonormal basis, the scalar expansion of the product u ⋅ û may be written into the form:

u ⋅ û = ((u ⋅ n) n + (u ⋅ t 1 ) t 1 + (u ⋅ t 2 ) t 2 ) ⋅ ((û ⋅ n)n + (û ⋅ t 1 )t 1 + (û ⋅ t 2 )t 2 ) (1.30)
Introducing first the following expressions:

u n = u ⋅ n (1.31) u 1 = u ⋅ t 1 (1.32) u 2 = u ⋅ t 2 (1.33)
Eq. (1.30) may then be recast as:

u ⋅ û = u n ûn + u 1 û1 + u 2 û2 (1.34)
Finally it comes:

u ⋅ û c û ⋅ n = u n c + u 1 c û1 ûn + u 2 c û2 ûn (1.35)
The hypothesis is made that the mean velocity vector u is always normal to the boundary, i.e. u 1 = u 2 = 0. After introduction of the Mach number vector M = u c, one then has :

u ⋅ û c û ⋅ n = M ⋅ n (1.36)
Finally, the proper relation between Z p,u (Eq. (1.24)) and Z J,m (Eq. (1.25)) takes the following form :

Z J,m = Z p,u + M ⋅ n 1 + Z p,u M ⋅ n (1.37)
Similarly, it comes:

Z p,u = M ⋅ n -Z J,m Z J,m M ⋅ n -1 (1.38)
The situation depicted in Fig. 1.1 where a Helmholtz domain (M = 0) is connected to a LEE domain (M ≠ 0) through an interface raises the question of the consistency between the two sub-domains. More precisely, let us assume that the impedance of the LEE domain is known at the interface, noted Z M ≠0 p,u (or Z M ≠0 J,m ). This impedance models the way acoustic waves propagate into and reflect from the LEE domain. Since the flow is not at rest in this region, Z M ≠0 p,u and Z M ≠0 J,m are obviously different (Eqs. (1.37)-(1.38)) although they represent the same LEE domain. When solving for the Helmholtz equation in the M = 0 domain, a boundary impedance Z M =0 must be imposed at the interface in order to represent the effect of the LEE domain. Note that in the Helmholtz domain the Mach number is assumed zero so that from Eqs. (1.37)-(1.38):

Z M =0 = Z M =0 p,u = Z M =0 J,m .
The next natural question is then: how should Z M =0 be defined from the knowledge of Z M ≠0 p,u (or Z M ≠0 J,m ) ? Is it possible to use a Z p,u impedance in the Helmholtz solver which would provide the correct energy flux through the boundary?

Flux mismatch introduced by the impedance

The purpose of the following is to assess the flux difference F M ≠0 -F M =0 at the interface between the Helmholtz domain and the outer domain where the baseline flow is not assumed at rest (see Fig. 1.1). The energy flux F M ≠0 relative to the LEE domain may be reformulated with the help of the impedance Z M ≠0 J,m defined at Eq. (1.25). Introducing first for clarity reasons the notation {ĝ} = R ĝe -jωt and considering the spatial integration over a uniform boundary of section S, Eq. (1.22) then becomes:

F M ≠0 = S T 0 { m ⋅ n} Ĵ dt = S T 0 c ρ { m ⋅ n} Z M ≠0 J,m m ⋅ n dt
The LEE solution is made non-dimensional by setting m = m0 at the considered boundary, leading to:

F M ≠0 = ρcS T 0 m0 ⋅ n ρ Z M ≠0 J,m m0 ⋅ n ρ dt (1.39)
The energy flux F M =0 relative to the zero-Mach number domain may be reformulated with the help of the impedance Z M =0 . Eq. (1.23) then becomes:

F M =0 = S T 0 {û ⋅ n} {p} dt = ρcS T 0 {û ⋅ n} Z M =0 û ⋅ n dt (1.40)
In the Helmholtz domain m = ρû so that the acoustic field is made non-dimensional by setting at the boundary û = m0 ρ. The energy flux F M =0 then becomes:

F M =0 = ρcS T 0 m0 ⋅ n ρ Z M =0 m0 ⋅ n ρ dt (1.41)
The error introduced by the modelling of Z M =0 can be assessed by forming the flux difference ∆F = F M ≠0 -F M =0 where the flux expression Eq. (1.39) is used together with Eq. (1.41). Different flux differences are obtained depending on the modelling choice made for Z M =0 :

• If the impedance Z M =0 is set to Z M ≠0 p,u :
Asymptotic expansion with respect to the Mach number leads to:

Z M =0 = Z M ≠0 p,u = M ⋅ n -Z M ≠0 J,m M ⋅ nZ M ≠0 J,m -1 = Z M ≠0 J,m + O(M ) (1.42)
Injection of Eq. (1.42) into Eq. (1.41) leads to the following energy flux formulation for the zero-Mach number domain:

F M =0 = ρcS T 0 m0 ⋅ n ρ Z M ≠0 J,m m0 ⋅ n ρ dt + O (M ) (1.43)
so that, taking the difference between Eq.(1.39) and Eq.(1.43) and keeping only lower order terms, the flux difference is then:

∆F = O (M ) (1.44) • If the impedance Z M =0 is set to Z M ≠0 J,m :
The equation (1.41) then becomes:

F M =0 = ρcS T 0 m0 ⋅ n ρ Z M ≠0 J,m m0 ⋅ n ρ dt (1.45)
so that the theoretical flux difference is simply ∆F = 0. However, neglecting the mean flow in the zero-Mach-number domain also introduces a mismatch in the frequency 1 T of oscillation of order M 2 .

Indeed, considering first a domain of length L under the zero-Mach number assumption, a simple reasoning shows that the first eigenmode is composed of a downstream T + = L c and an upstream T -= L c travelling acoustic wave which propagates at the speed of sound c. Hence the frequency of oscillation is:

f 0 = 1 T = 1 (T + + T -) = c 2L (1.46)
Following the same reasoning, the introduction of the mean flow velocity u leads to:

T + = L c + u = L c(1 + M ) and T - = L c -u = L c(1 -M ) (1.47)
This leads to the following estimation of the frequency of oscillation:

f ≈ c(1 -M 2 ) 2L (1.48)

Energy considerations

Hence it gives the following relation between the frequencies of oscillation with M ≠ 0 and the zero-Mach number assumption M = 0:

f ≈ (1 -M 2 )f 0 (1.49)
Consequently the flux difference is actually:

∆F = O M 2 (1.50)
The two approaches to model Z M =0 are summarised in Table 1.1, which also displays the asymptotic behaviours of the difference between the exact energy flux F M ≠0 (see Eq. (1.22)) and its zero-Mach number counterpart F M =0 (see Eq. (1.23)). A perfect model for the boundary impedance Z M =0 would produce zero difference between these two fluxes, meaning that the acoustic energy leaving the Helmholtz domain would exactly match the amount of energy entering the LEE domain. However, assuming the mean flow being at rest in the computational domain introduces a loss in acoustic energy conservation and a mismatch between F M =0 and F M ≠0 . The theoretical asymptotic analysis shows that the Z M ≠0 p,u approach produces an error F M ≠0 -F M =0 proportional to the Mach number between the Helmholtz and the LEE domains (noted O(M ) in Table 1.1). At the same time, the Z M ≠0 J,m approach produces a flux mismatch which is quadratic in the Mach number ( F M ≠0 -F M =0 = O(M 2 )), thus a smaller error in the acoustic energy conservation.

Z M =0 Z M ≠0 p,u Z M ≠0 J,m F M ≠0 -F M =0 O(M ) O(M 2 )
Table 1.1 -Summary of the flux mismatch error introduced in the zero-Mach number domain by the choice of boundary impedance Z M =0 .

The next section is devoted to the evaluation of the two approaches to model Z M =0 . The eigenmodes of a simple configuration are computed either by solving the full set of LEE, either by solving a truncated domain under the zero-Mach-number assumption. The influence of the modelling approach can then be assessed by comparison with results from the LEEs.

Testing impedance formulation in a simple case

The configuration considered to test the impedance expression of Table 1.1 is shown in Fig. 1.2.a. It consists of two connected tubes of section S 1 and S 2 . The subscripts 1 and 2 refer to parameters in the left and right tube, respectively. The boundary conditions are m1 = 0 and ŝ1 = 0 at the entrance of tube 1 and p2 = 0 at the exit of tube 2, respectively. This simple test case is chosen so as to mimic a burner where the inlet and outlet Mach numbers are different, but without entropy sources at all so as to keep the isentropic assumption.

The eigenmodes of this simple configuration are computed by three different methods:

• Method I is the analytical solution of the LEE in the whole domain: x ∈ [0, L] (see Fig. 1.2.a). This is considered as the reference solution and will be crossvalidated with a numerical LEE solver described at section 1.3.1.

• Method II is the solution of the LEE equations in a restricted region of the domain:

x ∈ [l z1 , l z2 ] (see Fig. 1.2.b). This solution shows that it is appropriate to solve the problem in a restricted part of the configuration if the proper impedances are set at x = l z1 and l z2 .

• Method III is the solution of the Helmholtz equation (i.e. assuming that M = 0) in a restricted region of the domain (see Fig. 1.2.b). This solution assesses the influence of the zero-Mach-number assumption and is used to test the two modelling approaches of Z M =0 summarised in Table 1.1. To this end, two versions of Method III are used in the following:

III pu when Z M =0 = Z M ≠0 pu and III Jm when Z M =0 = Z M ≠0 Jm .
Table 1.2 summarizes the analytical procedures described above.

Method Computed domain Mean flow assumption Impedance formulation (Z 1 and Z 2 )

I Full M ≠ 0 × II Truncated M ≠ 0 Z M ≠0 p,u III pu Truncated M = 0 Z M ≠0 p,u III Jm Truncated M = 0 Z M ≠0 J,m

Method I: Global domain at M ≠ 0

Method I num : Numerical LEE solver For an homogeneous reacting mixture with constant heat capacities C p and C v , the linearised harmonic form of conservation equations for mass (Eq. (1.3)), momentum (Eq. (1.4)) and energy (Eq. (1.5)) in a quasi-1D domain of cross section area S(x) read respectively [START_REF] Nicoud | About the zero mach number assumption in the calculation of thermoacoustic instabilitie[END_REF]:

û ∂ρ ∂x + ρ ∂ û ∂x + u ∂ ρ ∂x + ρ ∂u ∂x + ρû + ρu S ∂S ∂x = jω ρ (1.51) 1 ρ ∂ p ∂x + û ∂u ∂x + u ∂ û ∂x + uρ ρ ∂u ∂x = jω û (1.52) û ∂s 0 ∂x + u ∂ŝ ∂x + rq p 2 p - r p q = jωŝ (1.53)
with the linearised state equation and entropy expression (see [START_REF] Shapiro | The Dynamics and Thermodynamics of Compressible Fluid Flow[END_REF]; [START_REF] Anderson | Modern Compressible Flow[END_REF]:

p p - ρ ρ - T T = 0 and ŝ C v = p p -γ ρ ρ (1.54)
where γ and C v are the heat capacity ratio and heat capacity at constant volume, respectively.

Assuming that the unsteady heat release rate amplitude q is modelled as a linear operator, i.e. q = q ρ ρ + q û û + q ŝ ŝ, Eqs. (1.51) -(1.53) define an eigenvalue problem MW = jωW with

M = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∂u ∂x + u ∂ ∂x + u S ∂S ∂x ∂ρ ∂x + ρ ∂ ∂x + ρ S ∂S ∂x 0 1 ρ ∂c 2 ∂x + u ρ ∂u ∂x + c 2 ρ ∂ ∂x ∂u ∂x + u ∂ ∂x (γ -1) T 1 p ∂p ∂x + ∂ ∂x γrq ρp - r p q ρ ∂s ∂x - r p q û u ∂ ∂x + (γ -1) q p - r p q ŝ ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1.55)
where the eigenvector is W = (ρ, û, ŝ) T and ω denotes the angular frequency.

The system MW = jωW can be solved as an eigenvalue problem by a numerical discretisation of the operator M and with the help of linear algebra algorithms such as the Arnoldi method. It provides the complex eigenfrequencies and associated eigenvector W that contains the mode shapes. The period of oscillations as well as the growth/decay rate of the eigenmode are contained in the real and imaginary part of each eigenfrequency, respectively. In the present study, the code developed by [START_REF] Wieczorek | Numerical study of mach number effects on combustion instability[END_REF] during her thesis has been used. This code uses a discretisation of the acoustic quantities based on a central finite difference scheme and with a staggered grid method.

Testing impedance formulation in a simple case

Method I an : Global analytical acoustic model at M ≠ 0

In the case of a ducted system with constant-cross-section tubes, it is more practical to describe the fluctuations of pressure, velocity and entropy (respectively p, û and ŝ) with wave amplitudes: A + and A -for the forward and backward propagating acoustic waves, respectively and σ for the entropic wave. One then has:

p(x) = A + e jk + x + A -e -jk -x (1.56) û(x) = 1 ρc [A + e jk + x -A -e -jk -x ] (1.57) ŝ(x) = σe jksx (1.58)
where

k + = ω c + u = k 1 + M (1.59) k - = ω c -u = k 1 -M (1.60) k s = k M (1.61)
and k = ω c is the acoustic wave number. Note that for the sake of simplicity the normal is taken along the x-axis direction and the velocity variable is no longer written in a vector notation. After injection of Eq. The one-dimensional formalism can also be used to reformulate the impedances previously described at Eqs. (1.24) and (1.25) with the help of a reflection coefficient R:

R p,u = A - A + and R J,m = A -(1 -M ) A + (1 + M ) (1.63)
where the amplitudes A + and A -have been defined in Eqs. ( As shown in Fig. 1.2.a, the forward and backward waves (cf. Eqs. (1.56) and (1.57)) traveling in both tubes lead to four unknown: A + 1 , A - 1 , A + 2 and A - 2 . These wave amplitudes are solutions of a homogeneous linear system of equations obtained by requiring that the boundary conditions (x = 0 and x = L) and jump condition (x = l) are fulfilled. The two boundary conditions form a first set of relations between waves amplitudes:

B.C Z p,u R p,u Z J,m R J,m p = 0 0 -1 M - 1 -M 1 + M û = 0 ∞ 1 1 M 1 -M 1 + M Ĵ = 0 -M - 1 + M 1 -M 0 -1 m = 0 - 1 M 1 + M 1 -M ∞ 1
Table 1.3 -Impedances and reflection coefficients for canonical boundary conditions: constant pressure (p = 0), constant velocity (û = 0), fixed mass flow rate ( m = 0) and fixed total enthalpy ( Ĵ = 0).

• m1 = 0, at x = 0 yields

A + 1 (1 + M 1 ) -A - 1 (1 -M 1 ) = 0 (1.64) • p2 = 0, at x = L yields A + 2 e jk + 2 L + A - 2 e -jk - 2 L = 0 (1.65)
As explained by [START_REF] Davies | Practical flow duct acoustics[END_REF] the jump relations for the section change located at x = l can be expressed by integrating over an infinitesimal control volume the conservation of mass and energy as long as the specific stagnation enthalpy within the volume remains invariant. It should be noted that the physics at the section change involves an additional equation that introduces the entropy fluctuations. However it has been checked by solving numerically the full set of LEE that in the present configuration the generated entropy wave is negligible (see section 1.4.1). As pointed out by [START_REF] Davies | Practical flow duct acoustics[END_REF] in the case of an isentropic processes, the jump relations reduce to the following set of equations:

• Mass conservation at x = l:

(1 + M 1 )A + 1 e jk + 1 l -(1 -M 1 )A - 1 e -jk - 1 l S 1 c 1 = [(1 + M 2 )A + 2 e jk + 2 l -(1 -M 2 )A - 2 e -jk - 2 l S 2 c 2 (1.66)
• Energy conservation at x = l:

(1 + M 1 )A + 1 e jk + 1 l + (1 -M 1 )A - 1 e -jk - 1 l 1 ρ 1 = [(1 + M 2 )A + 2 e jk + 2 l + (1 -M 2 )A - 2 e -jk - 2 l 1 ρ 2 (1.67)
Equations (1.64) to (1.67) are recast in the form of a linear system MW = 0, with M:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (1 + M 1 ) -(1 -M 1 ) 0 0 0 0 e jk + 2 L e -jk - 2 L S 1 c 1 (1 + M 1 ) e jk + 1 l -S 1 c 1 (1 -M 1 ) e -jk - 1 l -S 2 c 2 (1 + M 2 ) e jk + 2 l S 2 c 2 (1 -M 2 ) e -jk - 2 l (1 + M 1 ) ρ 2 ρ 1 e jk + 1 l (1 -M 1 ) ρ 2 ρ 1 e -jk - 1 l -(1 + M 2 ) e jk + 2 l -(1 -M 2 ) e -jk - 2 l ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1.68) and W = (A + 1 , A - 1 , A + 2 , A - 2 )
T contains the four unknowns. The dispersion relation is then obtained by requiring the matrix M to be singular, producing the exact solution of the acoustical problem.

Method II : Truncated analytical acoustic model at M ≠ 0

As shown in Fig. 1.2.b the computational domain is reduced to the region between x = l z1 and x = l z2 , where specific boundary impedances Z M ≠0 p,u,1 and Z M ≠0 p,u,2 are imposed at these locations. The Mach number remains finite in the computational domain. Replacing Eq. (1.64) (x = 0) and Eq. (1.64

) (x = L) by impedance conditions p(x = l z1 ) = Z M ≠0 p,u,1 û(x = l z1 ) and p(x = l z2 ) = Z M ≠0
p,u,2 û(x = l z2 ) leads to a dispersion matrix:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (1 -Z M ≠0 p,u,1 )e jk + 1 l z1 (1 + Z M ≠0 p,u,1 )e -jk - 1 l z1 0 0 0 0 (1 -Z M ≠0 p,u,2 )e jk + 2 l z2 (1 + Z M ≠0 p,u,2 )e -jk - 2 l z2 S 1 c 1 (1 + M 1 ) e jk + 1 l -S 1 c 1 (1 -M 1 ) e -jk - 1 l -S 2 c 2 (1 + M 2 ) e jk + 2 l S 2 c 2 (1 -M 2 ) e -jk - 2 l (1 + M 1 ) ρ 2 ρ 1 e jk + 1 l (1 -M 1 ) ρ 2 ρ 1 e -jk - 1 l -(1 + M 2 ) e jk + 2 l -(1 -M 2 ) e -jk - 2 l ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (1.69)
Equations (1.68) and (1.69) should correspond to the same dispersion relation as long as Z M ≠0 p,u,1 and Z M ≠0 p,u,2 properly represent the truncated regions 0 ≤ x ≤ l z1 and l z2 ≤ x ≤ L. A simple translation of the impedance at a constant Mach number leads to:

Z M ≠0 p,u,1 = ζ 1 (M 1 -1) -(M 1 + 1) ζ 1 (M 1 -1) + (M 1 + 1) (1.70) Z M ≠0 p,u,2 = 1 -ζ 2 1 + ζ 2 (1.71)
where

ζ 1 = e 2jk 1 l z1 1-M 2 1 and ζ 2 = e 2jk 2 l z2 1-M 2 2
(1.72)

Method III : Truncated analytical acoustic model at M = 0

Method III is similar to method II in the sense that only the l z1 ≤ x ≤ l z2 domain is solved for. However, the Mach number is now assumed to be zero in this region, so that the Z M =0 approximations (see Table 1.1) can be tested. Methods III pu or III Jm are obtained by replacing

Z M =0 1 and Z M =0
2 by Z M ≠0 p,u,1 and Z M ≠0 p,u,2 , or Z M ≠0 J,m,1 and Z M ≠0 J,m,2 , respectively. The solution of this problem corresponds to the linear system (1.73) which is the system (1.69) where M 1 and M 2 are set to zero.

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (1 -Z M =0 1 )e jk 1 l z1 (1 + Z M =0 1 )e -jk 1 l z1 0 0 0 0 (1 -Z M =0 2 )e jk 2 l z2 (1 + Z M =0 2 )e -jk 2 l z2 S 1 c 1 e jk 1 l -S 1 c 1 e -jk 1 l -S 2 c 2 e jk 2 l S 2 c 2 e -jk 2 l ρ 2 ρ 1 e jk 1 l ρ 2 ρ 1 e -jk 1 l -e jk 2 l -e -jk 2 l ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A + 1 A - 1 A + 2 A - 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = 0 (1.73)
The expression of Z M ≠0 p,u,1 , Z M ≠0 p,u,2 , Z M ≠0 J,m,1 and Z M ≠0 J,m,2 are recalled as follow:

Z M ≠0 p,u,1 = ζ 1 (M 1 -1) -(M 1 + 1) ζ 1 (M 1 -1) + (M 1 + 1) (1.74) Z M ≠0 J,m,1 = ζ 1 + 1 ζ 1 -1 (1.75) Z M ≠0 p,u,2 = 1 -ζ 2 1 + ζ 2 (1.76) Z M ≠0 J,m,2 = - ζ 2 (M 2 -1) + (M 2 + 1) ζ 2 (M 2 -1) -(M 2 + 1) (1.77)
where

ζ 1 = e 2jk 1 l z1 1-M 2 1 and ζ 2 = e 2jk 2 l z2 1-M 2 2
(1.78)

1.4 Results

Results

Thermodynamic and geometric parameters for the computation of the eigenmodes for this test case are presented in Table 1.4. The temperature T 1 and pressure p 1 are imposed, while the values in tube 2 are deduced from the conservation equations of the steady state.

The relations are detailed in Appendix A.

S 1 [m 2 ] S 2 [m 2 ] T 1 [K] p 1 [Pa] γ 0.05 0.1 300 101325 1.4 r [SI] L [m] l [m] l z1 [m] l z2 [m] 287 1 0.5 0.1 0.9
Table 1.4 -Thermodynamic and geometric parameters for the computation of the configuration of figure 1.2.

Effect of the truncation of the domain

The objective of this section is:

• to validate the analytical method I an by comparison with the numerical LEE solver I num ,

• to demonstrate the validity of the computations in the truncated domain, by comparing methods I and II.

Several computations have been performed on the configuration with different inlet Mach numbers M 1 ranging from M 1 = 0.0005 to M 1 = 0.2. Results for the first eigenmode are depicted in Fig. 1.3 for the modes shapes and in Fig. 1.4 for the frequency and growth rate. The numerical LEE solver described at §1. for the first eigenmode with the numerical LEE solver I num ( ) and the reference method I an (×). First two rows: modulus and phase of the fluctuating total enthalpy Ĵ, respectively. Last two rows: modulus and phase of the fluctuating mass flow rate m, respectively. The eigenmode shapes were also verified to be identical although not reported here. This demonstrates that the differences between methods I and III discussed in the next section are due to the zero-Mach number assumption and the modelling of the boundary impedances rather than the truncation procedure.

Results

Effect of the zero-Mach number assumption

This section is devoted to the comparison of the methodologies, presented in Sec. 1.2.2, to prescribe the boundary impedances in method III. The reference is the computation using methods I. The frequency and growth rate of the first two eigenmodes are reported in Fig. 1.5 for various inlet Mach number M 1 .

Results with the zero-Mach number assumption methods III pu and III Jm exhibit the same behaviour for the two computed eigenmodes. The frequency of oscillations departs from the reference solution when M 1 increases. This result was expected since a simple speed argument in a constant cross section area tube shows that the frequency of oscillation at M ≠ 0 is approximatively (1 -M 2 )f 0 where f 0 is the frequency of oscillation under the zero-Mach number assumption. The growth rate is significantly underpredicted with method III pu but very well captured by method III Jm . This last result is coherent with the asymptotic behaviours of the flux mismatch presented in section 1.2.2: the method leading to the smallest mismatch (III Jm ) also leads to the best result in terms of damping rate.

In order to confirm the theoretical assessment of the flux mismatches introduced by methods III pu and III Jm , viz. O(M ) and O(M 2 ) respectively, the flux difference at the boundaries x = l z1 and x = l z2 is plotted in Fig. 1.6 for the two first eigenmodes (top and bottom graphs, respectively). As explained in 1.2.2, the energy flux are made nondimensional with respect to m0 , which is taken at x = l z1 and x = l z2 for the left and right ∆F, respectively. As expected, the difference between the exact energy flux F M ≠0 Eigenmode shapes computed at M 1 = 0.0005 and M 1 = 0.2 are depicted on the left and right columns of Fig. 1.7, respectively. The first two rows present the modulus and phase of the fluctuating total enthalpy Ĵ, respectively, while the last two rows present the modulus and phase of the fluctuating mass flow rate m, respectively. All data are made non-dimensional by setting max m = 1 and arg( m) = 0 where max m is reached. As expected, for the low Mach number case, method III Jm gives the same results as the reference methods I. For a larger Mach number M 1 = 0.2, method III Jm still gives satisfying results, despite small errors in the estimation of the phase shift.

M 1 = 0.0005 M 1 = 0.2 _-_ ʌ ʌ $UJ1-3 ʌ ʌ _P_ $UJsP .ʌ ʌ $[LVg>P@

Results

The present study gives a more precise insight about the use of impedances issued from a non zero mean flow framework. The theory of the energy flux conservation indicates that using a formulation which neglects the mean flow (method III pu ) leads to significant errors, because the error on the imposed energy flux is O(M ). This statement goes against the previous assumption of Lamarque & Poinsot (2008) that prescribed the use of a boundary condition formulated with purely acoustic state variables p and û. As recalled in the introduction of this chapter (see section 1.1), many authors in the aeroacoustic community prescribe the use of a so-called energetical correction of the reflection coefficient which is the (1 -M ) (1 + M ) term in Eq. (1.25). This term originates from the fact that forward and backward waves do not propagate at the same speed and is directly found in the formulations of acoustical reflection coefficients of ducted flow with area-expansion and/or outside radiation conditions [START_REF] Durrieu | Quasisteady aero-acoustic response of orifices[END_REF][START_REF] Boij | Flow effects on the acoustic end correction of a sudden in-duct area expansion[END_REF].

To the authors knowledge, although Method III Jm is equivalent to this well-known correction term, a detailed analysis of its effect when performing a zero-Mach number analysis of thermoacoustic modes has not been provided before. The present study shows that a non-zero-Mach number element can simply be taken into account in a zero-Mach number framework by using method III Jm . It has been shown that this latter method introduces an error of O(M 2 ) and is a good approximation to capture the growth rate and the shapes of eigenmodes at different frequencies. The following chapter is devoted to the modelling of the coupling between the acoustic and entropy waves through the extension of method III Jm .

Chapter 2

Accounting for entropy convection in a zero Mach mean flow In the previous chapter, a methodology has been presented to prescribe the acoustic impedance of the surrounding elements of a combustor by imposing the proper energy flux at the boundary. Thus, the mechanical work due to convective effects is taken into account in the zero-mach Number framework and the growth/decay rate of an instability is better estimated.

Entropy fluctuations accelerated in a mean flow generate acoustic waves (see [START_REF] Morfey | Amplification of aerodynamic noise by convected flow inhomogeneities[END_REF][START_REF] Ffowcs-Williams | The generation of sound by density inhomogeneities in low mach number nozzle flows[END_REF][START_REF] Howe | Contributions to the theory of aerodynamic sound, with application to excess jet noise and the theory of the flute[END_REF][START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF], which is the case when the combustor opens onto a high pressure distributor. This mechanism of indirect combustion noise is important for combustion instabilities, but is not accounted for in a Helmholtz solver since entropy fluctuations do not propagate if the mean flow is at rest. The aim of this present chapter is to propose a methodology to introduce the coupling between entropy and acoustic waves within a zero-Mach number framework. The methodology is decomposed into two parts: first, as the entropy variable does not exist in the computational domain, it must be modelled by a so-called transfer function. Then, the entropy variable is related to the acoustic variable via the imposition of a proper boundary condition, which has been named Delayed Entropy Coupled Boundary Condition (DECBC). The first two sections are devoted to a purely analytical study so as to validate the DECBC methodology. Section 2.3 is then devoted to the extension of DECBC to a numerical 3D Helmholtz solver named AVSP.

Introducing entropy-acoustic coupling in the Helmholtz framework

The configuration depicted in Fig. 2.1 has been chosen in order to exhibit the presence of a set of unstable mixed entropic/acoustic modes. As this kind of modes cannot be predicted by a Helmholtz equation, our objective is to propose a simple model which, used within the Helmholtz framework, allows predicting mixed entropic/acoustic modes. In a combustion chamber of constant cross section and length L, a compact flame (thickness δ f = 0.005 m) is located at x = x f . A choked nozzle is located at x = L while the sonic throat where the Mach number reaches unity is located at x = x th . Following the methodology of section 1.3, this configuration will be computed using three methods:

• Method I is the solution of the numerical LEE solver. The whole domain is computed: the combustion chamber as well as the nozzle.

• Method II is the solution of the LEEs in a restricted region of the domain: only the combustion chamber. The nozzle is modelled by its acoustic impedance and the flame is considered infinitely small. As this method gives rather the same solution as method I and is purely analytical, it is then considered as the reference solution.

• Method III Jm is the solution of the Helmholtz equation (i.e. assuming that M = 0) in the combustion chamber. This solution assesses the influence of the zero-Mach number assumption when

Z M =0 = Z M ≠0
Jm is imposed at the boundaries.

• Method III DECBC Jm is the solution of the Helmholtz equation (i.e. assuming that M = 0) coupled with the Delayed Entropy Coupled Boundary Condition (DECBC) method that introduces the coupling between entropy and acoustic waves at the nozzle. The acoustic boundary conditions are imposed with the formalism

Z M =0 = Z M ≠0
Jm .

Method I : Global thermoacoustic model at M ≠ 0

It is difficult to model the nozzle with the analytical method I an based on waves propagation. It would lead to a discretisation of the nozzle with a collection of several small elements modelling the progressive section change. Although it is theoretically possible, it would lead to a very big linear algebra problem to solve as an eigenvalue. Consequently, the numerical LEE solver I num described in Sec. 1.3.1 is preferred.

Method II : Truncated analytical thermoacoustic model at M ≠ 0

Assuming that the nozzle is properly represented by its acoustic and entropic response at x = L, solving the problem presented in figure 2.1 requires six equations for the six unknown wave amplitudes

(A + 1 , A - 1 , A + 2 , A - 1 , σ 1 and σ 2 ).
It is assumed that the flame is infinitely small, i.e. δ f → 0 and acts as an interface. Moreover the downstream nozzle is modelled with a compact hypothesis. Following the methodology presented by [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF], the following relations can be written:

• Inlet boundary condition at x = 0 : it is assumed that no entropy wave is injected at x = 0 therefore:

σ 1 = 0 (2.1)
Moreover, the mass flow rate is assumed constant, i.e. m1 = 0. Then using Eqs. (1.56) -(1.61) yields:

(1 + M 1 ) A + 1 -(1 -M 1 ) A - 1 = 0 (2.2)
• Mass conservation translates to:

ρ 1 û1 + u 1 ρ1 = ρ 2 û2 + u 2 ρ2 (2.3)
which is recast into:

(1 + M 1 )e jk + 1 x f A + 1 + (M 1 -1)e -jk - 1 x f A - 1 -(M 2 + 1) c 1 c 2 e jk + 2 x f A + 2 -(M 2 -1) c 1 c 2 e -jk - 2 x f A - 2 (2.4) + M 2 c 1 c 2 γp 2 C p σe jksx f = 0 • Momentum conservation yields: p1 + 2ρ 1 u 1 û1 + u 2 1 ρ1 = p2 + 2ρ 2 u 2 û2 + u 2 2 ρ2 (2.5)
which is recast into:

(1 + M 1 ) 2 e jk + 1 x f A + 1 + (M 1 -1) 2 e -jk - 1 x f A - 1 -(M 2 + 1) 2 e jk + 2 x f A + 2 -(M 2 -1) 2 e -jk - 2 x f A - 2 (2.6) + M 2 2 γp 2 C p σe jksx f = 0
• Energy conservation yields:

C p T 01 (ρ 1 û1 + ρ1 u 1 ) + ρ 1 u 1 C p T1 + u 1 û1 = (2.7) C p T 02 (ρ 2 û2 + ρ2 u 2 ) + ρ 2 u 2 C p T2 + u 2 û2
where T 0i = T i + 1 2 u2 i C p is the mean stagnation temperature in section i. In terms of wave amplitudes, the energy conservation takes the form:

M 31 A + 1 + M 32 A - 1 + M 33 A + 2 + M 34 A - 2 + M 3 2 2 γp 2 C p σe jksx f = 0 (2.8)
where

M 31 = c 1 c 2 (1 + M 1 ) M 1 + 1 (γ -1) + M 2 1 2 e jk + 1 x f M 32 = c 1 c 2 (1 -M 1 ) M 1 -1 (γ -1) -M 2 1 2 e -jk - 1 x f M 33 = -(1 + M 2 ) M 2 + 1 (γ -1) + M 2 2 2 e jk + 2 x f M 34 = -(1 -M 2 ) M 2 -1 (γ -1) -M 2 2 2 e -jk -
• Outlet boundary condition at x = L : [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] derived the transfer function of nozzles by assuming the compact assumption, which states that the wavelength of the acoustic and entropy waves are larger than the length of the nozzle. Hence, acoustic and entropy waves propagate quasi-steadily through the nozzle. Matching conditions between the entrance and the exit of the nozzle can be written so that the fluctuations of mass flow rate, total temperature and entropy are considered constant through the nozzle. The compact relation that models a supersonic choked nozzle at x = L then reads:

û2 (L) c 2 - γ -1 2 M 2 p2 (L) γp 2 - 1 2 M 2 ŝ2 (L) C p = 0 (2.9)
Using the wave amplitudes defined in Eqs. (1.56) -(1.61), Eq. (2.9) takes the following form:

(1 -Z M ≠0 p,u,2 )A + 2 e jk + 2 L + (1 + Z M ≠0 p,u,2 )A - 2 e -jk - 2 L -R sa 1 + Z M ≠0 p,u,2 e jksL γp 2 C p σe jksL = 0 (2.10)
where Z M ≠0 p,u,2 is the impedance of the nozzle submitted to an incident acoustic wave:

Z M ≠0 p,u,2 = 2 (γ -1)M 2 (2.11)
and R sa is its acoustic response when submitted to a purely entropic perturbation:

R sa = - M 2 2 + (γ -1)M 2
(2.12)

Equations (2.1) to (2.10) may be recast in the form of a linear system MW = 0, where M is the matrix:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ (1 + M 1 )e jk + 1 x f (M 1 -1)e -jk - 1 x f -(M 2 + 1) c 1 c 2 e jk + 2 x f -(M 2 -1) c 1 c 2 e -jk - 2 x f M 2 c 1 c 2 (1 + M 1 ) 2 e jk + 1 x f (M 1 -1) 2 e -jk - 1 x f -(M 2 + 1) 2 e jk + 2 x f -(M 2 -1) 2 e -jk - 2 x f M 2 2 M 31 M 32 M 33 M 34 M 3 2 2 (1 + M 1 ) -(1 -M 1 ) 0 0 0 0 0 M 53 M 54 M 55 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (2.13) and W = (A + 1 , A - 1 , A + 2 , A - 2 , Σ) T , with Σ = γp 2 σ C p e jksx f
and with :

M 53 = (1 -Z M ≠0 p,u,2 )e jk + 2 L
(2.14)

M 54 = (1 + Z M ≠0 p,u,2 )e -jk - 2 L
(2.15)

M 55 = -R sa 1 + Z M ≠0 p,u,2 e jks(L-x f ) (2.16)
The corresponding dispersion relation is obtained by zeroing det(M).

Method III Jm : Truncated analytical thermoacoustic model at M = 0

Assuming that M = 0 imposes that there is no entropic wave, i.e. σ 1 = σ 2 = 0. The 5-by-5 linear system in Eq. (2.13) reduces to a 4-by-4 system:

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ e jk 1 x f ρ 1 c 1 - e -jk 1 x f ρ 1 c 1 - e jk 2 x f ρ 2 c 2 e -jk 2 x f ρ 2 c 2 e jk 1 x f -e -jk 1 x f -e jk 2 x f -e -jk 2 x f 1 -1 0 0 0 0 (1 -Z M ≠0 J,m,2 ) e jk 2 L (1 + Z M ≠0 J,m,2 ) e -jk 2 L ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A + 1 A - 1 A + 2 A - 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = 0 (2.17)
Only method III Jm is used in this section given its better behaviour compared to method III pu (see section 1.4.2). Thus, at x = 0 the boundary condition m1 = 0 leads to Z M ≠0 J,m,1 = ∞ (see Table 1.1), which may be written into the relation

A + 1 -A - 1 = 0. At z = L the impedance Z M ≠0
J,m,2 may be expressed with Eq. (1.37) and Eq. (2.11), leading to:

Z M ≠0 J,m,2 = 1 M 2 1 + (γ -1)M 2 2 2 1 + (γ -1) 2 (2.18)
Of course since the conversion of the entropy wave into a reflected acoustic wave at x = L has been suppressed from the equations, the system (2.17) cannot capture modes where this mechanism is important. The next section shows how a simple model (called DECBC) can be used to re-introduce this mechanism in a zero-Mach formulation.

Method III DECBC

Jm

: Delayed Entropy Coupled Boundary Condition

As the Helmholtz framework relies on the zero-Mach number assumption, the effect of elements, such as the nozzle, where a mean flow occurs can simply be taken into account with method III Jm by imposing the right energy flux through the proper impedance formulation (see Sec. 1.2.2). However imposing only the acoustic impedance, Eq. (2.18), means neglecting the entropy-acoustics coupling and the subsequent sound wave being reflected by the nozzle when entropy spots flow through the exit. Indeed, ŝ2 (L) is not available in the Helmholtz domain. Thus the entropy fluctuations at the edge of the combustion chamber must be modelled before Eq. (2.10) can be applied. This is the purpose of the proposed DECBC (Delayed Entropy Coupled Boundary Condition) methodology described in the following of this section. Assuming that the entropy fluctuations entering the nozzle at x = L have first been generated in the flame region before being convected by the mean flow (see Fig. 2.2), one obtains:

ŝ2 (L) = G c e jωτc ŝ2 (x f ) (2.19)
where ŝ2 (x f ) is the amount of entropy generated by the flame, τ c is the convection time from the flame to the exit and G c a gain (smaller than unity) related to the dissipation of entropy waves during their convection, if any. Eq. (2. 19) is useful only if the generated entropy ŝ2 (x f ) can be related to some acoustic quantity. It is done here in a way similar to the n-τ model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF] for unsteady heat release, connecting the entropy fluctuation to the acoustic velocity taken at a reference point x ref located upstream of the flame region, here denoted by the subscript 1 :

ŝ2 (x f ) = G us e jωτus û1 (x ref ) (2.20)
where G us and τ us are the gain and the time delay of the entropy generation from an acoustic perturbation at a reference location x ref . The two transfer functions described at Eqs. (2. 19) and (2.20) are schematically depicted in Fig. 2.2.

In the present study, the quantities G us and τ us can be assessed analytically in the simple case of a 1D premixed flame [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF]. One obtains the following expression for the entropy produced by the flame at x = x f :

G us = - C p 2 T 02 -T 01 (γ -1)
u 2 c 2 2 ρ 1 ρ 2 ; τ us = 0 (2.21)
where T 0 = T + 1 2 u 2 C p is the mean stagnation temperature. Note that τ us = 0 in this case because the referential velocity û(x ref ) is considered at the flame location x ref = x f . Entropy fluctuations at the exit x = L may then be analytically expressed with the help of Eq. (2.19):

ŝ2 (L) = G us G c e jω(τc) û1 (x f ) (2.22)
where the coefficients related to the convection part are:

G c = 1 ; τ c = L -x f u 2 (2.23)
Note that no dissipation is taken into account in the present study. Injecting Eq. (2.22) into Eq. (2.10), and expressing û1 (x f ) in terms of waves A + 1 and A - 1 allows to write the following condition at x = L:

A + 1 e jk 1 x f βR sa (1 + Z M =0 2 ) -A - 1 e -jk 1 x f βR sa (1 + Z M =0 2 ) + A + 2 e jk 2 L 1 -Z M =0 2 + A + 2 1 + Z M =0 2 e -jk 2 L = 0 (2.24) with β = 1 2 C p (T 02 -T 01 )(γ -1)e jωτc u 2 c 1 (2.25)
The contribution of the DECBC method appears in the first two terms of Eq. (2.24), by coupling acoustic waves at the boundary condition (x = L) with the acoustic waves upstream of the flame. Finally, the 5×5 problem defined in Eq. (2.13) reduces to Eq. (2.26):

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ e jk 1 x f ρ 1 c 1 - e -jk 1 x f ρ 1 c 1 - e jk 2 x f ρ 2 c 2 e -jk 2 x f ρ 2 c 2 e jk 1 x f -e -jk 1 x f -e jk 2 x f -e -jk 2 x f 1 -1 0 0 M 41 M 42 M 43 M 44 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ A + 1 A - 1 A + 2 A - 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = 0 (2.26) with M 41 = βR sa 1 + Z M ≠0 J,m,2 e jk 1 x f (2.27) M 42 = -βR sa 1 + Z M ≠0 J,m,2 e -jk 1 x f (2.28) M 43 = 1 -Z M ≠0 J,m,2 e jk 2 L
(2.29)

M 44 = 1 + Z M ≠0 J,m,2 e -jk 2 L (2.30)

Results

Results

Thermodynamic parameters used in the computation of eigenmodes of the configuration depicted in Fig. 2.1 are presented in Table 2.1. For computations with the method III DECBC

Jm

, the parameters of the DECBC model are recalled in Table 2.2 for different inlet Mach number M 1 .

T 1 [K] p 1 [Pa] T 2 [K] γ r [SI] L [m] x f [m] x th [m] δ f [m] 300

Validations and effect of the truncation

Several computations have been performed on the configuration of Fig. 2.1 with different inlet Mach number M 1 ranging from M 1 = 0.00625 to M 1 = 0.05. Results for the first two eigenmodes are depicted in Fig. 2.3. The numerical LEE solver (method I num ) described at section 1.3.1 as well as the global analytical model (method II, see section 2.1.2) give very similar results in terms of frequency. The evaluation of the growth rate exhibits very small discrepancies, which can be attributed to the fact that in the numerical solver the flame is not infinitely thin and has a thickness δ f of 0.005 meters. Although not shown, the mode shapes are also very close. It confirms the robustness of the analytical method II and validates both the compactness assumption of the flame and the nozzle.

Frequency analysis

Results are presented in Fig. 2.4 for various inlet Mach numbers M 1 . The upper left graph presents the frequency of the first three eigenmodes predicted by the reference method II and zero-Mach number methods III Jm and III DECBC

Jm

. The associated growth rates are depicted on the upper right, lower left and lower right graphs, respectively.

The behaviour of the first eigenmode shows a strong dependency on the Mach number. The frequency is proportional to the convective time delay of the flow, suggesting that this first eigenmode is coupled with entropic waves. This eigenmode is not predicted by method III Jm (as expected). The second and third eigenmodes do not present the same behaviour and their frequencies do not follow a linear behaviour with respect to M 1 . At very small Mach number, the reference solution II and Helmholtz computations with method III Jm give the same estimation of the frequencies. It suggests that the second and third eigenmodes are mostly acoustic, but influenced by their interactions with the mean flow and the coupling with entropy waves at the exit nozzle. Note that when approaching Mach number M 1 = 0.05, the first two eigenmodes tend to merge into a unique unstable mode. Such behaviour has already been observed by Goh & Morgans (2013) who computed several academic combustors while taking into account the dissipation and dispersion of the entropy wave. It was shown that the entropy wave can destabilise a stable combustor and vice-versa, and that some configurations can exhibit a mode switch between unstable eigenmodes.

As the first eigenmode is purely dependent on the convective time delay, it can obviously not be captured by the Helmholtz solver when using method III Jm . In the same way, the two other eigenmodes computed with III Jm are invariant in frequency whichever the Mach number prescribed, because no convective effect is taken into account. The relative growth rates of the eigenmodes are also very close to zero, suggesting that the stability map of this kind of configuration cannot be estimated with a simple Helmholtz solver but that the DECBC extended model can. It is remarkable that the zero-Mach number method III DECBC Jm is able to capture the frequency of each eigenmode. By comparison with the reference method II, the growth rate is also correctly estimated, although slightly overestimated for M 1 > 0.025.

Mode shapes analysis

Eigenmode shapes computed at M 1 = 0.00625 and M 1 = 0.05 are depicted on left and right columns of Fig. 2.5 for the second eigenmode. The first two rows present the modulus and phase of the fluctuating total enthalpy Ĵ, respectively, while the last two rows present the modulus and phase of the fluctuating mass flow rate m, respectively. All data are made non-dimensional by setting max m = 1 and arg( m) = 0 where max m is reached.

As expected, the zero-Mach number method III Jm where entropy waves are not taken into account give good results at very low Mach number. At larger Mach numbers, the modes shapes are mispredicted because the frequency is not well predicted and corresponds to the zero-Mach number solution. However, modes shapes computed with method III DECBC Jm are very close to the reference method II. The enthalpy and mass-flow rate jumps across the flame at x = x f are well captured. The present results show that the use of the DECBC approach is able to capture mixed-modes based on coupling between entropy and acoustics, and to predict them correctly, despite the fact that the mean flow is neglected within the computational domain. ). First two rows: modulus and phase of the fluctuating total enthalpy Ĵ, respectively. Last two rows: modulus and phase of the fluctuating mass flow rate m, respectively.

Results

Extension of DECBC to 3D

The purpose of this section is to describe the implementation of the DECBC methodology within a 3D numerical Helmholtz solver. The derivation of the equations which leads to the inhomogeneous Helmholtz equation is presented first and completed by the reformulation into an eigenvalue problem discretised and solved by the code. The implementation of the boundary condition is then discussed.

The zero-Mach number Helmholtz equation

The starting point for the derivation of the Helmholtz equation is the set of linearised equations of conservation for a fluid (see Eqs. (1.6) -(1.8)). The Helmholtz equation relies on the zero-Mach number assumption, which means that u = 0. Application of this assumption to the momentum equation for the mean flow Eq. (1.10) also implies that ∇p = 0. The set of linearised equations then reduces to:

∂ρ ′ ∂t + u ′ ∇ρ + ρ∇ ⋅ u ′ = 0 (2.31) ρ ∂u ′ ∂t + ∇p ′ = 0 (2.32) ∂s ′ ∂t + u ′ ∇s = rq ′ p (2.33)
The next step is to express a link between entropy and the pressure and the density. First, the variation of entropy may be expressed with the following relation [START_REF] Shapiro | The Dynamics and Thermodynamics of Compressible Fluid Flow[END_REF] As the mean flow quantities are independent of time, Eq. (2.36) reduces to the following expression for the mean entropy gradient:

∇s = C v p ∇p - C p ρ ∇ρ (2.37)
As explained above the zero-Mach number assumption implies that ∇p = 0, so that Eq. (2.37) reduces to ∇s = -Cp ρ ∇ρ. Introducing the expression of ∇s in the linearised energy equation (2.33) and using the linearised state and entropy equations (1.54) to eliminate ρ ′ , the set of equations (2.31) -(2.33) can be simplified to:

1 γp ∂p ′ ∂t + ∇ ⋅ u ′ = rq ′ C p p (2.38) ∂u ′ ∂t + 1 ρ ∇p ′ = 0 (2.39)
Finally, taking the time derivative of Eq. (2.38) and adding the divergence of Eq. (2.39) allows to write the following inhogeneous wave equation:

1 γp ∂ 2 p ′ ∂t 2 -∇ ⋅ 1 ρ ∇p ′ = γ -1 γp ∂q ′ ∂t (2.40)
As recalled in section 1.2.2, this study is conducted in the frequency domain: a fluctuating quantity, g ′ (x, t), is defined by its complex amplitude, ĝ, and angular frequency, ω, through g ′ (x, t) = R ĝ(x)e -jωt . This assumption implies that a time derivative becomes a simple multiplication by -jω. Eq. (2.40) can then be reformulated in the frequency domain as the following Helmholtz equation:

ρ (x) c 2 (x) ∇ ⋅ 1 ρ(x) ∇p(x) + ω 2 p(x) = jω (γ(x) -1) q(x) (2.41)
where the sound speed c, the mass density ρ and the heat capacity ratio γ are local variables.

In order to close the problem, the flame that constitutes a forcing source is often modelled as a purely acoustic element thanks to a nτ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF]:

q(x) = n local (x)û(x ref )e jωτ local (x) ⋅ n ref (2.42)
Eq. (2.42) essentially relates the unsteady heat release q to acoustic quantities û at a reference location denoted x ref by an amplification factor n local (unit [J.m -4 ]) and a time delay τ local (unit [s]). The n local -τ local parameters can be estimated experimentally, by analytical models or directly extracted from the LES database by post-processing the heat release rate and velocity signals of an oscillating flame (forced or self excited). Appendix B presents more details about the active flame model and especially about the method employed to compute the nτ fields from the LES. Note that the Helmholtz equation only involves the acoustic pressure p and that the velocity vector û can be deduced from p by using the momentum equation: jωρû = ∇p. Finally Eq. (2.41) solved by the Helmholtz solver is rewritten as follows:

ρ (x) c 2 (x) ∇ ⋅ 1 ρ(x) ∇p(x) + ω 2 p(x) = (γ(x) -1) n local (x)e jωτ local (x) ∇p(x ref ) ⋅ n ref ρ(x ref ) (2.43)

Boundary conditions and DECBC implementation

In this section, the boundary condition that couples acoustic and entropy waves is derived first in a general manner with mean flow effects included. Then the transition to the zero-Mach number assumption is made to match the Helmholtz solver formalism. As far as purely acoustic aspects have been treated previously in chapter 1 and do not bring new difficulties, the boundary condition implemented in the solver must be reformulated with a new term that introduces the entropy contribution.

The acoustic impedance is imposed to the computational domain as a uniform boundary condition over a discretised surface. In this context the boundary condition relies on a quasi-1D formulation. Hence, the fluctuating pressure and velocity fields can be reformulated in a general manner with the following 1D relations [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF]:

p(ξ) γp = A out e jk out ξ + A in e -jk in ξ (2.44) û(ξ) ⋅ n c = -A out e jk out ξ -A in e -jk in ξ (2.45)
where k in = ω c(1-M ) and k out = ω c(1+M ) are the acoustic wave number of the acoustic waves A in and A out propagating along on the curvilinear axis ξ (see Fig. 2

.6).

A in is the wave entering into the computational domain while A out is the wave leaving the domain at c(1 + M ). Note that the minus sign in front of the RHS of Eq. (2.45) stems from the fact that the unit vector normal to the exit boundary is colinear with the ξ-axis, pointing toward the negative ξ. Following [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] an acoustic wave is generated by the entropy spots flowing through the exit nozzle or turbine. At the location ξ = 0 where the boundary condition is imposed, the wave A in represents the acoustics generated by the nozzle and entering the computational domain. It can be assumed linear in the two outgoing waves:

A in = R p,u A out + R p,s σ out (2.46)
where R p,u and R p,s are the reflection coefficients for an impacting acoustic and entropy wave, respectively. The term σ is the wave of the entropy fluctuations ŝ(ξ) C p = σ out e jksξ that is leaving the domain, C p is the heat capacity at constant pressure and k s = ω (cM ) the wave number associated to the convection of entropy spots by the mean flow. At ξ = 0, Eqs. (2.44) and (2.45) become:

p(ξ) γp = A out + A in (2.47) û(ξ) ⋅ n c = -A out + A in (2.48) while ŝ C p = σ out (2.49)
Eq. (2.46) may then be recast into the following relation:

p(R p,u -1) -ρ c û ⋅ n(R p,u + 1) + 2 γp C p R p,s ŝ = 0 (2.50)
with R p,u = A in A out and R p,s = A in σ out . Note that all the parameters γ, p, ρ, c and C p must be taken at the location where the boundary condition represented by Eq. (2.50) is imposed. After introducing the relation between the acoustic impedance and the acoustic reflection coefficient Z M =0 p,u = (R p,u + 1) (R p,u -1), the implementation of Eq. (2.50) into the Helmholtz solver AVSP leads finally to the following reformulation of the boundary condition (2.50): p

Z M =0 p,u -ρ c û ⋅ n + Z M =0 p,u -1 Z M =0 p,u γp C p R p,s ŝ = 0 (2.51)
As the Helmholtz solver relies on the zero-Mach number assumption, the effect of elements, such as the nozzle, where a mean flow occurs can simply be taken into account by imposing the right energy flux through the proper acoustic impedance formulation (see Sec. 1.2.2). Under the zero-Mach number assumption, the general boundary condition described at Eq. (2.51) takes then the following formulation:

p Z M =0 p,u -ρ c û ⋅ n = 0 (2.52)
The use of the momentum equation (jωρû ⋅ n = ∇p ⋅ n) allows to write Eq. (2.52) only as a function of the pressure variable as follow:

∇p (x) ⋅ n = jω c (x) Z M =0 p,u p (x)
(2.53) where (x) represents here the vector of nodes located at the boundary where the impedance is imposed.

However imposing only the acoustic impedance, Eq. (2.52), means neglecting the entropy-acoustics coupling and the subsequent sound being generated by the entropy spots flowing through the exit nozzle or turbine. Now, ŝ in Eq. (2.53) is not available in the Helmholtz domain and the entropy fluctuation at the edge of the combustion chamber must be modelled before Eq. (2.51) can be applied. This is done here by applying the DECBC modelling formalism described in Sec. 2.1.4. The entropy fluctuations at the boundary condition takes the following expression:

ŝ = G us G c e jω(τus+τc) û (x ref ) ⋅ n ref (2.
∇p (x) ⋅ n = jω c (x) Z M =0 p,u p (x) + γp C p R p,s Z M =0 p,u -1 Z M =0 p,u G us G c e jω(τus+τc) ∇p(x ref ) ⋅ n ref c (x) ρ ref (2.55)

Numerical resolution strategy

Eq. (2.41) corresponds to a non-linear eigenvalue problem which can be solved numerically by using appropriate algorithms [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF][START_REF] Salas | Aspects numériques et physiques des instabilités de combustion dans les chambres de combustion annulaires[END_REF]. First the Eq. (2.43) is discretised using a finite element method and turned into the following matrix form:

Ap + B(ω)p + ω 2 p = C(ω)p (2.56)
where the matrix A, B and C are defined as follows:

A = ρ (x) c 2 (x) ∇ ⋅ 1 ρ (x) ∇ B = c (x) ⎛ ⎝ jω Z M =0 p,u + γp C p R p,s Z M =0 p,u -1 Z M =0 p,u G us G c e jω(τus+τc) ∇ ref ρ(x ref ) ⎞ ⎠ (2.57) C = (γ (x) -1) n local (x) e jωτ local (x) ∇ ref ρ (x ref )
The discretisation procedures of the operator ∇ ⋅ 1 ρ(x) ∇ contained in the matrix A are described in Appendix C. B is the matrix containing the impedances at boundary points, and C the active flame operator. In both B and C the term ∇ ref should be understood as the gradient of the pressure taken at the reference point and along the reference direction n ref . Eq. 2.56 is a linear eigenvalue problem if B = 0 and C = 0 and can be solved via the Arnoldi algorithm. However this is not the case in the present study as complex boundary impedances as well as an active flame are accounted for. Eq. (2.56) is then linearised in ω 2 with the following iterative procedure: where ω n is known and ω * n is solution of the problem:

ω n+1 = αω * n + (1 -α)ω n (2.58)
[A + ω n B + C(ω n )] p + (ω * n ) 2 p = 0 (2.59)
Eq. (2.58) allows to obtain ω n+1 for the iteration n + 1. A relaxation coefficient α may be used to smooth the process in case of convergence difficulties. The algorithm is considered converged when = ω n -ω * n < 0.01. The present algorithm is known as the fixed-point method and is represented schematically in Fig. 2.7. Note that if the impedance Z M =0 p,u is not constant, it should be computed at each iteration to provide a value for the complex excitation frequency ω n . Finally, the real part of ω n is the eigenfrequency (unit [Hz]) while the imaginary part represents the damping/growth rate of the eigenmode (unit [s -1 ]), depending on its negative/positive sign, respectively.

Validations

The implementation of the Delayed Entropy Coupled Boundary Condition in the numerical Helmholtz solver AVSP is validated by computing the eigenmodes of the configuration described in Sec. 2.1 and depicted in Fig. 2.1. The results from AVSP/DECBC are compared to the analytical reference method II (see Sec. 2.1.2) and the zero-Mach number analytical method III DECBC Jm . For both methods using DECBC, the gain G and time delay τ of the transfer functions for the generation and convection of entropy fluctuations (see Eqs. (2.21) and (2.23)) are the same (see Table 2.2). Similarly, the imposed impedance and acoustic/entropy reflection coefficient are also the same (see Eqs. (2.18) and (2.12)). ).

Results are depicted in Fig. 2.8. All methods give exactly the same frequency. The growth rates still differ between the reference method and the zero-Mach number computations using DECBC. However computations with AVSP/DECBC are slightly closer to the reference method than the analytical method III DECBC Jm . Such discrepancies may be attributed to the fact that the location of the reference point x ref for the DECBC model is not exactly the same between the analytical and the numerical methods. Indeed, in the numerical solver the location of the reference point depends on the distribution of the mesh nodes. In the present case, x ref is actually located 0.001 m upstream of the flame.

Chapter 3

Investigation of a realistic gas turbine A combustion instability in a combustor terminated by a nozzle is analysed and modelled thanks to a low order Helmholtz solver. A Large Eddy Simulation (LES) of the corresponding turbulent, compressible and reacting flow is first performed and analysed thanks to Dynamic Mode Decomposition (DMD). The mode with the highest amplitude shares the same frequency of oscillation as the experiment (approx. 320 Hz) and it shows the presence of large entropy spots generated within the combustion chamber and convected down to the exit nozzle. The lowest thermo-acoustic mode being in the range 700 -750 Hz, it is postulated that the instability observed around 320 Hz stems from a mixed entropy/acoustic mode where the acoustic generation associated with entropy spots being convected throughout the choked nozzle plays a key role.

The objectives of this chapter are then as follows:

• build and analyse a large-eddy simulation data base to establish the presence of mixed modes in a realistic combustion chamber;

• develop and validate a methodology to estimate the entropy/acoustic transfer functions so as to feed the DECBC model in a Helmholtz solver;

• apply the obtained low-order method to the 3D industrial-like combustor.

First, the industrial configuration, the numerical set-up and the associated lowfrequency instability are presented in section 3.1, followed by preliminary results in section 3.2. In section 3.3 a Dynamic Mode Decomposition (DMD) is applied to the LES results in order to investigate the presence of an entropy-acoustic coupling via the reconstruction of a reduced-order model. An acoustic study is performed with the help of a zero-Mach number Helmholtz solver (see section 3.4.2), demonstrating that a classical thermo-acoustic approach cannot represent the observed low-frequency mode, even though the coupling between the heat release fluctuations from the flame and the acoustics is taken into account. In section 3.4.4 results from the zero-Mach number Helmholtz solver with the DECBC approach are presented to illustrate the benefit of the method. Finally, in section 3.5 a sensitivity study is performed to optimise the shape of the downstream nozzle. It is shown that for certain geometry of the nozzle, the growth rate of the unstable mode can be reduced. LES computations with an extended geometry of the nozzle confirm that the combustion instability can be significantly attenuated with the proposed methodology.

Experimental configuration and numerical set-up

Description of the experimental configuration

The case considered in this study is a combustor developed by the SAFRAN Group for aero-engine applications. Figure 3.1 shows the combustion chamber and the casing with primary and dilution holes although some parts of the geometry are not displayed for confidentiality reasons. The air inlet connected to the upstream compressor is also displayed. Note that the downstream high pressure distributor which connects the combustor to the turbine in the actual experiment is not represented here nor directly accounted for in the remaining of the study. Instead it is modelled by a choked nozzle with equivalent cross section area. For choked conditions, this ensures that the chamber sees an outlet condition close to the actual one. The whole mass-flow rate is not entering into the nozzle: approx. 16% is routed downstream of the surrounding casing (modelled in the simulation by suction walls). The fuel line is also visible, as well as the swirled injector used to mix fuel and air and to generate the recirculation zone which stabilises the flame. In the engine, several injectors are mounted all around the azimuthal combustion chamber although only one sector with one injector is displayed in the figure. Note that this geometry is representative of a typical aero-engine gas turbine combustor, as studied in other works (see [START_REF] Dowling | Acoustic Analysis of Gas Turbine Combustors[END_REF].

Experiments were performed by SAFRAN on the present configuration in order to study the effect of liquid water ingestion. Depending on the ratio between the mass flow rates of water and air, the overall stability of the engine evolves, the combustor becoming increasingly unstable as the mass flow rate of water increases. As the temperature of the ingested liquid water is approx. 300 K, it is postulated that all of the liquid water is vaporised when flowing through the compressor and that the main effect of its ingestion is to decrease the temperature of the air feeding the combustor. Although not detailed in the present thesis, several LES computations have been performed at SAFRAN for different inlet boundary condition:

• one which takes into account the mass fraction of vaporised water in the gas entering into the combustor,

• one where only pure air enters the combustion chamber.

No significant differences were observed on the flame dynamics or the frequency and amplitude of the instability. In consequence, the effect of liquid water is taken into account in the present study by only decreasing the air entering into the combustor. It is possible to link the inlet temperature T inlet to the mass flow rates of liquid water and air by a simple enthalpy balance. The engine is found stable for T inlet approx. 760 K but unstable for T inlet approx. 560 K. In unstable cases the measured pressure signal exhibit strong fluctuations at approximately 330 Hz and with a maximum amplitude of approx. 5 × 10 4 Pa.

Numerical set-up

To analyse this instability, large eddy simulations (LES) were performed using the AVBP code which embeds a cell-vertex formulation and a set of finite element / finite volume schemes for unstructured meshes. A two-step Taylor-Galerkin scheme (3 nd order in time and space) was retained [START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF]. Two regimes were computed by LES, corresponding to two representative operating conditions investigated experimentally at SAFRAN: one which contains a strong instability at nearly 330 Hz (T inlet = 557K; OP-A) and one which shows no instability (T inlet = 762K; OP-B). In each case, only one angular sector (out of 20) was computed, axiperiodic conditions being prescribed at the lateral boundary. It is consistent with the experimental data showing that the low frequency instability is not an azimuthal mode.

The main idea of the LES is to resolve the large scale eddies in the flow which contain most of the energy, while the small scales are represented by subgrid models. For any given quantity g, the filtered component g is resolved in the numerical simulation while the subgrid scales g ′ = gg are modelled. Hence, the filtering operation is directly related to the grid mesh resolution. As the density is variable, a mass-weighted Favre filtering is introduced according to ρg = ρg. One difficulty raised by LES is that the spatial filtering of the equations leads to unresolved terms such as the Reynolds stresses or the heat release and reaction rates. These terms need to be modelled. Another key issue is to represent the flame front properly. Several approaches are available in the literature to handle this problem and the present study uses the Dynamic Thickened Flame Model (Colin et al., 2000). Basically the idea is to artificially thicken the flame by a factor F. However in turbulent flows, the interaction between turbulence and chemistry is altered and the thickening of the flame reduces the ability of the vortices to wrinkle the flame front. In order to conserve the flame speed and to not underestimate the reaction rate, the diffusivity and the reaction rate is corrected by an efficiency function E. In the dynamic formulation of the thickened flame model, F is not constant and is driven by a sensor S depending on the local temperature and mass fractions. This sensor locates the flame zone so as to not apply the factor F in non reactive zones. The efficiency function describes effects of the sub-grid scale wrinkling as a function of the local sub-grid turbulence velocity u ′ and spatial filter size ∆. More details can be found in the literature (Poinsot & Veynante, 2011). In summary, the LES computed in the present thesis solve the following set of filtered modelled equations:

• Momentum ∂ρũ i ∂t + ∂ ∂x j (ρũ i ũj ) = - ∂ ∂x j pδ ij -2 (µ + µ t ) Sij - 1 3 Sll δ ij (3.1)
where ρ is the density, u the velocity, p the pressure, µ the dynamic viscosity of the gas, µ t the turbulent viscosity and S ij the strain rate tensor:

S ij = 1 2 ∂ũ i ∂x j + ∂ũ j ∂x i (3.2)
Subscripts i,j and l refer to the three spatial dimensions while µ t is the subgrid scale turbulent viscosity which must be modelled.

• Energy

∂ρ Ẽ ∂t + ∂ ∂x j ρ Ẽũ j = - ∂ ∂x j ũi pδ ij -2µũ i Sij - 1 3 Sll δ ij + C p EF µ P r + (1 -S) µ t P t r ∂ T ∂x j + ∂ ∂x j ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ N k=1 ⎛ ⎝ ⎡ ⎢ ⎢ ⎢ ⎣ EF µ S c,k + (1 -S) µ t S t c,k ⎤ ⎥ ⎥ ⎥ ⎦ W k W ∂ Xk ∂x j -ρ Ỹk Ṽ c j + Ṽ c,t j ⎞ ⎠ hs,k ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + E ωT F (3.3)
where E is the total energy per mass unit, T the temperature, C p the specific heat capacity, W is the mean molecular weight of the mixture while X and Y are the species molar and mass fractions, respectively. E is the efficiency function, F is the flame thickening factor and S the sensor. The vectors Ṽ c j + Ṽ c,t j are the corrected diffusion velocities, h s,k correspond to the sensible enthalpies per mass unit and finally ωT is the reaction rate modelled by an Arrhenius law.

• Species ∂ρ Ỹk ∂t + ∂ ∂x j ρ Ỹk ũj = ∂ ∂x j ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ ⎛ ⎝ EF µ S c,k + (1 -S) µ t S t c,k ⎞ ⎠ W k W ∂ Xk ∂x j -ρ Ỹk Ṽ c j + Ṽ c,t j ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ + E ωk F (3.4)
where ωk are the reaction rate of each species computed from the Arrhenius law.

A simple two-step kinetic scheme was used to represent kerosene-air chemistry [START_REF] Franzelli | A two-step chemical scheme for Large-Eddy Simulation of kerosene-air flames[END_REF]. The global power is approx. 1.4 × 10 6 W while the global equivalence ratio is approx. 0.3, which is characteristic of lean combustors.

Boundary conditions are imposed through the NSCBC formulation [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] to control spurious acoustic reflections. The solid walls are adiabatic and a classical logarithmic law is imposed to mimic the near-wall shear stress [START_REF] Schmitt | Large-eddy simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high-pressure burner[END_REF]. The coincident walls between the combustion chamber and the casing are actually multiperforated and modelled in the simulation by an adiabatic homogeneous suction/injection boundary condition (Mendez & Nicoud, 2008a,b). Only the hydrodynamics effects of the multiperforated walls are taken into account by the model, not the associated acoustic damping. For each LES computations data are recorded upon a period of 120 ms after the whole system reaches a well established limit cycle state. As the expected instability phenomenon has a frequency of approx. 330 Hz, this latter choice ensures that the averaging procedure of physical fields and signal processing like Fourier transforms cover at least 40 cycles. As mentioned in Sec. 3.1.1, the turbine is replaced by a choked nozzle with equivalent cross section area. The natural question that arises is whether the response of the nozzle to entropic (temperature) perturbations is sensitive to the three-dimensional nature of the temperature hot-spots. Two nozzle geometries were tested so as to ensure the robustness of the results to the shape of the nozzle: a flat one and a circular one, depicted in the top and bottom of Fig. 3.1, respectively. Two meshes were also tested for the geometry coupled with a circular nozzle: a coarse one (approx. 5.6 millions cells) and a fine one (approx. 15 millions cells). Note that in both meshes the veins inside the swirler as well as the dilution holes are refined to capture the flow details. An axial cut is depicted on the top and bottom of Fig. 3.2 for the coarse and fine meshes, respectively. Moreover, two subgrid scale models for the turbulent viscosity µ t in Eqs. (3.1) to (3.4) have been tested:

• the classical Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF]

µ t = (C S ∆) 2 2 Sij Sij (3.5)
where ∆ denotes the filter characteristic length and C S = 1.8 is a constant. The turbulent viscosity µ t is then proportional to the resolved filtered strain tensor Sij . The Smagorinsky model was developed in the 1960's and is extensively tested over a wide range of flow configurations. This model is easy to implement and presents good properties for flows away from solid walls at a low computational cost. Furthermore, it supplies the right amount of kinetic energy dissipation in homogeneous isotropic turbulent flows. The Smagorinsky model is also known as being too dissipative in wall bounded flows and to over predicts the fluid friction near solid boundaries. This problem can be avoided by the use of wall-law boundary conditions.

• the Sigma model (Nicoud et al., 2011)

µ t = (C σ ∆) 2 σ 3 (σ 1 -σ 2 )(σ 2 -σ 3 ) σ 2 1 (3.6)
where C σ = 1.5 is a model constant, σ 1 ≥ σ 2 ≥ σ 3 ≥ 0 are the singular values of the resolved velocity gradient tensor. The σ model was developed to overcome some drawbacks observed on most of the static models without using additional test filter. Indeed, it presents the interesting properties to vanish in various laminar flow configurations (pure shear, pure rotation, pure dilatation/contraction, ...) for which no subgrid scale activity is expected. It also has the property to vanish with the proper asymptotic behaviour near solid boundaries.

Table 3.1 presents a summary of the LES computations performed on the two operating points depending on the nozzle geometry and the mesh refinement. In order to discriminate either the effect of the nozzle or the inlet temperature, all numerical and physical parameters in the computations are strictly the same.

Mesh Nozzle

Run Operating Point Baseline Fine Circular Flat Regime

1 OP-A × × Unstable 2 OP-B × × Stable 3 OP-A × × Unstable 4 OP-B × × Stable 5 OP-A × × Unstable 6 OP-B × × Stable
Table 3.1 -Summary of the LES computations performed on the two operating points OP-A (T inlet = 557 K) and OP-B (T inlet = 762 K) depending on the nozzle geometry (circular or flat) and the mesh (coarse mesh approx. 5.6 millions cells, fine mesh approx. 15 millions cells). The last column denotes if the computed flow was found stable or unstable. Recall that the experiment is unstable for OP-A and stable for OP-B. ) are superimposed to the temperature field. Both cases exhibit a typical swirled flame with a central recirculation zone around the axis of the swirler. An interesting feature of the flame is that chemical reactions occur also near the cooling film, as well as along the dilution jet where fresh air flows into the combustor. Theoretically, the role of the dilution holes is to confine the flame in the upstream zone by suddenly reducing the fuel/air ratio. In both cases, this objective is reached and the temperature drops down after the dilution jets. However, in the case of the operating point OP-A (Fig. 3.3 top), a hot pocket of temperature is released after the dilution jets and the mixture is still reacting when convected to the downstream nozzle by the mean flow. A set of temperature snapshots is displayed in Fig. 3.4 showing how this hot pocket is created, convected and ejected by the nozzle. On the contrary in the case of the operating point OP-B, the dilution jets play their role and the flame is confined (Fig. 3.3 bottom). Only burned gas are convected to the downstream nozzle.

LES preliminary results

The pressure signal is recorded at a probe (denoted A, see Fig. 3.3) located at the entrance of the nozzle. Results for the operating point OP-A and OP-B are depicted in Fig. 3.5 by the solid and dotted line, respectively. In the case of OP-A, the pressure exhibits strong harmonic fluctuations at a fundamental low-frequency close to the experimental value (330 Hz). The mean pressure is approx. 1.65 × 10 6 Pa. The amplitude of the limit cycle is quite large and is in the range of pressure fluctuations amplitudes measured in the experiments (approx. 5 × 10 4 Pa). In the case of OP-B, the pressure signal is typical of broadband noise with a smaller fluctuation level. The mean pressure is also higher (approx. 1.76 × 10 6 Pa) than for the OP-A case.

In order to assess the influence of the shapes of the downstream nozzle, two computations have been performed using the baseline mesh for the operating point OP-A, both with the flat and circular geometry (Runs #1 and #3, see Table 3.1). Figure 3.6 presents the amplitude of the Fourier transform of each signal recorded at probe A. The first two peaks at 320 Hz and 640 Hz are virtually the same for both nozzle geometries while the frequency of the third and fourth peaks (approx. 960 Hz and approx. 1280 Hz, respectively) are slightly shifted. This shows that the estimation of the instability is not sensitive to the nozzle geometry for the low-frequency mode of interest. Only the circular nozzle geometry will be considered in the remaining of the study.

To make sure that the LES results are mesh independent, OP-A and OP-B were solved with the two levels of mesh refinement (Runs #5 and #6 in Table 3.1). Figure 3.7 presents the modulus of the Fourier transform of the pressure signal recorded at probe A for Runs #5 (fine mesh) and #1 (baseline mesh). Results for the first unstable mode are similar even if the amplitude of fluctuations computed with the baseline mesh is slightly underpredicted and their frequency slightly shifted. This can be explained by the fact that the baseline mesh is more dissipative and that entropy waves are dissipated too rapidly during their convection to the downstream nozzle. The modulus of the Fourier transform of the pressure signal recorded at probe A for the operating point OP-B is also displayed in Fig. 3.7. A peak is present at the frequency of approx. 435 Hz and has an amplitude of 2500 Pa. It is remarkable that pressure fluctuations for the operating point OP-B are very weak compared to OP-A, confirming the capacity of LES to capture and discriminate between stable and unstable regimes on both meshes of Table 3.1.

Additionnaly, to make sure that the LES results are independent of the subgrid scale model for the turbulent viscosity µ t , the Run #5 was computed with the classical Smagorinsky and the Sigma models as described at Sec. 3.1.2. Results for the modulus of the Fourier transform of the pressure signal recorded at probe A are depicted in Fig. 3.8. Results for the first unstable mode are similar both in frequency and amplitude. For the second unstable mode, results with the Sigma model slightly differ: the frequency is shifted of a few Hertz while the amplitude is under-estimated. However for the phenomena of our interest at approx 320 Hz, Fig. 3.8 demonstrates that the results are independent of the modelling of the turbulent viscosity µ t . Finally, as results are very similar even for mesh refinement and changes of subgrid scale model, the LES computations are considered robust.

Dynamic Mode Decomposition analysis

Dynamic Mode Decomposition analysis

Theoretical background

Experiments or simulations of flows often lead to complex, dynamic and non-linear systems which can not be analysed directly by simple visualisation or statistical techniques. Advanced methods have emerged during the last decades to provide a detailed analysis of many flow fields. In this context, the extraction of coherent structures from the flow as well as their stability analysis have become essential. The time-evolution of a flow is described by: dv

(t) dt = Av(t) (3.7)
where v is a vector containing the physical variables representative of the flow and A a linear operator. When integrated between t and t + ∆t while keeping A constant, Eq. (3.7) leads to:

v(t + ∆t) ≈ e ∆tA v(t) (3.8)
In practice, the flow is time-discretised and represented by a set of N vectors whose M components (M >> N ) contain any relevant quantity taken either from experimental or numerical data at N different instants. Mathematically speaking, the data are represented by the following sequence of vectors:

V N 1 = {v 1 , v 2 , v 3 , ⋯, v N } (3.9)
where N denotes the last snapshot. If the state vectors contain enough information, the flow dynamics can be approximated by a linear operator à = e ∆tA such that :

v n+1 = Ãv n , ∀n (3.10)
The à operator may be studied through its eigenvectors Φn and eigenvalues λn which provide the model structure and frequency/growth rate respectively. In practice, the eigenvalue problem: ÃΦ n = λn Φn (3.11) can not be solved since à is either unknown or too large. The Dynamic Mode Decomposition (DMD) developed by [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] is an algorithm to approximate the matrix Ã.

The DMD is an extension of a previous method presented by [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF] and based on the Arnoldi algorithm. The equations (3.9) and (3.10) may be expanded into the following relation:

ÃV N -1 1 = V N 2 = V N -1 1 S + re T N -1 (3.12)
where r is the residual vector which tends to 0 when the number of snapshots N increase, e N -1 is the (N -1) th unit vector and the matrix S is of companion type with:

S = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 a 1 1 0 a 2 ⋱ ⋱ ⋮ 1 0 a N -2 1 a N -1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (3.13)
where the coefficient {a 1 , a 2 , a 3 , ⋯, a N -1 } are the only unknowns. The matrix S shifts the snapshot index 1 to N -1 and is a low-dimensional representation of the full system. Then the main idea is that S and à share the same eigenvalues since they are similar from Eq. (3.12) if the residual vector in the RHS is indeed negligible. The advantage of considering S is its reduced size (N -1 × N -1) compared to the size of à (M × M ). One drawback of the DMD is that the operator S may be ill-conditioned in practice, especially when dealing with noisy signals. Several formulations of S are available in the literature, as for example by expanding V N -1 1 via a Singular Value Decomposition (SVD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] or a QR-decomposition [START_REF] Schmid | Applications of the dynamic mode decomposition[END_REF]. In the present work, the variant proposed by Chen et al. (2012) is implemented and described as follows. The operator S is reformulated by using Eq. (3.12):

S = (V N -1 1 ) -1 V N 2 (3.14)
Introducing first the following SVD: V N -1 1 = UΣW H where U and W are unitary and Σ is diagonal, the Moore-Penrose inversion of (V N -1 1) -1 leads to:

V N -1 1 -1 = WΣ -1 U H (3.15)
Introducing the substitution U = V N -1 1 WΣ -1 into Eq. (3.15) and replacing the whole expression into Eq. (3.14) gives:

S = WΣ -1 Σ -1 W H V N -1 1 H V N 2 (3.16)
The operators W and Σ may be then simply computed by forming the following diagonalization:

(V N -1 1 ) H V N -1 1 = WΣ 2 W H (3.17)
This construction of S is robust even when V N -1 1 is nearly rank-deficient, and computing the product

V N -1 1 H V N -1 1
requires less memory (see [START_REF] Chen | Variants of dynamic mode decomposition: Boundary condition, koopman, and fourier analyses[END_REF] than performing directly a SVD. The last step is to solve the eigenproblem Ss n = µ n s n where s n and µ n are the eigenvectors and eigenvalues of S. From Eq. (3.12) (neglecting the residual term), the relations between the eigenvectors Φn and eigenvalues λn of à and their counterparts s n and µ n are:

Φn = V N -1 1 s n (3.18) λn = µ n (3.19)
Recalling that à = e ∆tA , the eigenvectors and eigenvalues of A are:

Φ n = Φn = V N -1 1 s n (3.20) λ n = ln( λn ) ∆t = ln(µ n ) ∆t (3.21)
Thus studying the (N -1) × (N -1) matrix S gives information about the whole flow dynamics A.

Application of DMD to the SAFRAN combustor

Dynamic Mode Decomposition (DMD) was applied to the LES data in order to better understand the nature of the instability illustrated in Sec. 3.2. For this purpose, 200 snapshots were recorded over a time range corresponding to approx. 8 cycles of the instability, thus leading to a sampling of 28 snapshots by period. This amount of data is sufficient for the DMD to break down the reactive turbulent flow into dynamically relevant structures with periodic evolution over time [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. The input vectors (see Eq. (3.9)) for the DMD algorithm were built from the nodal values of pressure, static temperature, entropy, reaction rate and the three-components of the velocity at each grid point of the mesh used for the LES. Note that the entropy variable, denoted s, is not directly an output of the code and has been computed for each node of the domain from LES data with the following expression:

s(x) = C p (x) γ(x) ln p(x) ρ(x) γ(x) (3.22)
Although not detailed in the present study, several trials have been performed by changing the total number of snapshots or by composing the input vector with different sets of variables, showing that the number of variables taken into account in the DMD does not change the conclusions drawn in the remaining of this study. Moreover, variations in the total number of snapshots have a limited impact and only shift the frequency of the most unstable mode by a few Hertz.

The DMD results of the LES of OP-A and OP-B are depicted in the left and right images of Figs. 3.9 and 3.10, respectively. Figure 3.9 shows the L 2 norm Φ n L 2 = Φ H n Φ n of the pressure fluctuations associated to the different DMD modes. Figure 3.10 shows the amplification rate I(λ n ) of the pressure DMD mode. The DMD analysis for operating point OP-B does not exhibit a strong emerging mode. Although not evident, a peak is present at the frequency of 434 Hz and has an amplitude of approx. 3500 Pa, in agreement with experiences and the Fourier analysis. However results for the unstable operating point OP-A presented in the left of Fig. 3.9 exhibit a strong peak at 323 Hz, in good agreement with the Fourier analysis conducted at Sec. 3.2 (see Fig. 3.7). Additional peaks with a decreasing amplitude are also present at the frequencies of 648 Hz, 969 Hz and 1290 Hz. A close look at the spectrum of the amplification rate I(λ n ) depicted in Fig. 3.10 shows that each mode featuring a strong pressure activity has a null amplification. This result was expected because the DMD analysis is performed on LES data recorded when the thermoacoustic instability has reached a limit cycle. The transition from the stable to unstable operating points, i.e when the energy is growing inside the combustor, has not been considered in the present study. Note also that the eigenmodes extracted by the DMD are related by a simple relation f n = nf 0 , suggesting that the modes of higher frequency are harmonics of a fundamental mode f 0 .

The modulus of the fluctuating temperature field for the DMD mode with the highest amplitude at 323 Hz is depicted in the top of Fig. 3.11. The relative phase of the temperature and pressure mode is also recorded along the combustor, from the swirler to the exit nozzle. The path where the data are extracted is represented by a white line and the phase is plotted in the bottom of Fig. 3.11. Two different zones can be distinguished, each with a specific behaviour of the temperature fluctuations:

• Zone I corresponds to the primary zone between the swirler and the dilution holes.

The unsteady swirled flame naturally generates temperature fluctuations, but this latter mechanism is confined in the primary zone (I). Indeed, the flame is surrounded by a zero fluctuating temperature area, suggesting that no pockets of temperature generated by the flame are connected with the zone II downstream of the dilution holes. This statement is supported by the behaviour of the phase, which shifts from -π to π inside the zone I. Fluctuations of temperature are then rolling up inside the flame and no convection to the downstream zone II occurs.

• Zone II extends from the dilution holes to the exit nozzle. Strong fluctuations of temperature occur along the jets of fresh air that flow from the dilution holes. A large pocket of fluctuating temperature is formed at the top end of the jets, in the middle of the combustor. The phase unrolls monotonously from the formation area of the hot pockets to the downstream nozzle.

The analysis of the modulus and phase of the temperature for the first DMD mode at 323 Hz suggests that hot pockets of fluctuating temperature are released after the dilution holes and convected to the downstream nozzle. This statement confirms the scenario observed in Fig. 3.4 showing a series of temperature fields taken from the LES. However, the DMD gives more insight about the origin of these hot pockets: they are not created directly by the flame, which is confined in the primary region, but are generated by the hot gases interaction with the dilution jets.

In order to better understand the mechanism at the origin of the thermoacoustic instability, the fluctuating pressure, temperature and velocity fields from the first DMD mode at 323 Hz were reconstructed over one period. The four phases displayed in Fig. 3.12 support the idea that the unstable mode of interest relies, at least partly, on an entropyacoustic coupling. Note that the velocity is represented only by its component along the y-axis, i.e. the streaming axis of the dilution jets.

• At phase 0, the pressure is low everywhere within the combustion chamber and a pocket of cold gas is present downstream of the primary reaction zone, roughly at the middle of the combustion chamber. The velocity perturbations of the dilution jets are maximum: more fresh air is blown into the combustor,

• At phase π 2, the cold pocket is convected downstream and the unsteady pressure in the chamber is approximately zero. The velocity perturbations of the jets are very small,

• At phase π, the cold pocket interacts with the exit nozzle and a new pocket of hot gas is generated downstream of the primary zone, while the fluctuating pressure within the chamber is now positive and the dilutions jets have reached their minimum velocities, • At phase 3π 2, the pocket of hot gas is convected by the mean flow and the pressure within the chamber decreases, while the velocity perturbations of the dilution jets are now again very small.

As recalled in the introduction, each interaction between hot or cold pocket of gas and the nozzle generates acoustics which may propagate downstream (generating what is known as indirect noise) or upstream (generating another perturbation of the flame region and promoting the creation of a new entropy spot). This creates a pressure drop between the surrounding casing and the combustor. In the same phase relation, depending on the sign of the pressure drop, the dilution jets are either: (a) flowing fresh air into the combustor, decreasing then the combustion and creating the cold pocket, (b) or blowing fresh air out of the combustor, enhancing the formation of a hot pocket. When the pressure drop decreases to zero, the dilution jets are annihilated and the pocket of cold/hot temperature is no longer sustained and can be released and convected by the mean flow.

The acoustic response of the nozzle to an impacting temperature fluctuation can be characterised by the analytical theory based on the Euler equations of Duran & Moreau (2013). The amplitude of the reflected acoustic wave scales with the incident temperature fluctuation and the ratio p T is typically of the order of 100 Pa/K for most nozzles. As a consequence, for temperature fluctuations of approx. 350 K measured at 323 Hz from the DMD analysis of the LES data, the pressure fluctuations generated by the nozzle can then be computed with the help of the reflection coefficient and qualitatively estimated to be approx. 35000 Pa, which is in fair agreement with the pressure fluctuations (see Fig. 3.9) measured in the LES. The angle of p T is between 3π 4 and π, which is also in agreement with LES/DMD results (see Fig. 3.11). Hence, for this nozzle and flow regime, temperature pockets impacting the nozzle generate strong opposite-sign pressure fluctuations in the combustion chamber. Note that any acoustic pressure fluctuation generates an amount of temperature fluctuations, even in absence of nozzle or accelerated flow. Under the isentropic assumption one obtains:

T = T γ -1 ρ c 2 p (3.23)
Figure 3.13 shows that the isentropic temperature fluctuations are almost 25 times smaller that the temperature fluctuations measured in the LES via the DMD and are also out of phase. This demonstrates that the temperature fluctuations seen in zone II are not the thermal signature of the pressure (acoustics) fluctuations in the chamber; instead they are generated by the primary zone/dilution jets interaction. After interacting with the nozzle, they generate a strong pressure fluctuation of opposite sign. This demonstrates that the combustion instability does not rely on a purely acoustic mode. The idea that the unstable mode close to 320 Hz involves a coupling between entropy and acoustics is further supported by the numerical Helmholtz analysis of the combustor performed at Sec. 3.4.2, which shows that the smallest purely thermoacoustic frequency mode is close to 690 Hz, very far from the observed 320 Hz. Since Eq. (2.43) assumes no mean flow, it is necessary to restrict the study of thermoacoustic instabilities to only the combustion chamber (where the Mach number is always small). As depicted in Fig. 3.14, the domain is truncated at the outlet, before the nozzle, where the Mach number is M ≈ 0.1. Several computations have been performed with the AVSP Helmholtz solver (see Sec. 2.3) of the operating points OP-A and OP-B. In all cases, a very coarse mesh (approx. 10 5 nodes) was generated from the LES grid in order to represent the geometry properly. An axial cut of this latter mesh is depicted in Fig. 3.15. Additionally, a zero normal acoustic velocity condition was imposed over all boundaries except a) over the chamber exit where the nozzle impedance Z M ≠0 J,m (see Sec. 1.2.2) is imposed, and b) over the air inlet which is modeled by a constant mass flow rate m ⋅ n = 0, i.e. Z M ≠0 J,m = ∞.

The sound speed c is computed from the averaged temperature field given by LES and interpolated on the coarse acoustic mesh: the maps of sound speed are depicted in Fig. 3.16. In the cold region such as the casing where fresh air is flowing, the sound speed is approx. 300 m.s -1 for OP-A and approx. 500 m.s -1 for OP-B, while in hot regions where the flame takes place the sound speed is approx. 900 m.s -1 for both operating points. It is interesting to notice that in the combustor, the fields of sound speed are rather similar for both OP-A and OP-B. As explained in Sec. 3.3.2, one of the major effect of the combustion instability is to creates temperature heterogeneities that convect downstream to the nozzle with a certain frequency. In one period the temperature seen by the combustor after the dilution holes oscillates around a mean value. Hence, the time-averaged temperature field does not sensibly differ from OP-A and OP-B. However the effect of the combustion instability is evident when looking at the rate of heat release. For this purpose, the fields of n local and τ local are computed from the DMD mode to extract the fluctuating heat release rate q(x) at each point of the domain as well as the fluctuating velocity û(x ref ) ⋅ n ref at a reference point x ref located inside the swirler. The nτ fields are computed from the first DMD mode (see Appendix B). For Helmholtz computations of case OP-A, Fig. 3.17 presents the fields of n local and τ local computed from the DMD mode at 323 Hz. The reference point x ref for the velocity as well as the normal vector of reference n ref are also shown. For computations of the case OP-B, the fields of n local and τ local are extracted at 436 Hz which is the frequency (see Fig. 3.9) where most of the activity occurs and are depicted in Fig. 3.18. It is notable that for the OP-A unstable case, strong fluctuations of heat release rate occur in all the primary zone, from the edge of the swirler to the dilution jets. On the contrary, the OP-B stable case does not feature a significant pulsating flame. Most fluctuations of heat release rate are included in small compact zones located at the edge of the swirler or near dilutions holes.

Results

The main properties of the Helmholtz computations performed as well as the results for the first mode are gathered in Table 3.2 for OP-A and Table 3.3 for OP-B. The Active Flame column reports whether or not the flame-acoustic coupling (the RHS term of Eq. (2.43)) is included. The Outlet Impedance column reports whether the value of Z M ≠0 J,m (see Sec. 1.2.2) at the outlet is set to ∞ (simple m ⋅ n = 0 condition) or computed from the analytical nozzle theory of Duran & Moreau (2013). Computations of the operating points OP-A and OP-B feature the same behaviour. The first eigenmode obtained for OP-A is in the range of 690 -730 Hz, very far from the instability at 320 Hz encountered in the experiments and the LES computations. Similarly, the first mode computed for OP-B is found in the range of 817 -840 Hz, contrary to the (limited) thermoacoustic activity observed at approx. 435 Hz. Runs B show that the coupling with the flame only shifts the frequency by a few Hertz and has a damping effect on the acoustic modes. Runs C show that the use of more accurate boundary conditions has a limited effect on the frequency of the first eigenmode. Using an active flame model and proper boundary conditions within a Helmholtz solver framework is the most employed method to calculate thermoacoustic modes and to draw the stability map of a reactive system. However Runs D confirm that such a method totally fails to capture the desired unstable mode of interest (320 Hz) in this particular study. This result suggests that the instability encountered in the present combustor is not of the nature of a direct flame-acoustic coupling and supports the idea that it relies on an other mechanism. The next natural question is: how to evidence this mechanism if it exists and how to model it in the purely zero-Mach number acoustic framework? The next section is devoted to an advanced post-processing of the LES by using the Dynamic Mode Decomposition in order to compute the entropy/acoustic transfer functions required by the DECBC model. Consequently, the fluctuating velocity ûref for both the active flame and DECBC approach is chosen at the same reference point. Following the rationale developed in Sec. 2.1.4, the DMD analysis was also used to measure the convection time τ LES c from the flame region to the end of the combustion chamber (see Fig. Table 3.4 -Gain and time delay of the transfer functions modeling the entropy fluctuations for the operating points OP-A (T inlet = 557 K) and OP-B (T inlet = 762 K).

Helmholtz/DECBC results

Helmholtz computations are performed for each operating point OP-A and OP-B by using the respective fields of sound speed, active flame nτ fields and transfer functions from Table 3.4. The outlet impedance Z M ≠0 J,m is computed from the Duran & Moreau (2013) theory. Results for the first mode are gathered in Table 3.5. As already stated in Sec. 3.4.2, the first thermo-acoustic mode for the operating point OP-A oscillates at around 690 Hz, very far from the 320 Hz instability observed in the LES and experiment; this thermo-acoustic mode is slightly damped (decay rate of approximatively -40 s -1 ), even while taking into account the proper acoustic flux at the inlet/outlet as well as the coupling between acoustics and combustion through the active flame model. Introducing the DECBC approach corrects this behaviour: the Helmholtz solver now gives a strong unstable mode (growth rate 267 s -1 ) oscillating at 318 Hz, in fair agreement with the LES and experiment. Note that the frequency of oscillation obtained (318 Hz) is rather close to the first DMD mode (323 Hz), thus justifying a posteriori the fact that the transfer function gain and delay were supposed independent of frequency when performing the Helmholtz computations. The shapes of the first modes found either by the DMD analysis or the Helmholtz computations with the DECBC approach are depicted in Fig. 3.21. The fluctuating pressure p is made non dimensional by setting max(R(p)) = 1. It is remarkable that the eigenmode computed by the Helmholtz solver with the DECBC approach (bottom of Fig. 3.21) is rather similar to the DMD mode (top of Fig. 3.21). The maximum amplitude of fluctuation is located inside the whole combustion chamber, while the casing exhibits a structure close to a quarter wave mode.

Operating Point Frequency

[Hz] Growth/Decay rate [s -1 ] OP-A 318 +267 OP-B 429 +12
Table 3.5 -Helmholtz computations with the DECBC approach.

The methodology described previously is also applied to the operating point OP-B. Results exhibit the same behaviour (Table 3.5). A mode is found at 429 Hz with a small growth rate, consistent with the LES result which showed a small amplitude mode at 435 Hz. In [START_REF] Silva | Combining a helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF], authors mention that the growth rates assessed by a Helmholtz solver should be offset by a negative quantity which represents all the damping effects not included in the zero-Mach number framework (acoustic boundary layers, vortex generation in shear layers, etc.). Consequently the mode predicted for case OP-B at (+12 s -1 ) is consistent with the DMD analysis and experiments. Note also that the growth rate obtained for OP-A is much larger than for OP-B in full agreement with the LES/experimental results which display a strong instability in the first case only. The acoustic-entropy coupling allows, once included in a low-order Helmholtz solver, to recover the low frequency instability of the LES/experiments. This demonstrates that this instability is indeed driven by the acoustic waves generated by strong entropy fluctuations being convected through the nozzle. To the authors knowledge, this is the first time that such an instability is identified in an aeroengine configuration.

Application of DECBC: Optimisation of combustion instabilities

A novel strategy to take into account mean flow effects in a zero-Mach number Helmholtz solver has been developed in chapters 1 and 2. It relies on:

• a new formulation of the acoustic impedance to capture the energy damping due to the effect of flow convection on the propagation of acoustic waves;

• a new boundary condition (DECBC) that introduces the entropy/acoustic coupling with the help of transfer functions that model the generation and convection of entropy fluctuations to the downstream exit.

When applied to a realistic aero-engine configuration consisting of a choked nozzle mounted on a combustion chamber, such approach proves able to capture a low frequency unstable mode in full agreement with LES computations and DMD analysis (see Sec. 3.4). Moreover, it was shown that a purely thermo-acoustic analysis cannot reproduce such instability.

The purpose of the present section is to study the ability of a Helmholtz solver coupled with DECBC to design a new combustor so as to reduce the thermoacoustic instability. The present study focuses on the optimisation of the downstream nozzle via the modification of its geometry. In a first step, the ability of AVSP/DECBC to follow the impact of the changes of the nozzle length on the unstable modes is tested by computations of a simple 1D test case. Then, AVSP/DECBC computations of the realistic aero-engine combustor described in chapter 3 are resumed. The gradual increase of the nozzle length shows that the combustion instability is progressively damped. Finally, a new set of LES computations confirms that the new nozzle designed during the optimisation study with AVSP/DECBC makes the combustor more stable.

Influence of the nozzle length on thermoacoustic instabilities

Description of the test case

The configuration investigated in the present section is similar to the one studied in chapter 2 and depicted in Fig. 3.22. It consists of a combustion chamber of constant cross section and fixed length L c where a compact flame (thickness δ f = 0.005 m) is located at x = x f . The chamber is connected to an isentropic choked nozzle of length L-L c . The aim of the present section is to study the impact on thermoacoustic eigenmodes of a gradual increase of the length L of the nozzle. Thermodynamic parameters used in the computation of eigenmodes of the configuration are presented in Table 3.6. The eigenmodes of this simple configuration are computed by two different methods:

T 1 [K] p 1 [Pa] M 1 T 2 [K] γ r [SI] L c [m] x f [m] δ f [m]
• Method I num is the solution in the whole domain of the LEE numerical solver described at Sec. 1.3.1. This is considered as the reference solution. Note that the mean flow field inside the nozzle is described by an analytical distribution of the mean flow Mach number as follow:

M 2 (x) = M 2 (L c ) + (M out (L) -M 2 (L c )) x -L c L -L c 3 (3.26)
In this equation, M 2 (L c ) is the Mach number at the inlet of the nozzle and is deduced from the flow field parameters upstream of the flame. M out (L) is the Mach number imposed at the exit of the nozzle. In the present study, M out (L) = 1.5 so as to ensure that the flow is choked at the nozzle throat located at x th . As the total length L is a variable parameter, the throat location x th will then follow accordingly.

• Method AVSP/DECBC is the solution of the Helmholtz solver only in the combustion chamber region of the domain. The DECBC approach described at Sec. 2.1.4 is used as well as the impedance defined with the set variables Ĵ, m so as to ensure that the right acoustic energy flux is imposed. The impedance Z M ≠0 J,m is computed from the analytical method of Duran & Moreau (2013) by imposing the nozzle geometry taken from the the computation with method I num and which follow the Mach distribution from Eq. (3.26). The sound speed and density fields required to solves Eq. (2.43) are deduced from the imposed temperature field T 1 for x < x f and T 2 for x > xf . The gains and time delays required by the transfer function of the DECBC model (see Eq. (2.54)) are invariant and summarised in Table 3.7, as well as the location of the reference point x ref for the velocity fluctuations, taken slightly upstream of the flame.

G us τ us [s] G c τ c [s] x ref [m]

Results

In the present study, the length L -L c of the nozzle is gradually raised from 0 (compact conditions) to 0.5 m, which is the half of the combustor length. The frequencies and growth rates of the first four eigenmodes are depicted in Fig. 3.23. As expected, the main effect of extending the nozzle length is to create a much larger combustor. Hence the frequency decreases linearly as the nozzle length increase. The impact of the nozzle length on the growth rate is however less intuitive. Depending on the configuration, eigenmodes of the combustor may be stabilised or become unstable. Results with the AVSP/DECBC method are close to the reference method although some discrepancies occur. However the curves follow the same behaviour, demonstrating the capacity of the DECBC model to reproduce the evolution of both the frequency and growth rate depending on the nozzle length, which is virtually taken into account through the impedance.

Application to the realistic gas turbine AVSP/DECBC Results

Following the methodology employed in Sec. 3.4, computations of the OP-A unstable case are performed with the 3D Helmholtz solver with only the modification of the exit impedance Z M ≠0 J,m . For each computation, the geometry of the nozzle used for the computation of the impedance Z M ≠0 J,m is completed by an additional segment L noz located between the baseline nozzle and the combustor (see Fig. 3.24). Note that the actual dissipation of the entropy along the nozzle is not taken into account in the present model. The fields of sound speed (see Fig. 3.16) and nτ for the active flame (see Fig. 3.17), as well as the DECBC parameters (see Table 3.4) remain unchanged.

Results are depicted in Fig. 3.25. The frequency and the growth rate of the unstable mode decrease linearly as the length L noz increases. Although the frequency of the unstable mode is only shifted by a few Hz, it is noteworthy that the growth rate can be strongly damped by the extension of the nozzle. Indeed, the introduction of an additional length at the entrance of the nozzle has the effect to add an additional time lag that impact the relative phase between p and T . For the present configuration, the resulting effect is to lower the rate of acoustic pressure generated at the nozzle when impacted by temperature heterogeneities. The present study should be understood as qualitative. Indeed, as explained in Sec. 3.4.4, because all of the damping effects are not taken into account in the modelling, it is difficult to know exactly which value of the growth rate defines the limit between a stable or an unstable system. However such optimisation study provides a qualitative information about the evolution of the growth rate depending on the evolution of certain parameter, here L noz . According to Fig.

LES Results

The Run #3 (see Table 3.1) is chosen to test the geometry change of the nozzle. A new mesh is generated so as to handle the new geometry of the nozzle. Note that for the combustor as well as the geometrical details of the swirler, the grid refinement is similar to the baseline mesh used in Run #3. Fig. 3.26 presents the pressure signal recorded at the location of probe A (see Fig 3 .3). At 0.12s of Run #3, the LES solution is interpolated on the new mesh containing the extended nozzle. After a short transient, the flow reaches a new state. The pressure is still oscillating with a pronounced frequency, but the amplitude of the oscillation is severely reduced. When the combustor is mounted with the nozzle extended by a segment of length L noz = 0.03 m, the pressure amplitude of the first mode is approximatively five times smaller than in the baseline configuration. This result exhibit a stronger damping of pressure oscillations than in the AVSP/DECBC study, but one should note that in this latter methodology the dissipation of entropy waves is not taken into account during their convection in the modelled nozzle. Consequently the damping rate found by AVSP/DECBC is expected to be under-predicted. The present LES result is now in the range of the amplitude found with the analysis of the stable OP-B case (see right image of Figs. 3.9 and 3.10). Moreover, the frequency of this mode (323 Hz with the baseline nozzle) is shifted to 314 Hz, which is consistent with the results predicted by the AVSP/DECBC methodology (see Fig. The modulus of the pressure and temperature for the DMD mode at 314 Hz are depicted in the top and middle images of Fig. 3.30. The relative phase of the temperature and pressure mode is also recorded along the combustor, from the swirler to the exit nozzle. The path where the data are extracted is represented by a white line and the phase is plotted in the bottom of Fig. 3.30. Similarly with the DMD analysis of the unstable OP-A case with the baseline nozzle (see Fig. 3.11), the relative phase between p and T is still π. However in the case of the extended nozzle, the amplitude of the temperature heterogeneities are approximatively reduced by five times. Consequently the amplitude of the pressure generated at the nozzle is severely reduced with the same order of magnitude.

Conclusions

In the present section, the Helmholtz solver AVSP coupled with the DECBC approach has been used to perform a sensitivity analysis of the length of the downstream nozzle that ends the combustion chamber. Such methodology shows that the growth rate of the unstable mode can be significantly reduced when playing with the length of the nozzle. A new set of LES computations have been performed with an extended geometry of the nozzle. The combustion instability is severely reduced while the frequency is slightly shifted down, in fair agreement with the AVSP/DECBC prediction. It could be argued that the length extension of the nozzle is not applicable in an actual aero-engine configuration. This is partly true but one should keep in mind that the purpose of this present section is to demonstrate the potential of the methodology AVSP/DECBC to optimise the combustion instability. A novel approach to design nozzles has been presented by [START_REF] Giauque | Thermoacoustic shape optimization of a subsonic nozzle[END_REF]. It is based on genetic algorithms and allow to adjust several geometrical parameters at the same time, in addition to the length of the sonic throat.

Conclusions and perspectives

High-fidelity computations of a realistic aerocombustor that features a strong lowfrequency instability under certain operating conditions have been performed. Experiments show that the establishment of an instability depends to the first order on the fraction of liquid water ingested by the engine. Large-Eddy Simulations (LES) proved able to distinguish between the stable/unstable behaviours, the presence of water being accounted for by lowering the inlet temperature of the air entering the combustion chamber. Based upon two operating points representative of the stable/unstable states of the combustor, a Dynamic Mode Analysis (DMD) has been performed on the recorded LES data. For each case, a low frequency mode is found, featuring either strong or very weak fluctuations of all physical variables. The DMD analysis demonstrates that pockets of fluctuating entropy/temperature are created inside the combustor and convected to the downstream nozzle, creating a reflected acoustic perturbation that goes back to the flame.

A purely thermo-acoustic analysis has been performed with the Helmholtz solver AVSP. The coupling between acoustics with the combustion is taken into account by an active flame model based on a nτ formalism. It has been shown that when the computational domain is represented under the M = 0 assumption, a non-zero-Mach number element can simply be taken into account by imposing a proper acoustic impedance at the boundaries so as to mimic the mean flow effects in the outer, not computed flow domain. To this respect, it has been demonstrated that the best way to define the boundary impedance is the use of energetic state variables Ĵ and m. This latter choice introduces an acoustic energy flux mismatch at the M = 0 M ≠ 0 interface proportional to the square of the Mach number, instead of being proportional to M when relying on the acoustic variables p and û. The analytical asymptotic theory has been validated on a simple test case by computing eigenmodes either by solving the Linearized Euler Equation or the zero-Mach number Helmholtz equation. However, despite the fact that the high speed flow turbine is properly taken into account with the help of a formalism that reduces the error introduced in the acoustic domain, for both operating points results exhibit a damped eigenmode which frequency is very far from the low-frequency mode observed in the experiments and LES computations. This suggests that the mechanism leading to the combustion instability relies on a coupling between entropy and acoustic waves inside the nozzle, and as this latter mechanism is neglected in the Helmholtz analysis due to the zero-Mach number assumption, the low-frequency mode can not be reproduced with this approach.

A Delayed Entropy Coupled Boundary Condition (DECBC) model has been developed as a mean to recover some of the convective effects when representing a thermo-acoustic system under the zero-Mach number formalism. In this view, a simple model was first used in order to assess the entropy fluctuations at the exit of the combustion chamber. This model consists of two steps, one for the entropy generation in the flame region, a second one for the convection/dissipation of the entropy spots through the combustion chamber. This results in a transfer function which relates a reference acoustic velocity in the burner and entropy at the exit of the chamber and which was assessed by post-processing LES data via the DMD analysis. The acoustics generated by the convection of the entropy spots through the exit nozzle are then treated by applying a proper boundary condition which couples entropy and acoustic quantities.

Helmholtz computations coupled with the DECBC approach are conducted for both stable/unstable operating points. Results exhibit low-frequency eigenmodes with frequency, growth rate and shape in very good agreement with the experimental and LES data. This enlightens the major role played by the coupling between entropy and acoustic waves inside the nozzle on the establishment of this combustion instability, contrary to the coupling between the flame and the acoustics or the sound speed field which have a minor influence.

A limitation of the proposed approach is that the Helmholtz/DECBC solver must be fed by data coming from LES results. This may appear as a strong limitation since the cost of a LES must be paid before the DECBC methodology can be applied. Still, the Helmholtz/DECBC approach remains useful in situations where the same upstream velocity-exit entropy transfer function can be reused for multiple test conditions. The potential of this strategy has been demonstrated by conducting a set of Helmholtz/DECBC computations with different shapes of exit nozzle. It has been shown that the growth rate of the unstable mode can be significantly reduced for a certain length of the nozzle. A new set of LES computations have been performed with a selected nozzle geometry. The combustion instability is severely reduced while the frequency is slightly shifted down, in fair agreement with the AVSP/DECBC prediction.

Another potential field of application is the study of annular combustors; in this case, one can imagine that a simple one sector LES (containing only one burner) could be performed to feed a Helmholtz/DECBC computation of the whole annular combustor with several identical burners, thus avoiding the extremely CPU demanding LES of the full combustor. This strategy is justified if the flames behave independently and was already used for accounting for the flame response in annular combustors.

In conclusion, the present thesis proposes a practical methodology to reproduce and optimise low-frequency combustion instability at a very low computational cost. This opens the way to possible repeated optimisation studies for the design of aero-engines combustion chamber.

Publications

Articles in International Journals The global model described at Eq. (B.4) is only justified in the case of a compact flame (the thickness of the flame is smaller than the wavelength of the eigenmode), which is not the most general case. The alternative, more consistent with Eq. (2.41), is to use a local flame model that allows to write a flame transfer function for each point of the domain, like Eq. (B.2). The relation between the global and local formulations is given by the relation between the global heat release rate Q and the heat release rate by volume unit q(x ref ), as follow:

Q = V f q(x ref ) dV (B.5)
where V f is the volume of the flame.

Three possible scenario arise depending upon the origin of the data characterising the flame response:

B.1 From experimental measurements

In some cases, one should want to use experimental data. The experiments allow to measure the global characteristics of the flame response, i.e. n g and τ g , but do not give any informations on the local information. One must check with the experimenters that the given values of n g are consistent with Eq. (B.4). Indeed, experimenters often use another N g defined by Q = Q û u N g e jωτg where N g is made non-dimensional. As the time delay is generally not a problem, the interaction index must be checked. If experimenters give N g , hence n g = N g Q u where Q must be given by experimenters and corresponds to the total heat release of the burner (in [W]). Of course it exists an infinity of n local -τ local fields that satisfy this latter relation. One solution is to suppose that n local and τ local are 0 everywhere except in the flame zone where they are constant. This assumption leads to:

n local = n g V f (B.7) τ local = τ g if x ∈ V f (B.8)

B.2 From a theoretical model

The analytical model of [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF] for an infinitely thin 1D flame is:

Q = S f γp γ -1 n crocco û(x ref )e jωτcrocco ⋅ n ref (B.9)
where S f is the flame surface and p the mean pressure. Similarly to the procedure described for the experimental case, the use of Eqs. (B.2), (B.4) and (B.9) leads to:

S f γP 0 γ -
τ local = τ crocco if x ∈ V f (B.13)
where δ f is the flame effective thickness solved by the AVSP code and V f = S f δ f is the flame volume in the simulation.

B.3 From a LES computation

In this case a LES computation of either a forced flame or an auto-sustained flame can be used. The fields of n local and τ local can be directly extracted from a series of snapshots sufficiently time resolved. From Eq. (B.2) it comes: The case where a node j is located at a boundary differs to interior case in the sense that the dual cell is truncated. The node j is then surrounded by three elements only. The pressure gradient may be expressed with the help of the impedance. This is the so-called Robin boundary condition. The Eq. (2.53) is recalled here in the formalism of the node j:

F (x, ω) = q(x) û(x ref
∇p ⋅ n j = jω c j Z j pj (C.8)
Application of Eq. (C.4) over the three elements surrounding the node j leads to:

∇ ⋅ (1 ρ ∇) p j = 1 S N j [1 ρ Ω e1 ∇p Ω e1 ⋅ 1 2 L 1 n 1 + 1 ρ Ω e2 ∇p Ω e2 ⋅ 1 2 L 2 n 2 + 1 ρ Ω e3 ∇p Ω e3 ⋅ 1 2 L 3 n 3 (C.9) -1 ρ Ω e1 ∇p Ω e1 ⋅ l 13 n -1 ρ Ω e3 ∇p Ω e3 ⋅ l 33 n ]
From Eq. (C.8) it comes:

1 ρ Ω e1 ∇p Ω e1 ⋅ l 13 n -1 ρ Ω e3 ∇p Ω e3 ⋅ l 33 n = 1 ρ j c j

  c concerning the convection of entropy waves (DECBC model) us concerning the generation of entropy waves (DECBC model) p,s formulation with fluctuating pressure and entropy as state variables p,u formulation with fluctuating pressure and velocity as state variables J,m formulation with fluctuating enthalpy and mass as state variables f concerning the flame th concerning the throat of a nozzle i imaginary part of complex value r real part of complex value in wave entering the domain out wave leaving the domain ref reference state or position inlet concerning boundary condition at the inlet of the domain Superscripts ˆcomplex amplitude of a fluctuation mean value ̃Favre averaged value M =0 formulation under the zero-Mach number assumption (i.e. Helmholtz solver) M ≠0 formulation with Mach number effects (i.e. LEE) + concerning the downstream travelling acoustic wave-

Figure 1 -

 1 Figure 1 -Premises of the gas turbine. Left: The Aeolipile of Hero. Right: The Chimney Jack of Leonardo Da Vinci

  ; Ducruix et al., 2003; Lieuwen, 2003; Culick & Kuentzmann, 2006; Candel et al., 2009; Huang & Yang, 2009; DeGoey et al., 2011; Poinsot & Veynante, 2011; Swaminathan & Bray, 2011; Candel et al., 2013).

Figure 2 -

 2 Figure 2 -Combustor damaged by a combustion instability. Left: initial combustion system. Right: damaged combustion system. Reproduced from Huang & Yang (2009).

Figure 3 -

 3 Figure 3 -Typical pressure signal of the time evolution of a combustion instability. Reproduced from Poinsot & Veynante (2011).

  ; Schmitt et al., 2007; Wolf et al., 2012a; Gicquel et al., 2012).

Figure 4 -

 4 Figure 4 -Schematic view of the proposed modelling strategy. Figure (a): The whole domain is computed while taking the mean flow into account. Figure (b): The combustion chamber is solved under the zero-Mach number assumption. The acoustic environment from compressor and turbine is accounted for through impedances.

Figure 5 - 1 . 1

 511 Figure 5 -Schematic view of the organisation of the thesis. Some parts of the geometry are hidden for confidentiality reasons.

Fig. 4 .

 4 b).

  (2001) studied the acoustic response of a diaphragm which is crossed by a jet and constructed their analytical model with the fluctuating total enthalpy variable. Such a choice was also followed by[START_REF] Dupère | The absorption of sound near abrupt axisymmetric area expansions[END_REF] as well as[START_REF] Eldredge | The absorption of axial acoustic waves by a perforated liner with bias flow[END_REF] to model the absorption of sound near abrupt area expansions or close to perforated liners with a bias flow. In an experimental point of view, Karlsson & Åbom (2010) convert the pressure measured by microphones into fluctuating mass and total enthalpy so as to construct scattering matrix of T-junctions with the presence of a mean flow.

  25) Injecting Eqs. (1.20) and (1.21) into Eq. (1.25) leads to :

Figure 1 . 1 -

 11 Figure 1.1 -Schematic representation of the boundary between a domain computed with a Helmholtz solver and an acoustic element with the zero-Mach-number assumption.

Figure 1 . 2 -

 12 Figure 1.2 -Configuration investigated in section 1.3. (a): Complete geometry with its known acoustic boundary conditions. (b): Only a part of the domain is computed while the rest is modelled through impedances. (c): same as (b) but the computed flow domain is assumed at rest.

  (1.58) into Eq. (1.54), the entropy equation and perfect-gas equation of state yield: 56) -(1.58) and (1.62) describe the linear harmonic perturbations at finite Mach number.

3 . 1 Figure 1 . 3 -

 3113 Figure 1.3 -Mode shapes at M 1 = 0.0005 (left column) and M 1 = 0.2 (right column) for the first eigenmode with the numerical LEE solver I num () and the reference method I an (×). First two rows: modulus and phase of the fluctuating total enthalpy Ĵ, respectively. Last two rows: modulus and phase of the fluctuating mass flow rate m, respectively.

Figure 1 . 4 -

 14 Figure 1.4 -Frequency and growth rate of the first eigenmode for different inlet Mach number M 1 for the numerical LEE solver I num ( ), the analytical global method I an (×) and the truncated method (M ≠ 0) II (○).

Figure 1 . 5 -Figure 1 . 6 -

 1516 Figure 1.5 -Frequency and growth rate of the first (top) and second (bottom) eigenmode versus inlet Mach number M 1 with the reference methods I (×) and methods III pu ( ) and III Jm ( ).

Figure 1 . 7 -

 17 Figure 1.7 -Mode shapes at M 1 = 0.0005 (left column) and M 1 = 0.2 (right column) for the first eigenmode with the reference methods I (×) and methods III Jm ( ), and for the second eigenmode with the reference methods I (⊙) and method III Jm (). First two rows: modulus and phase of the fluctuating total enthalpy Ĵ, respectively. Last two rows: modulus and phase of the fluctuating mass flow rate m, respectively.

Contents 2. 1

 1 Introducing entropy-acoustic coupling in the Helmholtz framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.1.1 Method I : Global thermoacoustic model at M ≠ 0 . . . . . . . . . 47 2.1.2 Method II : Truncated analytical thermoacoustic model at M ≠ 0 47 2.1.3 Method III Jm : Truncated analytical thermoacoustic model at M = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.1.4 Method III DECBC Jm : Delayed Entropy Coupled Boundary Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 2.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 2.2.1 Validations and effect of the truncation . . . . . . . . . . . . . . . . 53 2.2.2 Frequency analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 2.2.3 Mode shapes analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.3 Extension of DECBC to 3D . . . . . . . . . . . . . . . . . . . . . . 58 2.3.1 The zero-Mach number Helmholtz equation . . . . . . . . . . . . . 58 2.3.2 Boundary conditions and DECBC implementation . . . . . . . . . 60 2.3.3 Numerical resolution strategy . . . . . . . . . . . . . . . . . . . . . . 62 2.3.4 Validations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

Figure 2 . 1 -

 21 Figure 2.1 -Configuration investigated in section 2.1.

Figure 2 . 2 -

 22 Figure 2.2 -Schematic representation of the transfer functions that models the generation of entropy fluctuations at the flame, and their convection to the exit boundary condition.

Figure 2 . 3 -

 23 Figure 2.3 -First two eigenmodes computed with methods I ( × ) and II ( ) versus inlet Mach numbers M 1 .

Figure 2 . 4 -

 24 Figure 2.4 -First three eigenmodes computed with methods II ( × ), III Jm ( ) and III DECBC Jm

M 1 =Figure 2 . 5 -

 125 Figure 2.5 -Mode shapes at M 1 = 0.0125 (left column) and M 1 = 0.05 (right column) for the second eigenmode with the reference method II (×), method III Jm ( ) and method III DECBC

Jm

(

Figure 2 . 6 -

 26 Figure 2.6 -Schematic representation of the curvilinear coordinates system of the computational domain. The nozzle boundary condition is located at ξ = 0. The normal n of the boundary is taken inward.

Figure 2 . 7 -

 27 Figure 2.7 -Schematic representation of the fixed-point method.

Figure 2 . 8 -

 28 Figure 2.8 -First three eigenmodes computed with the reference analytical method II ( × ), the zero-Mach number analytical method III DECBC Jm () and the numerical Helmholtz solver AVSP coupled with DECBC ().

Contents 3 . 1

 31 Experimental configuration and numerical set-up . . . . . . . . 67 3.1.1 Description of the experimental configuration . . . . . . . . . . . . 3.1.2 Numerical set-up . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 LES preliminary results . . . . . . . . . . . . . . . . . . . . . . . . . 73 3.3 Dynamic Mode Decomposition analysis . . . . . . . . . . . . . . . 79 3.3.1 Theoretical background . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Application of DMD to the SAFRAN combustor . . . . . . . . . . 3.4 Helmholtz analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 3.4.1 Numerical setup for Helmholtz solver . . . . . . . . . . . . . . . . . 3.4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.3 Application of the DECBC approach . . . . . . . . . . . . . . . . . 3.4.4 Helmholtz/DECBC results . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Application of DECBC: Optimisation of combustion instabilities 98 3.5.1 Influence of the nozzle length on thermoacoustic instabilities . . . Description of the test case . . . . . . . . . . . . . . . . . . . . . . . Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.5.2 Application to the realistic gas turbine . . . . . . . . . . . . . . . . 100 AVSP/DECBC Results . . . . . . . . . . . . . . . . . . . . . . . . . . 100 LES Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

Figure 3 . 1 -

 31 Figure 3.1 -Description of the configuration of interest. One sector of the azimuthal SAFRAN combustor is represented. The high pressure distributor downstream of the combustion chamber is modelled by a nozzle with equivalent cross section and which is either rectangular (top) or circular (bottom). Some parts are hidden for confidentiality reasons.

Figure 3 . 2 -

 32 Figure 3.2 -Axial cut of the meshes for LES computations. Top: baseline mesh (approx. 5.6 millions cells). Bottom: fine mesh (approx. 15 millions cells).

Figure 3 .

 3 Figure 3.3 displays typical snapshots (instantaneous temperature fields) of the LES for the operating points OP-A and OP-B. Contours of the heat release rate (1 × 10 9 W.m -3) are superimposed to the temperature field. Both cases exhibit a typical swirled flame with a central recirculation zone around the axis of the swirler. An interesting feature of the flame is that chemical reactions occur also near the cooling film, as well as along the dilution jet where fresh air flows into the combustor. Theoretically, the role of the dilution holes is to confine the flame in the upstream zone by suddenly reducing the fuel/air ratio. In both cases, this objective is reached and the temperature drops down after the dilution jets. However, in the case of the operating point OP-A (Fig.3.3 top), a hot pocket of temperature is released after the dilution jets and the mixture is still reacting when convected to the downstream nozzle by the mean flow. A set of temperature snapshots is displayed in Fig.3.4 showing how this hot pocket is created, convected and ejected by the nozzle. On the contrary in the case of the operating point OP-B, the dilution jets play their role and the flame is confined (Fig.3.3 bottom). Only burned gas are convected to the downstream nozzle.

Figure 3 . 3 -

 33 Figure 3.3 -Instantaneous fields of temperature with superimposition of one contour of Heat Release rate (white line, 1 × 10 9 W.m -3 ). Top: operating point OP-A. Bottom: operating point OP-B. The snapshots are from Runs #5 and #6, respectively (see table3.1). Some parts are hidden for confidentiality reasons.

  Figure 3.3 -Instantaneous fields of temperature with superimposition of one contour of Heat Release rate (white line, 1 × 10 9 W.m -3 ). Top: operating point OP-A. Bottom: operating point OP-B. The snapshots are from Runs #5 and #6, respectively (see table3.1). Some parts are hidden for confidentiality reasons.

Figure 3 . 4 -Figure 3 . 5 -

 3435 Figure 3.4 -Instantaneous fields of temperature with superimposition of one contour of Heat Release rate (white line, 1 × 10 9 W.m -3 ). Full cycle of the evolution of the temperature. Top left: the flame is confined in the primary zone. Top right: a hot pocket of temperature is created under the dilution holes. Bottom left: the hot pocket is convected downstream. Bottom right: the hot pocket is ejected through the exit nozzle.

Figure 3 . 6 -

 36 Figure 3.6 -Spectra of pressure signals from LES data at probe A performed on the coarse mesh and for the operating point OP-A. Solid line ( ): circular nozzle, Run #1 (see top of Fig. 3.1 and table 3.1). Dotted line ( ): flat nozzle, Run #3 (see bottom of Fig. 3.1 and table 3.1).

Figure 3 . 7 -

 37 Figure 3.7 -Spectra of pressure signals from LES data at probe A and for the circular nozzle geometry. Solid line ( ): OP-A and fine mesh (Run #5). Dotted line ( ): OP-A and coarse mesh (Run #1). Dot-dashed line with black circles ( • ): OP-B and fine mesh (Run #6).

Figure 3 . 8 -

 38 Figure 3.8 -Spectra of pressure signals from LES data at probe A and for the circular nozzle geometry and fine mesh (OP-A, Run #5). Solid line ( ): Smagorinsky model (Smagorinsky, 1963). Dotted line ( ): Sigma model (Nicoud et al., 2011).

Figure 3 . 9 -

 39 Figure 3.9 -Pressure amplitude (L 2 norm) spectrum from the DMD analysis of the LES computation. Left: OP-A. Right: OP-B. The frequency is R(λ n ).

Figure 3 .

 3 Figure 3.10 -Amplification spectrum (I(λ n )) from the DMD analysis of the LES computation. Left: OP-A. Right: OP-B. The frequency is R(λ n ).

Figure 3 . 11 -

 311 Figure 3.11 -Temperature fluctuations for the first DMD mode at 323 Hz of operating point OP-A. Top: Modulus over the medium plane. Bottom: phase of the temperature () and pressure ( ) along the path depicted by the white line.

Figure 3 . 12 -

 312 Figure 3.12 -Fluctuating pressure (left), temperature (middle) and velocity (right) from the DMD mode at 323 Hz for the operating point OP-A (T inlet = 557 K). From top to bottom, the four rows correspond to phases 0, π 2, π and 3π 2.

Figure 3 . 13 -

 313 Figure 3.13 -Solid line ( ): reconstructed temperature signal at the entrance of the nozzle measured from the first DMD mode at 323 Hz. Dotted line (): reconstructed acoustic temperature signal at the entrance of the nozzle computed with Eq. (3.23) and the pressure signal measured from the first DMD mode at 323 Hz. Note that the second curve is multiplied by a factor 25 to be comparable to the first one.

3. 4

 4 Helmholtz analysis 3.4.1 Numerical setup for Helmholtz solver

Figure 3 .

 3 Figure 3.14 -Isocontour of Mach number in the combustor for the computation at operating point OP-A. The black line represents the location where the domain is truncated (M ≈ 0.1).

Figure 3 .

 3 Figure 3.15 -Axial cut of the very coarse mesh (approx. 10 5 nodes) for AVSP computations.

Figure 3 .

 3 Figure 3.16 -Fields of sound speed c time-averaged from LES data. Top: OP-A, Run #5. Bottom: OP-B, Run #6.

Figure 3 .

 3 Figure 3.17 -Fields of n local and τ local computed from the DMD mode at 323 Hz and for the operating point OP-A, Run #5. The black dot indicates the location of the reference point x ref for the acoustic velocity.

Figure 3 . 18 -

 318 Figure 3.18 -Fields of n local and τ local computed from the DMD mode at 435 Hz and for the operating point OP-B, Run #6. The black dot indicates the location of the reference point x ref for the acoustic velocity.

Figure 3 . 19 -

 319 Figure 3.19 -Computational domain for the Helmholtz analysis. The point of reference and the zone of averaging used for the entropy generation modeling (Eq. 3.24) are displayed as well as the exit section where the entropy-acoustic boundary condition (Eq. 2.55) is applied.

Figure 3 . 20 -

 320 Figure 3.20 -Reconstructed signals from the first DMD mode at 323 Hz for the operating point OP-A (T inlet = 557 K). Dotted line ( ): reference velocity ûref ⋅ n ref . Solid line ( ): entropy < ŝf > averaged in the middle of the combustor. Dot-dashed line ( ): entropy < ŝ > averaged at the exit plane upstream of the nozzle.

Figure 3 . 21 -

 321 Figure 3.21 -Shapes of the first pressure eigenmode for the operating point OP-A (T inlet = 557 K). Top: DMD analysis (323 Hz). Bottom: Helmholtz computation with the DECBC approach (318 Hz). The fluctuating pressure p is made non dimensional by setting max(R(p)) = 1.

Figure 3 .

 3 Figure 3.22 -Configuration investigated in section 3.5.1.

Table 3 . 6 -

 36 Thermodynamic parameters for the computation of the configuration in Fig.3.22.

Table 3 . 7 -

 37 Parameters of the transfer functions used in the DECBC model to compute the configuration depicted in Fig 3.22.

Figure 3 . 23 -

 323 Figure 3.23 -First four eigenmodes (from top to bottom) computed with the M ≠ 0 method I num ( ) and the M = 0 method AVSP/DECBC (•), versus nozzle length L -L c .

Figure 3 .

 3 Figure 3.24 -Schematic representation of the extension of the nozzle geometry by a new segment L noz .

  (

  3.25), the growth rate predicted by the Helmholtz solver with DECBC is divided by two for a nozzle length L noz = 0.03 m. This value corresponds approximatively to the actual length of the baseline nozzle. In order to confirm that the present methodology is relevant, a new LES computation of the OP-A unstable case is performed with a nozzle extended by L noz = 0.03 m.

Figure 3 .

 3 Figure 3.25 -Results of AVSP/DECBC computation for different extension L noz of the nozzle geometry. Left: Frequency. Right: Growth rate.

Figure 3 .

 3 Figure 3.26 -Pressure signal from LES data at probe A.

Figure 3 .

 3 Figure 3.27 -Spectra of pressure signals from LES data at probe A and for the OP-A unstable case. Dotted line ( ): Baseline nozzle geometry (Run #3). Dot-dashed line with black circles ( • ): Extended nozzle (L noz = 0.03 m).

Figure 3 .

 3 Figure 3.27 presents the modulus of the Fourier transform of the pressure signal recorded at probe A for computations with the baseline nozzle (Run #3) and the extended nozzle. The mean value has been subtracted. A DMD analysis is also performed and results are gathered in Figs.[START_REF] Vervisch | Direct numerical simulation of non premixed turbulent flames[END_REF].28 and 3.29 for computations with the baseline (left images) and the extended nozzle (right images).

  3.25).

Figure 3 .

 3 Figure 3.28 -Pressure amplitude (volume weighted L 2 norm) spectrum of the DMD analysis of the LES computation of Run #3. Left: baseline nozzle. Right: extended nozzle (L noz = 0.03 m). The frequency is R(λ n ).

Figure 3 .

 3 Figure 3.29 -Amplification spectrum (I(λ n )) of the DMD analysis of the LES computation of Run #3. Left: baseline nozzle. Right: extended nozzle (L noz = 0.03 m). The frequency is R(λ n ).

Figure 3 .

 3 Figure 3.30 -Pressure and temperature fluctuations for the first DMD mode at 314 Hz of operating point OP-A with the optimised nozzle. Top: Pressure modulus over the medium plane. Middle: Temperature modulus over the medium plane Bottom: phase of the temperature () and pressure ( ) along the path depicted by the white line.

•

  E. Motheau, L. Selle and F. Nicoud. -Accounting for convective effects in zero-Mach number thermoacoustic models -Journal of Sound and Vibration, 333(1):246-262, 2014. • E. Motheau, Y. Méry , F. Nicoud and T. Poinsot. -Analysis and modelling of entropy modes in a realistic aeronautical gas turbine -Journal of Engineering for Gas Turbines and Power, 135(9):092602, 2013. • E. Motheau, F. Nicoud, and T. Poinsot. -Using boundary conditions to account for mean flow effects in a zero mach number acoustic solver -Journal of Engineering for Gas Turbines and Power, 134(11):111502, 2012.Peer-reviewed Papers • E. Motheau, L. Selle, Y. Méry, T. Poinsot, and F. Nicoud. -A mixed acousticentropy combustion instability in a realistic gas turbine -Center for Turbulence Research, Proceedings of the Summer Program, Stanford University, 2012 International Conferences (with scientific committee) • E. Motheau, Y. Méry , F. Nicoud and T. Poinsot. -Analysis and modelling of entropy modes in a realistic aeronautical gas turbine -ASME Turbo Expo 2013, pages GT2013-94224, San Antonio (Texas, USA), June 3-7 2013. • E. Motheau, F. Nicoud, and T. Poinsot. -Using boundary conditions to account for mean flow effects in a zero mach number acoustic solver -ASME Turbo Expo 2012, pages GT2012-68852, Copenhagen (Denmark), June 11-12 2012.

  The use of Eqs. (B.2), (B.4) and (B.5) allows to write the relation between the fields n local -τ local and n g and τ g : V f n local (x)e jωτ local (x) dx = n g e jωτg (B.6)

  ) = n local (x)e jωτ local (x) φ(x, ω) = arccos R (F (x, ω)) F (x, ω) B.3 From a LES computation is injected within stoichiometric proportions (Y C3H8 = 0.06, Y O2 = 0.219, Y N 2 = 0.72) at a temperature of T = 300 K. The inlet velocity is imposed with the following profile u inlet (y) = u 0 1 -y 2 y 2 0 with u 0 = 4 m/s and y 0 = 0.65 mm. The flame is forced by a pulsation of amplitude 0.4 m/s and frequency 1000 Hz is introduced at the inlet velocity. The time step between each recorded snapshot is 0.1e -4 ms. The DMD is performed over a base constructed with 300 snapshots, i.e. 3 periods. For the The Fourier transform method the signals are recorded over 30 periods. Results are depicted in Figs. B.2 and B.3. As expected the fields of n local and τ local are exactly the same.This simple test case proves that the DMD can be used for the computation of the fields of n local and τ local . Results are relevant for a snapshot database that contains a few periods. The frequency does not need to be known a priori and is a result of the DMD analysis.

Figure B. 2 -

 2 Figure B.2 -Fields of n local at 1000 Hz computed with the Fourier transform method (top) and the DMD (bottom).

Figure B. 3 -

 3 Figure B.3 -Fields of τ local at 1000 Hz computed with the Fourier transform method (top) and the DMD (bottom).

Figure C. 1 -

 1 Figure C.1 -Schematic representation of a dual cell (in grey) surrounding a node j.

Figure C. 2 -

 2 Figure C.2 -Schematic representation of a generic element Ω ei .

Figure C. 3 -

 3 Figure C.3 -Computation of ∇p Ω ei over a generic element Ω ei . The values of the pressure field is stocked at the nodes j, i and i + 1.

Figure C. 4 -

 4 Figure C.4 -Schematic representation of a dual cell (in grey) surrounding a node j located on a boundary.

jω Z j pj (l 13 +

 13 l 33 ) (C.10) Note that Eq. (C.10) is stocked into the diagonal matrix B in Eq. (2.56). For the Delayed Entropy Coupled Boundary Condition (DECBC), Eq. (C.10) takes the following form:1 ρ Ω e1 ∇p Ω e1 ⋅ l 13 n -1 ρ Ω e3 ∇p Ω e3 ⋅ l 33 n = 1 ρ j c j jω Z j pj + R SA Z j -1 Z j G us G c e jω(τus+τc) ∇p ref ρ j ref (l 13 + l 33 ) (C.11)

  Symbols

		c	speed of sound			[m s]
		f	frequency			[Hz]
		g	any quantity			[-]
		h	specific enthalpy			[J kg]
		j	standard imaginary unit			[-]
		k	complex wave number			[1 m]
		l	location of the section change			[m]
		l z	location of impedances			[m]
		m	mass flux			kg (s m 2 )
		n	local interaction index (flame model)			J m 4
		n	normal unit vector			[-]
		p	pressure			[Pa]
		q	heat release rate per unit volume			J (s m 3 )
		r	specific gas constant			[J (kg K)]
		s	entropy			[J (kg K)]
		t	time			[s]
		t	tangential unit vector			[-]
		u	velocity (1D)			[m s]
	A	u	complex amplitude of acoustic wave velocity vector	[Pa]	[m s]
	C s	x	constant of sub-grid scale Smagorinsky model coordinates of a local point 1D	[-]		[m]
	C σ	x	constant of sub-grid scale Sigma model coordinates of a local point	[-]		[m]
	C v		heat capacity at constant volume	[J (kg K)]
	C p E Greek Symbols heat capacity at constant pressure acoustic energy density	[J (kg K)] J m 3
	E		flame efficiency function	[-]
	F γ		acoustic energy flux adiabatic index (heat capacity ratio)	W m 2 [-]
	F φ, α	flame thickness factor phase	[-]	[rad]
	G δ f		gain (DECBC model) flame length	[J s (kg K m)] [m]
	J ∆		total enthalpy LES filter size	[J]	[m]
	L ∆F	location of the end of the tube energy flux difference	[m]	W m 2
	M λ		Mach number wavelength	[-]	[m]
	P ω		energy source term complex angular frequency	J m 3 [rad s]
	P r ωT		Prandtl number reaction rate	[-]		J (s m 3 )
	Q ωk		heat release rate mass reaction rate of species k	[J s] kg (s m 3 )
	R S ρ τ		complex reflection coefficient surface density time delay (flame model)	[-] m 2 [s] kg m 3
	S ij S τ i,j µ		strain rate tensor viscous stress tensor Thickened Flame Model sensor viscosity	[1 s] N m 2 [-] [kg (m s)]
	T σ		temperature complex amplitude of a fluctuating entropy wave	[K] [J (kg K)]
	T V ξ		period of oscillation volume curvilinear coordinates system of a combustion chamber in 3D [m] [s] m 3
	V c		corrected diffusion velocity	[m s]
	W		mean molecular wight	[kg mol]
	X		species molar fraction	[-]
	Y		species mass fraction	[-]
	Z		complex impedance	(N s) m 3

  At last, basic experiments studied in labs are usually operated at atmospheric pressure[START_REF] Lawn | A model for the thermoacoustic response of a premixed swirl burner, part i: Acoustic aspects[END_REF][START_REF] Kornilov | Experimental and numerical investigation of the acoustic response of multi-slit bunsen burners[END_REF][START_REF] Karimi | Linear and non-linear forced response of a conical, ducted, laminar premixed flame[END_REF][START_REF] Kim | Response of partially premixed flames to acoustic velocity and equivalence ratio perturbations[END_REF]

	•
	;
	Cumpsty & Marble, 1977; Bake et al., 2009; Howe, 2010; Miles, 2010; Leyko et al., 2011;
	Duran et al., 2013). However, the contribution of this mechanism to the unstable be-
	haviour of air-flames remains controversial.
	• One reason is that purely thermo-acoustic arguments usually allow the observed
	instabilities to be explained very nicely in most cases (McManus et al., 1993; Candel
	et al., 2009).
	• Another reason is that many interactions between chemistry, hydrodynamic and
	acoustics (Ducruix et al., 2003; DeGoey et al., 2011) are involved when dealing
	with combustors equipped with multi-staged, co-rotating swirlers (Wang et al., 2007;
	Huang & Yang, 2009) and/or effusion-cooled walls (Hughes & Dowling, 1990; Mendez
	& Eldredge, 2009; Scarpato et al., 2012) that couple the combustion chamber to a
	plenum and a casing. Identifying the entropy noise generation as the mechanism
	which feeds the acoustic field and leads to an unstable regime is thus challenging.

  ; Dowling & Stow, 2003; You et al., 2005) rely on the resolution of a wave equation governing the oscillatory field in the chamber. Complemented by experiments (Hield et al., 2009; Yu et al., 2010; Sisco et al., 2011) these low order models show that a new family of instability modes that depends on the mean flow velocity may exist. However, as pointed out by Eckstein et al.

Table 1 .

 1 2 -Summary of the analytical procedures (methods I, II, III) investigated to compute the two connected tubes configuration.

  1.56) -(1.57) and represent the propagation/reflection of acoustics waves through/from the LEE domain. The values of the impedances and reflection coefficients for some canonical boundary conditions are recalled in Table 1.3.

  54) where û (x ref ) ⋅ n ref is the fluctuating acoustic velocity taken at a reference point x ref and in direction n ref . Reformulation of Eq. (2.51) leads to the following Delayed Entropy Coupled Boundary Condition (DECBC):

  1 n crocco e jωτcrocco = V f n local e jωτ local dx (B.10) = S f δ f n local e jωτ local (B.11) Hence, it comes the following expressions for n local and τ local :

	n local =	γP 0 γ -1	1 δ f	n Crocco	(B.12)

x f

Acknowledgements Remerciements

Acknowledgements Merci aussi pour m'avoir permis d'aller plusieurs fois dans des conférences internationales. D'ailleurs je ne sais comment exprimer ma gratitude à Stéphane Moreau pour avoir assuré ma présentation un lundi matin 9h heures (et en ayant reçu les slides 20 minutes auparavant), alors que j'étais coincé à l'aéroport de New-York. Par contre je ne te remercie pas pour la soirée au Miyake, ce n'est pas un endroit pour acousticien ça !

Application of the DECBC approach

Recall that the configuration investigated in the present work exhibits a low frequency mode of oscillation at a frequency significantly smaller than any acoustic mode. Helmholtz analysis performed in Sec. 3.4.2 failed to capture this low frequency mode. It is thus natural to investigate whether the DECBC approach described in chapter 2 helps to recover this latter mode. The transfer function between the acoustic velocity upstream of the flame and the entropy generated downstream of the primary zone was assessed by post-processing the LES via the DMD analysis. More precisely, this entropy transfer function was defined as:

where n ref is a unitary vector of reference aligned with the main axis of the combustor, ûref is the acoustic velocity at the reference point depicted in Fig. 3.19 and < ŝf > is the entropy fluctuation averaged over a small volume located downstream of the primary combustion zone, in agreement with the mode structure observed from the DMD analysis of Sec. 3.3.2. Obtaining the acoustic velocity fluctuations ûref is challenging because the reference point is located inside the swirler where hydrodynamic fluctuations are present. However, as the present study focuses on a low-frequency instability, the hydrodynamic component is considered negligible. Note that this transfer function is similar to but different from the classical flame transfer function which relates the upstream acoustic velocity to the unsteady heat release thanks to a nτ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF][START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF].

National Conferences

Isentropic relations for quasi-one-dimensional flow

Following the pedagogy of [START_REF] Shapiro | The Dynamics and Thermodynamics of Compressible Fluid Flow[END_REF] or [START_REF] Anderson | Modern Compressible Flow[END_REF], before any derivation of equations it is convenient to introduce the concept of stagnation and critical states. If one consider a point in an arbitrary flow field, this point is moving with a certain velocity and Mach number M , and carry a certain static pressure p and temperature T . Now if one imagine that the velocity is slow down to zero, the fluid will stagnate, and for this particular state the static pressure and temperature will virtually change to their stagnation values p 0 and T 0 , respectively. In the same manner, is the flow is accelerated to M * = 1 from its initial state, the flow properties will take the critical values T * and p * . During this virtual change of the flow velocity, the process is assumed isentropic, i.e. the flow is both adiabatic and reversible. Consequently, no heat is added to the flow and there are no friction of dissipative effects.

A.1 Stagnation relations

Starting from the first law of the thermodynamic and the definition of the entropy (see [START_REF] Shapiro | The Dynamics and Thermodynamics of Compressible Fluid Flow[END_REF][START_REF] Anderson | Modern Compressible Flow[END_REF], one may write the following relation that express the difference between two entropy states:

The isentropic assumption assumes that ss ref = 0. From Eq. (A.1) it then comes the following useful isentropic relations:

Isentropic relations for quasi-one-dimensional flow

In the same manner, any change in pressure dp is accompanied by a corresponding isentropic change in density dρ, so that:

where c = γp ρ is the sound speed.

All states along a streamline lie on a line of constant entropy and have the same stagnation temperature T 0 . The static temperature T at any point in the flow is related to the stagnation temperature T 0 by the following relation that only involves the gas constant γ and the Mach number M = u c:

In the same way, introducing the relations of Eq. (A.2) into Eq. (A.4) leads to following relations between the stagnation and static states for the variables p and ρ:

A.2 Critical Mach number and area ratio

In the case of a quasi-one-dimensional flow, the cross-section of the streamtube that contains the flow may vary with the distance. From the law of the mass conservation formulated at Eq. (1.3), it comes the following relation:

where A is the cross-section area, ρ the density and u the flow velocity. Differentiation of this latter equation and introduction of Eq. (A.3) leads to the useful area-velocity relation:

Different states depending on the Mach number M can be deduced from Eq. (A.8):

• Incompressible flow M → 0: is the limit corresponding to incompressible flow. Eq. (A.8) reduces to Au = const.

• Subsonic flow 0 ≤ M < 1: a decrease in velocity is associated with a decrease in area. This means that the velocity increases in a converging duct and decreases in a diverging duct.

A.3 Relations at a section change

• Supersonic flow M > 1: an increase in velocity is associated with an increase of the area, and vice versa. Contrary to subsonic flow, the velocity increases in a diverging duct and decreases in a converging duct.

• Sonic flow M = 1: Eq. (A.8) yields dA A = 0. This mathematically corresponds to two possible cross-section area, but only the minimum value has a physical meaning: this is the sonic throat.

The sonic throat is interesting because it corresponds to the critical state M * and its associated critical area A * . From the mass conservation equation, this means that ρuA = ρ * u * A * . Consequently, it is possible to define a relation between the cross-section area and the Mach number:

The Mach number at any location of the duct is a function of the ratio of the local duct area to the sonic throat.

A.3 Relations at a section change

If a stream tube has a section change defined by an area ratio A 2 A 1 , all flow properties in the section 2 can be deduced from the properties of the section 1. From the isentropic assumption, critical and stagnation properties do not vary during the section change, so that for example:

Finally, one can express the direct relations between state variables in section 2 from their values in section 1. These relations are directly used in Sec. 1.4 for the computation of the acoustic eigenmodes in the configuration where two tubes with different section area are connected each other.

(A.12)

The Mach number M 2 in the section 2 can not directly be expressed as a function of M 1 . From Eq. (A.9) and using the relation

, it comes the following reformulation:

The Mach number M 2 can be found by solving the roots of Eq. (A.14) with the help of a Newton-Raphson method [START_REF] Moler | Numerical computing with MATLAB[END_REF]. Two peculiarities must be noted:

• For a given Mach number M 1 and a given area ratio A 2 A 1 , there are either two solutions for the final state of M 2 : one is subsonic, the other is supersonic. In the particular case studied in §1.4, there are no sonic throat and the flow is subsonic everywhere.

• When A 2 is smaller than A 1 , for certain values of M 1 there are no solutions for M 2 . Physically this signifies that for given conditions at section 1, there is a maximum contraction which is possible and corresponds to the sonic velocity at section 2. Hence in the subsonic flow there is a maximum initial Mach number that can be maintained. This phenomena is called choking.

Appendix B

About the active flame model and the evaluation of n localτ local fields

The AVSP code solves numerically as an eigenvalue problem the inhomogeneous Helmholtz equation described at Eq. (2.41) and reproduced here as follow:

where the sound speed c, the mass density ρ and the heat capacity ratio γ are local variables. In order to close the problem, the flame that constitutes a forcing source is often modelled as a purely acoustic element thanks to a nτ type of model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF]:

where n ref and x ref are the normal and the reference point, while û is the fluctuating velocity vector.

The fields n local [J/m 4 ] and τ local [s] are dependent on space and frequency. The modelling of the flame is the most important part as it directly has an influence on the results, especially on the stability of a thermoacoustic eigenmode. Most of the models rely on a global description of the flame, connecting the global rate of heat release to the fluctuations of velocity in the fresh gas upstream of the flame:

In the latter expression, the term Q [W] represents the fluctuations of heat release rate integrated over the total volume, n g and τ g are parameters representative of the amplitude and time delay of the flame response. Their unit are [J/m] and [s], respectively.

Note that the main difficulty is to compute the complex values of q(x) and û(x ref ). This can be done by performing a Fourier transform of the recorded signals of q ′ (x, t) and u ′ (x ref , t) at the frequency of interest. However the use of a Fourier transform rises two problems:

• the total length time of the recorded signals needs to be significant so as to provide enough unstable cycles. When a flame transfer function is sought for a low-frequency phenomena, this can lead to a high computational cost.

• the frequency at which the flame transfer function is computed needs to be known a priori. This is not a problem in the case of a flame that reacts linearly to an acoustic forcing, but in the case of an auto-sustained unstable flame the frequency is not necessarily known.

An alternative to the Fourier transform technique is to use the Dynamic Mode Decomposition (see Sec. 3.3). The DMD addresses the two problems cited above, because this method gives relevant results even when only a few period are considered (see [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF], and the frequency is a result and does not need to be known a priori.

Discretisation procedures

This appendix is devoted to the details of the discretisation of the operator ∇ ⋅ (1 ρ ∇) applied to any node of the mesh denoted j. The AVSP code relies on a volume finite formulation and the discrete values of the pressure field are stocked at the nodes of the mesh. The pressure gradients are computed at the center of each element. Two situations are considered:

• the node j is located inside the domain;

• the node j is located at a boundary of the domain.

Note that only the case in 2D is considered, the extension to 3D follows the same procedures.

C.1 Evaluation of ∇⋅(1 ρ ∇) p j for a node j inside the domain As depicted in Fig. C.1 the node j is surrounded by 6 elements Ω ei , with i = 1 . . . 6 the number of the element. A dual cell Ω N j with an associated volume S N j is virtually created by connecting the centroid of each element Ω ei , the node j being the center of Ω N j . In the element Ω ei the normals n i1 and n i2 are facing outward while l i1 and l i2 are the lengths of the vertices connecting the centroid of Ω ei and the midpoints. As the values of the pressure are stocked at the nodes of each elements, it comes the following relation: Then the pressure gradient takes the expression:

Note that the values of ∇ ⋅ (1 ρ ∇) p j are stocked in the matrix A in Eq. (2.56).