Instabilities in Lean

Keywords:

Thermo-acoustic Instablilities

By L. Gicquel F. Nicoud

& T. Poinsot

Thermo-acoustic instabilities arise from the coupling between acoustic waves and flames and can lead to high amplitude instabilities [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF][START_REF] Putnam | Combustion driven oscillations in industry, fuel and energy science series[END_REF][START_REF] Poinsot | Interactions between acoustics and combustion[END_REF][START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems-an overview[END_REF]. In general, these instabilities induce oscillations of all physical quantities (pressure, velocities, temperature, etc.); in the most extreme cases, they can destroy the burner by inducing large amplitude flame motion (flashback) or unsteady pressure (material fatigue). Since the equivalence ratio oscillates when instabilities are present, there is a general trend for combustors to be more unstable when operating in the lean regime. Also, due to new international constraints, pollutant emissions must be reduced and gas turbine manufacturers need to operate their systems under leaner and leaner conditions. Consequently, there is a need to understand combustion instabilities and to be able to predict them at the design level [START_REF] Lieuwen | [END_REF].

The objective of the following sections is to provide the reader with the relevant information regarding the description, modeling and computation of thermo-acoustic instabilities. The basic equations are first recalled in Section 1.1.1. Among the possible levels of description, two are discussed in more details in the subsequent Sections: the Large-Eddy Simulations approach in 1.1.2 and a 3D linear description based on the Helmholtz equation in 1.1.3. Since using appropriate acoustic boundary conditions is critical when analyzing combustion instabilities this issue is discussed in Section 1.1.4. At last, the different tools and approaches discussed are used in order to study the thermo-acoustic behavior of an industrial annular combustor in section 1.1.5.

Basic Equations and Levels of Description

Three types of numerical or semi-analytical methods have been considered so far to predict/describe these instabilities:

1. Large Eddy Simulation (LES) of all relevant scales of the reacting, turbulent, compressible flow where the instability develops. Many recent studies have demonstrated the ability of this method to represent the flame dynamics [6,[START_REF] Murota | Large-eddy simulations of combustion oscillation in premixed combustor[END_REF][START_REF] Caraeni | [END_REF]9,10,11,12], as well as the interaction between reaction zone and acoustic waves [13,14,15,16]. However, even when simulations confirm that a combustor is unstable, LES calculations do not say why and how to control the instability. Besides, because of its intrinsic nature (full three-dimensional resolution of the unsteady

THERMO-ACOUSTIC INSTABLILITIES

Navier-Stokes equations), LES remains very CPU demanding, even on today's computers, 2. Low-order methods where the geometry of the combustor is modelled by a network of homogeneous (constant density) 1D or 2D axisymmetric acoustic elements where the acoustic problem can be solved analytically [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF][START_REF] Stow | Modelling of circumferential modal coupling due to helmholtz resonators[END_REF][START_REF] Evesque | Low-order acoustic modelling for annular combustors: Validation and inclusion of modal coupling[END_REF][START_REF] Evesque | Spinning and azimuthally standing acoustic modes in annular combustors[END_REF][START_REF] Morgans | [END_REF]. Jump relations are used to connect all these elements, enforcing pressure continuity and mass conservation and accounting for the dilatation induced by infinitely thin flame, if any. The acoustic quantities in each segment are related to the amplitudes of the forward and backward acoustic waves which are determined such that all the jump relations and the boundary conditions are satisfied. This can only be achieved for a discrete set of frequencies, ω, which are the roots of a dispersion relation in the complex plane. The main advantage of low-order methods is that they allow the representation of a complex system with only a few parameters, thus allowing an extensive use for pre-design/optimization/control purposes. However, the geometrical details of the combustor cannot be accounted for and only the first "equivalent" longitudinal or orthoradial modes are sought.

3. As an intermediate step between LES and low-order methods, one may consider using a finite-element or finite-volume technique to solve for an equation (or a system of equations) describing the space-time evolution of small amplitude perturbations. A set of linear transport equations for the perturbations of velocity, temperature and density can be derived by linearizing the Navier-Stokes equations [23], where the local unsteady heat release appears as a forcing term. The resulting system of linear partial differential equations for the fluctuating quantities can be solved, for example in the time domain [START_REF] Pankiewitz | Time domain simulation of combustion instabilities in annular combustors[END_REF]. Depending on the coupling between the flame and acoustics, especially the phase between the pressure and heat release fluctuations, some modes present in the initial field can be amplified and grow exponentially; after a while, the unsteady field is dominated by the most amplified mode which can then be analyzed [START_REF] Pankiewitz | Time domain simulation of combustion instabilities in annular combustors[END_REF]. To facilitate the description of time delayed boundary conditions and also to obtain more information about the damped or less amplified mode, it is worth solving the set of linear equations in the frequency space, as proposed by [25] for the wave propagation through a complex baseline flow. If applied within the combustion instability framework, this would give rise to an eigenvalue problem, the eigenvalues being related to the (complex valued) frequencies of the thermo-acoustic modes. Combined with LES, this approach proved

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

useful to understand the structure and nature of the instabilities observed in academic or industrial burners [26,27,15,28].

Except when the thermo-acoustic analysis relies on LES, viscous contributions are generally neglected together with the mixture inhomogeneities. The latter assumption amounts to considering a gas mixture where all species share the same molar weight and heat capacity which is acceptable for typical practical flames. A direct consequence is that the difference in heat capacities r = C p -C v is constant even if C p , C v and their ratio, γ, may depend on temperature.

Under the above assumptions, the mass, momentum and entropy equations read respectively:

Dρ Dt = -ρ ∂u ∂x , (1.1)
ρ Du Dt = - ∂p ∂x , (1.2)
Ds Dt = r Ω p , (1.3)
where Ω is the heat release per unit volume. Together with the state equation and entropy expression

p ρ = rT and s -s st = T Tst C p (T ) T dT -r ln p p st , (1.4)
these transport equations describe the spatio-temporal evolutions of all relevant physical flow quantities.

Although thermo-acoustic instabilities can lead to high amplitude fluctuations, it is meaningful to consider the linear regime to analyze the conditions under which these instabilities appear. Eqs. (1.1) to (1.4) can be linearized by considering a simple case of large scale small amplitude fluctuations, denoted by , super-imposed on a zero Mach number mean flow, denoted by an over-bar, which depends only on space. The instantaneous pressure, density, temperature, entropy and velocity fields can then be written as p = p + p , ρ = ρ + ρ , T = T + T , s = s + s and u = u where the quantities p /p, ρ /ρ, T /T , s /s and u u /c are of order , where 1 and c = γp/ρ is the mean speed of sound. Note that the zero Mach number assumption implies that ∂p/∂x = 0, from Eq. (1.2), and Ω = 0, from Eq. (1.3), the latter condition being acceptable because only the fluctuating quantities are of interest in the linear analysis. The same assumption also implies that the approximation D/Dt ≈ ∂/∂t holds for any fluctuating quantity since u 0, the

THERMO-ACOUSTIC INSTABLILITIES

non linear convective terms are always of second order in . For simplicity, the temporal fluctuations of the heat capacities are often neglected. Injecting the above expansions for the instantaneous flow quantities into Eqs (1.1) to (1.4) and keeping only terms of order , one obtains the following set of linear equations for the fluctuating quantities ρ , u , s and p :

∂ρ ∂t + u ∂ρ ∂x + ρ ∂u ∂x = 0, (1.5)
ρ ∂u ∂t + ∂p ∂x = 0, (1.6
)

∂s ∂t + u ∂s ∂x = r Ω p 0 . (1.7)
The linearized state equation and entropy expression are:

p p - ρ ρ - T T = 0 and s = C p T T -r p p . (1.8)
In order to close the set of Eqs. ( 1.5) to (1.8), a model must be used to express the unsteady heat release Ω in terms of the other fluctuating quantities.

Flame response: Modeling the unsteady behavior of the flame is the most challenging part in the description of thermo-acoustic instabilities [29]. Several models have been proposed in the past to describe the response of conic or V-shape laminar flames [30], accounting for non-linear saturation effects [31,32] and equivalence ratio fluctuations [33,[START_REF] Sattelmayer | Influence of the combustor aerodynamics un combustion instabilities from equivalence ratio fluctuations[END_REF]. Most models describe the global (integrated over space) heat released in the whole flame zone. For premixed flames, the most natural way to proceed is to relate this global quantity to the acoustic velocity in the cold gas region upstream of the flame region. The idea behind this approach is that the heat release is mainly controlled by the fresh gas flow rate, if the flame speed is specified. The most classical model follows seminal ideas by Crocco [START_REF] Crocco | [END_REF]36] and is referred to as the n -τ model. This 1D formulation stipulates that the global heat release at time t is proportional to a time lagged version of the acoustic velocity at a reference upstream position x ref , usually taken at the burner mouth:

Ω tot = V Ω (t) dx = S ref γp γ -1 n u (x ref , t -τ ).
(1.9)

In this expression, the LHS term is the heat release fluctuations integrated over the flow domain V , S ref is the cross section area of the burner mouth,

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

u denotes the fluctuating velocity component in the direction x of the main flow which feeds the flame, the interaction index n controls the amplitude of the flame response to acoustic perturbations and τ is the time delay between the acoustic perturbation and the response of the flame. This latter parameter controls the phase between the acoustic pressure and the unsteady heat release in the flame zone, and thus the value of the Rayleigh index is

R = t V
p Ω dx dt.

(1.10)

According to the classical Rayleigh criterion, flame/acoustics coupling promotes the appearance of instabilities if R > 0, showing the importance of the parameter τ in the description and prediction of thermo-acoustic instabilities. Models for the global response of the flame are only justified for acoustically compact flames, where the typical length of the flame region L f is small compared to the characteristic acoustic wavelength L a . This condition is not always met. It is then natural to use a local flame model which relates the local unsteady heat release to a reference acoustic velocity in the injector mouth. The natural way to proceed is then to write:

Ω (x, t) Ωtot = n u (x) u (x ref , t -τ u (x)) n ref, U bulk , (1.11)
where n u (x) and τ u (x) are fields of interaction index and time lag and n ref, are the components of a fixed unitary vector defining the direction of the reference velocity. The scaling by the total heat release Ωtot and the bulk velocity U bulk have been used to make sure that n u (x) has no dimension. Obviously this modeling approach allows more degrees of freedom than any global model to represent the actual response of a typical industrial flame (two fields of parameters instead of two real numbers). However, a large amount of pointwise data is required to tune such models and for obvious technological reasons these data can hardly be obtained experimentally. As discussed in section 1.1.2, the alternative is then to use compressible reacting LES to investigate the response of a turbulent flame submitted to acoustic perturbations.

Using the local flame model given in Eq. (1.11), the transport equation for s , Eq. (1.7), can be re-written as:

∂s ∂t + u ∂s ∂x = r p Ωtot U bulk n u (x) u (x ref , t -τ u (x)) n ref, , (1.12)
and the set of Eqs. (1.5) to (1.12) can be solved to determine the thermoacoustic properties of the system.

LES of Compressible Reacting Flows

LES [START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF][START_REF] Pope | Turbulent flows[END_REF] is nowadays recognized as an intermediate approach to the more classical Reynolds Averaged Navier-Stokes (RANS) methodologies [START_REF] Pope | Turbulent flows[END_REF][START_REF] Chassaing | Turbulence en mécanique des fluides, analyse du phénomène en vue de sa modélisation à l'usage de l'ingénieur[END_REF].

Although conceptually very different these two approaches aim at providing new systems of governing equations to mimic the characteristics of turbulent flows. Recent studies using LES have shown the potential of this approach for reacting flows (see reviews in [12] or [START_REF] Poinsot | Theoretical and numerical combustion[END_REF]). LES is able to predict mixing [16,[START_REF] Prière | [END_REF]42,43,44], stable flame behaviour [45,46,47,48] and flame acoustic interaction [16,49,[START_REF] Staffelbach | Highly parallel large eddy simulation of multiburner configurations in industrial gas turbines[END_REF][START_REF] Staffelbach | Highly parallel Large Eddy Simulations of multiburner configurations in industrial gas turbines[END_REF]. It is also used for flame transfer function evaluation [27,52,53] needed for Helmholtz solvers (see subsection refsubsec:3DHelm). Although LES seems very promising for industrial applications, it remains computationally too intensive to integrate in the design cycle of the next generation of gas turbines. For example a typical single-sector LES computation as presented below usually costs of the order of 50, 000 CPU-hours. Helmholtz solvers, on the other hand, offer great flexibility and allow the prediction of combustion instabilities while designing new combustion chambers. The computional cost with this approach and for the complete combustion chamber is more of the order of 200 CPUhours It is also important to note that while most academic set ups used to study combustion instabilities [START_REF] Lieuwen | [END_REF][START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Candel | 24th Symp. (Int.) on Combustion[END_REF][START_REF] Crighton | Modern methods in analytical acoustics[END_REF] are limited to single burners and are subjected mainly to longitudinal acoustic modes, real gas turbines exhibit mostly azimuthal modes [START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF]56,[START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF] due to the annular shape of their chambers [START_REF] Lieuwen | [END_REF].

The governing equations for RANS and LES are respectively obtained by ensemble averaging [START_REF] Pope | Turbulent flows[END_REF][START_REF] Chassaing | Turbulence en mécanique des fluides, analyse du phénomène en vue de sa modélisation à l'usage de l'ingénieur[END_REF] and filtering the set of compressible Navier-Stokes equations. These operations yield unclosed terms which are to be modelled. In RANS, the unclosed terms are representative of the physics taking place over the entire range of frequencies present in the ensemble of realizations used for averaging. In LES, the operator is a spatially localized time independent filter of given size, , to be applied to a single realization of the studied flow. Resulting from this "spatial average" is a separation between the large (greater than the filter size) and small (smaller than the filter size) scales. The unclosed terms are representative of the physics associated with the small structures (with high frequencies) present in the flow. Figure 1.1 illustrates the conceptual differences between (a) RANS and (b) LES when applied to a homogeneous isotropic turbulent field.

Due to the filtering approach, LES allows a dynamic representation of the large scale motions whose contributions are critical in complex geometries. The LES predictions of complex turbulent flows are henceforth closer to the physics since large scale phenomena such as large vortex shedding and acoustic waves are embedded in the set of governing equations [START_REF] Poinsot | Theoretical and numerical combustion[END_REF].

For the reasons presented above, LES has a clear potential in predicting turbulent flows encountered in industrial applications especially in the context of thermo-acoustic instabilities. In particular and in conjunction with Helmholtz solvers, LES can provide the estimation and validation of the model used to represent the thermo-acoustic coupling: i.e. the Flame Transfer Function (FTF).

The LES Sub-Grid Scale (SGS) models: LES for reacting flows involves the spatial Favre filtering operation that reduces for spatially, temporally invariant and localised filter functions [58,59] to,

f (x, t) = 1 ρ(x, t) +∞ -∞ ρ(x , t) f (x , t) G(x -x) dx , (1.13)
where G denotes the filter function.

In the mathematical description of compressible turbulent flows with chemical reactions and species transport, the primary variables are the species densities ρ i (x, t), the velocity vector u (x, t), the total energy E(x, t) ≡ e s + 1/2 u u and the fluid density ρ(x, t) = N i=1 ρ i (x, t). The application of the filtering operation to the instantaneous set of compressible Navier-Stokes transport equations with chemical reactions yields the LES transport equations [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] which contain the so-called Sub-Grid Scale (SGS) quantities that need modelling [START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF][START_REF] Ferziger | Large eddy simulation: an introduction and perspective[END_REF].

THERMO-ACOUSTIC INSTABLILITIES

The SGS velocity stress tensor : The unresolved SGS stress tensor τ ij t , is modelled using the Boussinesq assumption [START_REF] Smagorinsky | [END_REF][START_REF] Pope | Turbulent flows[END_REF][START_REF] Chassaing | Turbulence en mécanique des fluides, analyse du phénomène en vue de sa modélisation à l'usage de l'ingénieur[END_REF]:

τ ij t - 1 3 τ t δ ij = -2 ρ ν t S ij , with, S ij = 1 2 ∂ u i ∂x j + ∂ u j ∂x i - 1 3 ∂ u ∂x δ ij . (1.14)
In Eq. (1.14), S ij is the resolved strain rate tensor and ν t is the SGS turbulent viscosity.

Most SGS turbulent viscosity model [START_REF] Smagorinsky | [END_REF]62,63,64] take on the generic form,

ν SGS = C m ∆ 2 OP (x, t), (1.15)
where C m is the constant of the model, ∆ is the subgrid characteristic length scale (in practice the size of the mesh) and OP is an operator of space and time, homogeneous to a frequency, and defined from the resolved fields.

The SGS species and energy flux models: The SGS species flux J i t and the SGS energy flux q t are, in most cases, respectively modelled by use of the species SGS turbulent diffusivity

D i t = ν t /Sc i t ,
where Sc i t is the turbulent Schmidt number (= 0.7 for all i). The SGS thermal conductivity for energy flux is also obtained from ν t by λ t = ρ ν t C p /P r t where P r t is a turbulent Prandtl number (= 0.7),

J i t = -ρ D i t W i W ∂ X i ∂x -Y i V c , with, q t = -λ t ∂ T ∂x + N i=1 J i t h i s . (1.16)
In Eq. (1.16), the mixture molecular weight W and the species molecular weight W i can be combined with the species mass fraction to yield the expression for the molar fraction of species i: X i = Y i W/W i . V c is the diffusion correction velocity resulting from the Hirschfelder and Curtiss approximation [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] and T is the Favre filtered temperature which satisfies the modified filtered state equation p = ρ r T [65,66]. Finally, h i s stands for the enthalpy of species i. Note that the performances of the closures could be improved using dynamic formulations described in [65,62,67,68,69].

The Dynamic Thickened Flame (DTF) model : LES of turbulent reacting flows imply the modelling of SGS combustion terms. One model employed in the context of thermo-acoustic instabilites is the Thickened Flame (TF) model [9]. Following the theory of laminar premixed flames [START_REF] Williams | Combustion theory[END_REF], the flame speed S o L and the flame thickness δ o L may be expressed as,

S o L ∝ √ α A , and, δ o L ∝ α S o L = α A , (1.17)
where α is the thermal diffusivity and A the pre-exponential constant of the reaction rate. Increasing the thermal diffusivity by a factor F , the flame speed is kept unchanged if the pre-exponential factor is decreased by the same factor [START_REF] Butler | A numerical method for two-dimensional unsteady reacting flows[END_REF]. This increases the flame thickness by factor F which is easily resolvable on a coarser mesh. Additional information needs however to be supplied to reproduce the effect of turbulence-chemistry interaction at the subgrid-scales [START_REF] Angelberger | Large eddy simulations of combustion instabilities in premixed flames[END_REF]13,[START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]. This is the intent of the so-called efficiency function, E [9]. When thickening is applied everywhere in the flow, the model is limited to fully premixed combustion. If mixing is present then thickening will strongly interfere with the physics by increasing the diffusion artificially everywhere in the flow. Thus, diffusion flames are inappropriate configurations for such an approach. To compute partially premixed or nonpremixed flames [START_REF] Poinsot | Theoretical and numerical combustion[END_REF], a modified version of the Thickened Flame model (TF) is used [START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]26,16,[START_REF] Boudier | [END_REF]: the Dynamic Thickened Flame model (DTF).

With the DTF model, the SGS fluxes are modified to become:

J i t = -(1 -S) ρD i t W i W ∂ X i ∂x + ρ Y i V c , with, q t = -(1 -S) λ t ∂ T ∂x + N i=1 J i t h i s , (1.18)
where S is a sensor detecting reaction zones: i.e. derived from an Arrhenius type of law for example. The local thickening factor depends on the local mesh size: typically thickening must ensure that enough points are present in the flame zone and the thickening factor F is given by,

F = 1 + (F max -1) S, and, F max = N c ∆ x δ o L , (1.19)
where N c is the number of grid points, typically 5 to 10, used to resolve the flame front and ∆ x is the local mesh spacing. Although this approach is still being developed and further validations are needed, its ease of implementation and it success in prior applications [26,16,[START_REF] Sengissen | Simulation aux grandes échelles des instabilités de combustion: vers le couplage fluide/structure[END_REF][START_REF] Boudier | [END_REF]28] suggest its suitability for the problem of thermo-acoustic instabilities such as presented in this chapter.

Numerical issues in LES solvers: In many cases, turbulence results from the development, amplification and saturation of unstable hydrodynamic modes of the main flow. Any numerical method used to compute such flow must therefore be able to represent the growth of these modes: i.e. it must not be too dissipative. In high-Reynolds number flows, the scale separation can be large (the integral to Kolmogorov length scale ratio is Re 3/4 , with Re = k 2 /(νε), where k is the turbulent kinetic energy, ν the kinematic

THERMO-ACOUSTIC INSTABLILITIES

viscosity and ε the rate of turbulent kinetic energy dissipation), and it is worth minimizing the number of grid point necessary to represent the smallest scales. Finally, the effective dissipation at the Kolmogorov scale must not be over estimated if the actual flow Reynolds number is to be accounted for (with Re = k 2 /(νε), any extra dissipation decreases Re). Numerics, especially in the context of LES, faces an important constraint: i.e. numerical dissipation must be as small as possible for all the length scales present in the flow. Such a constraint is the reason why spectral methods [START_REF] Orszag | [END_REF] have been considered, until the early 90's, as the only appropriate methods for performing DNS or LES of turbulent flows. However, in the case of complex geometries or boundary conditions, spectral methods cannot be used and the simulations must be based on either finite-volume, finite-element or finite-difference methods. The three methods can be used for unsteady simulations as long as appropriate spatial and temporal time stepping procedures are used.

As mentioned earlier, the numerical error must be controlled and minimized for all length scales present in the unsteady flow to be computed. This means that the accuracy of a numerical scheme cannot simply be reduced to its order of accuracy. As far as unsteady flow computations are concerned, it is necessary to perform a wavelength based numerical analysis where one considers a harmonic perturbation and compares how the discrete and the exact derivatives operate on this perturbation [START_REF] Vichnevetsky | Fourier analysis of numerical approximations of hyperbolic equations[END_REF]. The effective-to-exact wavenumber ratio, κ /κ, can be used to quantify the errors related to the numerical scheme. Considering the simple linear convection equation, this ratio can be interpreted as an error in the speed of propagation of a perturbation of wavelength κ. In a general case, κ /κ can be written in the form κ /κ = E(κ∆x) = E r (κ∆x) + i E i (κ∆x) and the effective equation solved numerically reads, ∂f ∂t

+ u 0 E(κ∆x) ∂f ∂x = 0. (1.20)
Assuming that the initial condition is f (x, t = 0) = exp(iκ), the exact solution of Eq. (1.20), or equivalently the solution for the linear convection equation provided by the numerical scheme (with perfect time advancement), is simply,

f (x, t > 0) = exp [κ E i (κ∆x) u 0 t] exp [iκ (x -E r (κ∆x) u 0 t)] . (1.21)
When κ∆x tends to zero or the number of grid point per wavelength tends to infinity, κ /κ = E(κ∆x) tends to unity and the exact solution is recovered: i.e. exp (i κ (x -u 0 t)). When the imaginary part of κ /κ is not zero, E i (κ∆x) = 0, the amplitude of the harmonic perturbation is not conserved; it

FD name E r [f i -f i-1 ]/∆x 1st order upwind [sin(κ∆x)]/κ∆x [f i+1 -f i-1 ]/2∆x 2nd order centered [sin(κ∆x)]/κ∆x [3f i -4f i-1 + f i-2 ]/2∆x 2nd order upwind [sin(κ∆x)][2 -cos(κ∆x)]/κ∆x [-f i+2 + 8f i+1 -8f i-1 + f i-2 ]/12∆x 4th order centered [sin(κ∆x)][4 -cos(κ∆x)]/3κ∆x
Table 1.1: Classical finite difference formula for the spatial first derivative and associated error.

0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 k ! x Re[E(k ! x)] 0 0.5 1 1.5 2 2.5 3 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 k ! x Im[E(k ! x)]
Figure 1.2: Effective-to-exact wavelength ratios for the schemes displayed in Table . 1.1. The imaginary part (right plot) is zero for centered schemes. The 1 st order upwind and the 2 nd order centered schemes share the same real part (left plot). These graphs can be interpreted as the effective-to-exact convection velocity ratio, or as the effective-to-exact first derivative ratio. The legend goes as: 2 nd order centered, continuous line; 4 th order centered, dash line; 2 nd order upwind, dot-dash line. is damped if E i (κ∆x) < 0 and unbounded if E i (κ∆x) > 0. The effective-toexact wavelength ratios of some classical finite difference schemes are reported in Table . 1.1 and plotted in Fig. 1.2.

Centered finite difference schemes have real valued κ /κ ratio and thus are non dissipative, whatever their order is. This property is not shared by the biased schemes which all introduce dissipation. Note also that a property shared by all finite difference schemes is that they cannot propagate wiggles accurately (E(π) = 0).

The same analysis can be performed for second derivative schemes for which κ 2 /κ 2 ratios are of interest. For LES, the effective energy dissipation issued by second-order derivatives at the small scales cannot be close to zero, whatever the SGS model is. Indeed, the energy flux from the largest to the

THERMO-ACOUSTIC INSTABLILITIES

13 smallest scales should be balanced to avoid an energy accumulation at the smallest resolved scales [START_REF] Nicoud | Towards direct and large eddy simulations of compressible flows in complex geometries[END_REF].

A key issue when performing LES of turbulent flows is the necessity to use virtually non-dissipative schemes to handle flow fields which contain a lot of energy at large wave numbers. Since numerical errors are large for the smallest scales, the risk for such computation to run unstable is high. As far as incompressible Navier-Stokes equations are concerned, experience shows that the kinetic energy must be conserved if a stable and dissipation-free numerical method is sought. Indeed, such property ensures that the sum of the square of the velocities cannot grow, even through non-linear interactions between modes exist: a numerical scheme which conserves kinetic energy cannot lead to unbounded growth of oscillations [START_REF] Morinishi | [END_REF]80,81,82]. It should be noted that the concept of kinetic energy conservation is only valid for incompressible or low-Mach number flows. Since the extension of this principle to compressible situations is quite difficult [83], one usually relies on numerical stabilization via artificial viscosity to stabilize LES of such flows [START_REF] Jameson | Numerical solution of the euler equations by finite volume methods using runge-kutta time stepping schemes[END_REF][START_REF] Colin | Simulations aux grandes échelles de la combustion turbulente prémélangée dans les statoréacteurs[END_REF].

To conclude on LES and the difficulties of the approach, modelling is clearly needed for the problem to be solved adequately from a purely mathematical and physical point of view. Numerics is however crucial and is also to be adequately addressed so as to properly qualify the modelling strategy at hand. Fundamental issues such as error propagation and cancellation are still not clear in LES codes. Parallelization and associated alogorithms are another difficulty which emphasize the need for careful validation and developments of such tools.

Flame transfer function based on LES: evaluation and validation:

To predict and avoid combustion instabilities [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Poinsot | [END_REF][START_REF] Candel | 24th Symp. (Int.) on Combustion[END_REF][START_REF] Crighton | Modern methods in analytical acoustics[END_REF], a well-known method is the identification of the combustion chamber response or Flame Transfer Function (FTF) to acoustic waves introduced in the combustor using loudspeakers or rotating valves. This identification of the FTF is usually performed either numerically or experimentally. Two methods may be found in the literature to analyze this response: the Purely Acoustic (PA) approach [87,[START_REF] Polifke | Determination of thermo-acoustic transfer matrices by experiment and computational fluid dynamics[END_REF][START_REF] Paschereit | Characterization of lean premixed gasturbine burners as acoustic multi-ports[END_REF][START_REF] Paschereit | Prediction of combustion oscillations in gas turbine combustors[END_REF][START_REF] Pankiewitz | Computation of transfer matrices for gas turbine combustors including acoustics/flame interaction[END_REF] and the n -τ approach [START_REF] Crocco | [END_REF]36,[START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Krebs | Thermoacoustic flame response of swirl flames[END_REF][START_REF] Bernier | [END_REF][START_REF] Cheung | Measurement and analysis of flame transfer function in a sector combustor under high pressure conditions[END_REF]. In the PA approach, the burner is considered as a "black box" and a two-ports formulation (based on acoustic pressure and velocity perturbations) is used to construct a transfer matrix linking acoustic fluctuations on both sides of the burner. In the n -τ approach, pressure measurements are replaced by a global heat release measurement (usually based on optical methods). The heat release fluctuations are then related to the flow velocity modulations at the combustor inlet. Both PA and n -τ methods can be used experimen-

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

tally or numerically. Numerical experiments comparing FTF results show that the n -τ approach often leads to an ill-defined problem where the measured transfer function depends on acoustic impedances upstream and downstream of the combustor.

The essential drawback of the original n -τ model comes from the fact that it tries to correlate heat release perturbations to velocity perturbations only. With the FTF model presented in [START_REF] Truffin | [END_REF] (called "extended n -τ "), a consistent formulation can be used for any location of the reference point by introducing the effects of pressure perturbation on heat release: The FTF model is formulated using the local unsteady pressure and velocity measured upstream of the flame,

Ω Ωtot = A u u (x a , t -τ u ) c + A p p (x a , t -τ p ) p , (1.22)
where the unsteady velocity, pressure and heat release are scaled respectively by the sound speed c, the mean pressure p and the mean integrated heat release Ωtot = V Ω(x) dx. 1 Assuming that all the fuel is burnt Ωtot may be estimated using the fuel mass flow rate ṁF and the heat of reaction Q: The concept of FTF and its modeling being clearly defined in the context of the PA or extended n -τ formalisms [START_REF] Truffin | [END_REF], two options are available for its estimation: laboratory measurements or numerical simulations. In the former, turbulent closure and combustion models are introduced and validation of the FTF is needed. For real turbulent flow configurations, two approaches are conceptually able to produce estimates of the FTF: RANS (or more suitably Unsteady RANS) and LES. Figure 1.6 compares numerical predictions [52] and experimental measurements of the FTF using the original n -τ model [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Crocco | Research on combustion instability in liquid propellant rockets[END_REF] in a gas turbine configuration with turbulent flow, shown in Fig. 1.5. Experimental results are given by the phase comparison between the signal given by a photomultiplier (with a OH* frequency filter) Figure 1.5: Forced turbulent premixed flame in a real gas turbine configuration: iso-surface of heat release colored by the local value of the delay as estimate with LES and the classical n -τ model [52].

Ωtot = Q ṁF (1.23) (a) (b) (c) (d)
which is directly linked to heat release, and a hot wire probe velocity signal at the inlet of the combustion chamber. Values of the LES phase φ ω = -τ ω * ω (where ω = 2π f is the the angular frequency, f the frequency and τ ω is the time delay) are in good agreement with experimental values (Fig. 1.6) and confirm the potential of LES when compared to RANS.

The differences between RANS and LES show that the heat release fluctuations are not solely linked to the convection of an initial perturbation but can result from a complex interaction between acoustics and the flow in the chamber. That observation supports the need for a fully spatially and temporally dependent numerical approach for a proper numerical estimation of the FTF. As of today, only LES seems able to provide such a comprehensive framework as illustrated on real applications.

3D Helmholtz Solver

Despite the potentials of LES, simpler and faster methods are often needed to design stable combustors. Solving the wave equation in reacting flow and identifying linearly unstable modes is such a technique and is found in Helmholtz solvers. These tools offer great flexibility and permit the prediction of combustion instabilities beforehand while designing the new combustion chambers. As noted earlier, while most academic set ups used to study combustion instabilities [START_REF] Lieuwen | [END_REF][START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Candel | 24th Symp. (Int.) on Combustion[END_REF][START_REF] Crighton | Modern methods in analytical acoustics[END_REF] are limited to single burners and are subjected mainly to longitudinal acoustic modes, real gas turbines exhibit mostly azimuthal modes [START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF]56,[START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF] due to the annular shape of their chambers [START_REF] Lieuwen | [END_REF] for which LES is very difficult and expensive to use while Helmholtz solver can be applied economically. This section presents the required wave equation and its solution methodology using linear algebra techniques.

Wave equation: Taking the time derivative of Eq. (1.5), adding the divergence of Eq. (1.6) and using Eqs. (1.7) and (1.8) to eliminate ρ yields the following wave equation for p , ∂ ∂x

1 ρ ∂p ∂x - 1 γp ∂ 2 p ∂t 2 = - γ -1 γp ∂ Ω ∂t , (1.24)

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

when the Mach number of the mean flow is zero (i.e., u /c = 0). An order of magnitude analysis suggests that this assumption (u 0) is valid when the characteristic Mach number M = √ u u /c of the mean flow is small compared to L f /L a where L f is the flame zone thickness and L a is the typical acoustic wavelength [START_REF] Nicoud | [END_REF]. However, the effect of the approximation M 0 on the shape, frequency of oscillation and stability of the thermo-acoustic modes is far from being well understood [98,99,[START_REF] Polifke | Suppression of combustion instabilities through destructive interference of acoustic and entropy waves[END_REF][START_REF] Polifke | [END_REF]102]. Recent studies suggest that the validity domain of the zero mean flow assumption might be rather small [103]. Nevertheless, this somewhat restrictive assumption is necessary to derive a wave equation for the thermo-acoustic perturbations. This situation is different for classical aeroacoustics where combustion is not present and where a wave equation for the perturbation potential can be derived if the baseline flow is assumed homentropic and irrotational [START_REF] Howe | Acoustics of Fluid-Structure Interaction[END_REF]. Since assuming the mean flow to be homentropic is not realistic when dealing with combustion instabilities, assuming that the mean flow is at rest is the most convenient way to simplify the formalism, the alternative being to deal with the complete set of Linearized Euler Equations [25,103].

Equation (1.24) being linear, it is natural to introduce harmonic variations at frequency f = ω/(2π) for pressure, velocity and local heat release perturbations, p = (p(x) exp(-j ωt)) , u = (û (x) exp(-j ωt)) ,

Ω = Ω(x) exp(-j ωt) .

Introducing Eq. (1.25) into Eq. (1.24) leads to the following Helmholtz equation, ∂ ∂x

1 ρ ∂ p ∂x + ω 2 γp p = j ω γ -1 γp 0 Ω(x), (1.26)
where ρ and γ depend on the space variable x and the unknown quantities are the complex amplitude p(x) of the pressure oscillation at frequency f and angular frequency ω. In the frequency space, the zero Mach number assumption leads to j ω û = (∂ p/∂x )/ρ and the flame model, Eq. (1.11), translates into,

Ω(x) = Ωtot j ω ρ(x ref ) U bulk n u (x) e j ω τu(x) ∂ p(x ref ) ∂x k n ref,k . (1.27)
Introducing Eq. (1.27) into Eq. (1.26) leads to, ∂ ∂x

1 ρ ∂ p ∂x + ω 2 γp p = γ -1 γp Ωtot ρ(x ref )U bulk n u (x)e j ω τu(x) ∂ p(x ref ) ∂x k n ref,k .(1.28) 1.1. THERMO-ACOUSTIC INSTABLILITIES 19
The applications discussed in references [26,27,[START_REF] Sensiau | [END_REF] as well as in section 1.1.5 are based on the solution of Eq. (1.28) but the methodologies developed can be applied to a more general case where the complex amplitude of the heat release is given by [START_REF] Nicoud | [END_REF]:

Ω(x) = Lu ∂ p(x) ∂x + Lp [p(x)] , (1.29)
where Lp and Lu are two linear operators acting on p and its gradient respectively. Of course the general formulation Eq. ( 1.29) has potential to include more physical effects than the local n -τ model described by Eqs. (1.11) and (1.27). Notably, it allows relating the unsteady heat release to the complete acoustic field at the reference position x ref instead of the velocity field only, consistently with the matrix identification approach for flame modeling [106,107] (see also Eq. (1.22)). Although the effects of the acoustic pressure are often neglected in flame transfer formulations, relating the unsteady heat release to the complete acoustic field (velocity and pressure) is highly desirable for cases where the flame is not compact or when its distance to the injector mouth is not small compared to the acoustic wavelength [START_REF] Truffin | [END_REF].

Boundary conditions: Denoting by n BC = (n BC,1 , n BC,2 , n BC,3 ) the outward unit normal vector to the boundary ∂Ω of the flow domain, three types of boundary conditions are usually used for acoustics:

• Zero pressure: this corresponds to fully reflecting outlets where the outside pressure is imposed strongly at the flow boundary, zeroing the pressure fluctuations:

p = 0, on boundary ∂Ω D , (1.30)
where the subscript D in ∂Ω D refers to the subset of ∂Ω where this Dirichlet boundary condition holds.

• Zero normal velocity, viz. ûk n BC,k = 0: this corresponds to fully rigid walls or reflecting inlets where the velocity of the incoming flow is imposed, zeroing the velocity fluctuations. Under the zero Mach number assumption, Eq. (1.6) can be used to re-write this condition as a Neumann condition for the acoustic pressure:

∂ p ∂x k n BC,k = 0, on boundary ∂Ω N , (1.31)
where the subscript N in ∂Ω N refers to the subset of ∂Ω where this Neumann boundary condition holds.

• Imposed reduced complex impedance Z = p/ρc ûk n BC,k . Under the zero Mach number assumption, this condition can be re-written as a linear relationship between the acoustic pressure and its gradient in the normal direction to the boundary:

cZ ∂ p ∂x n BC, -j ω p = 0, on boundary ∂Ω Z , (1.32)
where the subscript Z in ∂Ω Z refers to the subset of ∂Ω where the reduced impedance is imposed.

Associated with the homogeneous boundary conditions (1.30), (1.31) and (1.32) on ∂Ω = ∂Ω D ∂Ω N ∂Ω Z , Eq. (1.28) defines a non-linear eigenvalue problem whose solutions provide the shape, frequency and growing/damping rate of the relevant thermo-acoustic modes.

Assuming that the sound speed c and the density ρ distributions over space are known, Eq. (1.28) can be solved using a Galerkin finite-element method to transform this equation into a nonlinear eigenvalue problem of size N (the number of nodes in the finite element grid used to discretize the geometry, except those nodes belonging to ∂Ω D where p = 0 is known) of the form, ). Thus, Eq. (1.33) defines a nonlinear eigenvalue problem which must be solved iteratively, the k th iteration consisting in solving the quadratic eigenvalue problem in ω k defined as,

[A][P ] + ω[B(ω)][P ] + ω 2 [C][P ] = [D(ω)][P ], (1.33
([A] -[D(ω k-1 )]) [P ] + ω k [B(ω k-1 )][P ] + ω 2 k [C][P ] = 0. (1.34)
A natural initialization is to set [D](ω 0 ) = 0 so that the computation of the modes without acoustic/flame coupling is in fact the first step of the iteration loop. Usually, only a few (typically less than 5) iterations are enough to converge toward the complex frequency and associated mode. Note that a quadratic problem must be solved at each iteration Eq. (1.34). These problems are rather well known from a theoretical point of view; they can be efficiently solved numerically once converted into an equivalent linear problem of size 2 × N [108], for example by making use of a parallel implementation of the Arnoldi method [START_REF] Lehoucq | Arpack: Solution of large scale eigenvalue problems with implicitly restarted arnoldi methods[END_REF] available in the P-ARPACK library.

THERMO-ACOUSTIC INSTABLILITIES

21

Another option is to solve the quadratic eigenvalue problem directly without linearizing it; a specific algorithm must then be used instead of the Arnoldi approach. A good candidate is the Jacobi-Davidson method [START_REF] Sleijpen | [END_REF] which has recently been applied successfully to combustion instability problems [111]. Another way to proceed is to define the k th iteration in the following way:

([A] -[D(ω k-1 )] + ω k-1 [B(ω k-1 )]) [P ] + ω 2 k [C][P ] = 0, (1.35)
so that a linear eigenvalue problem must be solved at each sub-iteration and the classical Arnoldi iterative method [START_REF] Lehoucq | Arpack: Solution of large scale eigenvalue problems with implicitly restarted arnoldi methods[END_REF] can be used. This latter formulation showed good potential for large scale problems (N of order 10 6 ) arising from the thermo-acoustic analysis of annular combustors [START_REF] Gullaud | Damping effect of perforated plates on the acoustics of annular combustors[END_REF]. More details can be found in [START_REF] Nicoud | [END_REF].

Accounting for dissipative effects: The linear formulation described above is dissipation-free since no damping terms have been taken into account for its derivation (except for the acoustic radiation at boundaries which can be modeled via a complex valued impedance). However, damping effects should be included in some practical cases, for example when dealing with modern combustors for which perforated liners are increasingly used. Multiperforated plates (MP) are widely used in combustion chambers of turbofan engines to cool the chambers walls exposed to high temperatures [START_REF] Lefebvre | Gas Turbines Combustion[END_REF]. These plates consist of submillimeter apertures, across which the mean pressure jump forces a cold jet through the holes, from the casing into the combustion chamber. The micro-jets then coalesce to form a cooling film. Due to the tiny diameter of the perforations, the holes cannot be meshed for numerical computations and a model is required for the effect of perforated plates. This problem is encountered not only in CFD calculations [START_REF] Mendez | [END_REF]115], but also while computing acoustic modes of combustion chamber. Indeed, MP are known to have a damping effect on acoustics [116,117], which is enhanced by the presence of a mean bias flow [118]. Acoustic waves interact with the shear layer, creating a vortex breakdown at the rims of the apertures [START_REF] Dasse | Les of the acoustic response of a perforated plate[END_REF], which converts part of the acoustic energy into vortical energy (see Fig. 1.7).

To study the acoustic behavior of such a device, it is useful to introduce the Rayleigh conductivity K R [START_REF] Rayleigh | The Theory of Sound, Mac Millan[END_REF] of the aperture, relating the harmonic volume flux Q to the acoustic pressure jump across the plate:

K R = j ω ρ Q p+ -p-, (1.36)
where p+ and p-are the harmonic pressures upstream and downstream of the aperture respectively. Of course the following equation holds: where û± is the acoustic velocity on the plate, equal on both sides. Hence,

Q = d 2 û± , ( 1
K R = j ω ρd 2 û± p+ -p-. (1.38)
Howe expressed the Rayleigh conductivity for a circular aperture in an infinitely thin plate [117] as:

K R = 2a(Γ R -j∆ R ), (1.39)
where

Γ R -j∆ R = 1 + π 2 I 1 (St)e -St -jK 1 (St)sinh(St) St( π 2 I 1 (St)e -St + jK 1 (St)cosh(St)) , (1.40)
with St as the Strouhal number defined by ωa/U . Using the momentum equation and Eq. (1.38), one obtains:

∂ p ∂x n = K R d 2 [p + -p-] (1.41)
This simple analytical model showed good agreement with experiments [118] at least in the linear limit [START_REF] Tran | Damping combustion instabilities with perforates at the premixer inlet of a swirled burner[END_REF]. Improvements of this model have been made to include the plate thickness [122] and the interaction between the apertures [123].

THERMO-ACOUSTIC INSTABLILITIES

When using a linear Helmholtz solver to compute the modes of a thermoacoustic system, Eq. (1.41) can be used as a Neumann boundary condition at both sides of a multi-perforated plate present in the computational domain: once the geometrical properties of the plate (a, d) have been selected together with the bias flow velocity (U ), Eqs. (1.39), (1.40) and (1.41) can be used to express the pressure gradient normal to the multiperforated plate as a function of the angular frequency ω. A numerical procedure similar to the one used for accounting for complex valued impedance (see section 1.1.3) can then be used to solve the eigenvalue problem. This allows us to account for the acoustic damping related to the acoustic-to-vortical energy transfer while keeping an inviscid formulation based on the zero Mach number assumption [START_REF] Gullaud | Damping effect of perforated plates on the acoustics of annular combustors[END_REF].

Upstream/Downstream Acoustic Conditions

In practical applications, the combustion chamber where the zero Mach number thermo-acoustic analysis is relevant is surrounded by decelerated/accelerated regions. Thus, upstream and downstream boundary conditions must be prescribed in order to account for the acoustic impedance of the compressor and turbine stages. These complex valued impedances can be assessed analytically under the so-called compact assumption discussed next or numerically in the more general case.

Acoustic impedance under the compact assumption: In the low frequency limit, the acoustic wavelength is much larger than the characteristic length of the upstream and downstream devices which surround the combustor. Using the mass, energy and entropy conservations, Marble and Candel [124] established the relations linking the different perturbations in the case of planar waves traveling throughout quasi-1D devices (some analytical results can also be obtained in the case of circumferential modes in a choked nozzle [125] or for 2D baseline flows [START_REF] Cumpsty | Proc. R. Soc. Lond. A[END_REF]). For example, one can show analytically that the reflexion coefficient of a compact choked nozzle equals:

1

-(γ -1)M /2 1 + (γ -1)M /2 , (1.42)
where M is the Mach number of the flow entering the nozzle. More details regarding the analytical treatment can be found in the references cited above as well as in [START_REF] Leyko | Numerical and analytical investigation of the indirect noise in a nozzle[END_REF].

2A + exp jk + x in = p + ρcû,
where A + exp (jkx) stands for the forward propagating wave, k = ω/c is the acoustic wave number and A + is the associated pre-exponential factors which is set to any non zero value to ensure that the inlet condition is non homogeneous, 3. define the appropriate boundary condition to be prescribed at the outlet section depending on whether the mean flow is subsonic or supersonic, typically imposed pressure or no boundary condition respectively, 4. solve the corresponding linear system Eq. 1.62,

compute the acoustic equivalent impedance as

Z in = p/ρcû assessed at x = x in .
Note that when the mean flow is subsonic, there is a backward propagating wave entering the domain through the outlet section S out so that an outlet Acoustic impedance of non-compact elements: A numerical approach can also be used to compute the acoustic impedance of diffusers/nozzles under the isentropic mean flow assumption [START_REF] Lamarque | [END_REF]; it can be seen as a way to extend previous analytical results [124] to non compact nozzles. The appropriate equations to be considered are the quasi-1D linearized Euler equations written in the frequency space and under the constant mean entropy assumption. Once discretized, these equations can be converted into an linear algebraic system:

[A][V ] = [BT ], (1.43)
where [V ] is the discrete counterpart of the vector of acoustic unknowns V = (ρ, û) T , the matrix [A] depends on both ω and the details of the spatial discretization and the right-hand-side term comes from possibly non homogeneous boundary conditions. For any quasi-1D flow domain with inlet and outlet section S in and S out respectively (see Fig. 1.8), the following procedure is used to compute the equivalent acoustic impedance:

1. fix the frequency ω, 2. impose a non zero forward propagating acoustic wave at the inlet section S in . The corresponding boundary condition relates p and û at x = x in , viz.

2A + exp j k + x in = p + ρc û,
where A + exp (j kx) stands for the forward propagating wave, k = ω/c is the acoustic wave number and A + is the associated pre-exponential factors which are set to any non zero value to ensure that the inlet condition is non homogeneous, 3. define the appropriate boundary condition to be prescribed at the outlet section depending on whether the mean flow is subsonic or supersonic, Note that when the mean flow is subsonic, there is a backward propagating wave entering the domain through the outlet section S out so that an outlet boundary condition is required. In this case, the above procedure turns out to provide a way to transform a supposedly known acoustic boundary condition at S out to another condition at S in . For example, when the flow domain is a nozzle, this procedure allows us to displace an acoustic boundary condition at a high speed section to an upstream, low Mach number position. In the case where the nozzle is choked, no extra acoustic condition is required since no wave can enter the domain through the outlet section S out . In the particular case where the outlet section coincides with the location of the throat, the proper acoustic impedance to impose at S out is given by [124,[START_REF] Lieuwen | [END_REF],

Z th = 2du/dx -j ω (γ -1)du/dx -j ω , (1.44)
and the above procedure allows us to convert this impedance condition valid at the sonic throat to another condition valid at an upstream, low Mach number location.

Application to an Annular Combustor

The target configuration chosen to illustrate the proposed LES / Helmholtz solver strategy corresponds to an annular helicopter combustion chamber equipped with fifteen burners designed for a helicopter by Turbomeca, shown in Fig 1 .9. Each burner contains two co-annular counter-rotating swirlers. The fuel injectors are placed in the axis of the swirlers. To avoid uncertainties in boundary conditions the chamber's casing is also computed. The computational domain starts after the inlet diffuser and ends at the throat of the high pressure stator. In this subsection, the flow is choked allowing for an accurate acoustic representation of the outlet. The air and fuel inlets use non-reflective boundary conditions [129]. The air flowing at 578 K in the casing feeds the combustion chamber through the swirlers, films and dilution holes. To simplify the LES, fuel is supposed to be vaporized at the lips of the injector and no model is used to describe liquid kerosene injection, dispersion and vaporization. All LES's presented here use the Smagorinsky approach [START_REF] Smagorinsky | [END_REF] to model SGS stresses. Combustion is modeled using Arhenius type reaction rates: a reduced one-step scheme for JP10 / air flames fitted to match the full scheme's behavior for equivalence ratios ranging from 0.4 to 1.5 [START_REF] Boudier | [END_REF]130] is used. Five species explicitly solved are JP 10, O 2 , CO 2 , H 2 O and N 2 . Turbulence/flame interaction is modeled with the DTF model [9,16,[START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]131] described earlier.

A high-order spatial and temporal scheme (TTGC [132]) is used to propagate acoustic waves with precision. Computations are obtained for (a) the entire configuration (15 burners) and (b) a single sector (see Fig. 1.9 (b)) computational domain for FTF evaluations prior to (c) Helmholtz analysis.

Massively parallel LES of the full annular chamber: In the first computation [START_REF] Staffelbach | Proc. of the Combustion Institute[END_REF][START_REF] Boudier | [END_REF], shown in Fig. 1.10, the whole chamber is simulated from the diffuser outlet to the high pressure stator nozzle by a 9, 009, 065 nodes and 42, 287, 640 cell mesh 2 . The LES captures the self-excited instability and results (unsteady pressure RMS and phase fields) show that it is character-ized by two superimposed rotating modes with different amplitudes3 shown in Fig. 1.11. The turning modes4 are found to modulate the flow rate through the fifteen burners and the flames oscillate back and forth in front of each burner shown in Fig. 1.12, leading to local heat release fluctuations. Due to the rotating motion of the modes, all individual transfer functions of all the burners are the same: i.e., no mechanism of flame interactions between burners within the chamber is identified. Note that although such computations are very CPU demanding and not realistically possible in an industrial context, they allow validation of the hypotheses introduced in the use of a hierarchy of computational modelling based on single sector LES (b) and Helmholtz solvers (c).

The flow rate fluctuations in this computation are such that they impact the flame response and the heat release perturbations, Fig. 1.13. Swirler flow rate and pressure fluctuations are out-of-phase. The phase between global heat release and swirler flow rate fluctuations satisfies a Rayleigh criterion for the fifteen sectors.

It is possible to evaluate the response to the flow rate oscillations by computing the transfer function between inlet velocity fluctuations and mean single sector unsteady heat release [START_REF] Crocco | [END_REF]36]. Figure 1.14 shows the modulus n and the phase τ of the classical n -τ model for each of the fifteen burners. Amplitudes and delays are fairly constant suggesting a common response for all burners. This has important implications for steps (b) and (c) of the proposed computational hierarchy.

Single sector LES: FTF evaluation: Based on the previous observations, single sector forced LES should be sufficient to retrieve the FTF necessary for the 3D Helmholtz solver computations. For the case of interest, a single sector forced LES is obtained for an inlet acoustic modulation at 600 Hz. The local response amplitude coefficient, n(x) at 600 Hz is plotted in three different planes in Fig. 1.15 . One remarks that the flame response is spread in the primary zone, it is neither homogenous nor symetric. In the Figure 1.10: 3D view of the computational domain; temperature field on a cylindrical plane passing through all the swirlers with velocity magnitude isocontours. Black dots denote typical probe locations for which diagnostics are provided below.

following, the response amplitude coefficient n(x) is kept in its local formulation in order to take into account the flame response inhomogeneities. This flame response is assumed to be independent of the frequency of the acoustic forcing. This assumption is false if the forcing frequency varies in a broad band, but since in this particular case the known frequency of the first azimuthal mode varies between 550 Hz and 610 Hz, the flame response can be considered constant. To obtain the flame response over a broader frequency range, a white noise inlet forcing, together with a Whiener-Hopf inversion should be used [106]. At the oscillation frequency of 600 Hz, the delay and the phase obtained between peak flow rate oscillation and peak global heat release is 0.6761 ms or 2.548 rad respectively and they compare quite well Thermo-acoustic stability prediction: 3D Helmholtz solver: In this section, a detailed thermo-acoustic study of the annular helicopter combustor is provided by the use of a 3D Helmholtz solver. First, the influence of the geometry of the chamber in a non reactive flow is investigated. It is shown that the casing, the swirler and the primary holes have a strong influence on the acoustics of the azimuthal mode. A methodology to compute azimuthal instabilities in the context of Helmholtz solvers is then presented. The main advantages of this approach is its low CPU time cost. The coupling between acoustics and combustion is accounted for using a classical n -τ model. These n and τ parameters are computed by post processing the previous single sector LES and the corresponding fields are extended to multi-injection annular combustion chamber under the Independence Sector Assumption in Annular Combustor [START_REF] Sensiau | [END_REF] (ISAAC). This assumption is validated thanks to the LES of the full annular helicopter combustion chamber discussed above. Then the stability of the first azimuthal mode is investigated.

Helmholtz equations being essentially elliptic in their nature, the solution of the model equation strongly depends on the geometry and boundary conditions. This raises a simple but critical question in the context of industrial systems: what computational domain, technological devices, asperities and boundary conditions are needed to properly capture the right physics as observed in LES or on the engine ? To illustrate the impact of this critical step several 3D Helmholtz simulations are provided in Table 1.2. For these results combustion effects (FTF) are not taken into account and geometrical complexity is increased while keeping the boundary conditions identical. The vector P gives the pressure distribution in the domain, the real part of ω gives the eigen-frequency of the mode and the imaginary part of ω the growth rate or the damping of the eigen-mode due to acoustic flux at the boundaries ∂Ω Z . In the calculation presented in Table 1.2 there is no acoustic flux at the boundaries since Z = ∞ is imposed everywhere. The combustion chamber (CC) and the casing (C) are first computed separately. Then they are connected with the swirler (CC+C+S) and finally primary holes are taken into account in the calculation (CC+C+S+PH).

In this particular case the frequency of the first azimuthal mode remains within a narrow band between 500 Hz and 700 Hz which is typically the range of frequency involved in combustion instabilities. A first approximation of the frequency can be easily obtained using the formula: f approx = c mean /(2πR mean ), where c mean stands for the mean sound speed in the domain and R mean is the mean radius of the annular configuration. With c mean = 750 m/s, f approx = 680 Hz. This approximation is in the range of the results of the full calculation given in Table 1.2, as well as the frequency observed in the full annular LES. It is clear that the frequency of the first azimuthal mode is greater for CC+C+S than for CC and it is lower than for C. Note also that the computed modes are not purely azimuthal and that the azimuthal modes also have a longitudinal component. Adding the

THERMO-ACOUSTIC INSTABLILITIES

primary holes in the calculation changes this longitudinal component as well as the mode frequency. Those results emphasize the necessity to treat carefully all the geometrical details to model correctly the acoustics involved in combustion instabilities.

The coupling between acoustics and combustion is accounted for using a classical n -τ model. In its global formulation, it relates the fluctuations of the total heat release to the fluctuating velocity at a reference point, noted x ref . As underlined previously, this reference point must be chosen in the injection zone where acoustic fluctuations influence the flow rate or the equivalence ratio and it must be as close as possible to the combustion zone [START_REF] Truffin | [END_REF]. û1 (x ref ). n is the longitudinal velocity leaving the swirler. The time delay, τ , controls the stability of the configuration according to the Rayleigh criterion [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF]. The local n -τ approach allows us to account for the inhomogeneities of the flame response. Obviously, in this full annular chamber, the coupling between the fluctuating heat release and the acoustic fluctuation cannot be described by a flame transfer function involving only one point of reference. The local formulation is thus extended as follows, (1. 46) In this approach, the annular combustion chamber is split in 15 sectors and the Independence Sector Asumption in Annular Combustor (ISAAC) is assumed. The heat release fluctuations in a given sector are driven only by the fluctuating mass flow rates due to the velocity perturbations through its own swirler. The ISAAC assumption is implicitly used in most studies of annular combustors [START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF][START_REF] Morgans | [END_REF][START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF]. It is true only if flames issuing from neighboring burners do not interact, a property which is known to be false in certain cases [START_REF] Poinsot | Analyse des instabilités de combustion de foyers turbulents prémélangés[END_REF] but seems to be accceptable in this case, as presented above in the massively parallel LES of the full helicopter combustion chamber, Fig. 1.14. Based on ISAAC, FTF can be retrieved from the single-sector LES presented above and used in the context of the full annular Helmholtz solver as described above.

Note that in the following section, τ will be used in its global form and a sensitivity analysis of the combustor stability to this parameter is performed. The reaction index coefficient is used in its local form, n l (x), in order to take into account the space inhomogeneities of the flame response to a given excitation.

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

Under the ISAAC assumption, the local amplitude single sector response coefficient n l (x) is duplicated in all sectors and test values for τ are chosen between 0 and 2 ms. Because of the geomtry which azimuthally periodic, the stability behavior in this direction is also periodic in τ . Results are plotted in Fig. 1.16. When τ varies the frequency of the azimuthal eigenmode varies between 555 Hz and 610 Hz. When τ is in [0, T /2], with T the period, the azimuthal mode is damped. The delay measured in the full LES and in the single sector LES is ≈ 0.65 ms so that the burner is obviously operating very close to its stability limits. The same stability ranges are predicted analytically in a simple annular configuration with an infinitely thin flame [START_REF] Krebs | Thermoacoustic analysis of annular combustor[END_REF].

Conclusions

The prediction of thermo-acoustic instabilities in aeronautical gas turbine engines and more generally premixed and partially premixed burners is not an easy task. As of today, industrials are faced with such instabilities at the end of the design process. Recent developments in numercial applications allow to partially apprehend such problems. In the present document, a strategy based on the use of the fully unsteady LES approach and Helmoholtz solvers is proposed. Applications to a single sector and full annular gas turbine helicopter combustion chamber proves the methodology to be applicable. Technological details such as the swirler geometry, primary and secondary holes... are proven to be of importance in the determination of the eigen-modes and eigen-frequencies of the acoustics field obtained with the Helmholtz solver. With the help of the ISAAC hypothesis the prediction of the azi,uthal insta-bilities is also accessible and in very good agreement with full annular LES for the same configuration. Note however that current limitations, applicable to LES and Helnholtz solvers, rely on the proper determination of the acoustic impedances at the inlet and outlet of the computational domains. Simple approximations are possible but their extension to more complex systems is not so clear and further investigations seem needed. Finally, the estimation of the FTF, needed for the Helmholtz solver to determine the stability of the eigen-modes, remains a critical point although it can be retrieved from single sector LES for a given forcing frequency, determination of the FTF to the entire frequency range of interest inferring large compter costs.

Figure 1 . 1 :

 11 Figure 1.1: Conceptual representation of (a) RANS and (b) LES applied to a homogeneous isotropic turbulent field.

1

 1 Equation (1.22) contains four unknowns A u , A p , τ u and τ p which depend on the point where velocity and pressure fluctuations are measured. These parameters may be determined when adding a new state to have enough equations exactly as in the PA method.

Figure 1 . 3 :

 13 Figure 1.3: Schematic view of the laminar premixed flame configuration used to validate the PA and extended n -τ models.

Figure 1 . 4 :

 14 Figure 1.4: The four coefficients of the transfer matrix M of a laminar burner for different positions of the reference point x a . Solid line: absolute value; dashed line: phase [rad]. Line and symbols: PA; line: extended n -τ [95].

1. 1 . 17 Figure 1 . 6 :

 11716 Figure 1.6: Comparison of LES, RANS and experiment phases [52].

  ) where [P ] is the column vector containing the nodal values of the eigenmode at frequency ω and [A] and [C] are square matrices depending only on the discretized geometry of the combustor and mean flow fields c and ρ. Matrix [B] contains information related to the boundary conditions and thus depends on ω since in general Z is frequency dependent. Matrix [D] contains the unsteady contribution of the flame, i.e., Ω , and usually depends non-linearly on the mode frequency ω, see Eq. (1.27

. 37 )Figure 1 . 7 :

 3717 Figure 1.7: Array of circular apertures, of diameter 2a and aperture spacing d, with a bias flow of speed U .

CHAPTER 1 .Figure 1 . 10 :

 1110 Figure 1.10: Definition of the flow domain for the derivation of the quasi-1D approximation.

Figure 1 . 8 :

 18 Figure 1.8: Quasi-1D flow domain for the computation of impedance.

1. 1 .

 1 THERMO-ACOUSTIC INSTABLILITIES 25 4. solve the corresponding linear system Eq. (1.43), 5. compute the acoustic equivalent impedance as Z in = p/ρc û assessed at x = x in .

Figure 1 . 9 :

 19 Figure 1.9: The annular helicopter combustor: (a) full annular view of the computational domain, (b) identification of one sector composing the full annular combustor, (c) detailed view of one single sector and (c) the expected flow distribution within the sector.

Figure 1 .

 1 Figure 1.11: (a) Angular variation of P rms P mean on a ring passing through probes located in the lower part of the chamber casing and (b) pressure signal phase for probes located in front of all swirlers ( Fig. 1.10). Note that sector 1 is used as reference: -model, • LES.

Figure 1 . 12 :

 112 Figure 1.12: Detailed view of half of the burners. Top : Temperature field with temperature isocontour. Bottom: Pressure fluctuations with p = 0 isoline.

Figure 1 .

 1 Figure 1.13: (a): -Maximum flame position (along the burner axis) and Average heat release; (b) • Heat release and -Pressure fluctuation average. All results are obtained for sector 1.

Figure 1 . 14 :

 114 Figure 1.14: Sector by sector flame response to the self-exited mode: • Response amplitude • Response delay (ms).

Figure 1 . 15 :

 115 Figure 1.15: Inhomogeneities of the amplitude response coefficient, n(x).

  A cut of the first azimuthal mode Table1.2: Influence of the geometry on the first azimuthal eigenmode calculated by Helmholtz solver (light grey denotes a pressure anti-node and black gray a pressure node).

  ω)e iωτ (x) û (x ref 1 ). n 1 U bulk for x ∈ Sector 1, n(x, ω)e iωτ (x) û (x ref 2 ).n 2 U bulk for x ∈ Sector 2, ... n(x, ω)e iωτ (x) û (x ref 15 ). n 15 U bulk for x ∈ Sector 15.

Figure 1 . 16 :

 116 Figure 1.16: Influence of the parameter τ on the frequency and the stability of the first azimuthal mode: (a) eigenfrequency, • with and --without combustion; (b) growth rate, with and --without combustion.

Resolution effects on LES of real configurations have been addressed in[130,[START_REF] Boudier | Mesh dependency of turbulent reacting large-eddy simulations in a gas turbine combustion chamber[END_REF].1.1. THERMO-ACOUSTIC INSTABLILITIES

Paschereit et al[START_REF] Schuermans | Thermoacoustic modeling and control of multiburner combustion systems[END_REF][START_REF] Schuermans | Non-linear combustion instabilities in annular gas-turbine combustors[END_REF] propose a non-linear theoretical approach showing that standing wave modes can be found at low oscillation amplitudes but that only one rotating mode is found for large amplitude limit cycles.

[START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems-an overview[END_REF] A simple model can be used to recover the amplitude and phase of the pressure signals found in the LES for the different burners by considering two counter-rotating acoustic modes P + and P -such that:

CHAPTER 1. INSTABILITIES IN LEAN FLAMES

Bibliography