Romain Garby
email: romain.garby@imft.fr

Etude

M Sébastien

Ducruix M Mauro

Valorani M Thierry

Poinsot M Laurent

M Gérard

NICOUD Ordonneau M Franck

Franck Nicoud

Hélène, Catherine, Véronique Alain Thomas

Gilles Dominique

Benoît, Eric, Franck Jean-Luc Pascal

Moïse Rudy Bernard

Nico Florian

Zafer, Marco, mais aussi Viviane Enrica Ali

T A Thermoacoustic

Etude numérique de la stabilisation de flamme et des instabilités de combustions dans les systèmes de propulsion

Keywords: combustion instability, Large-Eddy-Simulation (LES), Helmholtz solver, rocket engine, nozzle, p.19 DTF Dynamically Thickened Flame, p.42

teaching and research institutions in France or abroad, or from public or private research centers.

5

Résumé

Cette thèse se focalise sur la compréhension et la prédiction des instabilités de combustion dans les systèmes à haute pression. Elle s'oriente autour de la simulation numérique d'un banc d'essai, opéré à l'université de Purdue, comprenant un injecteur caractéristique des moteurs-fusées et dont les propriétés acoustiques peuvent varier à l'aide d'un tube d'injection mobile. Une méthode d'initialisation et d'allumage pour les calculs LES de chambres de combustions terminées par une tuyère est présentée. Un point de fonctionnement instable est choisi pour étudier le mécanisme de l'instabilité. Les simulations sont comparées aux résultats expérimentaux en terme de fréquence et structure du mode instable. La fonction de transfert de flamme est calculée à l'aide du modèle nτ puis implémentée dans un solveur acoustique (ne résolvant que les perturbations acoustiques à partir de l'équation de Helmholtz en écoulement réactif). Différents modèles d'impédance de tuyère, extraits de la littérature, sont comparés et leurs impacts sur les résultats de stabilité sont analysés. Le théorème d'impédance translatée est implémenté dans le solveur acoustique pour analyser, à faible coût de calcul, l'influence de la variation de la longueur du tube d'injection. Des écarts entre les fréquences prédites et celles trouvées expérimentalement subsistent mais la carte de stabilité de l'expérience est bien reproduite.

Mots clés : instabilité de combustion, Simulations aux Grandes Echelles (LES), solveur de Helmholtz, moteur-fusée, tuyère.

Nomenclature

Introduction

The framework of this thesis is combustion instabilities in high-speed flows and more specifically in combustors terminated by a choked nozzle. This class of combustors is found in multiple systems used for propulsion (helicopter [1], aircraft [2,3,4,5,6], rockets [7]) or energy production [8,9,10,11,12]. Choked nozzles are commonly used in all these combustors (Fig. 1) to increase the operating pressure and the efficiency. It is interesting to note that most combustors developed in laboratories on the other hand operate at atmospheric pressure using an open ended outlet (Fig. 2). Combustion instabilities are observed in all combustors and they represent a major danger in multiple present programs because they cannot be predicted at the design stage and are very difficult to eliminate at later times [9,14,15,16]. We will not describe the general phenomenology of combustion instabilities here: this has been done in numerous recent PhDs focusing on this question [17,18,19,20,21,22]. In the present work, we will focus on the specificities of high-speed flows in choked combustors which constitute a subset of unstable combustors but an important one in terms of applications even though they have not been studied as much as atmospheric burners.

There are multiple reasons for the fact that laboratories prefer using atmospheric experiments: the cost of experiments varies typically like pressure to a large exponent, diagnostics become difficult, flow rates are much higher. In general the fundamental understanding gained on mechanisms is often not very different at atmospheric or at high pressure. This is not always the case however and for certain problems, operating at high pressure and using a choked nozzle becomes mandatory because the mechanisms at play are simply too different with and without choked nozzles. Combustion instabilities are one of these examples.

When a combustor becomes unstable, the mean pressure at which it operates is not a major parameter but the impedances which terminate the combustor (inlet and outlet) are [14]. We will discuss nozzle impedances in more details in this work but to first order it is interesting to compare the two dump combustors sketched in Fig. 3. The first one (left) is typical of an academic experimental setup and its outlet is directly blowing into free atmosphere. The corresponding impedance (defined by the normalized ratio of pressure (P ) to velocity (u ) perturbations: Z = P /(ρcu ) ) is zero since pressure is essentially imposed by the atmosphere. On the other hand, in a combustor terminated by a nozzle (Fig. 2 right), the nozzle condition controls the impedance.

The impedances of nozzles have been studied for a long time [23,24,25]. To first order, at low frequencies and when the Mach number upstream of the nozzle is low, the nozzle essentially behaves like a rigid wall: u = 0 and Z = ∞. As shown in Fig. 3, this means that for an atmospheric burner, the first acoustic mode which can be forced by combustion is the quarterwave mode of the combustion chamber while it is the half-wave mode for the choked combustor. An additional complexity introduced by a choked outlet is the fact that the mean pressure in the chamber and the outlet flow rate are continuously linked by the unity Mach number condition at the throat. When the combustion rate in the chamber oscillates, however, this operating condition also changes because the stagnation temperature varies, introducing an additional constraint on the unstable mode. In certain cases, the pressure in the combustor can actually become low enough for the throat to become subsonic again, opening the path for an additional instability mechanism where the throat would oscillate between choked and unchoked conditions.

This implies that the two combustors of Fig. 3 behave very differently from an acoustic point of view and that studying combustion instabilities in an atmospheric system may be an irrelevant model when considering choked systems. This is a disturbing observation since most of the academic knowledge gained on combustion instabilities in the last fifty years has been obtained on atmospheric systems for obvious reasons of cost and simplicity.

The importance of impedances when studying instabilities has become clear in multiple fields recently: it was obvious for rocket experts who have been studying choked nozzles impedances for a long time [14,26]. It has become clearer for the gas turbine community too when they realized that a given combustion chamber, unstable in the real engine (where the chamber is choked), could often become stable when it was operated outside the engine in a test bench without outlet nozzle. These observations are not surprising for acoustic experts who fully expect to see a determinant effect of impedances on instabilities but they also demonstrate the need for studies dedicated to the behavior of combustors terminated by choked nozzles (and therefore necessarily operating at high pressure). Performing such a study numerically is the first objective of this work.

Finally a major difference between atmospheric burners and choked combustors is the level of pressure oscillations. In atmospheric burners, because of the outlet condition which imposes an almost zero pressure oscillation, the typical levels of oscillations remain low and rarely reach more than 0.1 bar (which already corresponds to a very noisy output). In a choked combustor however, the outlet losses are much smaller, the mean pressure is higher and pressure oscillations of the order of a few bars are commonly observed. Obviously, such levels of oscillations generate mechanisms which can be quite different from those observed at atmospheric pressure. The most Introduction obvious one is the non linearity of all mechanisms involved in the unstable loop.

The goal of this study is to address these issues on one specific example and it is based on three main points:

• A typical target experiment was chosen, corresponding to a recent set of data obtained at Purdue University [27,28], on a system called CVRC which has two interesting characteristics: (1) as required, it operates at high pressure, is terminated by a choked nozzle but (2) in addition, the inlet length can be continuously adjusted, allowing to explore a wide range of acoustic characteristics and therefore providing interesting data on the chamber stability. Adjustable lengths are useful in combustion instability studies [29] because they allow to change the acoustic eigenmodes of the chamber continuously and to excite the flame over a wide range of acoustic conditions with a reduced modification of the geometry.

• For this target experiment, a Large Eddy Simulation (LES) numerical experiment was developed, based on the AVBP code of CERFACS [30,31,32,33,34]. This work has required to develop specific chemical schemes adapted to the combustion of CH 4 with O 2 and H 2 O as used in the CVRC and to adapt the models to high pressure conditions. However, the first originality of this part is the methodology used to compute choked reacting and non-reacting flows. A method based on the NSCBC technique [35,15] was used to handle outlet flow conditions in a combustor which can be either choked or not. The method also allows to transition smoothly from subsonic conditions at the outlet throat to choked conditions. LES was used to study one specific operating point with fixed inlet length to analyze the steady flow and its unstable modes. An additional issue identified during this work is the importance of heat losses. While most atmospheric burners operate almost in adiabatic conditions, choked combustors generate much more heat and losses at walls become important. In the CVRC experiment, a simple preliminary model (based on the operating conditions of the nozzle) will show that one third of the power created by combustion is lost to the walls demonstrating that a wall heat loss model should be integrated in the LES to mimic the operating conditions of the combustor (Sec. 4.6). As shown by previous LES [32,36], introducing heat losses not only changes the mean flow in the combustion chamber but also the instabilities. We will study these mechanisms in details.

• Even though LES is a powerful investigation method, it is also expensive. It is well adapted to the investigation of one case for fixed geometry and operating conditions, but not to studies of the effects of variable inlet length. Investigating all possible lengths for the post injection system used in the CVRC experiment would require too many simulations. Instead of following this path, the last part of the thesis uses a different numerical tool to address this issue: an acoustic solver called AVSP [37,[START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF] was adapted to the CVRC configuration and Flame Transfer Functions measured in the LES were used to introduce the flame effects in the acoustic equations. When the inlet length is varied, this flame response does not change so that the effects of variable inlet length are correctly captured by the acoustic solver. This allows to study the effects of variable inlet lengths on the stability map. In the CVRC, unstable modes are observed for intermediate values of the inlet post length while the burner is stable for short or long post lengths. Hysteresis is also observed. AVSP is used here to compute the frequency, the growth rate and the modal structure of all modes appearing in the CVRC for all post lengths and results are compared to the experimental data. Similarly to the study of choked nozzle which proved necessary for LES, the acoustic study with AVSP has required to address the specific question of nozzle impedances because they control the eigenmodes. The impedance of nozzle can be computed using several approaches, ranging from compact approximation methods proposed by Marble and Candel or Zinn [23,24] to more recent analytical or numerical approaches proposed by I. Durán at CERFACS [START_REF] Duran | Solution of the quasi one-dimensional linearized Euler equations using flow invariants and the Magnus expansion[END_REF][START_REF] Duran | Analytical and numerical study of direct and indirect combustion noise through a subsonic nozzle[END_REF] or A. Giauque at ONERA. These methods will be presented and their application to the CVRC nozzle discussed. Their impact on the results of AVSP will also be discussed. A major difficulty in acoustic approaches is the evaluation of the acoustic losses which are created at the inlet, the outlet but also the walls and within the domain. The discussion of which energy should be used to characterize the acoustic activity in a combustor and which terms should be used to measure losses remains controversial [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][START_REF] Martin | Analysis of acoustic energy and modes in a turbulent swirled combustor[END_REF][START_REF] Giauque | Budget of disturbance energy in gaseous reacting flows[END_REF] and was not pursued here: the usual Rayleigh criterion [15] will be used to study stability. Another difficulty for AVSP is that this solver works at zero Mach number while the Mach number in the CVRC is not small. This effect can be compensated by manipulating the inlet and outlet impedances as shown by Motheau et al [START_REF] Motheau | Using boundary conditions to account for mean flow effects in a zero mach number acoustic solver[END_REF] and we will also test these methods here.

The thesis is organized as follows: Chapter 1 first describes the experimental setup which will be simulated in this study, Chapter 2 describes the numerical tools, Chapter 3 the numerical setup for this particular application, Chapter 4 presents the results for a fixed unstable case and finally Chapter 5 presents the results for variable length configuration.

Introduction

Chapter 1

Description of the experimental setup

In order to discriminate the capacity of numerical simulation to predict combustion instability, an experimental setup which exhibits high-pressure oscillations in a high-speed flow (Mach 0.4 in the combustion chamber and 1 at the exit nozzle) is chosen as a test case. This burner has a typical rocket-engine-injector geometry. A variable inlet length allows to explore a wide range of acoustic characteristics and identify self excited modes by tuning the frequency of the acoustic modes to the chamber eigenmodes.

The following sections present its design (Sec. 1.1.1) and operating point (Sec. 1.1.2). Then a zero-dimensional model is derived to relate mean combustion parameters (Sec. 1.3): in a choked chamber, flow rates, mean pressure and mean heat losses can be linked by a simple model which will be useful to analyze the full LES of the following chapters. The stability map of this experiment is presented in Sec. 1.4 and finally Sec. 1.5 describes recent numerical studies of this device.

Configuration

Description

The burner studied here was designed and operated at Purdue University (Indiana, USA) [28,27]. The system is called Continuously Variable Resonance Chamber (CVRC) because the length of the oxidizer injector can be varied continuously, allowing for a detailed investigation of the coupling between acoustics and combustion in the chamber. The CVRC, presented in Fig. 1.1 and 1.2, is an axisymmetric burner with a L = 15-inch long (38.1 cm) combustion chamber closed by a choked nozzle. The reactants are fed into the chamber by a coaxial injector: oxidizer at the center and fuel at the periphery. For a continuous tuning of the acoustic response of the CVRC, the oxidizer-post length (L op ) can vary from L op = 3.5 in. (8.89 cm) to L op = 7.5 in (19.05 cm), with respect to the inlet of the combustion chamber. The inlet of the oxidizer post is also choked so that the acoustic boundary conditions are well defined and will correspond in the LES to fixed flow rates. Six pressure transducers are installed along the burner to measure pressure fluctuations: the first one is located 3 in upstream the injector face, and the others respectively at 0.5, 1.5, 2.5, 3.5 and 14.5 in downstream the injector face. These locations are respectively called P1 to P6 in the remainder of the thesis.

Operating point

The fuel is pure gaseous methane injected at a rate of ṁf = 0.0272 kg.s -1 and temperature T f = 300 K. The oxidizer is a mixture of Y H 2 O 2 = 90 % hydrogen peroxide (H 2 O 2 ) and Y H 2 O = 10 % water (where Y k is the mass fraction of the species k).

The hydrogen peroxide is fully decomposed following the reaction:

H 2 O 2 -→ 1 2 O 2 + H 2 O (1.1)
For an oxidizer which is a mixture of hydrogen peroxide and water, this leads to the following decomposition reaction for the mixture:

X 0 H 2 O 2 H 2 O 2 + X 0 H 2 O H 2 O -→ X 0 H 2 O 2 2 O 2 + X 0 H 2 O 2 + X 0 H 2 O H 2 O (1.2)
where X k is the molar fraction of the species k and the exponent 0 refers to the initial state.

Relation between mass and molar fraction is:

Y k = W k W X k and W = 1 Y k W k = X k W k (1.3)
where W is the global molar mass of the mixture and W k is the molar mass of the specie k.

Thus the molar fraction of the decomposed oxidizer is:

X 1 O 2 = 1 2 1 + X 0 H 2 O X 0 H 2 O 2 + 1 2
and X 1

H 2 O = 1 + X 0 H 2 O X 0 H 2 O 2 1 + X 0 H 2 O X 0 H 2 O 2 + 1 2 (1.4)
this leads to the fully decomposed oxidizer composition :

Y 1 O 2 = 42 % and Y 1 H 2 O = 58 %.
The oxidizer is injected in the chamber at T ox = 1030 K (both water and oxygen are in gaseous phase) and ṁox = 0.3180 kg.s -1 .

The single-step global reaction is written:

CH 4 + 2O 2 + H 2 O ←→ CO 2 + 3H 2 O (1.5)
For this reaction the stoichiometric ratio is

s = 2W O 2 /W CH 4 = 2 × 32/16 = 4.
Consequently, the global equivalence ratio is φ = s ṁf / ṁox = 0.813 because the fuel and the oxidizer flow rates are fixed for all experiments in the CVRC.

Computation of adiabatic flame temperature

Fresh gases are composed of fully decomposed hydrogen peroxide (basically oxygen and water) injected at T 0 ox = 1030 K and methane at T 0 CH 4 = 300 K (Sec. 1.1.2). It is useful to compute the temperature of a mixture of hydrogen peroxide and methane before and after complete combustion. To do this, the properties of the mixture of fuel and oxidizer have to be calculated versus mixture fraction z (Fig. 1.3). To find the temperature of the mixture T 0 mix the enthalpy of the mixture is calculated as a composition of the enthalpies of the different components at their temperature h k T 0 k :

h mix T 0 mix = Y CH 4 h CH 4 T 0 CH 4 + Y O 2 h O 2 T 0 O 2 + Y H 2 O h H 2 O T 0 H 2 O (1.6)
where h k is tabulated for each species for temperature from 0 to 5000 K. Then T 0 mix is found searching in the enthalpy tables the temperature verifying:

h mix T 0 mix = Y CH 4 h CH 4 T 0 mix + Y O 2 h O 2 T 0 mix + Y H 2 O h H 2 O T 0 mix (1.7)
The mass fractions Y CH 4 , Y O 2 and Y H 2 O are given by the mixture fraction z as follows (mixing lines):

Y CH 4 = z (1.8) Y O 2 = Y 0 O 2 (1 -z) (1.9) Y H 2 O = Y 0 H 2 O (1 -z) (1.10)
Two values of z are interesting: the average z of the burner z xp which corresponds to the mean outlet temperature, and the stoichiometric value z st which corresponds to the maximum temperature reached in the burner. Tab. 1.1 presents the fresh-gases properties at the stoichiometric mixture fraction z st = 0.095 as well as the average mixture fraction of the experimental operating point z xp = 0.079 (calculated from the mass flow rates ṁf and ṁox , z xp = s ṁf / ṁox ). A complete mechanism is used (GRI-Mech 3.0 [START_REF] Frenklach | GRI-mech: an optimized detailed chemical reaction mechanism for methane combustion[END_REF]) which account for all species (52 species), to compute the adiabatic burnt-gases temperature T ad . Eq. ratio

Eq.ratio z T 0 mix [K] Y CH 4 Y O 2 Y H 2 O Y CO 2 Y CO Y N 2 1 z st = 0.
z T ad b [K] Y CH 4 Y O 2 Y H 2 O Y CO 2 Y CO Y N 2 1
z st = 0.095 2661 < 10 -4 2.49.10 -2 7.12.10 -1 2.22.10 -1 2.45.10 -2 0 0.813 z xp = 0.079 2512 < 10 -4 7.50.10 -2 6.99.10 -1 2.06.10 -1 6.17.10 -3 0 Table 1.2: Burnt-gases properties in thermodynamic conditions of the CVRC calculated with GRI-Mech 3.0 mechanism.

Zero-dimensional model for mean combustion parameters 1.3.1 Prediction of mean chamber pressure

Before performing any detailed CFD of the experimental setup, it is useful to do a simple analysis of the performances of the burner. Since the outlet nozzle is choked, the mean operating pressure P c depends on the mass flow rate ṁ, the combustion efficiency η and the heat losses f HL :

P c = f ( ṁ, η, f HL ) (1.11)
This functional dependance can be obtained by first computing the stagnation pressure P 0 and temperature T 0 in the chamber (Sec. 1.2) and then using the chocked nozzle relation (Sec. 1.3.1). Since the total flow rate going through the choked nozzle is known ( ṁ = ṁox + ṁf = 1.3 Zero-dimensional model for mean combustion parameters 0.3452 kg.s -1 ) assuming that the chamber temperature, T c , is equal to the adiabatic burnt-gases temperature at equilibrium T ad b , the chamber pressure P c can be determined using:

ṁ = c d AP 0 √ rT 0 γ 1 2 2 γ + 1 γ+1 2(γ-1) (1.12)
where the stagnation pressure and temperature, respectively P 0 and T 0 are given by:

P 0 = P c 1 + γ -1 2 M 2 γ γ-1 ; T 0 = T c 1 + γ -1 2 M 2 (1.13)
where M is the Mach number of the burnt gases in the chamber, A the area of the choked nozzle, γ the ratio of heat capacities, r the perfect gas constant and c d the discharge coefficient.

Using

c d = 1, T c = T ad b , γ = 1.
182 and r = 389.2 (computed from the burnt-gases composition at stoichiometry), one predicts a theoretical mean chamber pressure P th c = 1.543 MPa. The mean pressure actually measured in the experiment is P xp c = 1.34 MPa, which implies that either heat losses or incomplete combustion lead to a lower burnt-gases temperature: T c < T ad b . Using Eq. 1.12 backwards in order to obtain P c = P xp c , yields T xp c = 1904 K. Because the simulations presented in chapter. 3 show that combustion is complete (there is no fuel leakage at the nozzle), this reduced temperature is attributed to thermal losses at the combustor walls. Since the regime is lean, the total power of the burner can be evaluated directly using the fuel flow rate ṁf and the mass heat of the reaction between oxygen and methane Q = 50, 1 MJ/kg, which leads to P b = ṁf Q = 1.36 MW. Thus, a first estimation of heat losses flux through the wall f HL which must be imposed to have P c = P xp c using the burnt-gases heat capacity at T ad (C p = 2529 J/kg/K) :

f HL = ṁC p T ad -T xp c = 530 kW (1.14)
which represents 39% of the total power of the burner. This is a large value compared to the global power of 1.36 MW : 39% of the heat released by combustion is lost to the walls.

Implication for LES modeling

The previous section has revealed that heat losses were significant for the operating point we plan to compute for the CVRC. This has a direct and important implication for the LES: the simulations should not be performed using adiabatic assumptions. Heat losses at walls should be included. Furthermore it also implies that the chemical kinetics model which will be used in the LES can not be based on the adiabatic assumption, thereby eliminating the possibility of a simple tabulation based on the mixture fraction z and the progress variable c [START_REF] Fiorina | Approximating the chemical structure of partially-premixed and diffusion counter-flow flames using FPI flamelet tabulation[END_REF][START_REF] Van Oijen | Modeling of premixed laminar flames using flamelet generated manifolds[END_REF]. Extensions of these FPI or FGM techniques to situations where heat losses must be included are possible but require a tabulation using three parameters (z, c, h) where h is the total enthalpy. Since the details of chemistry (emissions levels for example) are not the topic of this work, a different route will be chosen here, based on reduced chemistry as described in Sec. 2.1.2.

A second implication of the previous section is that wall treatments will be important. In many LES of reacting flows, especially at atmospheric pressure, crude models can be used at walls because they have a negligible influence on the flow. Here, since the walls will control heat losses, they will also control the mean pressure and overall organization of the flow. Therefore, they require specific attention: a law-of-the-wall, initially proposed by Schmitt et al [32] will be applied as discussed in Sec. 4.6. We will see that heat losses control both the stabilization (the place where the flame is anchored) and stability (the occurrence of combustion instabilities).

Chapter 1 : Description of the experimental setup

Experimental results on stability

The experimental stability map (amplitude of pressure fluctuations versus L op ) obtained by Pr. Anderson's team is presented in Fig. 1.5 for a number of runs where the oxidizer-post is either translated forward (L op increases), backward (L op decreases) or kept constant during the whole test [START_REF] Yu | Spontaneous longitudinal combustion instability in a continuously-variable resonance combustor[END_REF]. Three different behaviors of the system can be distinguished :

1. Stable regions for L op < 3.75 in (9.5 cm) and L op > 7.1 in (18 cm).

2. An unstable range for 3.75 < L op < 6.3 in (9.5 to 16 cm).

3. A zone of hysteresis for 6.3 < L op < 7.1 in (16 to 18 cm).

Value of L op for LES of Chapter 3 The LES simulations will focus on an unstable operating point : they are conducted for L op = 4.75 in (∼ 12 cm). This value of L op is also chosen to minimize potential hysteresis effects in the instability as all tests show strong pressure fluctuations under these conditions. The typical levels of pressure fluctuations in the CVRC is from 0.8 to 1.9 MPa. Figure 1.6 shows pressure signal versus time for a given evolution of the post length. Pressure oscillations vary from 10 to 70% of the mean pressure and the maximum pressure oscillation can reach 1 MPa (peak to peak). These very large oscillations make instabilities in choked systems specific and are rarely obtained in classical academic configurations operating at atmospheric pressure. Two operating points are studied, one fuel rich regime with a mixture ratio z rich = 0.13 (not presented here) and one lean regime z lean = 0.079 (the stoichiometric mixture ratio is z st = 0.095). They use two kinetic schemes:

Experimental results on stability

• a single step global reaction:

CH 4 + 2O 2 -→ CO 2 + 2H 2 O (1.15)
described by the Arrhenius law. This single step kinetic model has a low computational cost but is not accurate in lean conditions.

• a four-step model with six species: Comparison between 2D and 3D simulation: This comparison has also been performed by the authors of this thesis and published in [START_REF] Garby | Large-eddy simulation of combustion instabilities in a variable-length combustor[END_REF] (see Appendix A). In both studies, the 2D axi-symmetric calculations are wrong for the mean properties of the flow (such as velocity, flame length, temperature on the axis, etc,...) due to the axis condition. But both exhibit a combustion instability for which the amplitude is half of the amplitude of the 3D case or the experimental value, but with a good mode structure and an over estimated frequency which stays within 15% of the experimental one. Both authors agree on the interest of 2D simulations which are less reliable but allow very efficient preliminary studies due to their low computational cost to determine the causes of acoustic instabilities.

2CH 4 + O 2 -→ 2CO + 4H 2 (1.16) CH 4 + H 2 O -→ CO + 3H 2 (1.17) 2H 2 + O 2 -→ H 2 O (1.18) CO + H 2 O -→ CO 2 + H 2 (1.19) (1.
Chapter 1 : Description of the experimental setup

Impact of kinetic model:

The choice of a kinetic model has an effect on the location of the heat release rate. For the four-step case, the peak heat release is located at a pressure anti-node, for the single-step case the peak is further downstream, the pressure is not amplified at its peak value, leading to an overall lower level of instability. Thus the kinetic model has a drastic impact on the precision of the frequency and magnitude of the instability. Extrapolating to what a four-step threedimensional simulation might predict gives a value of around 1390 Hz which is very close to the experimentally measured frequency.

Study of the acoustic instabilities:

Time-averaged fields for the fully three dimensional simulation are presented in Fig. 1.9. The mean pressure is P mean = 1.55 MPa, much higher than the experimental value (P xp = 1.34 MPa), but consistent with the theoretical analysis, indicating complete combustion. The flame length is around 10 cm, the simulation operates in fuel-lean regime as there is oxygen exiting through the nozzle. There are strong perturbations in the first three centimeters of the oxidizer post in the Mach number field after the inlet nozzle. Instantaneous fields of pressure, temperature, Mach number and vorticity are presented in Fig. 1.10. There is a large pressure wave in the area surrounding the combustion zone. Time resolved images show that the wave is propagating upstream in the oxidizer post and downstream in the combustor. The temperature field presents multiple cold pockets convecting downstream until the nozzle. The pressure wave moving upstream in the oxidizer post is also visible in the Mach number contour. The pressure wave moving upstream in the oxidizer post is also visible in the Mach number field. The high Mach flow leaving the oxidizer slots is mixed and travels through the later two-thirds of the combustor at a relatively uniform Mach number. The area near the slots is also the source of high vorticity, which begins to dissipate downstream. Vorticity then increases in the area of the backstep where shedding takes place.

The shape of the first mode is plotted on Fig. 1.11. In the combustor, there are two pressure anti-nodes, one just after the nozzle and one at the injector face (x = 0 m). A pressure node lies in the middle of the combustor around x = 0.2 m . The velocity has two nodes, one near the injector face and one near the nozzle throat. A velocity anti-node lies in the center of the domain, aligned with the pressure anti-node. The velocity mode shape just past the injector face is distorted due to the recirculation region. The negative velocity present in the recirculation region lowers the recorded average velocity in this region giving rise to an elongated node. The dip in the pressure data just past the backstep is a result of the heat release. Chapter 2

Description of the numerical tools All the physics required to predict combustion instabilities is contained in the 'full' Navier Stokes equations. 'Full' equations indicate compressible, reacting, multi-species conservation equations including two-phase flows (for liquid fuels), radiation and heat transfer through combustor walls. Unfortunately, solving these equations in a turbulent flow is and will remain impossible for a long time. Simplifications are required. The last ten years have shown that two methods can be used to study combustion instabilities using high-fidelity large-scale simulation codes running on massively parallel systems (see Fig. 2.1). These two classes depend on the simplifications used to solve the compressible unsteady Navier-Stokes equations in a turbulent reacting flow:

1. Brute force LES (Large Eddy Simulation): recent progress in LES of reacting flows has allowed to simulate full combustion chambers, first for laboratory-scale burners [START_REF] Di Mare | Large eddy simulation of a model gas turbine combustor[END_REF]30] and recently for chambers corresponding to real industrial configurations [33,32,[START_REF] Fureby | LES of a multi-burner annular gas turbine combustor[END_REF]. In LES the equations are solved in time on a discretized space (using unstructured grids in most combustor cases). It involves only a limited set of assumptions (mainly a filtering operation of the smallest spatial scales of the flow), leading to the most precise method today but also the most computationally intensive. The image of a 'virtual burner' is appropriate here too to describe the power of this approach: a brute force LES of a 36

Chapter 2 : Description of the numerical tools combustor exhibiting CIs will provide a 'virtual unstable combustor', easier to explore than an experiment.

2. Thermoacoustics codes (TA codes) are the natural companion method for brute force LES. In most combustion instabilities, acoustics is the dominant resonant mechanism and multiple studies performed in the last sixty years (most of them inspired by the pioneering work of L. Crocco [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors[END_REF]) have shown that a proper method to decompose the physics of combustion instabilities was to use 'acoustics' codes which can track the propagation of waves in the combustor and the growth of eigenmodes [START_REF] Crighton | Modern methods in analytical acoustics[END_REF]14]. In this approach, the mean flow is frozen and the solver only tracks the acoustic modes of the system. The flames are replaced by active components (which can be compared to complex loud speakers). If the action of these active elements can be properly represented, the global stability of the combustor can be predicted. A major asset of these methods is to isolate the elements leading to combustion instabilities into different blocks (something a brute force LES cannot do): a) the acoustics of the combustor, b) the outlets and inlets impedances and c) the response of the flame which is quantified by a function called the Flame Transfer Function (FTF) describing how much unsteady heat release is produced by a flame when it is submitted to an acoustic velocity fluctuation. TA codes can operate in the frequency domain or in the time domain but in all cases, TA codes are faster than brute force LES and their information more useful to mitigate combustion instabilities. In the present thesis, both methods were used: AVBP was utilized to compute one regime with LES and AVSP was used to investigate changes due to post length (L op ) variations where the chamber response remains the same and only the acoustic properties of the system are varied. This chapter describes first the LES code AVBP and then the Helmholtz solver AVSP . For both solvers, specific tests are presented, related to the nature of the present work: for AVBP , the capacity to propagate waves correctly is tested. For AVSP , a method to investigate simply the effect of a variable post length is tested on simple cases where an analytical solution is available.

Large Eddy Simulation

The numerical code used in this work is AVBP . The reader is referred to the AVBP handbook1 and to references cited therein. AVBP [START_REF] Selle | Compressible large-eddy simulation of turbulent combustion in complex geometry on unstructured meshes[END_REF]30,[START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] uses the LES approach (Sec. 2.1.1.a) to solve the compressible Navier-Stokes equations (Sec. 2.1.1.b) on unstructured grids. The numerical scheme used here is a second-order in space and time Lax-Wendroff formulation [START_REF] Lax | Systems of conservation laws[END_REF]. Turbulence is modeled using the classical Smagorinski model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF] (Sec. 2.1.1.c). Flame / turbulence interaction is accounted for by the Dynamically Thickened Flame (DTF) model [START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]. In this approach, the flame is artificially thickened to be resolved on the LES mesh and subgrid reaction rate is recovered through an efficiency function [START_REF] Charlette | A power-law wrinkling model for LES of premixed turbulent combustion: Part I -non-dynamic formulation and initial tests[END_REF] (Sec. 2.1.1.d).

Description of the AVBP code

2.1.1.a The LES concept

Large Eddy Simulation (LES) [START_REF] Pope | Turbulent flows[END_REF][START_REF] Sagaut | Introduction à la simulation des grandes échelles[END_REF] is now well used in the research community and is slowly transfered to industry [31,[START_REF] Gourdain | High performance parallel computing of flows in complex geometries: II -applications[END_REF] to study cases where the RANS approach is insufficient. The founding idea comes from the turbulence theory of Kolmogorov [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large reynolds numbers[END_REF][START_REF] Kolmogorov | A refinement of hypothesis concerning the local structure of turbulence in a viscous incompressible fluid at high reynolds number[END_REF] which states that for a larger Reynolds number there is a scale separation between the energetic structure of the flow and the structure which disperse this energy. In the spectral plane the energy spectra have a structure in Re -5/3 between the large energetic structures and the small dissipative vortices. The idea of LES consists in applying a spatially localized time independent filter, of a given size ∆, on the governing equations. All the "large" (greater than the filter size) scales are explicitly resolved and the "small" scales (smaller than the filter size) are modeled (Fig. 2.2). Due to the filtering approach, LES allows a dynamic representation of the large scale motions whose contributions are critical in complex geometry. The LES predictions of complex turbulent flows are henceforth closer to the physics since large scale phenomena such as large vortices shedding and acoustic waves are embedded in the set of governing equations [15].

2.1.1.b Governing equations

Considering a quantity f , it is decomposed in a filtered part f and a subgrid scale part f due to the unresolved flow motion :

f = f + f (2.1)
For the compressible Navier-Stokes equations, in order to simplify the terms to be modeled, the equations are not written for the filtered variable f but for the mass-weighted Favre filtered f defined by :

ρ f = ρf (2.2)
The balance equations for Favre averaged quantities are obtained by filtering the instanta-Chapter 2 : Description of the numerical tools neous balance equations:

∂ρ u i ∂t + ∂ ∂x j (ρ u i u j ) = - ∂ ∂x j P δ ij -τ ij -τ ij t (2.3) ∂ρ E ∂t + ∂ ∂x j ρ E u j = - ∂ ∂x j u i P δ ij -τ ij + q j + q j t + ωT + Q r (2.4) ∂ρ Ỹk ∂t + ∂ ∂x j ρ Ỹk u j = - ∂ ∂x j J j,k + J j,k t + ωk (2.5)
where a repeated index implies summation over this index (Einstein's rule of summation) and the index k is reserved to refer to the k th species and does not follow the summation rule (unless specifically mentioned). The cut-off scale corresponds to the mesh size (implicit filtering). As usually done, we assume that the filter operator and the partial derivative commute [START_REF] Sagaut | Large Eddy Simulation for incompressible flows[END_REF][START_REF] Vasilyev | A general class of commutative filters for les in complex geometries[END_REF].

The flux tensor can be divided in three parts: the inviscid part, the viscous part and the subgrid scale turbulent part.

Inviscid terms:

  ρ u i + P δ ij ρ E u j + P u j δ ij ρ k u j   (2.6)
Viscous terms:

  -τ ij -u i τ ij + q j J j,k   (2.7)
Filtering the balance equations leads to unclosed quantities which need to be modeled.

Subgrid scale turbulent terms:

  -τ ij t q j t J j,k t   (2.8)
The filtered diffusion terms are [15]:

• the laminar filtered stress tensor τ ij is given by:

τ ij = 2µ S ij - 1 3 δ ij S ll (2.9) ≈ 2µ S ij - 1 3 δ ij S ll (2.10)
with:

S ij = 1 2 ∂ u i ∂x j + ∂ u j ∂x i . (2.11)
τ ij may then be written:

τ xx ≈ 2µ 3 2 ∂ u ∂x - ∂ v ∂y - ∂ w ∂z (2.
τ yz ≈ µ ∂ v ∂z + ∂ w ∂y (2.17)
• the filtered diffusive species flux vector is:

J i,k = -ρ D k W k W ∂X k ∂x i -Y k V c i (2.18) ≈ -ρ D k W k W ∂ X k ∂x i -Y k V i c (2.19)
where higher order correlations between the different variables of the expression are assumed negligible.

• the filtered heat flux is:

q i = -λ ∂T ∂x i + N k=1 J i,h h h,k (2.20) ≈ -λ ∂T ∂x i + N k=1 J i,h h h,k (2.21)
The forms assume that the spatial variations of molecular diffusion fluxes are negligible and can be modeled through simple gradient assumptions.

Subgrid scale turbulent terms for non-reacting LES Filtering the transport equations leads to a closure problem on the subgrid scale turbulent fluxes. Details on the forms and models available in AVBP to allow to solve the system numerically are given in this subsection.

• the Reynolds tensor is:

τ ij t = -ρ ( u i u j -u i u j ) (2.22)
and is modeled as:

τ ij t = 2ρν t S ij - 1 3 δ ij S ll (2.23)
The modelisation of ν t is explained in Sec. 2.1.1.c.

• the subgrid scale diffusive species flux vector:

J i,k t = ρ u i Y k -u i Y k (2.24)
is modeled as:

J i,k t = -ρ D t k W k W ∂ X k ∂x i -Y k V i c,t (2.25)
with

D t k = ν t S t c,k (2.26)
The turbulent Schmidt number is the same for all species and is fixed in the source code (like P t r ), S t c,k = 0.7. Note also that having one turbulent Schmidt number for all species does not imply, V c,t = 0 because of the W k /W term in Eq. 2.25.

• the subgrid scale heat flux vector:

q i t = ρ u i E -u i E (2.27)
is modeled as:

q i t = -λ t ∂ T ∂x i + N k=1 J i,k t h s,k (2.28)
with

λ t = µ t C p P t r . (2.29)
The turbulent Prandtl number is fixed in AVBP , usually P t r = 0.9. The correction diffusion velocities are henceforth obtained in AVBP using Eq. 2.26 from:

V c i + V c,t i = N k=1 µ ρS c,k + µ t ρS t c,k W k W . ∂ X k ∂x i (2.30)
2.1.1.c Subgrid models for the Reynolds tensor τ t ij

In AVBP three models are available:

• the Smagorinsky model,

• the Wall Adapting Linear Eddy (WALE) model,

• the Filtered Smagorinsky model.

The LES models are derivated on the theoretical ground that the LES filter is spatially and temporally invariant. Variations in the filter size due to non-uniform meshes or moving meshes are not directly accounted for in the LES models. Change of cell topology is only accounted for through the use of the local cell volume; ∆ = V cell 1/3 . Most of the subgrid scale models are based on the Boussinesq's assumption which allow to link the unresolved Reynolds tensor τ t ij to the resolved deformation tensor S through a turbulent viscosity ν t . Such an approach assumes the effect of the subgrid scale field on the resolve field to be purely dissipative, this hypothesis is essentially valid within the cascade theory of turbulence. This leads to a general model for the subgrid scale:

τ ij t = -ρ ( u i u j -u i u j ) (2.31) = 2ρν t S ij - 1 3 ∂ u k ∂x j δ ij , (2.32)
with

S ij = 1 2 ∂ u i ∂x j + ∂ u j ∂x i - 1 3 ∂ u k ∂x k δ ij .
(2.33)

Smagorinsky model

In the original Smagorinsky model, the turbulent viscosity ν t is written:

ν t = (C S ∆) 2 2 S ij S ij , (2.34)
where ∆ denotes the filter characteristic length (cube-root of the cell volume), C s is the model constant set to 0.18 but can vary between 0.1 and 0.18 depending on the flow configuration. The Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. the basic experiment[END_REF] was developed in the sixties and heavily tested for multiple flow configurations. This closure has the particularity of supplying the right amount of dissipation of kinetic energy in homogeneous isentropic turbulent flows. Locality is however lost and only global quantities are maintained. It is known as being too dissipative and not suited to transitioning flows [START_REF] Sagaut | Introduction à la simulation des grandes échelles[END_REF].

Chapter 2 : Description of the numerical tools

WALE model

The WALE model [START_REF] Ducros | Wall-adapating local eddy-viscosity models for simulations in complex geometries[END_REF] uses a different expression for ν t :

ν t = (C w ∆) 2 s d ij s d ij 3/2 S ij S ij 5/2 + s d ij s d ij 5/4
(2.35) with

s d ij = 1 2 g 2 ij + g 2 ij - 1 3 g 2 kk δ ij (2.36)
where ∆ is the filter characteristic length (cube-root of the cell volume), C w = 0.4929 is the model constant and g ij denotes the resolved velocity gradient. The WALE model was developed for wall bounded flows in an attempt to recover the scaling laws of the wall. Similarly to the Smagorinsky model locality is lost and only global quantities are to be trusted.

Filtered Smagorinsky model

The filtered Smagorinsky model was introduced to limit the strong dissipation induced by the Smagorinsky model in the sheared zones. It uses:

ν t = (C S F ∆) 2 2HP S ij HP S ij (2.37)
where ∆ is the filter characteristic length (cube-root of the cell volume), C S F = 0.37 is the model constant and HP S ij denotes the resolved strain rate tensor obtained from a highpass filtered velocity field. This field does not contain the mean shear anymore. This model was developed in order to allow a better representation of local phenomena typical of complex turbulent flows [START_REF] Ducros | Large-eddy simulation of transition to turbulence in a boundary layer developing spatially over a flat plate[END_REF]. With the Filtered Smagorinsky model transition is better predicted and locality is in general better preserved.

2.1.1.d The Thickened Flame model

A typical flame thickness is around 0.1 mm ; in order to resolve it correctly, ten grid points are needed in the flame front. So, for the simulation of industrial configurations of several dozen of centimeters, this leads to a numerical grid resolution out of the capacity of the current calculators. The thickened flame model (TFLES) thickens the flame, allowing to solve it on coarse grid. In effect the flame speed s 0 L and the flame thickness δ 0 L may be expressed as:

s 0 L ∝ D th B and δ 0 L ∝ D th s 0 L = D th B (2.38)
where D th is the thermal diffusivity and B the pre-exponential constant [15]. Increasing the thermal diffusivity by a factor F while decreasing the pre-exponential constant by the same factor F leads to a flame thickness multiplied by F while the flame speed remains the same.

The thickening factor F has to be chosen in order to have enough grid point in the thickened flame. In turbulent flows, the interaction between turbulence and chemistry is modified, eddies smaller than F δ 0 L no longer interact with the flame. As a result, the thickening of the flame reduces the ability of the vortices to wrinkle the flame front and as the flame surface is reduced, the reaction rate is underestimated [START_REF] Angelberger | Large eddy simulations of combustion instabilities in premixed flames[END_REF] (Fig. 2.3). An efficiency function corresponding to a sub-grid scale wrinkling factor, has been derived to account for this effect [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF].

The model used in this study is an improvement of the TFLES called Dynamic Thickened Flame Model (DTFLES) where the thickened factor F is not constant: it is unity outside of the reaction zone and reaches its maximum value inside the reaction zone. In this zone, F depends on the ratio of the mesh size to the flame thickness and is adjusted to have five points within the reaction thickness.

CH 4 + 3 2 O 2 -→ CO + 2H 2 O (2.39) CO + 1 2 O 2 ←→ CO 2 (2.40)
The first equation is irreversible while the second one is reversible and leads to an equilibrium between CO and CO 2 in the burnt gases. This scheme was developed for CH 4 /air flames and used here without modification for

CH 4 /O 2 /H 2 O combustion.
A common shortcoming of two-step mechanisms is that burning rates are over-predicted for rich mixtures. Following the procedure presented in [START_REF] Franzelli | A two-step chemical scheme for Large-Eddy Simulation of kerosene-air flames[END_REF], the present scheme uses a preexponential constant of the first reaction that is a function of the local equivalence ratio so as to recover accurate burning rates in rich mixtures. This scheme was validated for the thermodynamic conditions of the CVRC versus the Gri-Mech 3.0 mechanism.

First, equilibrium states were compared to ensure that no important species would be neglected. Tab. 2.1 presents the burnt-gases properties for both schemes and Fig. 2.4 shows the Eq.ratio equilibrium adiabatic temperature of the burnt gases T ad b versus mixture fraction z. Differences between the complete mechanism and the six species equilibrium are observed for rich condition (z > z st ), but for lean and stoichiometric combustion (z ≤ z st ). Figure 2.5 is a scatter plot of heat release rate and z in all the computational domain, extracted from the LES (Chapter 4). It shows that the combustion occurs mainly at low mixture fraction (z < 0.2) where considering only six species for the chemical scheme gives the right adiabatic temperature. Because it is expected that the combustion regime will be partially premixed and lean, the premixed laminar flame speeds at stoichiometry were also compared (for P = 14 bar, T = T 0 mix (z st ) = 881 K): s L = 1.59 m.s -1 for Gri-Mech 3.0 and s L = 1.52 m.s -1 for the present two-step scheme. It was also checked that auto-ignition times for Gri-Mech like for the 2-step scheme (0.122 s at z = z st ) are much larger than the residence time in the chamber (5.10 -3 s). This is consistent with the fact that under the conditions of the CVRC experiment, combustion is not controlled by auto-ignition. Overall, these calculations validate the approximations implied by the use of the reduced scheme. Of course, the main interest of the two-step scheme is that it mimics the behavior of the full scheme at a low computational cost, a very useful property for our study.

z Scheme T ad b [K] Y CH4 Y O2 Y H2O Y CO2 Y CO Y N2 1 z st = 0.095 Gri-

Wall treatment

No slip walls (AVBP WALL NOSLIP boundary condition) are used for all walls except the outlet nozzle (Sec. 3.1.2). They simply impose a zero-velocity condition at the walls. Nozzle's walls are slip walls (AVBP WALL SLIP boundary condition) where a zero normal velocity is imposed. Both adiabatic ( ADIAB) and heat losses ( LOSS) conditions will be used. For the heat losses, only convective losses at the walls are accounted for (no radiation) with a simplest Conjugate Heat Transfer (CHT) model. The code solves an implicit equation:

φ = T w2 -T out R = φ LES (T w2 ) (2.41)
where φ is the heat flux through the wall, T w2 the temperature at the wall (Fig. 2.6). T out and R are the two input parameters of the wall heat transfer model corresponding respectively to the external temperature (typically the water temperature if the system is water cooled) and the heat resistance defined by :

R = 1 h out + e λ (2.42)
where h out is the convective heat transfer coefficient on the wall, outside of the chamber itself.

The value of T out and R chosen for the simulation of the CVRC are given in Sec. 3.1.2.

Validation of the LES solver for wave computation

Since the focus of this work was combustion instabilities, an additional work was performed to verify the code performances for wave computation which are critical for LES of combustion instabilities: vorticity and acoustic waves must be propagated at the right speeds and amplitudes to capture unstable modes correctly.

In order to check the accuracy of AVBP for wave propagation, simple academic tests were performed for both the acoustic and the hydrodynamic mode resolution and compared to analytical solutions. Two phenomena (at least) must be correctly captured in a simulation of acoustically driven combustion instability;

• acoustic waves must propagate accurately because they control the eigen-modes present in the chamber (dispersive effects must be minimum to capture the right frequencies and a limited numerical dissipation is required to obtain the correct amplitudes)

• vortices created in shear zones must be predicted accurately because they control the formation of vortices in the chamber and therefore unsteady combustion. AVBP has already been tested in details2 to verify that it was able to propagate vortices with very limited dissipation and dispersion effects [START_REF] Moreau | A priori testing of subgrid stress models for eulereuler two-phase LES from euler-lagrange simulations of gas-particle turbulent flow[END_REF]. Here, we will focus on the code's capacities to predict the birth and growth of vortices in a shear zone. The hydrodynamic mechanism [START_REF] Lieuwen | Unsteady combustor physics[END_REF][START_REF] Drazin | Hydrodynamic stability[END_REF] which leads to the first unstable modes is essential and we will check if AVBP captures it and how this depends on the resolution and on the numerical scheme.

2.1.4.a Acoustic eigenmodes in a closed 1D box

For this test, the eigenmode of a 1D box with walls on both sides is initialised in AVBP (Fig. 2.7). All variables (for example the pressure P) are decomposed in a mean part (P ) and a perturbation (P ). The analytical expressions for pressure (P), velocity (u), density (ρ) and temperature (T) are (in the linear regime):

P = P + P with P = P 0 cos (kx) cos (ωt) (2.43) u = 0 + u with u = u 0 sin (kx) sin (ωt) (2.44) ρ = ρ + ρ with ρ = P c 2 (2.45)
u 0 = P 0 ρc (2.46)
T = P ρr (2.47)
where L is the length of the box, L = 0.1 m. The initial solution is set at time t = 0 using Eq. 2.43 to 2.47. The amplitude of the pressure wave is P 0 = 150 Pa, the mean pressure is 15.10 5 Pa (P 0 /P = 10 4 ). The gas has the same composition as the burnt gases for the global mixture fraction of the experimental operating point (z xp = 0.079) (Tab. 2.1) which will be the the main component in the combustion chamber studied in all the following chapters, the sound speed in this gas is c = 1096 m/s. Two numerical schemes are compared: Lax-Wendroff and TTGC [START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF]. Figure 2.8 shows pressure and velocity perturbation in the duct at the initial time and after 15, 30 and 45 periods. The mode shape is perfectly conserved for TTGC. A limited dissipation is observed for Lax-Wendroff. Even with a small number of points per wave length (N = λ/∆x = 20), the accuracy of the code is quite sufficient. To quantify the dissipation more accurately, Fig. 2.9(a) shows the cycle by cycle damping D of the signal pressure on the left wall of the chamber versus the number of points per wave length. The cycle by cycle damping rate is calculated by:

D = |P last cycle | |P f irst cycle | 1 number of cycles (2.48)
A good energy conservation is observed except for very poorly resolved acoustic wave lengths with less than five mesh points per spatial period for the Lax-Wendroff scheme (λ/∆x < 5). Fig. 2.9(b) shows the ratio between the frequency of the acoustic mode measured in AVBP (by taking the frequency of the highest peak of the FFT (Fast Fourier Transform) of the pressure signal), and its theoretical value f th = c/2L versus the number of mesh points per wave length. When more than twenty points per wave length are used, the shift of the frequency is negligible. As expected from numerical analysis [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF][START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF], TTGC is more precise than Lax-Wendroff.

In the case of the CVRC experiment, typical frequencies range from 1200 Hz to 1500 Hz so that the expected acoustic wave lengths are of the order of 1 m. For these meshes, the number of points per wave length will be much larger than 100 so that all AVBP schemes will propagate acoustic waves accurately with negligible errors. We will first recall the analytical solution for the Orr-Sommerfeld equation in a parallel flow with the velocity profile of Fig. 2.10(c) [START_REF] Lieuwen | Unsteady combustor physics[END_REF][START_REF] Drazin | Hydrodynamic stability[END_REF]. Using the notations of Sec. 2.1.4.a with the x-velocity is u(x, y) = U (y) + u and the vertical y-velocity is v(x, y) = 0 + v . Linearizing the equations of momentum and mass conservation in two dimensions for an incompressible flow leads to the following system:

2.1.4.b Linear stability of shear layers

               ρ ∂u ∂t + ρU ∂u ∂x + ρu ∂U ∂y = - ∂P ∂x + µ ∂ 2 ∂x 2 + ∂ 2 ∂y 2 u ρ ∂v ∂t + ρU ∂v ∂x = - ∂P ∂y + µ ∂ 2 ∂x 2 + ∂ 2 ∂y 2 v ∂u ∂x + ∂v ∂y = 0 (2.49)
Under the hypothesis of zero-viscosity, the current function ψ, can be introduced. ψ is defined by

u = ∂ψ ∂y et v = - ∂ψ ∂x
. Solutions of ψ can be searched as [START_REF] Drazin | Hydrodynamic stability[END_REF]: -ct) where

ψ (x, y, t) = φ (y) e iα(x
√ i = -1 and α = 2π λ (2.50)
where λ is the wave length of the perturbation and c its speed. With the same method, the perturbation of pressure can be searched as P = P e iα(x-ct) . Eq. 2.49 can be written as:

ρ -iαcφ + U iαφ -iαφU = -iα P ρ -α 2 cφ + α 2 U φ = -P (2.51)
where the exponent indicates the derivative a single-variable function. If U corresponds to an infinite shear (Fig. 2.10(b)) or is piecewise linear (Fig. 2.10(c)), a solution for φ can be expressed as φ (y) = Ae ±αy where A is a real constant. Jump relations are also needed where the U profile is not C 1 : at y = 0 for Fig. 2.10(b) or at y = ±d/2 for Fig. 2.10(c). The first one is based on the conservation of pressure in x-planes, ∆ P = 0 which leads to (with Eq. 2.51):

∆ U -c φ -U φ = 0 (2.52)
The second jump relation comes from the conservation of momentum [START_REF] Drazin | Hydrodynamic stability[END_REF]:

∆ φ U -c = 0 (2.53)

2.1.4.d Equation of the perturbation for the infinite shear case

In this section, expressions of the perturbation of velocity and pressure for the infinite-shear case (Fig. 2.10(b)) will be first demonstrated. Assuming the following shape for φ: φ (y) = Ae -αy if y 0 Be αy if y < 0 where A and B are real constants (2.54)

and introducing there expressions in the jump relations Eq. 2.52 and 2.53, gives the dispersion relation:

c = iU (2.55)
and leads to the analytical expression of the velocity perturbation for a chosen wave length λ:

u 1 = -iαBe -αy e iα(x-ct) v 1 = αBe -αy e iα(x-ct) and u 2 = αBe αy e iα(x-ct) v 2 = -iαBe αy e iα(x-ct) (2.56)
where the index 1 or 2 designates the y > 0 and y < 0 domains respectively. The perturbation of pressure can be found using Eq. 2.51:

P = ρ U φ -U -c φ (2.57)
Expressions Eq. 2.56 and 2.57 can be used to initialized an AVBP by setting t = 0 in the expression of u 1 , u 2 and P . The growth rate of the hydrodynamic instability for the infinite shear flow is σ ∞ = αU . The initial fields of u , v and P used to start AVBP at t = 0 are displayed in Fig. 2.11.

2.1.4.e Relation of dispersion for the linear shear case

For the linear case (Fig. 2.10(c)), the mean velocity has the following expression:

U (y) =            U 1 if y d 2 U 1 -U 2 d y + U 1 + U 2 2 if d 2 > y > - d 2 U 2 if y - d 2 (2.58)
Solutions of φ are searched under the shape:

φ (y) =            Ae -α(y-d 2 ) if y d 2 Be -α(y-d 2 ) + Ce α(y+ d 2 ) if d 2 > y > - d 2 De α(y+ d 2 ) if y - d 2 (2.59)
where A, B, C and D are real constants. Applying jump relations Eq. 2.52 and 2.53 on both transition points y = d/2 and y = -d/2 leads to the following expression written in matrix form:

    -α (U 1 -c) α (U 1 -c) + a [-α (U 1 -c) + a] e αd 0 0 [α (U 2 -c) + a] e αd -α (U 2 -c) + a α (U 2 -c) 1 -1 -e αd 0 0 -e αd -1 1         A B C D     = 0 (2.60)
The determinant of this matrix must be zero, leading to the dispersion relation:

c = U 1 + U 2 2 ± U 1 -U 2 2αd (αd -1) 2 -e -2αd 1 2
(2.61)

Chapter 2 : Description of the numerical tools Thus the theoretical growth rate of the mode of wavelength λ = 2π/α is given by taking the imaginary part of the term αc : σ lin = Im (αc). The growth rate is measured by observing the variation of the integral of kinetic energy of vertical velocity on the domain Ec(v):

Ec(v) = 1 2 S v 2 dxdy (2.63)
Since v varies like e -iαct , its time variation at all points scales like e αIm(c)t and Ec (v) change likes e 2αIm(c)t . Therefore plotting 0.5 log Ec(v) versus time must provide the growth rate Im (αc). Fig. 2.13 shows the evolution of 0.5 log Ec(v) versus time for cases described in Tab. 2.2. Linear interpolation is used to measure the growth rate σ. For most cases, a clear linear region is identified where the growth rate can be measured as shown in Fig. 2.13(a). For other cases, the growth rate can not be measured or is very small as shown in Fig. 2.13(b). For such cases, some errors on σ must be expected but in most cases, they correspond to small σ values. Figure 2.14 shows the evolution of vorticity and velocity at several times for the case λ = 0.8 m in Tab. 2.2. When vortices appear (typically after the linear phase, around t = 0.35 s), the linear growth hypothesis is not valid anymore. The comparison between AVBP and linear theory is plotted on Fig. 2.15, for a case where the resolution is given by λ/∆x = 100 and λ/∆y ≥ 45. For this resolution, the code recovers the theory for all wave lengths. The normalized growth rate σd/(U 1 -U 2 ) plotted versus the reduced wave length of the perturbation αd, given by AVBP matches the theoretical result for all numerical schemes which were tested ; LW, TTGC, TTG4A.

L x = λ (m) L y (m) d (m) N x N y λ ∆x λ ∆y λ d U 1 (m/s) U 2 (m/

2.1.4.g Effects of resolution

The example of Sec. 2.1.4.f shows that AVBP captures the growth of hydrodynamic modes in a shear layer accurately for all schemes (LW, TTGC, TTG4A) when the resolution is sufficient : λ/∆x = 100 and λ/∆y = 45 were used for these tests. 10 points were used to resolve the mean velocity gradient (d/∆y = 10) and 100 points were used to resolve the spatial wavelength. In practice, however, such a resolution is rarely obtained in combustor simulations where mean gradients are resolved over a few mesh points and vortices can have sizes of the order of a few cells. In the LES, actually, there are always vortices which have dimensions very close to the mesh size because, by definition, some part of the spectrum is unresolved and the filter is of the order of the mesh size. Therefore it is important to see how the accuracy of the code evolves when the resolution becomes poorer (or equivalently when the vortices become smaller on a given grid). Twelve simulations for each scheme were performed for a fixed shear layer thickness d = 0.1 m and a fixed wave length λ = 0.8 m but with varying resolution in the horizontal or in the vertical direction. For each simulation, the growth rate σ is measured on the growth of the kinetic energy of vertical velocity Ec(v) (Eq. 2.63). Fig. 2.16 shows isolines of σ/σ lin for the different numerical schemes (σ lin is the exact solution given by Eq. 2.61). This ratio is 1 for a perfectly accurate simulation.

Figure 2.16 displays the σ/σ lin ratio plotted versus the number of points used to resolve the wavelength of the most amplified mode in the axial direction (λ/∆x) and the number of points to resolve the mean velocity gradient in the transverse direction (d/∆y). Fig. 2. 16 shows that TTGC (Fig. 2.16(b)) and TTG4A (Fig. 2.16(c)) exhibit a robust behavior. They are precise (σ/σ lin > 0.8) as soon as the resolution is sufficient, typically as soon as d/∆y > 5 and λ/∆x > 50, TTGC being better than TTG4A as expected from their theoretical dispersion and dissipation properties [START_REF] Colin | Development of high-order taylor-galerkin schemes for unsteady calculations[END_REF]. For TTGC like TTG4A, when the resolution decreases (ie when λ/∆x or λ/∆y or d/∆y decrease), the growth rate of hydrodynamic perturbation decreases. In other words, these schemes damp the growth of vortices when the shear layers are underresolved. On the other hand, Lax-Wendroff behaves differently. At low resolution, probably because of its insufficient dispersive characteristics, it amplifies vortices growth more than they are in the real world (σ/σ lin > 1). When the transverse resolution is low, for small values of d/∆y (under-resolved velocity gradients), σ/σ th can reach values much larger than unity, suggesting that performing a LES with Lax-Wendroff on a poorly resolved grid can lead to the formation of non physical vortices.

Computation of choked flows with LES

Calculating configurations with a choked nozzle in LES is delicate and some attention must be paid, especially for the initialization of the simulation. This section recalls the basic feature of flows in nozzles and links them to their computation in AVBP . In a 1D nozzle, the flow parameters can be calculated analytically, using the conservation of mass flux and enthalpy under isentropic assumption allow to link the ratio of area A to throat area A * as a function of a Mach number M [START_REF] Candel | Mécanique des Fluides[END_REF]:

A A * = 1 M 2 γ + 1 1 + γ -1 2 M 2 γ + 1 2 (γ -1) .
(2.64)

The quantities to be imposed for the initial solution are then given by:

T 0 T = 1 + γ -1 2 M 2 , (2.65)
P 0 P = 1 + γ -1 2 M 2 1 γ -1 ,
(2.66)

ρ 0 ρ = 1 + γ -1 2 M 2 γ γ -1 , (2.67) (2.68)
where the subscript 0 corresponds to total values of otherwise static variable. But, as the geometry of real nozzles is complex, using Eq. 2.64 to calculate the Mach number and the other flow parameters is not easy and often, to initialize the simulation with a zero-velocity solution is preferred. Two kinds of nozzles are used in LES:

• simply converging nozzles (Fig. 2.17(a))

• converging-diverging nozzles (Fig. 2.17(b))

For the following explanation in both cases, the Mach number in the chamber upstream of the nozzle is negligible, therefore, the static pressure in the chamber P 1 is equal to the total pressure P 0 and is considered as constant.

Converging nozzle: in order to choke the flow, there must be a large difference in pressure between the chamber and the outlet. The flow will choke as soon as P 1 = r C P 2 where r C is the critical pressure ratio:

r C = P 0 P M =1 = γ + 1 2 γ γ -1 .
(2.69)

Converging-diverging nozzle: starting from a zero-velocity state in which P 2 = P 0 , reducing P 2 leads to the following events:

1. the flow begins to move, the highest Mach number is at the smallest section, the throat, but still smaller than one. 2. the Mach number reach one at the throat, the flow in the diverging section remains subsonic. This occurs at P 2 = P choke function of P throat and A 2 /A * with P throat = P 0 /r C .

3. P 2 is further decreased, the static pressure at the throat stays constant but the flow accelerate after the throat and the pressure continues to drop under supersonic conditions. The pressure of the flow exiting the nozzle is lower than the ambient pressure P 2 . In this case there is a two-dimensionnal oblique shock to rise the pressure of the jet to the ambient pressure (Fig. 2.18(a)). This flow is called over-expanded.

4. there is a lower P 2 = P adapt for which no shock is needed in the diverging section. The flow is fully supersonic in the diverging section, the Mach number is determined by A/A * (Eq. 2.64). It is called an adapted flow.

5. for P 2 < P adapt the flow is under-expanded. To match the ambient pressure, the flow must still expand and there are Prandtl-Meyer expansion waves outside the nozzle (Fig. 2.18(b)).

In combustion, the computational domain usually stops before the end of the nozzle where the flow is exactly the same than for P 2 = P adapt . In practice, the LES can be started from an initial solution which is choked at the throat or not (zero velocity everywhere for example). Starting from a zero-velocity initial solution to compute a choked flow is tough because, at one instant, the code has to transition from a subsonic outlet (where the pressure must be imposed) to a supersonic outlet (where nothing can be imposed anymore). This criterion is achieved in AVBP with the OUTLET RELAX P boundary condition which measures the velocity and the Mach number at each point on the outlet patch. This initialization method remains very violent for the code: in the first iterations, when the flow is subsonic, the pressure is highest at the outlet of the nozzle where a very low pressure is imposed, this leads to a large pressure gradient. Therefore, the relaxation coefficient of this boundary condition has to be chosen wisely: if it is too high, the initial pressure gradient is interpreted as a very strong wave hitting the boundary condition and will be then reflected into the domain as a strong acoustic wave. If it is too low, the strong gradient may cause an important pressure drift at the boundary which might choke the nozzle during the establishment and then keep it choked even if the flow is not supposed to choke. It is suggested [START_REF] Selle | The actual impedance of non-reflecting boundary conditions: implications for the computation of resonators[END_REF] to use a relaxation coefficient K = 4πc 0 /2L where L is the total nozzle length and c 0 the sound speed in the chamber, to cut off all frequencies higher than the first longitudinal harmonic.

For acoustic studies, it is important to notice that since the nozzle is choked, it acts as a wall at the throat and reflect all acoustic waves. In Chapter 3, we will describe a tool developed at Cerfacs (by T. Livebardon and C. Lapeyre) which allows to track the flow state in a (chamber pressure vs outlet flow rate) diagram to analyze if the LES is and remains choked or not at all instants of the simulation.

Acoustic solver for combustion instabilities

In the CVRC experiment, the post length is continuously varied during the experiment and the resulting unstable modes are recorded as a function of its length L op . Performing an LES where the length of the post changes continuously or even performing a set of 10 to 20 LES for variable post lengths is not a good solution. It is too expensive and interestingly in this case, it is probably not smart: in the CVRC experiment, the flame mean field and the FTF of the flame are very similar during the whole run: only the post length L op varies, modifying only the acoustic impedance of the combustor inlet. Therefore studying how the modes of the CVRC change with post length L op can be done by obtaining FTF information from one single LES run and assuming that it does not change 3 . Then a Helmholtz solver can be used to study the stability of the combustor when the inlet post length is changed. Such Helmholtz solvers are much cheaper than LES. Furthermore they allow also to gain more insight into the mode structure and into possible control methods. In this work the Helmholtz solver AVSP of CERFACS was used [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF][START_REF] Sensiau | A tool to study azimuthal and spinning modes in annular combustors[END_REF][START_REF] Gullaud | Impact des plaques multiperforées sur l' acoustique des chambres de combustion aéronautiques[END_REF].

To mimic the behavior of the CVRC with variable post length in AVSP , two solutions may be used (Fig. 2.20): • Approach 1 (Fig. 2.21): mesh the entire geometry (post + chamber) for various post lengths L op and compute the eigenmodes for each post length using a new mesh for each post length value.

• Approach 2 (Fig. 2.22): mesh only the combustion chamber and replace the acoustic effect of the entire post (length L op ) by its equivalent impedance. Indeed, since the post is a simple tube at constant temperature, it can be acoustically replaced by an impedance imposed at the combustor inlet (x = 0) (impedance translation theorem [23]). This impedance changes when L changes and can be determined analytically. The flame response (FTF, see Sec. 4.4) remains the same. In AVSP , this means that one single mesh and one single set of FTF can be used to explore the whole range of post lengths simply by changing the impedance at the combustor inlet Z(x = 0).

Approach 2 is more elegant than Approach 1 and it was chosen here. Note however that it raises additional questions: for each post length the impedance Z(x = 0) will depend on the frequency and require specific convergence methods in the Helmholtz solver. To validate this methodology, a simple non-reacting example corresponding to a case where the chamber is replaced by a straight tube and for which a fully analytical solution is easily derived is presented in Sec. 2.2.3 to verify that the method works. Before this test, the Helmholtz solver AVSP is first described (Sec. 2.2.1). Results for the CVRC will be presented in Chapter 5. The objective is to reproduce the experimental procedure where the burner is operated at one fixed point in terms of flow rates and mean pressure and where the post length L op varies continuously from 8.9 cm to 19.1 cm.

Description of the AVSP code

If the mean flow is supposed to be frozen and acoustic perturbations are small, the full Navier-Stokes equations can be linearized, leading to a wave equation for the pressure perturbation P which is solved in the frequential domain, writing P = P e -iωt . The code AVSP solves the eigenvalue problem associated to this Helmholtz equation [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF]:

∇. c 2 0 ∇p + ω 2 p = iω (γ -1) q (2.70)
where q is the amplitude of the unsteady heat release q ( x, t) = q ( x) e -iωt , c 0 is the speed of sound, and ω is the complex pulsation. The problem is closed by modeling the flame response with the nτ model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. part II[END_REF] which amplifies and delays the flame response q by reference to the acoustic perturbation at a given point. Eq. 2.70 is only valid under the zero Mach number assumption, assuming that the mean velocity of the flow is very small compared to the speed of sound. Furthermore, Eq. 2.70 does not account for entropy waves. Entropy waves passing through a nozzle can be an important source of noise [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF]25] which is ignored in AVSP .

Validation case : variable length tube

Considering a one-dimensional tube with a diameter of 5 cm filled with air at 300 K corresponding to a uniform sound speed c = 347 m/s. The right boundary condition (x = L) imposes a zero velocity fluctuation, and on the left boundary an acoustic impedance Z (x = -l) is imposed. For such a simple case, Sec. 2.2.2.a presents an analytical solution for all values of Z (x = 0). AVSP can then be run on the cut configuration (a duct from x = 0 to x = L) with a variable impedance Z (x = 0) mimicking a variable post length l to check that AVSP recovers the theoretical results.

2.2.2.a Analytical expression

This section presents the analytical solution of the problem of Fig. 2.23 between x = -l (fixed impedance Z (x = -l)) and x = L (wall). It is a useful exercise to validate our method in AVSP .

Under the hypothesis of:

• no combustion

• zero volume force

• zero volume heat source

• negligible viscous forces

• linear acoustics: acoustic variables (noted ) are supposed to be small compared to reference quantities • isentropic variations: with the assumptions of no heat release and no viscous terms, the flow remains isentropic if it is homogeneous and

• low speed mean flow: u 0 = 0
The acoustic solution in 1D is:

P = A + e ikx + A -e -ikx and u = 1 ρc A + e ikx -A -e -ikx (2.71)
The acoustic impedance is defined as:

Z = 1 ρc P u (2.72)
At x = L, the condition u = 0 leads to:

A + = A -e -2ikL (2.73)
At x = -l the impedance Z -l is imposed so that:

Z -l = A + e -ikl + A -e ikl A + e -ikl -A -e ikl (2.74)
which with Eq. 2.73 leads to:

e -2ik(l+L) = Z -l + 1 Z -l -1 (2.75)
To determine the wave number k two cases are possible depending on Z -l :

2.2 Acoustic solver for combustion instabilities 65

• if |Z -l | < 1, Eq. 2.
75 is equivalent to:

   e -2ikr(l+L) = -1 e 2k i (l+L) = Z -l + 1 Z -l -1 (2.76)
were k r is the real part of k and k i its imaginary part. Then:

     k r = (2n -1) π (l + L) with n ∈ N k i = 1 2 (l + L) ln Z -l + 1 Z -l -1 (2.77)
which corresponds to a quarter wave mode when n = 1.

• if |Z -l | > 1, Eq. 2.
75 is equivalent to:

   e -2ikr(l+L) = 1 e 2k i (l+L) = Z -l + 1 Z -l -1 (2.78)
And then:

     k r = nπ (l + L) with n ∈ N k i = 1 2 (l + L) ln Z -l + 1 Z -l -1 (2.79)
which corresponds to a half wave mode when n = 1.

The pressure and velocity perturbations can be expressed using Eq. 2.73:

P = A -e ikx e -2ikL
+ e -ikx e -iωrt e ω i t (2.80)

u = A - ρc e ikx e -2ikL + e -ikx e -iωrt e ω i t (2.81)
where ω i is the growth rate of the mode [23]:

ω i = k i c = c 2 (l + L) ln Z -l + 1 Z -l -1 (2.82)
In Sec. 3.2.2 we will study the impedance of nozzles. This expression will be useful to evaluate the losses due to the mean flow through the nozzle.

2.2.2.b Numerical configuration

The domain is smaller than the one pictured in Fig. 2.23. The duct from x = -l to x = 0 is not included and will be replaced by a translated impedance imposed at x = 0. The geometry computed with AVSP is a cylindrical tube of one half meter length and a diameter of five centimeter. An impedance Z 0 is chosen on the left face (x = 0) and all other walls have a zero velocity boundary condition (Fig. 2.24). The mesh is a very coarse mesh composed of 5800 tetrahedral unstructured cells with a size of 1 cm (Fig. 2.25). As a first test, AVSP was used for a fixed inlet length (l = 0) and a fixed inlet impedance Z 0 to verify that the code works for a fixed impedance case. The prediction of AVSP for the modulus and phase of P is compared with the theory for l = 0 and L = 0.5 m, using Z 0 = 0.5 for the fist mode (Fig. 2.26) and the second mode (Fig. : analytical expression, × : AVSP .

Impedance of a variable length tube: the impedance translation theorem

The negative x part of the tube (-l < x < 0) can be modelled by an impedance Z 0 at x = 0 which depends on l, on Z -l and on ω = kc. Impedances at x = 0 and x = -l are:

     Z 0 = A + + A - A + -A - Z -l = A + e -ikl + A -e ikl A + e -ikl -A -e ikl (2.83)
This leads to the expression of Z 0 in function of Z -l and l :

Z 0 = -Z -l e -2ikl + 1 + e -2ikl -1 Z -l (e -2ikl -1) -e -2ikl -1 (2.84)
which is known as the impedance translation theorem [START_REF] Pierce | Acoustics: an introduction to its physical principles and applications[END_REF].

2.2.3.a Validation on the tube test case

AVSP is now launched with only the x positive part of the tube meshed while the x negative part is modeled for different tube lengths l with a variable impedance Z 0 given by Eq. 2.84, see Fig. 2.24. The variation of the frequencies of the first and second mode regarding l is plotted for Z -l = 0.5 (Fig. 2.30) and Z -l = 2.0 (Fig. 2.31). Both cases agree with the theory for the two first modes frequency and amplification. : theoretical 1 st mode, × : AVSP for the 1 st mode, : theoretical 2 nd mode, • : AVSP for the 2 nd mode.

structure for l = 0.5 m with Z -l = 0.5 and Fig. 2.33 shows this one for Z -l = 2.0. These modes structures are compared to the theoretical structure of the full tube of length l + L: again AVSP recovers perfectly the theory. This method will be applied on the real case of the CVRC experiment in Sec. 5.

Impedance of nozzles

Several methods are available to model the impedance of the nozzle in an Helmholtz solver. Of course the Helmholtz solver which works under the Mach= 0 assumption can not be used in a nozzle itself where the flow can be choked. The usual procedure is to replace the whole nozzle by an equivalent impedance which must be computed by another model. The computation of nozzle impedances is a classical topic in acoustics and we will simply recall here the main results [23,24,[START_REF] Duran | Analytical and numerical study of direct and indirect combustion noise through a subsonic nozzle[END_REF]. From the simplest to more complex approaches, one can:

• impose a zero velocity perturbation at the nozzle throat (Fig. 2.34). Considering that most of the duct section is closed when the nozzle cross section is small compared to the duct section (high blockage). Assuming u = 0 at the throat is reasonable and leads to an infinite impedance. The inconvenient is that, near the throat, the code calculates acoustic perturbation on a high velocity flow where the Helmholtz hypothesis can not apply: in the converging nozzle, this may be a strong limitation of the method accuracy.

• cut the geometry just before the nozzle (Fig. 2.35) and impose a constant impedance equivalent to the whole nozzle. When the charactericic nozzle length is small in comparison to the characteristic wave length of the acoustic wave, the wave reflection at the nozzle can be considered instantaneous and the nozzle flow may be assumed to be quasi-steady [23,[START_REF] Bell | Experimental and theoretical determination of the admittances of a family of nozzles subjected to axial instabilities[END_REF][START_REF] Bell | the prediction of three-dimensional liquid-propellant rocket nozzles admittances[END_REF]. Solving the steady-state equation leads to the following relation developed by Marble and Candel [24]:

Z out = 2 (γ -1) M . n (2.85)
where M = u/c is oriented with a normal n. This expression is only valid for low frequencies (with large wave length). For high frequencies the compact assumption is not valid.

• cut the geometry before the nozzle, but impose a variable impedance Z out (f ). Duràn et al [START_REF] Duran | Solution of the quasi one-dimensional linearized Euler equations using flow invariants and the Magnus expansion[END_REF] have analytically extended the method of Marble and Candel to non-zero frequencies removing the compact nozzle assumption. They assumed however that the nozzle is long compared to its characteristic transversal wave-length so that the flow can be considered as one-dimensional. Then they use the Euler equations, neglecting viscous terms and develop an analytical method using the Magnus expansion [START_REF] Magnus | On the exponential solution of differential equations for a linear operator[END_REF][START_REF] Blanes | The magnus expansion and some of its applications[END_REF] to solve the differential equations systems for all frequencies and all geometries under these two assumptions. With this method, extracting the evolution of Mach number and radial area along the nozzle axis from the LES, allows to generate a tabulated complex impedance Z out (f ) for all frequencies.

Effects of Mach number

The Mach number is supposed to be zero in Eq. 2.70 but in practice, even if non-zero Mach effects are negligible within the flow, they are quite important in inlets and outlets where they change acoustic fluxes and the mode growth rates. In other words, it is possible to assume that M ≈ 0 in the domain but Mach effects should appear as inlet and outlet impedances. This can be achieved using the recent work of Motheau et al in [START_REF] Motheau | Using boundary conditions to account for mean flow effects in a zero mach number acoustic solver[END_REF] who compare acoustic formulation at zero Mach (as done in AVSP ) and more general formulation at non-zero Mach numbers based on Linearized Euler Equations (LEE). Acoustic impedances for Helmholtz solvers are expressed using the pressure and velocity perturbations (P and u ) while, working with the LEE approach, the linearized isentropic fluctuating parts of the total enthalpy (J = C p T + U 2 /2) and the mass flow rate (m) which include the mean flow velocity are more attune. Motheau et al have shown that these two sets of independent variables are equivalent and one may write the LEE or Helmholtz equation with any set (P , u ) or (J , m ) defined by:

J = P ρ 0 + u 0 u , (2.86)
m = ρ 0 u + u 0 c 0 2 P . (2.87)
Then an impedance Z and reflection coefficient R may be defined for each set of variables:

Z (P ,u ) = P ρ 0 c 0 u (2.88)
R (P ,u ) = A - A + (2.89)
and

Z (J ,m ) = ρ 0 c 0 J m (2.90) R (J ,m ) = A -(1 -M ) A + (1 -M ) . ( 2

Numerical setup for LES simulations

This chapter presents the LES study of a specific case for the CVRC combustor corresponding to an oxidizer post length L op = 4.75 in (Fig. 1.5). This is an unstable operating point.

Geometry and mesh

The geometry used for the simulation is the exact geometry of the experiment extracted from CAD files, the simulation domain starts in the oxidizer tube just past the movable nozzle. All LES simulations were performed for an oxidizer post length L op = 4.75 in. The mesh is composed of 2,500,000 nodes (14,000,000 cells) and is presented in Fig. 3.1. Cells sizes vary from 0.1 mm (in the straight fuel post in order to have ten cells on the radial length of the methane injector) to 3 mm at the end of the chamber were there is no combustion to solve and a coarser mesh is sufficient. The mixing zone between oxygen and fuel is refined with cell sizes varying from 0.1 to 0.5 mm. Preliminary simulations have shown that the walls of the chamber have to be well resolved in order to capture the burnt-gases recirculation which has a major impact on the flame stabilization (Sec. 4.2).

The numerical scheme used for these simulation is the Lax-Wendroff scheme. The thickening of the DTFLES model (Sec. 2.1.1.d) F varies from 1 to 100 with a mean thickening of 2.3.

Boundary conditions

The boundary conditions used for LES are presented in Fig. 3.2: they use the NSCBC method [35]. Walls are adiabatic with a slip condition in the nozzle and a no-slip condition in the rest of the chamber, injector tubes and in particular the fuel post are very thin, imposing zero velocity at wall has a strong impact on the flow in these posts. For the inlet the mass flow rate, temperature and mass fraction are imposed with a relaxation coefficient of 10 6 s -1 . At the outlet the pressure is imposed at 1 bar, with a low relaxation coefficient (10 3 s -1 ) imposed in the characteristic boundary condition, but as soon as the flow becomes supersonic through this outlet, no conditions are imposed. A study of the influence of the heat losses on the stability of the CVRC will be done in Sec. 4.6. For this part, convective heat losses are imposed only on the chamber wall, the other walls having a negligible surface. This heat losses boundary condition has been described in Sec. 2.1.3. The two parameters of the model (the external temperature T out and the heat resistance of the wall R) are set to:

T out = 300 K R = 1.5 × 10 -4 K•m 2 /W
because they allow the mean pressure to adapt to the measured value.

Initialization procedures for LES of choked flows

The initialization of the LES in a choked domain requires different steps and may change the final results as discussed below.

3.1.3.a Initialization of non-reacting choked flows

At initialization the chamber is filled with oxidizer (Y O 2 = 0.4235 and Y H 2 O = 0.5765) at 13.5 bar with zero velocity. There is no chemical reaction. At t = 0 the injection of oxidizer and fuel begins. 3.4 shows that pressure and outlet flow rate both vary strongly with time and demonstrates the need to follow the flow state in a (Pressure -Flow rate) diagram and to be able to determine when the nozzle is choked or not. Equation 1.12 gives the relation between the mean pressure in the chamber and the flow rate through a choked nozzle. Evolutions of chamber pressure versus flow rate are plotted on Fig. 3.5(a) 1 . The pressure first decreases slowly while the flow rate hardly increases, over-passing the theoretical value of the choked nozzle before going back to it at t = 3.8 × 10 -4 s. The next portion of this path is presented in Fig. 3.5(b). After a loop at high pressure (t = 7.3 × 10 -4 s), the path joins the theoretical line, and from here, pressure and flow rate decrease following the theoretical relation. From Eq. 1.12, for the injection flow rate at the mixing temperature T 0 mix = 924 K (Sec. 1.2), one can predict that a simultaneous decrease of both P and ṁout will take place until the chamber reaches a pressure of 8 bar (Fig. 3.4(b)). In practice, there is no need to wait for this stabilization and ignition is activated at t = 4 ms.

Figure 3.6 presents instantaneous fields of this cold flow before ignition. The methane has propagated in the chamber (Fig. 3.6(a)), velocity is not evolving anymore, the jet from the oxidizer post breaks into the chamber and the exit nozzle is choked (Fig. 3.6(b)). The pressure in the chamber is homogeneous (Fig. 3.6(c)).

3.1.3.b Transition from cold choked flow to reacting choked flow

Two methods have been tested to initialize the chemical reaction:

• ignition of the whole chamber (IGNIT-A): all gas after the entrance of the chamber (x > 0) are replaced by burnt-gases (with burnt gases composition at adiabatic burnt gases temperature T ad b ) • ignition depending on the local mixing ratio (IGNIT-Z): only the zones of the chamber where the mixture fraction z is between 0.08 and 0.11 (the stoichiometric z z st = 0.095) are replaced by burnt-gases.

In both methods the velocity field obtained from the cold flow LES is conserved. It is useful to analyze both methods because they might lead to different flow states if the flow exhibits hysteresis (as shown experimentally). In IGNIT-Z, the combustor contains much more unburnt premixed gases and it is expected to lead to a more violent ignition process and therefore to instabilities if these instabilities are non linear and sensitive to the forcing level. : pressure traces in the chamber versus outlet flow rate, Eq. 1.12 without combustion, Eq. 1.12 for ignited cases.

(b) Zoom at the end of ignition, from t = 0.0190 s to t = 0.020 s. : pressure traces in the chamber versus outlet flow rate, Eq. 1.12 for ignited cases. Ignition in mixing ratio (IGNIT-Z) Figure 3.11 shows the field of temperature used to ignite the LES using the mixing ratio method. The red pockets (where 0.08 ≤ z ≤ 0.11) are burnt-gases with the burnt-gases composition at T = T ad b . Here the hot volume is small compared to the chamber volume. The evolutions of outlet mass flow rate ṁout and mean chamber pressure are plotted on Fig. 3.12. The outlet flow rate ṁout is continuous because the gases in the nozzle are not replaced by burnt gases. On both curves, the perturbation is large at the beginning of the ignition: the pressure goes up to more than 20 bar. Even though only a small pocket of gases is replaced by burnt-gases, the whole chamber is premixed and the flame propagates through it to burn all the available fuel.

.13 presents the evolution of heat release rate in the chamber from the first iteration (this why it does not start at zero). The increase of heat release rate is initially important: at the first time step the combustion releases 6.5 × 10 9 W/m 3 and goes up to 11 × 10 9 W/m 3 (corresponding to a chamber power of 74 MW for the first instants). This type of "violent" ignition is something which experimentalists try to avoid by limiting the amount of premixed gas in the chamber at sparking time. Numerically, of course, this is not a difficulty. Once the flow stabilizes one can observe oscillations of pressure, flow rate and heat release (Fig. 3.12) corresponding to the acoustic instabilities discussed in Chapter 4.

Starting from the point where ignition starts (t = 4 × 10 -3 s), pressure and flow rate increase simultaneously, and decrease following the "cold" theoretical line. This is due to the fact that cold gas are still flowing through the nozzle for a certain time (t < 4.8 × 10 burns all fresh gases in the chamber. When burnt gases reach the nozzle (t = 4.8 × 10 -3 s) the trajectory in (P , ṁout ) plane, moves to the "hot" choking line. Later, pressure and flow rate oscillate, due to the instability, but stay on the theoretical line for ignited flow. Figure 3.14(b) shows how pressure and flow rate oscillate during one instability cycle: variations are very large; ṁout goes from 0.3 to 0.42 kg/s while P goes from 13.8 to 19.5 bar. However, both quantities remain very close to the "hot" choking line suggesting that combustion remains complete during the oscillations.

(a) From ignition to t = 0.020 s.

: pressure traces in the chamber versus outlet flow rate, Eq. 1.12 without combustion, Eq. 1.12 for ignited cases.

(b) Zoom on one cycle of the instability. : pressure traces in the chamber versus outlet flow rate, Eq. 1.12 for ignited cases. parameters of the flow such as the mean pressure in the chamber (Fig. 3.16(a)), the mean heat release in the chamber (Fig. 3.16(b)), and the temperature at the end of the chamber (Fig. 3.17). Both simulations converge to the same mean flow parameters but lead to different oscillation levels.

The ignition in mixing ratio (IGNIT-Z) is violent for the flow, and initiates the acoustic instability (at a frequency of 1510 Hz which is close to the experimental value 1400 Hz [START_REF] Yu | Spontaneous longitudinal combustion instability in a continuously-variable resonance combustor[END_REF]) while the method IGNIT-A does not. This confirms that this flow is dominated by non linearities which are already evident in the fact that the steady state depends on the initial condition, something which is consistent with the existence of hysteresis in the experiment. Moreover it suggests that this flame may also have a non-linear response to acoustic perturbations (as described by Noiray et al. [22]). In the next section, this hysteresis phenomenon is investigated in more details.

3.1.3.c Analysis of non-linearities and hysteresis in the CVRC flame

The previous section has shown that initialization of the LES could lead either to a stable case (IGNIT-A) or to an unstable behavior (IGNIT-Z) at 1510 Hz. This is an interesting feature which requires, however, further investigation. For example, if the flow is sensitive to the level of perturbation, it should be possible to start from the stable case and force it to instability by a strong perturbation. This is done in this section.

Figure 3.18 shows the local pressure at probe P2 and its FFT for the "stable" regime reached after ignition Ignit-A. Even though this regime was called "stable", it exhibits a limited unsteady activity at high frequency; these fluctuations are due to a transverse mode at 11500 Hz with a small amplitude. The IGNIT-A "stable" flow is perturbed by stopping the combustion reaction for 1 ms at t = 31 ms and reactivating it at t = 32 ms. The objective is to produce a strong perturbation of the flow and to see whether it comes back to the "stable" IGNIT-A case or evolves towards another state.

The corresponding evolutions of the mean pressure in the chamber and heat release rate are plotted on Fig. 3.19. The flow does not restabilize and does not come back to the stable IGNIT-A regime: after a transition period, it becomes unstable and presents high pressure and heat release rate oscillations. Figure 3.20 shows the pressure perturbation at probe P2 and its FFT. There is a strong mode at 1500 Hz. There is also still a high frequency mode now at 13,000 Hz and with a smaller amplitude than the 11,500 Hz mode of the stable regime. The limit cycle for this state corresponds to the limit cycle obtained after an ignition IGNIT-Z (3.21).

This procedure (schemed on Fig. 3.22) demonstrates that the flame of the CVRC presents a non-linear response to perturbations [2]: a smooth ignition like IGNIT-A can lead to a stable flow. But if this flow is perturbed, it becomes unstable. Similarly, a violent ignition like IGNIT-Z directly leads to an unstable regime (Fig.

Numerical setup for the Helmholtz simulations

The previous section has shown how LES was set up to study instabilities in the CVRC. The LES study of the self-excited mode will be described in Chapter 4. In this section, we discuss how our second tool, the acoustic solver, is applied to the CVRC. The full acoustic results will be given in Chapter 5.

Boundary conditions for the acoustic solver

The geometry of the burner for the Helmholtz AVSP calculation can be simplified (Fig. 3.23) because only acoustic waves have to be tracked. The outlet nozzle will not be calculated, and will be modeled by an equivalent impedance Z out (Sec. 3.2.2). In the same way, the inlet nozzle is modeled by an impedance Z in . All other boundaries are considered as walls, even the fuel inlet because its very small volume has no impact on the acoustic modes.

Two geometries are used for the AVSP calculation, one with an oxidizer post length L op = 4.75 in (Fig. 3.23(a)), and another one with an oxidizer post length L op = 3.50 in (Fig. 3.23(b)) where the fuel lines are removed to simplify the mesh. The calculation time with AVSP is very sensitive to the number of nodes of the mesh and to the variation of cell size: removing the fuel line allows to change the calculation time from 18h (for the L op = 4.75 in mesh) to 5h (for the L op = 3.50 in mesh). While only one calculation per mode is made on the L op = 4.75 in case, seventeen calculations are made for one mode on the L op = 3.50 in case in order to explore all the different oxidizer post length through the impedance translation theorem and thus, the gain in calculation time becomes significant. The independence of the acoustic results to the presence of the fuel lines will be demonstrated in Sec. 5.3. The corresponding meshes are presented in Fig. 3.24: both meshes are unstructured and a particular attention has been paid to have minimum changes in cell size (∆ = 3 mm). For the L op = 4.75 in configuration, the mesh (Fig. 3.24(a)) is composed of 40,000 nodes (200,000 cells), the small fuel posts force to have very small cells. This is not the case for the L op = 3.50 in geometry (Fig. 3.24(b)) where the mesh (composed of 27,000 nodes) is more regular.

Impedance of the CVRC boundary conditions

All models described in Sec. 2.2.4 will be tested and compared for the calculation of the CVRC both with the (P , u ) and (J , m ) formulation and are summarized in Tab. 5.1 (Chapter 5).

3.2.2.a Inlet impedance

Figure 3.25 shows the piston located at the combustor inlet and its open surface: there is not one nozzle but several holes of various sizes and there is not a throat but a large zone with a constant area, leading to a over-expanded flow with small conical shocks outside each hole in the injector tube as shown in Fig. 3.26. In AVBP , the inlet boundary condition is an imposed mass flow rate (ρu = constant): it is this boundary condition which will be replaced by an impedance in AVSP . Assuming isentropic hypothesis, a constant mass flow rate leads to a theoretical impedance Z in = -1/ M . n . As, in AVSP , all normal vectors are entering, for this inlet the impedance is: This value is a singular point for the relation which transform (P , u ) variables in (J , m ) ones (Eq. 2.92): with the oriented normal entering the domain by the inlet. It leads to:

Z in (P ,u ) = - 1 M = -3.13 (3.1)
Z (J ,m ) = Z (P ,u ) + M M Z (P ,u ) + 1 (3.2)
Z in (J ,m ) = +∞. (3.3)
which is logical since m = 0 for a choked nozzle by definition. An infinite impedance can be set imposing zero velocity fluctuation at the boundary, but in practice for the CVRC this method can not be applied due to the impedance translation method to model the stretching of the oxidizer post. So an impedance Z in (J ,m ) = 10 6 will be imposed. Several tests have been performed to insure that this value is large enough to mimic an infinite impedance.

Chapter 3 : Application to the CVRC burner

3.2.2.b Outlet nozzle impedance

For the outlet impedance, the three different models presented in Sec. 2.2.4 will be tested in (P , u ) formalism and in (J , m ) formalism to account for Mach number effects as discussed by Motheau et al [START_REF] Motheau | Using boundary conditions to account for mean flow effects in a zero mach number acoustic solver[END_REF], using Eq. 2.92. All cases are summarized in Tab. 5.1.

The infinite (P , u ) impedance case will be imposed by setting a zero velocity fluctuation at the outlet nozzle throat.

The typical wave length of the acoustic mode found experimentally is 71 cm corresponding to a frequency of 1400 Hz. The length of the nozzle is 5 cm and can be assumed to be compact compared to the acoustic wave length following the definition of Marble and Candel [24]. Accounting for the Mach number in the outlet section of the chamber, before the nozzle, with a normal vector entering the domain (in the negative x direction), the expression from Marble and Candel becomes:

Z out (P ,u ) = -2 (γ -1) M = -60. (3.4)
Eq. 2.92 is now:

Z out (J ,m ) = Z out (P ,u ) -M -M Z out (P ,u ) + 1 = -5.02 (3.5)
Figure 3.27(a) displays the reflection coefficient R versus the frequency, calculated with Duràn's method [START_REF] Duran | Solution of the quasi one-dimensional linearized Euler equations using flow invariants and the Magnus expansion[END_REF] (Sec. 2.2.4) which is valid at all (longitudinal) frequencies. For f = 0 it recovers the Marble and Candel value. What is imposed in the Helmholtz solver is the impedance Z, linked to R by (for an outlet):

Z out = R + 1 R -1 (3.6)
As R is around one, the variation of Z out is very important (Fig. 3.27(b)) and may have an important effect on the calculation of the acoustic mode. The (J , m ) impedance is calculated using Eq. 2.92 oriented for an outlet (entering normal vector): This chapter focuses on the analysis of the unsteady reacting regime reached in the full 3D LES after the flow in the CVRC has reached a stationary state using an initialization method IGNIT-Z (Sec. 3.1.3). This state corresponds to an oscillating flow when L op = 4.75". Mean flow results are analyzed first (Sec. 4.1), before investigating flame structure (Sec. 4.2), unsteady regimes (Sec. 4.3), flame transfer function (Sec. 4.4), and studying unsteady modes (Sec. 4.5). The importance of heat losses is demonstrated in Sec. 4.6 by comparing an adiabatic case and one with heat losses.

Z out (J ,m ) (f ) = Z out (P ,u ) (f ) -M -M Z out (P ,u ) (f ) + 1 . ( 3

Mean flow results

Once a permanent oscillating regime has been reached, the flow is averaged to compute its mean properties. The averaging time is 18 ms, which corresponds to more than 25 cycles of 99 the instability at 1510 Hz (Sec. 4.3) and about 10 flow-through times of the chamber. The axial velocity field is presented in Fig. 4.1. The black iso-contour is the zero axial velocity, which highlights the recirculation zones. These recirculation zones are important for flame stabilization and their prediction requires a mesh refinement along the combustor walls. Small recirculations are observed in the wake of the lip between oxidizer and methane but as shown in Fig. 4.2, they are not strong enough to allow flame stabilization. There is a high jet stream in the chamber, mainly due to the high oxidizer injection speed. After the middle of the chamber the flow becomes homogeneous, before reaccelerating through the nozzle.

The mean temperature field is plotted in Fig. 4.2 together with a black iso-contour at T = 2000 K in order to represent the mean position of the turbulent flame brush. The flame is stabilized inside the combustion chamber at the location of the section change. The mean temperature field shows no diffusion flame between the fuel and oxidizer streams in the wake of the fuel injector. Based on the T = 2000 K iso-surface, the average flame length is L f = 15 cm. It is expanding outwards for the first 5 cm and then inwards with a conical shape.

Flame structure

To analyze the stabilization of the flame in the combustion chamber, an instantaneous longitudinal cut of the heat-release rate is presented in Fig. 4.3. Two iso-lines are superimposed: temperature T = 2000 K (in purple), which serves as a delimitation between fresh and burnt gases and mixture fraction z = z st (in black) to highlight the regions where fuel and oxidizer are in stoichiometric proportions. As expected, the stoichiometric line starts at the fuel injector lip and continues into the chamber. Meanwhile, as suggested by the averaged-temperature field (Fig. 4.2), the temperature line is anchored at the junction between the chamber and the injector. The field is distorted by turbulent structures.

This peculiar flame structure can be interpreted as a special triple flame [START_REF] Kioni | Experimental and numerical study of a triple flame[END_REF], produced by a triple layer: fuel, oxidizer and burnt gases (Fig. 4.4(b)). It indicates that chemistry is not sufficiently fast in the CVRC to stabilize the flame at the lips where fuel is injected. At these high velocities ; stabilization can be produced only by the burnt gases recirculating at the dump plane. this mixture is heated by the recirculating hot gases and the flame reaches a maximum intensity when the high temperature line crosses the stoichiometric line, creating a structure very similar to the BGS (Burnt Gas Supported) triple flame computed by Jiménez and Cuenot [START_REF] Jimenez | DNS study of stabilisation of turbulent triple flames by hot gases[END_REF]. This observation also explains why a high resolution is required to predict flame stabilization in this combustor: the mixing zone between fuel and oxidizer must be properly meshed but the whole chamber must also be well resolved to capture the recirculation zones that enable the stabilization by bringing back burnt gases near the injection plane. shows the pressure oscillation at the entrance of the chamber. Large pressure oscillations are visible with a peak to peak amplitude of 8 bar, which is quite similar to the oscillations observed in the experiment at the same operating point (L op = 4.75 in). The signal also contains a high frequency perturbation (larger than 10 kHz). Pressure perturbations are also recorded at the end of the chamber and presented in Fig. 4.5(b). The amplitude has decreased to 5 bar and is dominated by the low frequency component at 1510 Hz. The variations of total heat release rate (heat release integrated over the whole computational domain) are plotted on Fig. 4.6 and show strong oscillations. The heat release signal is non-linear, which is to be expected given the large amplitude of pressure oscillations. Even the pressure signal (Fig. 4.5(a)) exhibits harmonics and a slightly non-linear behavior, something which is rarely observed usually in atmospheric burners.

Unsteady activity 4.3.1 Temporal evolutions

A Fast Fourier Transform (FFT) is performed on the LES pressure signal and plotted on Fig. 4.7. The frequency of the fundamental mode is f = 1510 Hz. The frequency of the other peak are summed up in Tab. 4.1, their normalized frequencies show that the three first subharmonics are well reproduced, and the fifth and sixth subharmonics have a lower frequency than expected. The experimental investigation of the operating point L op = 4.75 in shows that a strong unstable acoustic mode occurs at a frequency f xp = 1400 Hz with an amplitude (power spectral density at the frequency of the instability measured at probe P2) between 1200 and 1500 kPa 2 .Hz -1 . As discussed in Sec. 1.3, the present simulation does not account for heat losses, resulting in a higher mean temperature in the burnt gases and a higher mean pres-sure upstream of the choked nozzle: in this simulation, the mean pressure is 15.1 bar, versus 13.4 bar in the experiment. Further simulations with heat losses matching the mean pressure of the experiment will be performed and presented in Sec. 4.6.

A set of pressure probes located every cm on the axis of the whole computational domain are used to plot the spatial structure of the unstable mode. By performing a FTF at each point, one can build the amplitude and phase of pressure fluctuations at the frequency of the instability. These results are presented in Fig. 4.8, where the origin of the x-coordinate is at the dump plane so that x < 0 corresponds to the oxidizer post and x > 0 to the combustion chamber.

Looking at the shape in the chamber (x > 0), the instability has a half-wave mode shape, with a shape in "V" for the pressure modulus (Fig. 4.8(a)): one pressure node is observed at the middle of the chamber (x = 0.21 m) and two maximums on the chamber extremities (x = 0 and x = 0.38 m), with a difference of π/2 for the phase between the two sides of the chamber (Fig. 4.8(b)). In the oxidizer tube (x < 0), the pressure modulus varies progressively and the phase variation is linear suggesting that only a propagative mode takes place in this tube.

Instantaneous fields

Instantaneous fields of pressure, temperature and heat release rate are recorded on one period of the instability ( 6.6 × 10 -4 s). Figure 4.9 presents the evolution of pressure at the probe P2 on three cycles with a marker for the times when an instantaneous field is recorded (one screenshot is taken every 6 × 10 -5 s). Figure 4.10 shows these instantaneous fields of pressure. In the combustion chamber, the snapshots of pressure (Fig. 4.10) allow to recognize the half-wave mode described previously. On frame 2, 8 and 13 the two extremities of the chamber have a high pressure variation and are opposite in phase. Around frames 5 and 10, there is no pressure variation along the chamber. At every instant, the middle of the chamber (x = 0.21 m) stays at the same pressure: it is a pressure node as shown by Fig. 4.8(a). In the oxidizer post, frames 1 to 10 show the negative pressure perturbation moving upstream in the post and confirm that this is a propagative mode which goes upstream. All variations occur along the longitudinal axis, with no radial perturbation. Figure 4.11 presents instantaneous fields of temperature with a black iso-line at T = 2000 K. The cold methane line (blue) injected at 300 K mixes with the oxidizer which is warmer ( 1000 K). The recirculation zones, very important for flame stabilization (Sec. 4.2) are cooled. The iso-line of temperature highlights the flame position: as expected from Sec. 4.1, the flame is always attached at the entrance of the chamber, mainly on the wall on the left of the chamber, but sometimes (frame 8), the flame flashes back and enters the injector for a few millimeters. The flame has a large flapping movement and its length varies a lot, from 5 cm on frame 8 to 10 cm on frame 11. Its width also changes: it is very thin on the first frame and occupies all the width of the chamber on frame 5. Cold pockets are released (for example between frames 3 and 4) and convected downstream, but they are burnt before the middle of the chamber. The temperature is homogeneous at the outlet nozzle: this is an important observation because Marble and Candel [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] have shown that variations in temperature (and thus entropy) passing through a nozzle can be an important source of acoustics. This does not seem to be the case here so that an instability scenario based on entropy spots convection is not likely in this configuration.

The instantaneous heat release rate is shown on Fig. 4.121 . The most intense flame is situated in the first half of the chamber. The maximum of reaction (in red) changes in position and intensity along the unstable cycle and the flame even goes upstream in the injector. On the first frame, the reaction zone spreads radially in the chamber, while on the 9 th frame it is more intense and located on the axis. This positive perturbation of heat release rate occurs at the same time as the positive perturbation of pressure which, through the Rayleigh criterion, could be a source of acoustic instabilities [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF].

Unsteady activity 109

The coupling between the pressure perturbation P and the heat release rate perturbation ω can be observed through the Rayleigh criterion which results from the mean acoustic energy E a over a cycle of the instability, defined as:

E a = 1 τ t+τ t V 1 2 ρu 2 + P 2 ρc 2 dv dt (4.1)
where u , P , and ρ are respectively the acoustic velocity, the pressure fluctuation and the mean density, c is the speed of sound and τ is the period of one unstable cycle. A conservation equation for the acoustic energy may then be derived [15,[START_REF] Nicoud | Thermoacoustic instabilities: should the rayleigh criterion be extended to include entropy changes?[END_REF]:

dE a dt = R -F (4.2)
where the Rayleigh source term, R, and the flux of acoustic energy through the boundaries, F , are defined by:

R = 1 τ t+τ t V γ -1 γP P ω dv dt (4.3)
F = 1 τ t+τ t Σ P u .n dσ dt (4.4)
where V is the volume of the domain and Σ its boundary; ω is the fluctuation of heat release rate and n the outward normal of Σ2 . The source term of acoustic energy due to the unsteady combustion, R, is then compared to the mean power output, P b , of the burner and reported in Table 4. and volume averaged Rayleigh term, R, is a measure of the overall contribution of the flame to the acoustic energy. In order to assess the impact of the various flame elements, a local time-averaged Rayleigh term, r, is defined as :

r = 1 τ t+τ t γ -1 γP P ω dt (4.5)
A longitudinal cut of r through the combustion chamber of the CVRC is presented in Fig. 4.13. The Rayleigh term is overwhelmingly positive showing that the flame is compact: all points in the combustor feed the instability. In gas turbines, where the flow velocity is much smaller, alternative positive and negative patterns usually appear on fields of Rayleigh criterion.

Measurement of flame response

Flame transfer functions are usually measured by forcing flames [15,[START_REF] Tay Wo Chong | Identification of flame transfer functions using les of turbulent reacting flows[END_REF][START_REF] Durox | Experimental analysis of nonlinear flame transfer functions for different flame geometries[END_REF][START_REF] Ducruix | Theoretical and experimental determinations of the transfer function of a laminar premixed flame[END_REF][START_REF] Balachandran | Experimental investigation of the nonlinear response of turbulent premixed flames to imposed inlet velocity oscillations[END_REF]. However, they can also be evaluated on a self-excited regime by correlating inlet velocity fluctuations and heat release rate oscillations. Of course, this result is valid only at the instability frequency (1510 Hz here). The FTF is calculated by recording heat release fields every 2 × 10 -5 s. The point [0,0,0] at the entrance of the chamber is chosen for the reference velocity, reference normal is oriented as the x-axis, and in each point of the domain the two parameters n (x, y, z) and τ (x, y, z) of the FTF are calculated as:

DFT [q (x, y, z)] DFT [u (x ref , y ref , z ref )] = n (x, y, z) e iωτ (x,y,z) (4.6)
were DFT is the direct Fourier transform operator. The result of this calculation is presented in Fig. 4.14. Both n (Fig. 4.14(a)) and τ (Fig. 4.14(b)) effect are localised on the zone of heat release, they represent only the effect of the acoustic field on the heat release rate of the flame. One can introduce a global n (n g ) and a global τ (τ g ) to characterize the whole flame. These two parameters have to verify:

n g e iωτg .V = V n (x, y, z) e iωτ (x,y,z) dV (4.7)

Eq. 4.7 can be used to evaluate τ g and n g , the global FTF parameters. Note that τ g can also be estimated measuring the mean interval between the velocity fluctuation at the reference point and the global heat release fluctuation (Fig. 4.15). For the CVRC flame, the global FTF parameters are: It is assumed that the flame response will not change with L op , this is a reasonable hypothesis, because neither the fuel (or oxidizer) flow rates nor the equivalence ratio change for the different lengths of L op , so that we can assume that the combustion regime remains the same.

n g =

Analysis of unsteady modes

Since LES and experiments predict that an unstable mode exists at 1510 Hz, it is interesting to see what the acoustic analysis predicts for this regime and AVSP was used for the 4"75 in mesh. The frequency of the unstable mode predicted by AVSP is f avsp = 1621 Hz, higher than the frequency found for the LES (f LES = 1510 Hz) and the experiment (f xp = 1400 Hz). This difference can be explained by two factors:

• the sound speed field used for the calculation in AVSP comes from the LES results which have a higher pressure and thus, a higher sound speed (Sec. 1.3),

• AVSP is a Helmholtz solver which works under low-Mach assumption. The real frequency (f M =0 ) must be lower than the frequency calculated under the low-Mach assumption (f M =0 ), the difference can be estimated by:

f M =0 = 1 -M 2 f M =0 (4.10)
Here the mean Mach number in the chamber is M = 0.18, which leads to a corrected frequency f M =0 avsp = 1568 Hz (which is still higher than the experimental frequency).

Even though the AVSP and AVBP frequencies do not match perfectly, a better method to check whether both codes actually capture the same acoustic mode is to compare the P model shapes obtained by both solvers (knowing that only the model shapes can be compared, amplitude cannot because AVSP is working in the linear regime). The structure of the unstable mode is compared for AVBP and AVSP on Fig. 4.16. Both mode structures have a similar shape, but the variation of pressure modulus in the oxidizer post (x < 0) which was small in the combustion chamber for the LES, is larger in AVSP . In the combustion chamber (x > 0), the "V" shape of the mode (typical of a half-wave mode) is reproduced by both codes: the pressure node is at the same position and the extrema are quite equal. Obviously all codes confirm that the 1500 Hz mode corresponds roughly to a half-wave mode of the chamber. Some differences exist however for the phase variation of the pressure (Fig. 4.16(b)) which is calculated with reference to the point x = 0 which has been set to zero. In the oxidizer post (x < 0), the global level of phase variation is the same for both codes, but where the phase evolves progressively in the LES for the propagative mode, AVSP presents a strong variation of π/2 corresponding to a half-wave mode. In the combustion chamber (x > 0), the phase jump of the half-mode is well recovered for both calculations.

Effects of heat losses

The previous sections have focused on the adiabatic LES of the 4.75 in case. This LES has one deficiency: the mean pressure is 15 bar while the experimental value is close to 13.5 bar. This indicates either incomplete combustion or non adiabatic flow (see Sec. 1.3.1). In this section, a simulation with heat losses is performed to investigate their effects on the flow. Heat losses are introduced with the method presented in Sec. 2.1.3 in order to match the mean chamber pressure of the experiment. Figure 4.17(a) presents the evolution of the mean chamber pressure after the introduction of heat losses at time t = 0.0117 s. The pressure decreases as expected and stabilizes to the target value P = P xp = 13.5 bar. The amount of heat losses needed to reach this pressure relative to the global flame power is plotted on Fig. 4.17(b). Heat losses have been tuned to bring the pressure back to 13.5 bar. They represent around 37% of the flame power: this a high value but not unrealistic. It also corresponds to the simple estimation performed in Sec. 1.3.1 (39%). : adiabatic, : with heat losses.

Figure 4.18(a) presents an average velocity field where the burnt gases recirculation is highlighted, and Fig. 4.18(b) shows the average temperature field. The adiabatic case is presented on the top figures, the case with heat losses on the bottom. There is no major change in the velocity field of the flow but the introduction of heat losses on the wall induces a cooling of the recirculating burnt gases: the flame lifts off and is no longer attached at the entrance of the chamber. The instantaneous fields of heat release rate for both simulations are presented in Fig. 4.19. This confirms the observation made on the temperature field: with heat losses, there is no combustion in the first 5 cm of the chamber. The flame does not move downstream and ends at the same abscissa than the adiabatic case but is shorter. The flame still has the triple flame structure described in Sec. 4.2. This flame modification has an important effect on the local Rayleigh term displayed on Fig. 4.20. As there is no more combustion in the beginning of the chamber ( ω = ω + ω = 0), the local Rayleigh term becomes null. The volume in the first five centimeters of the chamber where the flame amplifies the acoustic instability (red) in the adiabatic case has disappeared after the introduction of heat losses. Moreover, the damping zone (blue) has increased. This will have a major impact on the instability. Top: adiabatic case, bottom: case with heat losses.

Chapter 4 : LES and acoustic study of the 4.75" unstable case (dashed line) and the case with heat losses (solid line) after the introduction of the heat losses.

After six cycles the oscillations of heat release rate have vanished. It is important to notice that the value of the time-average mean heat release rate has not changed: even if the flame shape is different, the same quantity of fuel is burnt, the combustion regime is the same with or without heat losses. The pressure signal at probe P2 (Fig. 4.21(b)), confirms that the instability is damped after a few cycles.

Conclusion on heat losses :

Introducing heat losses has allowed to match the experimental mean pressure, but has also modified the flame/acoustic coupling from an unstable mode at 15 bar to a detached flame and a vanishing unstable mode at 13.5 bar. Some points in the modelling may be improved. For example only convective heat losses have been taken into account and implementation of radiative losses may be necessary. Of course, this is a disturbing result because it means that the present LES does not capture the operating pressure and the limit cycle observed experimentally at the same time. Considering the other results obtained on the CVRC (hysteresis, observed experimentally for example or dependance of modes on initial condition, observed numerically in Sec. 3.1.3), this is another manifestation of the difficulty of LES in these high pressure devices.

There is no doubt that a modification of other parameters in the LES with heat losses can lead to an unstable mode again but this path was not followed here. It was decided to investigate effects of post length on instability rather than exploring LES at a fixed post dimension.

In the following of this work, all references to the LES will be to the adiabatic unstable case.

We will not discuss this issue here, but we only point out that even if one remains in the framework of Eq. 4.1, losses by molecular effects should be included in a proper acoustic analysis to evaluate the complete acoustic energy budget. This is a difficult task.

Experimentally, losses can be measured (see work at EM2C for example [29]). In the present work, an evaluation of these losses was made by measuring them during the limit cycle of the self-excited mode studied in Chapter 4. The expression of the additional loss term F mol due to molecular effects is given for example by Nicoud and Poinsot [START_REF] Nicoud | Thermoacoustic instabilities: should the rayleigh criterion be extended to include entropy changes?[END_REF]:

F mol = - Ω λ γ -1 γP 0 ∇P 1 .∇T 1 dΩ - Ω - → τ 1 : - → ∇.u 1 dΩ (5.1)
The evolution of acoustic energy and F mol are plotted on Fig. 5.1 for six cycles of the instability.

The damping during the unstable cycle is α mol = F mol /E ac = -10 s -1 . Assuming that the losses are the same for all limit cycles for the different L op , this value of α mod can be subtracted to the amplification rate of the mode given by AVSP to have at least an estimation of the internal losses.

Boundary impedances

This configuration has two boundary conditions: the inlet (Z in ) and the outlet (Z out ) impedances. Their values will have a major impact on the stability. The different possible values of Z in and Z out have been described in Sec. The range of oxidizer-post length L op , tested experimentally by [27], will be explored using the impedance-translation theorem developed in Sec. 2.2.3 for an active flame. In order to validate this method, the same case L op = 4.75 in is calculated with AVSP on both a complete geometry with L comp op = 4.75 in and an inlet impedance Z in (Fig. 5.2(a)), and a cut geometry with L cut op = 3.50 in and an inlet impedance Z 0 calculated by the impedance-translation theorem applied to Z in . The inlet and outlet impedances are Z in = Z out = -10, corresponding to case 6 in Tab. 5.1.

The frequencies and amplification rates of the acoustic modes found by both calculations are presented in Tab. 5.2. For the first three modes, the frequencies and amplification rates are very similar, while significant differences appear for modes L4 to L6. It seems that when an acoustic node is close to the boundary where the impedance-translation theorem is applied, the methodology fails to adequately capture both the amplification rate and the frequency of the eigen-mode. The structure of the modes are compared in Fig. 5.3 with the reference results on the complete geometry in dashed blue line and on the cut geometry in solid black line. All the calculated modes have a longitudinal structure. The first three modes can be identified as:

Mode L1 is a combination of a 1/2-wave mode in the chamber and a 1/4-wave mode in the oxidizer post, without phase shift at the junction.

Mode L2 has the same structure as mode L1, except that there is a phase shift of π at the junction.

Mode L3 is the combination of a 1/2-wave in the oxidizer post and a full wave in the chamber.

The calculations give different results for the following modes but in the range of interest (0 < f < 4000 Hz), both AVSP configurations lead to the same conclusions: modes L1 and L3 are damped and can be discarded while Mode L2 is amplified and has a frequency of 1657 Hz close to the experimental value (1400 Hz) and is the mode which we need to track as done in the next sections.

Computation of acoustic modes for variable post length

All computations are now performed on the L op = 3.50 in case mesh (Fig. 3.24(b)) to be able to change L op continuously.

Without active flame

First, calculations are performed without active flame for which the nτ field is not given to AVSP . Figure 5.4 presents the frequency of the second mode and its amplification rate. For all cases the frequency of the second mode varies in the same way within 40 Hz. Cases 1 and 2 give exactly the same frequency variation, such as cases 4 and 6 which give the higher estimation of the frequency. As expected for a case without active flame, the second mode is never amplified (Fig. 5.4(b)). Cases with (P , u ) formalism (cases 1, 2 and 3) have the same evolution for the amplification rate increasing with L op . Case 6 evolves in the same way but in a shorter range, and the cases with (J, m) formalism (cases 4 and 5) evolve in the opposite way; the amplification rate decreases with L op .

With active flame

The frequency of the second mode and its amplification rate are plotted on Fig. 5.5. Introducing the active flame model does not change the frequency of the second mode (Fig. 5.5(a)) but leads to amplified modes over a range of post lengths L op . The frequency has the same variations for all cases. The amplification rates (Fig. 5.5(b)), exhibit complex variations with L op . All curves have shapes which correspond to the observed experimental behavior; amplification is lower on the extrema values of L op and higher for the central value (around L op = 5 in = 12.7 cm) see Fig. 1.5 which shows that the pressure oscillations are maximum for 12 < L op < 14 cm. The slope of the curve is higher for low values of L op than for large ones: this could be an explanation to the large hysteresis zone observed on Fig. 1.5.

Cases 1 to 3 are very similar. Low values of L op lead to damped modes. The amplification rate increases with L op to a maximum for L op = 5.5 in but after it not decreases a lot and remains positive for large L op , even when accounting for acoustic losses F mol (Sec. 5.1). Case 4 and 5 have a better shape with an amplification rate increasing and decreasing after, but the small L op are too amplified. With theses impedances, the second mode is only damped for L op > 6.00 in. Case 6 has been tuned to have an amplification rate superior than the molecular losses for 4.00 < L op < 6.50 in and negative for L op < 4.00 in and L op > 6.50 in. The values of the impedances Z in and Z out do not correspond to theoretical values, but as discussed in Sec. 3.2.2, the shape of the inlet and outlet nozzle are particular and their impedances can not be determined easily.

Results show that this problem is highly dependent on the two parameters Z in and Z out : small variations of their values lead to different stable or unstable states. This is why the question of uncertainty quantification is so a important topic in the actual research projects. To tend to reliable CFD study, it is necessary to be able of quantifying the sensitivity of numerical codes to margin of error in the input parameters. In the present case, impedances are obviously the most uncertain parameters. : -α mol , × : case 1, : case 2, : case 3, + : case 4, : case 5, • : case 6.

Influence of the oxidizer-post length

In this section, the variation of the three acoustic modes L1 to L3 with L op is studied for case 6 which is the most representative of the experiment. For these variations of L op , the modes frequencies do not cross, eliminating possibilities of mode "switching" as discussed by Boudy et al [2] for example.

The amplification rates of these three modes are plotted on Fig. 5.7. The active flame has no influence on mode L3, while it damps the mode L1 specially for low L op and amplifies the mode L2.

Figure 5.7(b) confirms that only the second mode L2 can lead to acoustic instabilities with the shape described previously for P : modes L1 and L3 are always damped. It is important to point out that the flame transfer function has been calculated for only one frequency, close to the frequency of the two first modes, and thus the precision of the impact of this FTF on the third mode at very different frequency is questionable. The variation of L op has no impact on the amplification of L3 while the damping of L1 decreases with increasing L op . • : Mode L1, × : Mode L2 and : Mode L3.

The modulus of P for various L op for the three modes are plotted on Fig. 5.8 and their phases are presented on Fig. 5.9. All the modulus are set to be equal to one at x = 0 and the phases are set to π at the end of the chamber. Mode L2 presents no evolution at all when L op increases while modes L1 and L3 are impacted by the variation of L op . For L2, the structure of the mode in the chamber remains exactly the same for all L op . In the oxidizer post, the maximum value of |P | is at the entrance of the tube, as L op increases, this inlet moves away from the boundary of the computational domain and only the end of this quarter-wave mode shape remains visible. For L1, the value of the modulus of P at the inlet increases with L op while the position of the pressure node moves upstream. For L3, nothing changes while L op remains inferior at 5.00 in. At L op = 6.00 in, the structure of the mode in the oxidizer post changes: the phase at the entrance jumps of π/2, a maximum of pressure appears and increases until L op = 7.00 in, and decreases afterwards. In the chamber, the amplitude of the modulus of pressure evolves in the same way and the pressure nodes moves upstream.

Conclusion:

The impedance translation theorem has been applied on the CVRC's geometry and the first six longitudinal acoustic modes are identified. The only mode which is not damped and can lead to instability is the second one called L2. Various inlet and outlet impedances corresponding to different theories for modelling the inlet and outlet nozzle are tested. The best reproduction of the experimental curve is obtained with tuned impedance of whom values are between the different theoretical values. Introducing the molecular loss, the experimental stability map is well mimicked, but the frequency of the mode is still over estimated.

0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 4.00" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 5.00" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 6.00" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 6.50" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 7.00" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] L op = 7.50" 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m] 0 0.1 0.2 0.3 -π -π/2 0 π/2 π x [m]

Conclusions

The general objective of this thesis was to extend the understanding of combustion instabilities in high pressure devices with numerical simulation and propose a method to analyze the influence of geometric variation at low computational cost.

A reference test case has been identified in the unstable experiment (called CVRC) operated at Purdue University. It comprises two nozzles at the inlet and at the outlet in order to have a high pressure operating point more representative of the industrial configuration. A sliding oxidizer post allows to change the acoustic eigenmode of the combustion chamber and leads to stable and unstable regimes.

After an adaptation of the kinetic scheme for the extreme thermodynamic conditions of this experiment, the numerical code AVBP was validated for both acoustic and hydrodynamic modes resolution.

Different methods for initialization and ignition of the simulation for such nozzle ended devices were studied. Starting from a zero velocity state, the acceleration of the flow until the outlet nozzle choked was compared to theory. The ignition scenarios showed that the CVRC flame has a non-linear behavior. Two ignition methods were presented: the ignition of the whole chamber lead to a stable state while the ignition in limited close-to-stoichiometric zones led to an unstable state. Adding a perturbation to a stable flow made it become unstable.

Combustion instabilities were then studied with LES for a fixed oxidizer post length. A particular triple flame structure stabilized by burnt gases recirculation was observed. Variation of pressure and heat release rate were studied on a unstable cycle. The LES gave a good reproduction of the acoustic mode present in the experiment but with a slightly higher frequency. Introducing heat losses led to a better reproduction of the experimental flow conditions, but also induced a damping of the combustion instability. The mechanism of the instability was studied through the Rayleigh criterion. Its structure was analyzed and the flame transfer function was computed to feed the acoustic solver with an nτ model.

The Helmholtz solver AVSP was then used to analyze the eigenmode of the combustor. The stability results are highly dependent on the inlet and outlet impedances and several cases corresponding to different theories for nozzle impedances were tested. Good results were obtained for a tuned impedance which was a compromise between the different theoretical values. The impedance translation theorem was implemented in the acoustic solver AVSP to mimic at low computation cost the variation of the oxidizer post length. The experimental stability map was well reproduced.

The combination of the LES with the Helmholtz solver provides a good low cost tools to study the stability of a complex configuration, but the prediction in term of frequency and amplitude of the instability still lack precision. This article presents a simulation of a model rocket combustor with continuously variable acoustic properties thanks to a variable-length injector tube. Fully compressible Large-Eddy Simulations are conducted using the AVBP code. An original flame stabilization mechanism is uncovered where the recirculation of hot gases in the corner recirculation zone creates a triple flame structure. An unstable operating point is then chosen to investigate the mechanism of the instability. The simulations are compared to experimental results in terms of frequency and mode structure. Two-dimensional axi-symmetric computations are compared to full 3D simulations in order to assess the validity of the axi-symmetry assumption for the prediction of mean and unsteady features of this flow. Despite the inaccuracies inherent to the 2D description of a turbulent flow, for this configuration and the particular operating point investigated, the axi-symmetric simulation qualitatively reproduces some features of the instability.

© 2012 Académie des sciences. Published by Elsevier Masson SAS. All rights reserved.

r é s u m é

Cet article présente la simulation d'un injecteur expérimental de type moteur-fusée dont les propriétés acoustiques peuvent varier continûment grâce à un tube d'injection de longueur variable. Des simulations aux grandes échelles (LES) sont réalisées à l'aide du code AVBP. Un mécanisme de stabilisation de la flamme original montrant le rôle important de la recirculation de gaz chauds en entrée de chambre est observé. Un point de fonctionnement instable est ensuite choisi pour étudier le mécanisme de l'instabilité. Les simulations numériques sont comparées aux mesures du banc d'essai en termes de fréquence de l'instabilité et de structure du mode instable. Une simulation 2D axisymétrique est comparée à une simulation 3D avec pour objectif de tester l'hypothèse de symétrie pour la prédiction de l'écoulement moyen et instationnaire de cette configuration. Bien que la description en 2D d'un coulement turbulent soit une simplification importante, pour cette configuration et dans le cas du point de fonctionnement choisi,

Introduction

Most high-performance propulsion systems exhibit high-amplitude pressure oscillations at some stage in their design process. Unsteady heat release from combustion can couple with acoustic eigenmodes of the combustion chamber: this phenomenon is called combustion instability [1][2][3]. These instabilities can lead to important stresses and heat fluxes on the chamber walls and injectors, sometimes even resulting in the destruction of the engine [4,5]. Combustion instabilities are very difficult to predict because of the wide range of physical phenomena involved (acoustics, turbulence, chemical reactions, etc.) and they are often discovered at an advanced phase of the conception program.

A trial-and-error methodology with incremental modifications of the engine is usually necessary to recover stability but the associated cost is prohibitive. Therefore, the prediction of combustion instabilities through numerical simulation is a major issue. For gas turbine or aeronautic propellers, several studies have already been conducted [6][7][8][9][10][11][12] but for rocket engines the number of unsteady simulations is smaller [13][14][15][16][17][18]. This could be explained by the extreme thermodynamic conditions in rocket engines and the wide range of length scales: the characteristic length of the injector is a fraction of a millimeter while the flame length can be several dozens of centimeters.

In terms of numerical framework, the best way to study combustion instabilities and all the phenomena involved would be compressible Direct Numerical Simulation (DNS) with the full resolution of the Navier-Stokes equations. However, this method is inapplicable because of the prohibitive number of mesh points necessary to resolve the complex geometry and the high Reynolds-number flow of a rocket combustion chamber. Large-Eddy Simulation (LES) is a good candidate to study these instabilities and it has been used successfully in the field of gas turbines [19][20][21][22] for example to capture one mode of instability. What is more difficult to do is to predict transition from stable to unstable flow which is the ultimate objective in many practical cases: it is unclear today whether LES can predict stability limits. Moreover, the cost of full 3D LES for multiple operating points remains high. For geometries which are axi-symmetric, the question of performing only 2D axisymmetric LES to explore stability maps much more rapidly becomes interesting. The objective of this paper is to explore this point with a model rocket engine by performing both 2D axi-symmetric and full 3D LES and comparing results to experimental data. The CVRC setup operated at Purdue University [23] is axi-symmetric and exhibits very strong unstable modes. This system works at high pressure and has a variable injector duct length which allows to explore stability limits in a continuous manner. The experimental setup and operating point are described in Section 2. The numerical method is then presented in Section 3 together with the computational domain and boundary conditions. Finally, the results of the numerical simulations are presented in Section 4.

Experimental setup

Description

The burner studied in this paper was designed and operated at Purdue University (Indiana, USA) [23][24][25]. The system is called Continuously Variable Resonance Chamber (CVRC) because the length of the oxidizer injector can be varied continuously, allowing for a detailed investigation of the coupling between the acoustics and combustion in the chamber. The CVRC, presented in Fig. 1, is an axi-symmetric burner with a 15-inch long (38.1 cm) combustion chamber closed by a choked nozzle. The reactants are fed into the chamber by a coaxial injector: oxidizer at the center and fuel at the periphery.

For a continuous tuning of the acoustic response of the CVRC, the oxidizer-post length (L op ) can vary from L op = 3.5 in. (8.89 cm) to L op = 7.5 in. (19.05 cm), with respect to the inlet of the combustion chamber. The inlet of the oxidizer post is also choked so that the acoustic boundary conditions are well defined. Six pressure transducers are installed along the burner to measure high-frequency pressure fluctuations: the first one is located 3 in upstream of the injector face, and the remaining sensors are placed in the combustion chamber at 0.5, 1.5, 2.5, 3.5 and 14.5 in downstream of the injector face. These locations are respectively called P 1 to P 6 in the remainder of the paper.

Operating point

The fuel is pure gaseous methane injected at a rate of ṁ f = 0.0272 kg s -1 and temperature T f = 300 K. The oxidizer is a mixture of 90% hydrogen peroxide and 10% water (per mass unit), which is fully decomposed: this is equivalent to 42% oxygen and 58% water. The oxidizer is injected in the chamber at T ox = 1030 K (both water and oxygen are in gaseous phase) and ṁox = 0.3180 kg s -1 . Consequently, the global equivalence ratio is φ = 0.813.

Accounting for all the species in the GRI-Mech 3.0 mechanism (52 species), the adiabatic burnt-gases temperature at equilibrium is T ad b = 2512 K. Since the total flow rate going through the choked nozzle is known ( ṁ = ṁox + ṁfuel = 0.3452 kg s -1 ) assuming the chamber temperature, T c , is equal to the adiabatic burnt-gases temperature at equilibrium T ad b , the chamber pressure P c can be determined using:

ṁ = c d A P 0 √ rT 0 γ 1 2 2 γ + 1 γ +1 2(γ -1) (1)
The stagnation pressure and temperature, respectively P 0 and T 0 are given by:

P 0 = P c 1 + γ -1 2 M 2 γ γ -1 ; T 0 = T c 1 + γ -1 2 M 2 ( 2
)
where M is the Mach number of the burnt gases in the chamber, A the area of the choked nozzle, γ the ratio of heat capacities, r the perfect gas constant and c d the discharge coefficient. Using the simulations presented in Section 4 show that the combustion is complete (there is no unburnt fuel exiting the nozzle), this reduced temperature is attributed to thermal losses at the combustor walls. This approximate evaluation corresponds to a total heat loss at the walls equal to 25% of the available power (Table 3). Because the present work does not include heat losses, in the LES, the burnt-gases temperature is T ad b and the chamber pressure is 1.54 MPa (Table 1). The resulting increase in the speed of sound yields a small shift of the unstable mode frequency, as discussed in Section 4.3.

Stability

The experimental stability map (amplitude of pressure fluctuations versus L op ) is presented in Fig. 2 for a number of runs where the oxidizer-post is either translated forward (L op increases), backward or kept constant during the whole test. Three different behaviors of the system can be distinguished: (i) Stable regions for L op < 3.75 in (9.5 cm) and L op > 7.1 in (18 cm). (ii) An unstable range for 3.75 < L op < 6.3 in (9.5 to 16 cm). (iii) A zone of hysteresis for 6.3 < L op < 7.1 in (16 to 18 cm).

Because one the objectives of the present work is to compare 2D axi-symmetric and 3D computations to compute an unstable operating point, the simulations are conducted for L op = 4.75 in (∼ 12 cm). This value of L op is also chosen to minimize potential hysteresis effects in the instability as all tests show strong pressure fluctuations under these conditions.

Numerical method

Solver and models

The simulations are carried out with the AVBP code developed at CERFACS [6,7,26] that solves the compressible Navier-Stokes equations on unstructured grids. The numerical scheme is the second-order in space and time Lax-Wendroff formulation [27]. Turbulence is modeled using the classical Smagorinsky model [28]. The flame/turbulence interaction is accounted for by the Dynamically Thickened Flame (DTF) model [29]. In this approach, the flame is artificially thickened to be resolved on the LES mesh and subgrid reaction rate is recovered through an efficiency function [30].

The CH 4 /O 2 kinetics are modeled using a reduced two-step mechanism (BFER [31]) that accounts for six species (CH 4 , O 2 , CO 2 , CO, H 2 O and N 2 ) and two reactions:

CH 4 + 3 2 O 2 -→ CO + 2H 2 O CO + 1 2 O 2 ←→ CO 2 (3)
The first equation is irreversible while the second one is reversible and leads to an equilibrium between CO and CO 2 in the burnt gases. This scheme was developed for CH 4 /air flame and used here without modification.

A common shortcoming of two-step mechanisms is that burning rates are over-predicted for rich mixtures. Following the procedure presented in [31], the present scheme uses a pre-exponential constant of the first reaction that is a function of the local equivalence ratio so as to recover accurate burning rates in rich mixtures. This scheme was validated for the thermodynamic conditions of the CVRC versus the Gri-Mech 3.0 mechanism. First, equilibrium states were compared to ensure that no important species would be neglected. Table 2 presents the burnt-gases properties for both schemes at the stoichiometric mixture fraction z st = 0.095 as well as the average mixture fraction of the experimental operating point z xp = 0.079 (calculated from the mass flow rates ṁ f and ṁox ). Even though it is expected that the combustion regime is a diffusion flame, the premixed laminar flame speeds at stoichiometry are compared: s L = 1.59 m s -1 for Gri-Mech 3.0 and s L = 1.52 m s -1 for the present two-step scheme. It was also checked that auto-ignition times for Gri-Mech like for the 2-step scheme are much larger than the residence time in the chamber. This is consistent with the fact that under the conditions of the CVRC experiment, combustion is not controlled by auto-ignition. Overall, these calculations validate the approximations implied by the use of the reduced scheme.

Computational domain and boundary conditions

With the advent of massively parallel computers and codes that scale on thousands of processors, three-dimensional (3D) simulations of turbulent combustion in complex geometries are now quite common in the research community. The AVBP solver is particularly attune for such computations [21,32,33]. However, because the background idea of the use of numerical simulations is to provide a quick and relatively inexpensive evaluation of the stability of a given configuration, we opted for a different approach: it was decided to compare the accuracy of a 3D computation using significant resources and a 'light' axi-symmetric computation feasible on lost-cost computer clusters. This study is performed for a single unstable operating point. The typical shortcoming of axi-symmetric computations is that turbulence being essentially a 3D phenomenon, the 2D assumption may yield incorrect vortex dynamics and the symmetry condition that imposes zero turbulent fluxes on the axis is inadequate. The consequences of these approximations will be addressed in Section 4.

Axi-symmetric computations were performed on a relatively fine mesh using 7 × 10 5 cells while three-dimensional computations use a mesh of 14 × 10 6 cells (Fig. 3). Because of the variety of length scales in the geometry (oxidizer tube diameter is 4 cm whereas fuel injector thickness is 0.4 mm) unstructured meshes are required. Cell sizes were chosen to ensure a correct resolution of the dynamic of the flow in the mixing zone and of the chemical and heat transfer of the flame downstream the injector. The computational time for 1 ms of physical time (that is more than one cycle of the instability) is 5 CPU hours on 32 processors of the SGI Altix ICE 8200 of CINES (Intel Quad-Core E5472) for the axi-symmetric mesh and 2 CPU hours on 512 processors of the same machine for the 3D mesh (6.4 times more expensive). The computational domain starts slightly downstream of the choked inlet of the oxidizer post and ends after the exit nozzle, which allows for well defined acoustic boundary conditions. Fuel is injected through an annular slit corresponding to the experiment but the tubes feeding the slit are not included. The boundary conditions are prescribed using the NSCBC technique [34]: mass flow rate and temperature are imposed at the oxidizer ( ṁox = 0.32 kg s -1 and T ox = 1030 K) and fuel ( ṁ f = 0.027 kg s -1 and T f = 300 K) inlets while nothing needs to be imposed at the supersonic outlet. Walls are treated as no-slip and adiabatic.

Results and discussion

This section is devoted to the comparison of the 2D axi-symmetric and 3D computations as well as the validation against available experimental data. Mean flow fields are compared in Section 4.1 and the flame-stabilization and dynamics are presented in Section 4.2. Finally Section 4.3 presents the analysis of the acoustic mode of this operating point in terms of structure and mechanism of the instability.

Mean flow field

Once a permanent regime has been reached in both computations, the flow is averaged to compute the mean properties. For both simulations, the averaging time is 18 ms, which corresponds to more than 25 cycles of the instability (cf. Section 4.3) and about 10 flow-through times of the chamber. The axial velocity field is presented in Fig. 4 where the upper and lower parts show respectively the 2D and 3D simulation. The black iso-contour is the zero axial velocity, which highlights the recirculation zones. The flow fields are quite different: first, the axi-symmetric simulation presents a longer corner-recirculation zone but it also maintains a larger axial velocity close to the axis. These differences are not due to the grid resolution as a significant effort was done to resolve the recirculation zones (Fig. 3). Even though the no-slip boundary is clearly a rough condition that could be replaced by an adequate wall-law model, the 3D simulation with larger cells at the wall does a better job in the second half of the chamber. The mean temperature field is plotted in Fig. 5 together with a black iso-contour at T = 2000 K in order to represent the mean position of the turbulent flame brush. There is a striking difference between the two simulations as the flame length is over-estimated by a factor of three in the 2D computation. This feature is typical of axi-symmetric computations where the boundary condition on the axis imposes zero turbulent fluctuations. The flapping movement of the flame on the axis is impossible resulting in a longer average flame. This elongated flame could be the origin of the distorted axial velocity field. The 3D computation also has a lower mean temperature in the corner recirculation zones that is attributed to a wider flapping motion of the flame (cf. Section 4.2).

Flame stabilization and dynamics

For the analysis of the stabilization of the flame in the combustion chamber, an instantaneous longitudinal cut of the heat-release rate is presented in Fig. 6. Two iso-lines are superimposed: temperature T = 2000 K (in white) as a delimitation between fresh and burnt gases and mixture fraction z = z st (in black) to highlight the regions where the fuel and oxidizer are in stoichiometric proportions. In both simulations, the stoichiometric iso-line (black lines) starts at the fuel injector and continues into the chamber while the temperature iso-line (white lines) is anchored at the junction between the chamber and the injector. Both fields are distorted by turbulent structures but the convolution is more intense in the 2D case, which is typical of such simulations where the dynamics of the vorticity is massively affected by the 2D assumption. As inferred from the mean fields (Fig. 5), the lack of turbulent fluctuations on the axis of the 2D computation affects the dynamics of the flame front resulting in a much longer flame. In the 2D case, the mixture fraction and temperature isocontours are much closer but the structure of the reaction layer is similar in both computations: most of the reaction rate is sitting on the temperature iso-contour and extremely intense heat release rate occurs where the two contours cross. This peculiar flame structure can be seen as a special triple flame [35], produced by a triple layer: fuel, oxidizer and burnt gases (Fig. 7(b)). It shows that chemistry is not sufficiently fast in the CVRC to stabilize the flame at the lips where fuel is injected. Stabilization is produced by the burnt gases recirculating at the dump plane. Fig. 7 illustrates the mechanism suggested by these results: oxidizer and fuel start mixing at the lips where methane is injected; no combustion exists in this zone. Downstream, this mixture is heated by the recirculating hot gases and the flame reaches a maximum intensity when the high temperature line crosses the stoichiometric line, creating a structure very similar to the BGS (Burnt Gas Supported) triple flame computed by Jiménez and Cuenot [36]. This observation also explains why a very high resolution is required to predict flame stabilization here: the mixing zone between fuel and oxidizer must be properly meshed but the whole chamber must also be well resolved to capture the recirculation zones which are controlling the stabilization.

Analysis of the unstable acoustic mode

The experimental investigation of the operating point L op = 4.75 in shows that a strong unstable acoustic mode occurs at a frequency F xp = 1400 Hz with an amplitude (power spectral density at the frequency of the instability measured at probe P 2 ) between 1200 and 1500 kPa 2 Hz -1 . We will now test the ability of the numerical simulations to recover these characteristics. The temporal evolution of the pressure signal at probe P 2 is plotted in Fig. 8(a). The signal seems mildly non-linear and high-frequency harmonics are visible on top of the main unstable mode. The frequency of the fundamental is F 2D = 1572 Hz for the 2D axi-symmetric computation and F 3D = 1510 Hz for the 3D case. The discrepancy between the simulations and the experiment (F xp ≈ 1400 Hz) is caused by the different values of the mean pressure, resulting in a slight shift of the speed of sound. In terms of amplitude, the root mean square (rms) of the 3D signal is twice as much as that of the 2D computation. This difference will be explained later in this section with the analysis of the acoustic energy equation. Then the space-averaged heat release rate versus time is presented in Fig. 8(b). Unlike the pressure signal, the heat release response is non-linear, which is to be expected for finite amplitude combustion instabilities. The harmonic content of the 2D simulation is wider than that of the 3D case, which may be related to the overall higher small-scale wrinkling of the flame surface shown in Fig. 6. The spatial structure of the unstable mode (Fig. 9) was obtained by performing a power spectral density analysis of the whole field at the frequency of the instability. The longitudinal variation of the modulus and phase of the pressure fluctuations are plotted in Figs. 9(a) and 9(b), respectively. It appears that in the combustion chamber (x > 0), the instability has a half-wave mode shape, while in the oxidizer tube, the pressure fluctuations are mainly propagative. Both simulations predict nearly identical structures, with a minor deviation in the modulus around the junction between the injector and the chamber. The mean acoustic energy E a over a cycle of the instability is defined as:

E a = 1 τ t+τ t V 1 2 ρu 2 + p 2 ρc 2 dv dt (4)
where u , p , and ρ are respectively the acoustic velocity, the pressure fluctuation and the mean density, c is the speed of sound and τ is the period of one unstable cycle. A conservation equation for the acoustic energy may then be derived [20,37]:

∂ E a ∂t = R -F ( 5
)
where the Rayleigh source term, R, and the flux of acoustic energy through the boundaries, F , are defined by: R =

) 6
where V is the volume of the domain and Σ its boundary; ω is the fluctuation of heat release rate and n the outward normal of Σ .

The source term of acoustic energy due to the unsteady combustion, R, is then compared to the mean power output, P b , of the burner. The values reported in Table 3 indicate that the Rayleigh term is much larger in the 3D case, as expected from the amplitude of the pressure fluctuations presented in Fig. 8(a).

The value of the global Rayleigh term, R, is a measure of the overall contribution of the flame to the acoustic energy. In order to assess the impact of the various portions of the flame, a local time-averaged Rayleigh term, r, is defined as: A longitudinal cut of r through the combustion chamber of the CVRC is presented in Fig. 10. Despite being similar in size and shape, there are two major differences between the 2D and 3D simulations. First, in the 3D case, the Rayleigh term is overwhelmingly positive while a significant zone of negative r is observed in the recirculation zone of the 2D computation. Then, on the axis of the 2D computation, the lack of turbulent fluctuations drives r toward zero, which is not the case in 3D.

Finally, it seems that the longer tail of the flame seen in 2D (for x > 10 cm) does not contribute significantly to the acoustic energy. This may explain why in term of prediction of the instability, the 2D computation yields results similar to the 3D case.

Conclusion

An unstable operating point of a model rocket combustor with continuously variable acoustic properties is simulated using both low-cost 2D axi-symmetric and detailed 3D calculations. The axi-symmetric simulation does not predict the correct mean properties of the flow (velocity, the flame length, etc.) or the flapping movement of the flame due to the boundary condition on the axis that imposes zero turbulent fluctuations. An unusual structure of the stabilization flame is observed: the flame looks like a distorted triple flame with a strong rich premixed branch and a weak diffusion tail, stabilized by the recirculation of burnt gases at the entrance of the combustion chamber. Nevertheless, both simulations exhibit a combustion instability with virtually identical frequency and structure of the acoustic mode. The frequency of the instability compares favorably with the experimental data. The contribution of the flame to acoustic energy is larger for the 3D geometry than for the 2D one, resulting in a larger amplitude of the pressure fluctuation (7.4 bars instead of 4.7 bars). Despite the differences in the flow field, for this configuration and this operating point, the 2D axi-symmetric simulation qualitatively reproduces some features of the instability.

( a )Figure 1 :

 a1 Figure 1: Vulcain rocket engine ended by a nozzle (SNECMA).

Figure 2 :

 2 Figure 2: Experimental bench with atmospheric outlet for the study of ignition in swirl flows and acoustic instabilities (EM2C [13]).

  (a) Atmospheric combustor. Quarter wave mode. (b) Choked outlet combustor. Half wave mode.

Figure 3 :

 3 Figure 3: Two different types of combustor with the modulus of pressure perturbation.

20 Chapter 1 :

 201 Description of the experimental setup (a) Experimental setup installed at Purdue University: the CVRC. Figure courtesy of Pr. W. Anderson, Purdue University. (b) Scheme of the CVRC.

Figure 1 . 1 :

 11 Figure 1.1: Geometry of the experimental setup CVRC.

Figure 1 . 2 :

 12 Figure 1.2: Details of the CVRC experiment. Figure courtesy of Pr. W. Anderson, Purdue University.

Figure 1 . 3 :

 13 Figure 1.3: Computation of adiabatic flame temperature

Figure 1 . 1 : 1 Figure 1 . 4 :

 11114 Figure 1.4: Equilibrium adiabatic temperature of the burnt gases T ad b versus mixture fraction z calculated with complete GRI-Mech mechanism.: mixing temperature, : Z xp , : Z st .

Figure 1 . 5 :

 15 Figure 1.5: Stability map of the CVRC combustor (experimental results). Power spectral density of the pressure fluctuation at the frequency of the instability measured by probe P 2 (x=3.81 cm). : forward translation of the oxidizer post, : backward translation and • : fixed oxidizer-post length. Figure courtesy of Pr. W. Anderson, Purdue University.

Figure 1 . 6 :Chapter 1 :

 161 Figure 1.6: Temporal evolution of: : pressure at the 6 th pressure probe and : oxidizer post length L op .

Figure 1 . 7 :

 17 Figure 1.7: Geometry of the CVRC in the numerical work of Harvazinski et al. [50]

  = 5.50" in fuel lean regime in the DES of Harvazinski et al. [50] power spectral density (PSD) are plotted on Fig. 1.8.

Figure 1 . 8 :

 18 Figure 1.8: PSD analysis for the case L op = 5.50" in fuel lean regime in the DES of Harvazinski et al. [50]

Chapter 1 :

 1 Description of the experimental setup

Figure 1 . 10 : 1 . 5 33 (

 1101533 Figure 1.10: Instantaneous fields in the 3D DES of Harvazinski et al. [50]

Figure 1 . 11 :

 111 Figure 1.11: Mode shape of the first longitudinal mode in the 3D DES of Harvazinski et al. [50]

Figure 2 . 1 :

 21 Figure 2.1: Numerical methods to study combustion instabilities.

Figure 2 . 2 :

 22 Figure 2.2: Conceptual representation of LES applied to a homogeneous isotropic turbulent field.

5 Figure 2 . 3 :

 523 Figure 2.3: Direct Numerical Simulation of fame/turbulence interactions with iso-level of reaction rate and vorticity fields [15].

Figure 2 . 4 :

 24 Figure 2.4: Equilibrium adiabatic temperature T ad b vs mixture fraction z. × : 2S-CM2, • : GRI-Mech.

Figure 2 . 5 :

 25 Figure 2.5: Scatter plot of heat release versus mixture fraction z in the LES.

Figure 2 . 6 :

 26 Figure 2.6: Scheme of the heat losses boundary condition in AVBP .

Figure 2 . 7 :

 27 Figure 2.7: 1D box for study of resonant modes with a representation of the initialized pressure and velocity.

Figure 2 . 8 :

 28 Figure 2.8: Instantaneous fields of perturbation of pressure and velocity for a mesh with 10 points per wave length. • : Lax-Wendroff, × : TTGC, : t=0, : t = 15T , : t = 30T and : t = 45T where T = 2L/c is the period of the acoustic wave. The mode corresponds to the first acoustic eigenmode of the duct. (λ/∆x = 20).

Figure 2 . 9 :

 29 Figure 2.9: Damping of acoustic eigenmode in a 1D closed box as a function of the number of points per wave length (λ/∆x). o : Lax-Wendroff, x : TTGC.

Figure 2 .

 2 Figure 2.10(a) is a schematic representation of a general 2D shear layer flow, which is a good prototype of the shear zones in a combustor where vortices are created. The growth of vortices in a sheared flow depends strongly on the assumption used for the mean flow profile[START_REF] Betchov | Stability of parallel flows[END_REF][START_REF] Drazin | Hydrodynamic stability[END_REF][START_REF] Okong | Real-gas effects on mean flow and temporal stability of binaryspecies mixing layers[END_REF]. This profile can be simplified as an infinite shear flow (Fig.2.10(b)), or less roughly by a linear velocity profile with a transition thickness d (Fig.2.10(c)). It is that case which is calculated here with AVBP . The simulations are initialized by introducing a perturbation corresponding to the solution of the infinite shear case (Sec. 2.1.4.d) which has a simple analytical expression ; the growth rate of the instability is then compared to its analytical value.

Figure 2 . 10 :

 210 Figure 2.10: Configurations of shear flow

Figure 2 . 11 :

 211 Figure 2.11: Initial fields of horizontal velocity, vertical velocity and pressure perturbation for hydrodynamic mode study.

2. 1

 1 .4.f Comparison between AVBP and linear theory at fixed resolution Simulations are performed in a two-dimensional rectangular domain where left and right boundary conditions are periodic. The upper and lower limits of the domain correspond to a symmetry condition. The gas has the same composition as in Sec. 2.1.4.a, with a upper velocity U 1 = 5 m/s and a lower velocity U 2 = -5 m/s. The transition between these two velocities is linear with a thickness d =0.1 m (Fig. 2.10(c)). A perturbation is added to study the linear growth of the instability. This perturbation corresponds to the analytical solutions of a shear flow with infinite shear (d = 0) (Fig.2.11 and 2.12). It would have been more precise to use the eigenfunction corresponding to the linear profile of Fig.2.10(c) but they are more complex and play a limited role on the results. The growth rate for example becomes rapidly independent of the initial conditions:

Figure 2 . 12 :

 212 Figure 2.12: Initial field of horizontal velocity u for hydrodynamic mode study.

  (v)) (b) λ = 0.45 Lax-Wendroff

Figure 2 . 13 :

 213 Figure 2.13: Temporal evolution of 1/2 log E c (v), the two vertical lines represent the limits of the zone where linear interpolation is used.

Figure 2 . 15 :

 215 Figure 2.15: Growth rate of the hydrodynamic instability of a shear layer versus reduced wave length αd.

Figure 2 . 16 :

 216 Figure 2.16: Resolution of hydrodynamics. Variation of σ/σ lin as a function of axial (λ/∆x) and transverse resolution (d/∆y), with the ratio λ/d = 8.

  (a) Simply converging nozzle (b) Converging-diverging nozzle

Figure 2 . 17 :

 217 Figure 2.17: Schemes of the two types of nozzle used in LES.

  (a) Over-expanded flow (b) Under-expanded flow

Figure 2 . 18 :

 218 Figure 2.18: Schemes of the two types of nozzle used in LES.The variation of pressure along the nozzle axis for various outlet pressure P 2 is summarized in Fig.2.19.

Figure 2 . 19 :

 219 Figure 2.19: Static pressure through a nozzle for different imposed outlet pressure P 2

Figure 2 . 20 :

 220 Figure 2.20: Scheme of the CVRC geometry to be solved with an acoustic solver.

Figure 2 . 21 :

 221 Figure 2.21: N AVSP computations on N different meshes. Boundary conditions are imposed at x = -L op and x = l.

Figure 2 . 22 :

 222 Figure 2.22: N AVSP computations on one single mesh with a variable impedance imposed at x = -l. Boundary conditions are imposed at x = -l (its impedance Z (x = -l) changes with L op ) and x = L.

Figure 2 . 23 :

 223 Figure 2.23: Scheme of an academic one-dimensional tube with imposed inlet impedance and zero velocity fluctuation condition at outlet. The modulus of acoustic pressure of a half-wave mode is represented.

Figure 2 . 24 :

 224 Figure 2.24: Geometry of the test case for mode computing with a variable impedance at inlet.

Figure 2 . 25 : 2 . 2

 22522 Figure 2.25: Longitudinal and radial cut of the test case tube mesh.

2 0

 2 2.27). The same comparison is made for Z 0 =

Figure 2 . 26 :

 226 Figure 2.26: First acoustic mode shape for Z 0 = 0.5.: analytical expression, × : AVSP .

Figure 2 . 27 :Chapter 2 :

 2272 Figure 2.27: Second acoustic mode shape for Z 0 = 0.5. : analytical expression, × : AVSP .and presented in Fig.2.28 and 2.29. A good agreement is found between AVSP and the theory.We can now model the first part of the tube by the method of impedance translation.

Figure 2 . 28 :

 228 Figure 2.28: First acoustic mode shape for Z 0 = 2.0.: analytical expression, × : AVSP .

Figure 2 . 29 :

 229 Figure 2.29: Second acoustic mode shape for Z 0 = 2.0.: analytical expression, × : AVSP .

Fig. 2 .

 2 32 shows an example of the first mode's

Figure 2 . 30 :

 230 Figure 2.30: Variation of frequency f = kc/ (2π) in function of l for Z -l = 0.5.: theoretical 1 st mode, × : AVSP for the 1 st mode, : theoretical 2 nd mode, • : AVSP for the 2 nd mode.

Figure 2 . 31 :

 231 Figure 2.31: Variation of frequency f = kc/ (2π) in function of l for Z -l = 2.0.: theoretical 1 st mode, × : AVSP for the 1 st mode, : theoretical 2 nd mode, • : AVSP for the 2 nd mode.

Figure 2 . 32 : 2 . 2 71 -

 2322271 Figure 2.32: First acoustic mode shape for Z -l = 0.5 and l = 0.5.: analytical expression, × : AVSP .

Figure 2 . 33 :

 233 Figure 2.33: First acoustic mode shape for Z -l = 2.0 and l = 0.5.

Figure 2 . 34 :

 234 Figure 2.34: Imposed impedance at the nozzle's throat.

Figure 2 . 35 :

 235 Figure 2.35: Imposed impedance before the nozzle

Figure 3 . 1 :

 31 Figure 3.1: Detailed view of the mesh.

Figure 3 . 2 :

 32 Figure 3.2: Boundary condition for AVBP calculation.

Figure 3 .

 3 3 presents instantaneous fields 2 × 10 -4 s after the initialization. Methane has entered and begins to diffuse in the chamber (Fig. 3.3(a)). The velocity increases in the oxidizer post due to the injection of oxidizer and in the nozzle due to the low pressure condition at the outlet which drain the chamber (Fig. 3.3(b)). The nozzle is not choked yet. The pressure highly increases in the oxidizer post and decreases in the nozzle (Fig. 3.3(c)). Figure 3.4 shows (a) Mass fraction of CH4 (b) Velocity (c) Pressure

Figure 3 . 3 :

 33 Figure 3.3: Instantaneous field in a longitudinal cut of the CVRC at t = 2 × 10 -4 s after initialization.the evolution of flow rate through the nozzle (Fig.3.4(a)) and mean pressure in the chamber (Fig.3.4(b)). The nozzle first chokes at t = 1.8 × 10 -4 s and the mean pressure stops increasing.

Figure 3 . 4 :

 34 Figure 3.4: Initialization at t = 0 to ignition at t = 4.0 × 10 -3 s (no combustion).

Figure

  Figure 3.4 shows that pressure and outlet flow rate both vary strongly with time and demonstrates the need to follow the flow state in a (Pressure -Flow rate) diagram and to be able to determine when the nozzle is choked or not. Equation 1.12 gives the relation between the mean pressure in the chamber and the flow rate through a choked nozzle. Evolutions of chamber pressure versus flow rate are plotted on Fig.3.5(a)1 . The pressure first decreases slowly while the flow rate hardly increases, over-passing the theoretical value of the choked nozzle before going back to it at t = 3.8 × 10 -4 s. The next portion of this path is presented in Fig.3.5(b). After a loop at high pressure (t = 7.3 × 10 -4 s), the path joins the theoretical line, and from here, pressure and flow rate decrease following the theoretical relation. From Eq. 1.12, for the injection flow rate at the mixing temperature T 0 mix = 924 K (Sec. 1.2), one can predict that a simultaneous decrease of both P and ṁout will take place until the chamber reaches a pressure of 8 bar (Fig.3.4(b)). In practice, there is no need to wait for this stabilization and ignition is activated at t = 4 ms.Figure3.6 presents instantaneous fields of this cold flow before ignition. The methane has propagated in the chamber (Fig.3.6(a)), velocity is not evolving anymore, the jet from the oxidizer post breaks into the chamber and the exit nozzle is choked (Fig.3.6(b)). The pressure in the chamber is homogeneous (Fig.3.6(c)).

( a )

 a Figure 3.5:: pressure traces in the chamber versus outlet flow rate, theoretical relation (Eq. 1.12) for non reacting case.

  (a) Mass fraction of CH4 with two isolines z = 0.08 and z = 0.11 (b) Velocity (c) Pressure

Figure 3 . 6 :

 36 Figure 3.6: Instantaneous field in a longitudinal cut of the CVRC before ignition, at t = 4 × 10 -3 s (see Fig. 3.4).

Figure 3 . 9 :

 39 Figure 3.9: INGIT-A: average heat release rate in the chamber

  (a) From ignition to t = 0.020 s.

Figure 3 . 10 :

 310 Figure 3.10: Trajectories of flow state in a (chamber pressure vs outlet flow rate) diagram for IGNIT-A.

Figure 3 . 11 :

 311 Figure 3.11: Longitudinal cut of temperature field used for ignition in mixing ratio (IGNIT-Z).

Figure 3 . 12 :

 312 Figure 3.12: IGNIT-Z

  Figure3.13 presents the evolution of heat release rate in the chamber from the first iteration (this why it does not start at zero). The increase of heat release rate is initially important: at the first time step the combustion releases 6.5 × 10 9 W/m 3 and goes up to 11 × 10 9 W/m 3 (corresponding to a chamber power of 74 MW for the first instants). This type of "violent" ignition is something which experimentalists try to avoid by limiting the amount of premixed gas in the chamber at sparking time. Numerically, of course, this is not a difficulty. Once the flow stabilizes one can observe oscillations of pressure, flow rate and heat release (Fig.3.12) corresponding to the acoustic instabilities discussed in Chapter 4.Starting from the point where ignition starts (t = 4 × 10 -3 s), pressure and flow rate increase simultaneously, and decrease following the "cold" theoretical line. This is due to the fact that cold gas are still flowing through the nozzle for a certain time (t < 4.8 × 10 -3 s) until the flame

Figure 3 . 13 :

 313 Figure 3.13: IGNIT-Z: average heat release in the chamber

Figure 3 . 14 :Chapter 3 :Figure 3 . 15 :

 3143315 Figure 3.14: Trajectories of flow state in a (chamber pressure vs outlet flow rate) diagram for IGNIT-Z

3 ]

 3 (b) Mean heat release in the chamber

Figure 3

 3 Figure 3.16: : cold flow, : IGNIT-A, : IGNIT-Z

Figure 3 . 17 :Chapter 3 :

 3173 Figure 3.17: Temperature at the end of the chamber. : cold flow, : IGNIT-A, : IGNIT-Z

fFigure 3 . 18 :

 318 Figure 3.18: Stable regime reached after IGNIT-A

1

 1 Numerical setup for LES simulations 91 0.03 0.032 0.034 0.036 0.038 0.04 0Mean pressure in the chamber 0.03 0.032 0.034 0.036 0.038 0.04 0Mean heat release in the chamber

Figure 3 . 19 :

 319 Figure 3.19: Perturbation of stable regime by switching of the combustion reaction.

  FFT of the pressure signal from probe P2

Figure 3 . 20 :fFigure 3 . 21 :Figure 3 . 22 : 2

 3203213222 Figure 3.20: Steady state reached for IGNIT-A after a perturbation (Fig. 3.19)

Figure 3 . 23 :

 323 Figure 3.23: Boundary condition for AVSP calculation.

Figure 3 . 24 :

 324 Figure 3.24: Longitudinal cut of AVSP meshes.

Figure 3 . 25 :

 325 Figure 3.25: View of the inlet piston.

Figure 3 . 26 :

 326 Figure 3.26: Shock surfaces (purple) downstream the oxidizer slots and streamlines colored with pressure (from Harvazinski et al. [50]).

Figure 3 . 27 :

 327 Figure 3.27: Outlet boundary condition with Duràn's method [39].

Figure 4 . 1 :

 41 Figure 4.1: Average field of axial velocity in a longitudinal cut of the CVRC with u = 0 iso-line highlighting the recirculation zones.

Figure 4 . 2 :

 42 Figure 4.2: Average temperature field with T = 2000 K iso-line to highlight the mean position of the turbulent flame brush.

4. 2 Figure 4 . 3 :

 243 Figure 4.3: Instantaneous field of heat release rate [W.m -3 ] during the unstable cycle, : stoichiometric mixture fraction (z = z st ); : temperature T = 2000 K.

Figure 4 .

 4 4 illustrates the mechanism suggested by these results: oxidizer and fuel start mixing at the lips where methane is injected; no combustion exists in this zone. Downstream, (a) Flame stabilization and structure for the CVRC experiment. (b) BGS (Burnt Gas Supported triple flame).

Figure 4 . 4 :

 44 Figure 4.4: Schematic representation of the flame stabilization and structure for the CVRC experiment and a BGS (Burnt Gas Supported) triple flame.

Figure 4 .

 4 Figure 4.5 shows the evolution of the pressure signal recorded on the axis at the entrance of the chamber (1 cm downstream the injector, Fig.4.5(a)), and at the end of the chamber just before the outlet nozzle (Fig.4.5(b)). The time step in this simulation is around 3.5 × 10 -8 s and the signal is recorded every 100 time steps, leading to an acquisition frequency of 285 kHz ensuring a very high temporal resolution of the instability expected at a frequency around 1400 Hz.

Figure 4 . 5 :

 45 Figure 4.5: Evolution of P at the entrance and at the end of the chamber (just upstream of nozzle).

Figure 4 .

 4 Figure 4.5(a) shows the pressure oscillation at the entrance of the chamber. Large pressure oscillations are visible with a peak to peak amplitude of 8 bar, which is quite similar to the oscillations observed in the experiment at the same operating point (L op = 4.75 in). The signal also contains a high frequency perturbation (larger than 10 kHz). Pressure perturbations are also recorded at the end of the chamber and presented in Fig.4.5(b). The amplitude has decreased to 5 bar and is dominated by the low frequency component at 1510 Hz. The variations of total heat release rate (heat release integrated over the whole computational domain) are plotted on Fig.4.6 and show strong oscillations. The heat release signal is non-linear, which is to be expected given the large amplitude of pressure oscillations. Even the pressure signal (Fig.4.5(a)) exhibits harmonics and a slightly non-linear behavior, something which is rarely observed usually in atmospheric burners.A Fast Fourier Transform (FFT) is performed on the LES pressure signal and plotted on Fig.4.7. The frequency of the fundamental mode is f = 1510 Hz. The frequency of the other peak are summed up in Tab. 4.1, their normalized frequencies show that the three first subharmonics are well reproduced, and the fifth and sixth subharmonics have a lower frequency than expected. The experimental investigation of the operating point L op = 4.75 in shows that a strong unstable acoustic mode occurs at a frequency f xp = 1400 Hz with an amplitude (power spectral density at the frequency of the instability measured at probe P2) between 1200 and 1500 kPa 2 .Hz -1 . As discussed in Sec. 1.3, the present simulation does not account for heat

3 ]Figure 4 . 6 :Figure 4 . 7 :

 34647 Figure 4.6: Evolution of heat release rate in the chamber.

  Phase of P .

Figure 4 . 8 :

 48 Figure 4.8: Structure of the acoustic field obtained from the LES of the self-excited mode at 1510 Hz.

Figure 4 . 9 :

 49 Figure 4.9: Evolution of pressure at probe P2, • : saved snapshots of Fig. 4.10 to 4.12.

Figure 4 . 10 :

 410 Figure 4.10: Instantaneous field of pressure in [Pa] with a time step of 60 µs.: P = 15.5 bar iso-line.

Chapter 4 :

 4 LES and acoustic study of the 4.75" unstable case

Figure 4 . 12 :

 412 Figure 4.12: Instantaneous field of heat release rate [W.m -3 ] with a time step of 60 µs.

2 . 6 Table 4 . 2 :

 2642 The value of the time R [W] P b [W] 8.23.10 3 1.34.10 Time-averaged value of the Rayleigh source term, R, of acoustic energy compared to total burner power output, P b .

Figure 4 . 13 :

 413 Figure 4.13: Time-averaged field of local Rayleigh term, r.

  10, 525 J/m (calculated by Eq. 4.7) (4.8) τ g = 2.89 × 10 -4 s (estimated on Fig. 4.15) (4.9)

•

  (a) Field of n [J/m 4 ]. • (b) Field of τ [s].

Figure 4 . 14 :

 414 Figure 4.14: Flame transfer function of the CVRC calculated from instantaneous LES fields. • : reference point.

Figure 4 . 15 :

 415 Figure 4.15: Estimation of the global phase τ g for the FTF of the CVRC. : velocity fluctuation at reference point, : global heat release fluctuation.

4. 5 -

 5 Analysis of unsteady modes113 Phase of P .

Figure 4 . 16 :

 416 Figure 4.16: Structure of the acoustic field. • : AVBP , × : AVSP

  Heat losses / Flame power (b) Temporal evolution of the ratio between heat losses and flame power.

Figure 4

 4 Figure 4.17:: adiabatic, : with heat losses.

Figure 4 .

 4 Figure 4.18(a) presents an average velocity field where the burnt gases recirculation is highlighted, and Fig.4.18(b) shows the average temperature field. The adiabatic case is presented on the top figures, the case with heat losses on the bottom. There is no major change in the velocity field of the flow but the introduction of heat losses on the wall induces a cooling of the recirculating burnt gases: the flame lifts off and is no longer attached at the entrance of the chamber. The instantaneous fields of heat release rate for both simulations are presented in Fig.4.19. This confirms the observation made on the temperature field: with heat losses, there is no combustion in the first 5 cm of the chamber. The flame does not move downstream and ends at the same abscissa than the adiabatic case but is shorter. The flame still has the triple flame structure described in Sec. 4.2. This flame modification has an important effect on the local Rayleigh term displayed on Fig.4.20. As there is no more combustion in the beginning of the chamber ( ω = ω + ω = 0), the local Rayleigh term becomes null. The volume in the first five centimeters of the chamber where the flame amplifies the acoustic instability (red) in the adiabatic case has disappeared after the introduction of heat losses. Moreover, the damping zone (blue) has increased. This will have a major impact on the instability. Figure 4.21(a) compares the mean heat release rate in the chamber for the adiabatic case

  (a) Time-average field of axial velocity in a longitudinal cut of the CVRC with u = 0 iso-line to highlight the recirculation zones. (b) Time-average field of temperature with T = 2000 K iso-line to highlight the mean position of the turbulent flame brush.

Figure 4 . 18 :

 418 Figure 4.18: Top: adiabatic case, bottom: case with heat losses.

Figure 4 . 19 :

 419 Figure 4.19: Instantaneous field of heat release rate [W.m -3 ] during the unstable cycle, : stoichiometric mixture fraction (z = z st ); : temperature T = 2000 K. Top: adiabatic case, bottom: case with heat losses.

Figure 4 . 20 :

 420 Figure 4.20: Time-average field of local Rayleigh term.Top: adiabatic case, bottom: case with heat losses.

  Temporal evolution of mean heat release rate. Temporal evolution of pressure at the probe P2 (cf. Fig.1.1)

Figure 4

 4 Figure 4.21:: adiabatic, : with heat losses.

Figure 5 . 1 :

 51 Figure 5.1: Acoustic losses measured in the LES of Chapter 4.

3 . 2 . 2 , 4 Z 8 (Z out = - 10 Table 5 . 1 :

 322481051 they are summed up in Tab. 5.1. The impedances of cases 1 to 3 come from the theoretical value in (p, u) formulation, starting from the simplest model of no fluctuation (Case 1), then using the Marble and Candel model (Case 2), and the Duràn model (Case 3) discussed in Sec. 3.2.2. In cases 4 and 5, the (J, m) formulation is used for an outlet model with model of Marble and Candel or with model of Duràn[START_REF] Duran | Solution of the quasi one-dimensional linearized Euler equations using flow invariants and the Magnus expansion[END_REF]. Case 6 is an arbitrary set of values allowing a good reproduction of the experiment. The reflection coefficient R is calculated as:R = Z + 1 Z -1 (5.2)The stability maps for each case are presented without active flame in Sec. 5.4.1 and with active flame in Sec. 5.4out (f ) = Z Duran (f, M (x) , A (x))ex: Z P ,u out (1600 Hz) = -1.7 -i8.8 (at 1600 Hz)Case out (f ) -M o -M o Z P ,u out (f ) + 1 ex: Z J,mout (1600 Hz) = 4.4 -i2.List of inlet (Z in ) and outlet (Z out ) impedances for the AVSP calculation and their corresponding reflexion coefficient (R in and R out ). The Mach number at the inlet is M i = 0.32, and at the outlet, M o = 0.18.5.3Validation of the method on the case L op = 4.75 in

  (a) Complete geometry. (b) Cut geometry with impedance translation theorem.

Figure 5 . 2 :

 52 Figure 5.2: Geometries of computational domain for the validation of the impedance translation theorem on the CVRC.

Chapter 5 :

 5 Influence of the oxidizer-post length on the stability

  Phase of P .

Figure 5 . 3 : 5 . 4

 5354 Figure 5.3: Structure of the acoustic field.× : AVSP on mesh with L op = 4.75 in, : AVSP on short mesh with the method of translation of impedance for L op = 4.75 in with Z in = -10 and Z out = -10 (case 6 in Tab. 5.1).

  Frequency of mode 2 versus oxidizer post length Lop.

  Amplification rate of mode 2 versus oxidizer post length Lop.

Figure 5 . 4 :

 54 Figure 5.4: Mode 2 without active flame. × ; case 1, : case 2, : case 3, + : case 4, : case 5, • : case 6.

  Amplification of mode 2 versus oxidizer post length Lop.

Figure 5 . 5 :

 55 Figure 5.5: Mode 2 with active flame.: -α mol , × : case 1, : case 2, : case 3, + : case 4, : case 5, • : case 6.

Figure 5 .

 5 6 presents the evolution of frequency of these three modes without (Fig. 5.6(a)) and with active flame (Fig. 5.6(b)). Accounting for the active flame has no impact on the frequency of the three modes. The frequency (Fig. 5.6) of mode L1 evolves linearly from 1313 Hz for L op = 3.50 in to 875 Hz for L op = 7.50 in. Mode L2 varies from 1952 Hz to 1481 Hz and mode L3 from 2926 Hz to 2613 Hz.

  (a) Without active flame.

L

  op [in] Re(f ) [Hz] (b) With active flame.

Figure 5 . 6 :

 56 Figure 5.6: Frequency versus oxidizer post length L op without active flame for case 6 in Tab. 5.1.• : Mode L1, × : Mode L2 and : Mode L3.

5. 5 50 L

 550 Influence of the oxidizer-post length op [in] Im(f ) [s -1 ] (b) With active flame, : -α mol .

Figure 5 . 7 :

 57 Figure 5.7: Amplification rate versus oxidizer post length L op with active flame for case 6 in Tab. 5.1.• : Mode L1, × : Mode L2 and : Mode L3.

Chapter 5 :

 5 Influence of the oxidizer-post length on the stability

Figure 5 . 8 :

 58 Figure 5.8: Variation of modulus of P with L op for modes L1, L2 and L3.

5. 5

 5 Influence of the oxidizer-post length

Figure 5 . 9 :

 59 Figure 5.9: Phase variation of P with L op for modes L1, L2 and L3.

Fig. 1 .

 1 Fig. 1. Experimental setup installed at Purdue University: the CVRC [25]. Figure courtesy of Pr. W. Anderson, Purdue University.

  c d = 1, T c = T ad b , γ = 1.182 and r = 389.2 (computed from the burnt-gases composition), one predicts a theoretical mean chamber pressure P th c = 1.543 MPa. The mean pressure actually measured in the experiment is P xp c = 1.34 MPa, which implies that heat losses or incomplete combustion lead to a lower burnt-gases temperature: T c < T ad b . Using Eq. (1) backwards in order to obtain P c = P xp c , yields T xp c = 1904 K. Because

Fig. 2 .

 2 Fig. 2. Stability map of the CVRC combustor (experimental results [25]). Power spectral density of the pressure fluctuation at the frequency of the instability measured by probe P 2 (x = 3.81 cm). : forward translation of the oxidizer post, : backward translation and : fixed oxidizer-post length. Figure courtesy of Pr. W. Anderson, Purdue University.

Fig. 3 .

 3 Fig. 3. Detailed view of the two meshes.

Fig. 4 .

 4 Fig. 4. Average field of axial velocity in a longitudinal cut of the CVRC with u = 0 iso-line highlighting the recirculation zones. Top: 2D axi-symmetric, bottom: 3D.

Fig. 5 .

 5 Fig. 5. Average temperature field with T = 2000 K iso-line to highlight the mean position of the turbulent flame brush. Top: 2D axi-symmetric, bottom: 3D.

Fig. 6 .

 6 Fig. 6. Instantaneous field of heat release rate (W m -3 ). Black iso-line: stoichiometric mixture fraction (z = z st ); white iso-line: temperature T = 2000 K. Top: 2D axi-symmetric, bottom: 3D.

Fig. 7 .

 7 Fig. 7. Schematic representation of the flame stabilization and structure for the CVRC experiment and a BGS triple flame.

Fig. 8 .

 8 Fig. 8. ---3D, --axi-symmetric.

Fig. 9 .

 9 Fig. 9. Structure of the acoustic field, ×: 3D, •: 2D axi-symmetric.

Fig. 10 .

 10 Fig. 10. Time-averaged field of local Rayleigh term, r. Top: 2D axi-symmetric, bottom: 3D.

Table 1 . 1 :

 11 Fresh-gases properties in thermodynamic conditions of the CVRC.

		095	881.3	0.095 0.380 0.525	0	0	0
	0.813	z xp = 0.079	924.3	0.079 0.387 0.534	0	0	0

Table 1 . 3

 13 Previous numerical studies of the CVRC setup 29 which gives a good reproduction of the complete kinetic GRI-MECH model in lean conditions but drastically increases the computational cost. 3D simulation were not affordable with the four-step model, so it is only performed with the one-step model. 2D simulations are performed with and without the four-step model. Table1.3 presents the frequency of the first six modes in 2D axi-symmetric with both kinetic schemes and 3D with the 1-step kinetic scheme compared to experimental results. The corresponding

	Mode Experiment 3D 1-step 2D 1-step 2D 4-step
	1	1350 Hz	1550 Hz	1650 Hz	1460 Hz
	2	2679 Hz	3100 Hz	2810 Hz	2650 Hz
	3	4045 Hz	4650 Hz	4280 Hz	4170 Hz
	4	5390 Hz	6250 Hz	5980 Hz	5555Hz
	5	6750 Hz	7800 Hz	7100 Hz	-
	6	8090 Hz	9350 Hz	-	-

20)

1.5 : Frequency identification for the case L op

Table 2 . 1 :

 21 Mech 2661 < 10 -4 2.49.10 -2 7.12.10 -1 2.22.10 -1 2.45.10 -2 0 2-step 2781 < 10 -4 2.25.10 -2 7.39.10 -1 1.99.10 -1 3.94.10 -2 0 -4 7.50.10 -2 6.99.10 -1 2.06.10 -1 6.17.10 -3 0 2-step 2586 < 10 -4 7.81.10 -2 7.11.10 -1 2.02.10 -1 8.62.10 -3 0 Burnt-gases properties in thermodynamic conditions of the CVRC.

	0.813	z xp = 0.079	Gri-Mech	2512	< 10

Table 2 . 2 :

 22 Parameters used for the tests of Fig.2.13 and 2.14.

	s)

  .91) Chapter 3 : Application to the CVRC burner 2.2.1 Description of the AVSP code . . . . . . . . . . . . . . . . . . . . . . . 63 2.2.2 Validation case : variable length tube . . . . . . . . . . . . . . . . . . . 63 2.2.2.a Analytical expression . . . . . . . . . . . . . . . . . . . . . . . 63 2.2.2.b Numerical configuration . . . . . . . . . . . . . . . . . . . . . . 65 2.2.2.c Comparison with theory for a fixed impedance case Z (x = 0) = Z 0 = [0.5, 2.0] . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 2.2.3 Impedance of a variable length tube: the impedance translation theorem 69 2.2.3.a Validation on the tube test case . . . . . . . . . . . . . . . . . 69 2.2.4 Impedance of nozzles . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 2.2.5 Effects of Mach number . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

  .7)

	Chapter 4
	LES and acoustic study of the 4.75"
	unstable case

Contents 3.1 Numerical setup for LES simulations . . . . . . . . . . . . . . . . . . 76 3.1.1 Geometry and mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 3.1.2 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 3.1.3 Initialization procedures for LES of choked flows . . . . . . . . . . . . . 77 3.1.3.a Initialization of non-reacting choked flows . . . . . . . . . . . . 77 3.1.3.b Transition from cold choked flow to reacting choked flow . . . 79 Ignition of the whole chamber (IGNIT-A) . . . . . . . . . . . . . 82 Ignition in mixing ratio (IGNIT-Z) . . . . . . . . . . . . . . . . . 85 Conclusion on ignition methods: . . . . . . . . . . . . . . . . . . 88 3.1.3.c Analysis of non-linearities and hysteresis in the CVRC flame . 90 3.2 Numerical setup for the Helmholtz simulations . . . . . . . . . . . . 93 3.2.1 Boundary conditions for the acoustic solver . . . . . . . . . . . . . . . . 93 3.2.2 Impedance of the CVRC boundary conditions . . . . . . . . . . . . . . 93 3.2.2.a Inlet impedance . . . . . . . . . . . . . . . . . . . . . . . . . . 93 3.2.2.b Outlet nozzle impedance . . . . . . . . . . . . . . . . . . . . . 96

Table 4 . 1 :

 41 Frequency identification for the case L op = 4.75" with AVBP .

Table 5 . 2 :

 52 Frequency identification for the case L op = 4.75" with active flame.

  Contents lists available at SciVerse ScienceDirectComptes Rendus Mecaniquewww.sciencedirect.com Combustion, flow and spray dynamics for aerospace propulsion Large-Eddy Simulation of combustion instabilities in a variable-length combustor Simulation aux grandes échelles des instabilités de combustion dans un brûleur de longueur variable Romain Garby a,b, * , Laurent Selle a,b , Thierry Poinsot a,b a Institut de Mécanique des Fluides, UMR CNRS/INP-UPS 5502, Allée du Professeur Camille Soula, 31400 Toulouse, France

	b CNRS; IMFT; 31400 Toulouse, France
	a r t i c l e i n f o	a b s t r a c t
	Article history:
	Available online 10 January 2013
	Keywords:
	Combustion instabilities
	Large-Eddy Simulation
	Rocket propulsion
	Mots-clés :
	Instabilités de combustion
	Simulation aux Grandes Echelles
	Moteur-fusée

Table 1

 1 Operating point.

	Experiment	LES

Table 2

 2 Burnt-gases properties in thermodynamic conditions of the CVRC.

	Eq. ratio	z	Scheme	T ad b [K]	Y CH4	Y O2	Y CO2	Y CO	Y H2 O	Y N2
	1	z st = 0.095	Gri-Mech	2661						0
			2-step	2781	< 10 -4	2.25 × 10 -2	1.99 × 10 -1	3.94 × 10 -2	7.39 × 10 -1	0
	0.813	z xp = 0.079	Gri-Mech	2512	< 10 -4	7.50 × 10 -2	2.06 × 10 -1	6.17 × 10 -3	6.99 × 10 -1	0
			2-step	2586	< 10 -4	7.81 × 10 -2	2.02 × 10 -1	8.62 × 10 -3	7.11 × 10 -1	0

< 10 -4

2.49 × 10 -2 2.22 × 10 -1 2.45 × 10 -2 7.12 × 10 -1

Table 3

 3 Time-averaged value of the Rayleigh source term, R, of acoustic energy and burner power output, P b , for both simulations.

	Simulation	R [W]	P b [W]
	2D axi	3.08 • 10 3	1.37 • 10 6
	3D	8.23 • 10 3	1.34 • 10 6

Figure 1.9: Time-averaged fields in the 3D DES of Harvazinski et al. [50]

http://www.cerfacs.fr/4-26334-The-AVBP-code.php

http://elearning.cerfacs.fr/numerical/benchmarks/vortex2d/

This assumption might be wrong if the response depends on the amplitude as suggested by recent work on FDF[22]. Nevertheless, it simplifies considerably the problem and it was used here.

This plot is obtained with the Xnozzle tool in AVBP . It displays P in the chamber versus the outlet flow rate ṁ. The dotted line corresponds to the theoretical formula for a choked nozzle Eq. 1.12. All points on the left side of this curve are unchoked flows while all points on right side the line correspond to choked conditions

A very small level of heat release rate can be observed through the exit nozzle: this is due to a change of the equilibrium between CO and CO2 in the second reaction of the chemical scheme (Sec.

2.1.2) because of the change in pressure and temperature through the nozzle.

Figure 4.11: Instantaneous field of temperature with a time step of 60 µs, : T= 2000 K iso-line.

Note that viscous dissipative terms are neglected in Eq. 4.2. They will be discussed in Sec. 5.1

R.Garby et al. / C. R. Mecanique 341 (2013) 220-229

Acknowledgment

Acknowledgement

The first author is supported by a scholarship intended for PhD students funded by the DGA (Direction Générale de l'Armement).

This work was granted access to the high-performance computing resources of CINES under the allocation c2011_026625 made by Grand Equipement National de Calcul Intensif.

Chapter A : Publication in Comptes Rendus Mécanique Bibliography

Chapter 2 : Description of the numerical tools

The relation between the two impedances Z (P ,u ) and Z (J ,m ) is:

Z (P ,u ) + M . n M . n Z (P ,u ) + 1 .

(2.92)

Motheau et al show that, when the Mach number at the boundary of the domain becomes non-negligible (which is often the case in the combustion chamber which are fed by an inlet nozzle or a compressor and followed by a nozzle or a turbine), using a (J , m ) formulation in an Helmholtz solver ensures a more consistent description of the acoustic fluxes at these boundaries and a proper evaluation of the acoustic energy. This process consists in simply replacing Z (P ,u ) in AVSP by a Z (J ,m ) given by Eq. 2.92 in which the mean Mach number is an additional parameter. Motheau et al show that this simple change allows to have the correct energy fluxes and mode growth rates in AVSP accounting for non-zero Mach effect at the inlet and outlet. This will be used in Chapter 5 where the outlet Mach number of the CVRC is large enough to require such a model.

Chapter 3

Application to the CVRC burner Contents Chapter 3 : Application to the CVRC burner

Ignition of the whole chamber (IGNIT-A)

Gases after the entrance of the chamber (x > 0) are replaced by burnt-gases, the transition is linear between x = -1 cm through x = 0. Figure 3.7 presents the new field of temperature. The evolution of the mass flow rate exiting the nozzle is presented on Fig. 3.8(a). At ignition, as the temperature in the nozzle has been increased by a large factor, the density has decreased, so that the mass flow rate is also initially smaller. After ignition, it increases back to the value of the inlet flow rate. The pressure (Fig. 3.8(b)) grows up quickly and stabilizes around 15 bar. This value is higher than for the cold flow because the stagnation T in Eq. 1.12 has been increased due to reaction. The heat release rate (Fig. 3.9) increases rapidly and stabilizes at 2 × 10 9 W/m 3 . Figure 3.10 shows the path of the pressure and outlet mass flow rate after Chapter 5

Influence of the oxidizer-post length on the stability In this chapter, the method presented in Sec. 2.2 will be applied to the CVRC. The objective is to study the stability of all modes in the CVRC as a function of the oxidizer-post length, L op , which will vary between 3.5 and 7.5 in. For this study, the flame transfer function (obtained in Sec. 4.4) is supposed to be independent of the post length and the Helmholtz solver AVSP is used to predict linear stability using the setup discussed in Sec. 3.2.

Modeling acoustic losses

In the description of the acoustic solver used for combustion oscillations in Sec. 2.2, an important phenomenon was left out for simplicity: viscous effects (including molecular heat transfer). This makes the analysis much simpler and leads to a Rayleigh criterion (Eq. 4.2) which states that a mode is linearly unstable if the driving mechanism of the instability (the Rayleigh term) p ω averaged over time and over the whole volume is larger than the flux of acoustic energy through the boundary (u p averaged over the area of the inlets and outlets). This is true only in the absence of internal losses (for example by conversion of acoustic energy into vorticity).

If internal losses are accounted for, an additional term must be introduced in the acoustic energy equation. At this point, there is even a controversy in the literature, based on the work of Chu [START_REF] Chu | On the energy transfer to small disturbances in fluid flow (part i)[END_REF], who suggested that the acoustic energy defined by Eq. 4.1 may not be the proper one to use in a reacting flow where three independent variables (p , u and s where s is the entropy) should be used instead of only two (p and u ) in Eq. 4.1.

Conclusions

Further work may focus on the determination of the impedance of a complex nozzle like the one in the inlet of the CVRC in order to provide better model for the acoustic solver. Studies of the influence of the heat losses on the instability are also necessary to gain in precision.

Appendix A

Publication in Comptes Rendus Mécanique