F Nicoud
email: nicoud@math.univ-montp2.fr

L Benoit

C Sensiau

T Poinsot

Lieuwen

Zinn

Poinsot Veynante

Stow Dowling

Benoit

Duchamp De La Geneste

Huang Yang

Pierce Moin

Acoustic modes in combustors with complex impedances and multidimensional active flames

This paper presents a method for computing the thermo-acoustic modes in combustors. In the case of a non-isothermal reacting medium, the wave equation for the pressure fluctuations contains a forcing term related to the unsteady heat release. Depending on the phase relationship between the acoustics and the flame, certain linear modes may become unstable, leading to thermo-acoustic instabilities. The relevant Helmholtz equation is derived and two approaches for solving the corresponding non-linear eigen-

value problem are proposed. The first one is based on an asymptotic expansion of the solution, the baseline being the acoustic modes and frequencies for a steady (or passive) flame and appropriate boundary conditions. This method allows a quick assessment of any acoustic mode stability but is valid only for cases where the coupling between the flame and the acoustic waves is small in amplitude. The second approach is based on an iterative algorithm where a quadratic eigenvalue problem is solved at each sub-iteration.

It is more CPU demanding but remains valid even in cases where the response of the flame to acoustic perturbations is large. Frequency dependent boundary impedances are accounted for in both cases. A parallel implementation of the Arnoldi iterative method is used to solve the large eigenvalue problem that arises from the space discretization of the Helmholtz equation. Several academic and industrial test cases are considered to illustrate the potential of the method.

I. Introduction

Combustion oscillations are frequently encountered during the development of many combustion chambers for gas turbines. [START_REF] Poinsot | Vortex driven acoustically coupled combustion instabilities[END_REF][START_REF] Candel | Combustion instabilities coupled by pressure waves and their active control[END_REF][START_REF] Crighton | Modern methods in analytical acoustics[END_REF][START_REF] Poinsot | Theoretical and numerical combustion[END_REF] To date, these oscillations can not be predicted at the design stage and correcting actions at later stages proved to be extremely costly. Testing burners in simplified combustion chambers is a common method to verify their stability but is also an ambiguous approach because a given burner can produce unstable combustion in one chamber and not in another. Gas turbine manufacturers have devoted large efforts in order to handle combustion instabilities issues and a huge amount of expertise has been gathered by the major companies. [START_REF] Lieuwen | Combustion instabilities in gas turbine engines[END_REF] Still, the mechanisms are not fully understood yet and predictive methods providing stability analysis at the design level are requested. Table 1 gathers the main characteristics of the methods developed in the past.

A current practice is to model the geometry of the combustor by a network of homogeneous (constant density) 1D or 2D axisymmetric acoustic elements where the acoustic problem can be solved [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Stow | Thermoacoustic oscillations in an annular combustor[END_REF][START_REF] Polifke | Constructive and Destructive Interference of Acoustic and Entropy Waves in a Premixed Combustor with a Choked Exit[END_REF] analytically. The flame is supposed infinitely thin and only appears at interfaces between the low and large temperature segments. Jump relations are used to connect all these elements, enforcing pressure continuity and flow rate conservation. The amplitudes of the forward and backward acoustic waves in each segment are determined so that all the jump relations and the boundary conditions are satisfied. This can only be achieved for a discrete set of pulsations ω which are the roots of a dispersion relation in the complex plane. The main advantage of this approach is that it allows to describe a complex system with a few lumped elements only, leading to a very low order model that can be extensively used for pre-design purpose. Extension to azimuthal modes 9, 10 and 1.5D networks [START_REF] Benoit | Prédiction des instabilités thermo-acoustiques dans les turbines à gaz[END_REF] have been proposed recently. The main drawback of this type of approach is that the geometrical details of the combustor cannot be accounted for and only the first "equivalent" longitudinal or orthoradial modes are sought for.

An alternative is to perform Large Eddy Simulations (LES) which are well suited to study the dynamics of turbulent flames (see recent books on turbulent combustion [START_REF] Poinsot | Theoretical and numerical combustion[END_REF][START_REF] Peters | Turbulent combustion[END_REF] ). Multiple recent papers have demonstrated the power of these methods. [START_REF] Murota | Large-eddy simulation of combustion oscillation in premixed combustor[END_REF][START_REF] Desjardins | Two dimensional Large Eddy Simulation of soot formation in the near field of a strongly radiating nonpremixed acetylene-air jet flame[END_REF][START_REF] Angelberger | Large Eddy Simulations of chemical and acoustic effects on combustion instabilities[END_REF][START_REF] Caraeni | Modeling of Liquid Fuel Injection, Evaporation and Mixing in a Gas Turbine Burner Using Large Eddy Simulation[END_REF][START_REF] Colin | A thickened flame model for large eddy simulations of turbulent premixed combustion[END_REF][START_REF] Pitsch | Large Eddy Simulation of Premixed Turbulent Combustion using a level-set approach[END_REF][START_REF] Huang | Bifurcation of flame structure in a leanpremixed swirl-stabilized combustor: transition from stable to unstable flame[END_REF][START_REF] Pierce | Progress-variable approach for large eddy simulation of non-premixed turbulent combustion[END_REF][START_REF] Selle | Compressible Large-Eddy Simulation of turbulent combustion in complex geometry on unstructured meshes[END_REF] However, an important limitation of LES is its cost: the intrinsic nature of LES (full three-dimensional resolution of the unsteady Navier-Stokes equations) makes it very expensive, even on today's computers. Moreover, even when they confirm that a combustor is unstable, LES does not say why and how to control it. Therefore, tools are needed to analyze LES results and to provide capacities for optimization and control of thermo acoustic oscillations in chambers.

A set of linear transport equations for the perturbations of velocity, temperature and density can be derived by linearizing the Navier-Stokes equations. [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] The local unsteady heat release appears as a forcing term in the linearized energy equation and is responsible for combustion noise and thermo-acoustic instabilities. Assuming that an appropriate model for the unsteady heat release is available, the system of linear partial differential equations (PDE) for the fluctuating quantities is closed and can be solved. The most natural way to proceed is then to use a finite element based method to discretize and solve it in physical space. [START_REF] Pankiewitz | Time domain simulation of combustion instabilities in annular combustors[END_REF] Starting from initial random fields of pressure, velocity and density or temperature, each PDE is solved in the time domain, discretizing the time derivatives by a Runge-Kutta algorithm. 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 fields can be amplified and grow exponentially. After a while, the fluctuations are dominated by the most amplified one (or the less damped one if the flame/acoustics coupling stabilizes all the fluctuations) and it is possible to examine the spatial structure, frequency of oscillation and growing factor of this particular thermo-acoustic mode. This intuitive approach provides a natural way to account for the coupling with the flame as well as non-linear saturation effects. [START_REF] Pankiewitz | Time domain simulation of combustion instabilities in annular combustors[END_REF] However, accounting for frequency-dependent impedances at boundaries is not straightforward when performing the time domain integration. Moreover, from an engineering point of view, it is certainly not very convenient to gain knowledge only about the most unstable modes: experience shows that several modes are most often present or at least dangerous when thermo-acoustic instabilities appear.

A proper framework to analyze combustion stability is the wave equation in a reacting flow. [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] A natural way of gaining information about the whole set of thermo-acoustic modes is then to consider the Helmholtz equation, the frequency domain version of the wave equation. This equation can be derived by combining the linearized version of the equation of state with the set of linear PDE's for each fluctuating quantities. Under the low-Mach number assumption for the mean velocity field, an approximate equation controlling the propagation of pressure perturbations in a reacting flow is then obtained and, considering the associated Helmholtz equation, the initial value problem in the physical domain is transformed into an eigenvalue problem. The shape of the thermo-acoustic modes is related to the eigenfunctions and their frequency and growing rate are determined by the eigenvalues. In the case of classical acoustics in a homogeneous non-reacting medium, finite element based methods are often used to perform the spatial discretization and transform the Helmholtz equation into a finite dimension eigenvalue problem. 31 When combustion occurs, the flame/acoustic coupling term makes the eigenvalue problem non-linear with respect to the pulsation of the mode and classical methods cannot be applied to solve the Helmholtz equation. The "linear Galerkin" method proposed by Culick 22 consists in writing the solution of the thermo-acoustic problem as a linear combination of the acoustic modes of the homogeneous (without combustion forcing) Helmholtz equation. Assuming that a) the eigenvectors of the homogeneous wave equation form an orthogonal basis of the linear space of the square integrable functions defined over the flow domain and b) that the effect of the flame coupling is to slightly perturb the acoustic modes, Culick showed that one can assess the frequency shift induced by the unsteady heat release. The stability of the modes is then obtained by studying the sign of the imaginary part of the frequency shift. This method suffers from two major drawbacks related to the two assumptions given above. Assumption a) is in fact only valid for particular choices of the boundary impedance 11 (see appendix B for more details). Classical examples of such boundary conditions are pressure imposed outlet and impermeable wall. However, in the case of a general complex boundary impedance, this condition is not fulfilled and the Galerkin method can hardly be generalized. Assumption b) is only valid when the amplitude of the flame response to acoustic perturbations is small. However, considering the limiting case of small frequency oscillations, it can be shown that this assumption is most likely not valid for practical cases. [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF] In this paper, two methods to solve the thermo-acoustic problem written in the frequency space are described in section IV. The first one is based on an asymptotic expansion of the solution around the passive flame case and still requires assumption b) but not assumption a). The second is based on an iterative algorithm and is valid even if neither a) nor b) is fulfilled. The basic hypothesis are first recalled in section II where the derivation of the inhomogeneous Helmholtz equation is also given. The spatial discretization and linear algebra method are then described in section III while basic academic test cases are considering in section V. Typical results obtained in the case of an industrial annular combustor are also described to illustrate the 3D capabilities of the method.

II. Basic assumptions and equations A. Basic equations

To simplify the derivation, we will consider a gas mixture where all species share the same molar weight and heat capacity. This assumption is not necessary to derive the generalized acoustic wave equation although it makes the algebra simpler. It is valid for air flames but must be revisited for the case of H 2 -O 2 mixtures for example. A direct consequence is that the heat capacities difference r = C p -C v is constant even if C p , C v and γ may depend on temperature. Viscous terms (molecular diffusion of momentum and heat) will also be neglected in the present analysis as it is usual for the analysis of low frequency acoustic perturbations.

Under the above assumptions, the mass, momentum and entropy equations read respec-tively:

Dρ Dt = -ρ∇ • u, (1)
ρ Du Dt = -∇p, (2)
Ds Dt = rq p . (3)
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 , (4)
these transport equations describe the spatio-temporal evolutions of all relevant physical flow quantities.

Eq. 1 to 4 can be linearized by considering the simple case of large scale small amplitude fluctuations (index 1) super-imposed to a zero Mach number mean flow (index 0) which depends only on space. The instantaneous pressure, density, temperature, entropy and velocity fields can then be written as p = p 0 + p 1 , ρ = ρ 0 + ρ 1 , s = s 0 + s 1 and u = u 1 where the quantities p 1 /p 0 , ρ 1 /ρ 0 , T 1 /T 0 , s 1 /s 0 and √ u 1 • u 1 /c 0 are of order ǫ, where ǫ ≪ 1 and c 0 = γp 0 /ρ 0 is the mean speed of sound. For simplicity, the temporal fluctuations of the heat capacities are neglected. The zero Mach number assumption (u 0 ≃ 0) is valid as soon as the characteristic Mach number M = √ u 0 • u 0 /c 0 of the mean flow is small compared to L f /L a where L f is the thickness of the reaction zone and L a is the typical acoustic wavelength (see appendix A for more details). This result makes the approximation M ≃ 0 reasonable for many applications. Moreover the zero Mach number assumption implies that ∇p 0 = 0 (from Eq. 2) and q 0 = 0 (from Eq. 3 ), the latter condition being acceptable because only the fluctuating quantities are of interest in the present analysis. Finally this assumption implies that the approximation D/Dt ≈ ∂/∂t holds for any fluctuating quantity since, with u 0 ≃ 0, the non linear convective terms are always of second order in ǫ.

B. Linearization

Injecting the above expansions for the instantaneous flow quantities into Eqs 1 to 4 and keeping only terms of order ǫ, one obtains the following set of linear equations for the fluctuating quantities ρ 1 , u 1 , s 1 and p 1 :

∂ρ 1 ∂t + u 1 • ∇ρ 0 + ρ 0 ∇ • u 1 = 0, (5)
ρ 0 ∂u 1 ∂t + ∇p 1 = 0, (6)
∂s 1 ∂t + u 1 • ∇s 0 = rq 1 p 0 . (7)
The linearized state equation and entropy expression are:

p 1 p 0 - ρ 1 ρ 0 - T 1 T 0 = 0 and s 1 = C p T 1 T 0 -r p 1 p 0 . (8)
Eq. 7 shows that the fluctuating flow is not isentropic as soon as either:

• the mean entropy gradient ∇s 0 is not orthogonal to the velocity fluctuations or

• the heat release has an unsteady component q 1 = 0.

In practical applications the geometry of the flow domain is fully 3D and highly complex so that the first condition is most probably matched as soon as ∇s 0 = 0. Under the zero Mach number assumption the mean pressure gradient is null and from Eq. 4 the mean entropy gradient is approximately (neglecting the gradient of C p and C v ): ∇s 0 ≃ -C p ∇ρ 0 /ρ 0 = C p ∇T 0 /T 0 . This shows that the u 1 • ∇s 0 term in Eq. 7 is non zero in the reacting zone, even if there is no heat release, viz. q 1 = 0 (passive or steady flame). The RHS term of the same equation is non zero only if the flame has an unsteady behavior (active or unsteady flame).

In order to close the set of equations 5, 6, 7 and 8, a model must be used to express the unsteady heat release q 1 in terms of the other fluctuating quantities.

C. Flame response

Modeling the unsteady behavior of the flame is the most challenging part in the description of thermo-acoustic instabilities. [START_REF] Lieuwen | Modeling Premixed Combustion-Acoustic Wave Interactions: A Review[END_REF] Several models have been proposed in the past to describe the response of conic or V-shape laminar flames, [START_REF] Schuller | A unified model for the prediction of laminar flame transfer functions: comparisons between conical and V-flames dynamics[END_REF] accounting for non-linear saturation effects [START_REF] Dowling | Nonlinear self-exited oscillations of a ducted flame[END_REF] ] and equivalence ratio fluctuations. [START_REF] Lieuwen | The Role of Equivalence Ratio Oscillations In Driving Combustion Instabilities In Low NOx Gas Turbines[END_REF][START_REF] Sattelmayer | Influence of the combustor aerodynamics on combustion instabilities from equivalence ratio fluctuations[END_REF] Most models available so far (if not all of them) 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 the flame region. The idea behind this modeling approach is that heat release depends on the flame surface which, the flame speed being given, is mainly controlled by the fresh gas flow rate. The most classical model follows seminal ideas by Crocco 38,[START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. Part II[END_REF] and is referred to as the n -τ model. This is essentially a 1D formulation which 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:

Q(t) = Ω q 1 (t)dΩ = S ref γp 0 γ -1 × n × u 1 (x ref , t -τ ). (9)
In this expression, Q(t) is the heat release integrated over the flow domain Ω, S ref is the cross section area of the burner mouth, u 1 denotes the velocity component along 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:

R = t Ω p 1 q 1 dΩdt. ( 10
)
According to the classical Rayleigh criterion, flame/acoustics coupling promotes the appearance of instabilities in cases where 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, viz. cases 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 matched for modern combustors with high efficiency. It is then natural to use a local flame model which would relate the local unsteady heat release to a reference acoustic velocity in the injector mouth. The natural way to proceed is then to write:

q 1 (x, t) q tot = n u (x) u 1 (x ref , t -τ u (x)) • n ref U bulk , (11)
where n u (x) and τ u (x) are fields of interaction index and time lag and n ref is a fixed unitary vector defining the direction of the reference velocity. The scaling by the total heat release q 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 local data is required to tune this kind of model and for obvious technological reasons these data can hardly be obtained experimentally. The alternative is then to use compressible reacting LES to investigate the response of a turbulent flame submitted to acoustic perturbations. By performing a spectral analysis of the unsteady fields of heat release and velocity, it is then possible to determine the optimal parameter fields n u (x) and τ u (x) to match the actual flame response by using Eq. 11. [START_REF] Varoquie | Experimental analysis and large eddy simulation to determine the response of non-premixed flame submitted to acoustic forcing[END_REF][START_REF] Truffin | Measurements of transfer functions in reacting flows using Large Eddy Simulations[END_REF] Using the local flame model given in Eq. 11, Eq. 7 can then be re-written as:

∂s 1 ∂t + u 1 • ∇s 0 = r p 0 q tot U bulk n u (x)u 1 (x ref , t -τ u (x)) • n ref , (12)
and the set of equations 5, 6, 8 and 12 can be solved to determine the thermo-acoustic properties of the system.

D. Wave equation

Taking the time derivative of Eq. 5, adding the divergence of Eq. 6 over ρ 0 and using Eqs. 8 and 7 to eliminate ρ 1 yields the following wave equation for the pressure fluctuations p 1 :

∇ • 1 ρ 0 ∇p 1 - 1 γp 0 ∂ 2 p 1 ∂t 2 = - γ -1 γp 0 ∂q 1 ∂t . ( 13
)
If the zero Mach number assumption is not made, it is not possible to manipulate the set of equations 5, 6, 7 and 8 in order to obtain a scalar wave equation. It is then necessary to consider the original set of equations, either written in the time domain or in the frequency space, in order to solve the thermo-acoustic problem. Note also that no assumption has been made about the spatial evolution of the isentropic coefficient γ to derive Eq. 13. If γ is constant over space, one recovers the classical equation for inhomogeneous medium with the elliptic term being ∇ • (c 2 0 ∇p 1 ). Eq. 13 being linear, it is natural to introduce harmonic variations at frequency f = ω/(2π) for pressure, velocity and local heat release perturbations:

p 1 = ℜ (p(x) exp(-iωt)) , u 1 = ℜ (û(x) exp(-iωt)) , (14)
q 1 = ℜ (q(x) exp(-iωt)) .
Introducing Eq. 14 into Eq. 13 lead to the following Helmholtz equation :

∇ • 1 ρ 0 ∇p + ω 2 γp 0 p = iω γ -1 γp 0 q(x), (15)
where ρ 0 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 pulsation ω. In the frequency space, the zero Mach number assumption induces iω û = ∇p/ρ 0 and the flame model Eq. 11 translates into : q

(x) = q tot iωρ 0 (x ref )U bulk n u (x)e iωτu(x) ∇p(x ref ) • n ref . (16)
Introducing Eq. 16 into Eq. 15 leads to :

∇ • 1 ρ 0 ∇p + ω 2 γp 0 p = γ -1 γp 0 q tot ρ 0 (x ref )U bulk n u (x)e iωτu(x) ∇p(x ref ) • n ref . (17)
This paper will focus on the resolution of Eq. 17 but the methodologies developed can be applied to the more general case where the complex amplitude of the heat release is given by: q

(x) = Lu [∇p(x)] + Lp [p(x)] , (18)
where Lu and Lp are two linear operators acting on ∇p and p respectively. The wave equation for the complex pressure amplitude is then

∇ • 1 ρ 0 ∇p + ω 2 γp 0 p = iω γ -1 γp 0 Lu [∇p] + Lp [p] , (19)
which reduces to Eq. 17 for the particular choice:

Lp : ĝ(x) -→ 0, ( 20
) Lu : v(x) -→ q tot iωρ 0 (x ref )U bulk n u (x)e iωτu(x) v(x ref ) • n ref .
Of course the general formulation Eq. 18 has the potential to include more general effects than the local n -τ model described by Eqs. 11 and 16. For example, defining Lu as in Eq. 20 and Lp as:

Lp : ĝ(x) -→ q tot iωp 0 n p (x)e iωτp(x) ĝ(x ref ), (21)
would allow 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. [START_REF] Polifke | Reconstruction of acoustic transfer matrices by instationnary computational fluid dynamics[END_REF][START_REF] Paschereit | Measurement of transfer matrices and source terms of premixed flames[END_REF] In Eq. 21, n p (x) and τ p (x) are fields of interaction index and time delay describing the effect of acoustic pressure on the unsteady heat release. 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 | Comparison and extension of methods for acoustic identification of burners[END_REF] Note finally that the fields of interaction index n u (x), n p (x) and time delay τ u (x), τ p (x) may also depend on the frequency of the perturbations acting on the flame, although this dependency has not been explicitly written for simplicity.

E. Boundary conditions

Denoting by n BC 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 outer pressure is imposed strongly to the flow domain, zeroing the pressure fluctuations:

p = 0, on boundary ∂Ω D , (22)
• Zero normal velocity, viz. û • n BC = 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. 6 can be used to re-write this condition as a Neumann condition for the acoustic pressure:

∇p • n BC = 0, on boundary ∂Ω N , (23)
• Imposed reduced complex impedance Z = p/ρ 0 c 0 û • n BC . Under the zero Mach number assumption, this condition can be re-written as a Robin condition for the acoustic pressure:

c 0 Z∇p • n BC -iω p = 0, on boundary ∂Ω Z , (24)
Associated with the homogeneous boundary conditions 22, 23 and 24 on ∂Ω = ∂Ω D ∂Ω N ∂Ω Z , equation 17 defines a non-linear eigenvalue problem which is the basis of the numerical tool described in the present paper.

III. Numerical methods for the passive flame problem

The major difficulty in solving the thermo-acoustic problem given in Eqs. 19 or 17 comes from the strong non-linearity of the RHS term with respect to the pulsation ω. As a first step, it is then natural to consider the purely acoustic problem obtained by neglecting the unsteady flame effect. This means taking n u (x) = n p (x) = 0 (or more generally Lu = Lp = 0) and finding p(x) and ω satisfying:

∇ • 1 ρ 0 ∇p + ω 2 γp 0 p = 0, (25)
with boundary conditions 22, 23 and 24. Since complex geometries must be handled properly, a finite element based method is first used to discretize Eq. 25 and obtain a finite-dimension eigenvalue problem which can then be solved by using an appropriate linear algebra method.

A. Spatial discretization

A finite element strategy is used to discretize the exact geometry of the combustor so that no assumption is made a priori regarding the shape of the modes. This feature gives the Helmholtz solver the potential to test the effect of geometrical changes on the stability of the whole system. Denoting the flow domain and its boundary by Ω and ∂Ω respectively, lets consider a finite element mesh made of linear elements (triangles in 2D, tetrahedra in 3D). The mesh consists of N e elements Ω j (j ∈ 1, ..., N e ), N v vertices v k (k ∈ 1, ..., N n ) of which N d are located on the boundary ∂Ω D . We call S v the set of vertices of the mesh which do not belong to ∂Ω D and S Ω the set of all the elements. Also, for each node v j ∈ S v , we denote by φ j the piecewise linear shape function which equals 1 at node v j and 0 for all v k =j . The complex unknown function p(x) is then approximated by p(x) = j:v j ∈Sv pj φ j , the homogeneous Dirichlet boundary condition on ∂Ω D being then automatically satisfied. Finding p(x) for every x ∈ Ω is then equivalent to determine the N = N v -N d complex coefficients pj . To do so, we follow the classical Galerkin finite element method: start from Eq. 25, replace p(x) by its approximation, multiply by the test function φ k and integrate over the flow domain to obtain:

∀k : v k ∈ S v , Ω φ k ∇• 1 ρ 0 ∇ j:v j ∈Sv pj φ j dx + ω 2 Ω φ k γp 0 j:v j ∈Sv pj φ j dx = 0. ( 26
)
Integrating the first term by parts, inverting the integration and summation operations and making use of Eq. 22-24, one obtains :

∀k : v k ∈ S v , - j:v j ∈Sv Ω 1 ρ 0 ∇φ k • ∇φ j dx pj +iω j:v j ∈Sv ∂Ω 1 ρ 0 c 0 Z φ k φ j dσ pj + ω 2 j:v j ∈Sv Ω 1 γp 0 φ k φ j dx pj = 0. ( 27
)
Note that only ∂Ω Z contributes to the boundary term because φ k (x) = 0, ∀k : v k ∈ S v , ∀x ∈ ∂Ω D and because ∇p • n BC = 0 on ∂Ω N . Introducing the symmetric matrices A, B and C of size N and of generic element:

A kj = - Ω 1 ρ 0 ∇φ k ∇φ j dx, B kj = ∂Ω Z i ρ 0 c 0 Z φ k φ j dσ, (28)
C kj = Ω 1 γp 0 φ k φ j dx,
the following discrete form of the acoustic problem defined by Eqs. 25, 22, 23 and 24 is obtained:

AP + ωB(ω)P + ω 2 CP = 0, ( 29
)
where P is the column vector containing the N complex coefficients pj . Note that A and C are real matrices while B is complex except in the particular case where the reduced impedance Z is purely imaginary.

B. Linear algebra

In general, the boundary impedance Z is a function of the pulsation and the matrix B depends on ω. Equation 29 then constitutes a non-linear eigenvalue problem of size N which can be formally written as:

Find ω ∈ C and P = 0 such that R(ω

)P = 0, ( 30
)
where R is a square matrix of size N which depends on ω. No general algorithm exists for solving this type of eigenvalue problems except for particular dependencies of the operator R on ω. For example, when R is polynomial (of degree d) in ω, the problem 30 can be re-written as a linear eigenvalue problem of the form SQ = ωQ where Q is a matrix of size d × N which does not depend on ω and Q is a column vector of length d × N whose N first rows contains P.

In the particular cases where Z does not depend on ω, Eq. 29 is a quadratic eigenvalue problem since this equation can be re-written as R(ω)P = 0 where R = A + ωB + ω 2 C and A, B and C are constant matrices. Quadratic eigenvalue problem are rather well known from a theoretical point of view. [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF] One efficient way of solving them numerically is to convert them into an equivalent linear problem of size 2 × N . This can be done for example by introducing the new column vector

P ω = ωP, (31)
and re-writing Eq. 29 as

AP + BP ω + ωCP ω = 0. ( 32
)
Eqs. 31 and 32 can then be re-written under the following matrix form

0 -I A B P P ω + ω I 0 0 C P P ω = 0, ( 33
)
where I is the unit matrix of size N . Thus the solutions (ω, P) of the quadratic problem obtained from 29 when B is constant can be computed from the linear eigenvalue problem of size 2 × N defined by Eq. 33.

Several numerical methods can then be used to assess the eigenmodes. Direct methods are exact and have the advantage to provide all the eigenmodes. For example, the QR algorithm is the most efficient way of producing the Shur decomposition of the matrix defining the linear eigenvalue problem and is the method of choice when the size of the problem is typically less than 1000. Since only the first few frequencies are usually of interest from a physical point of view, it is more appropriate to use an iterative method which can be applied for large problems (N > 10 5 ). For example, Krylov-based algorithms performed a partial Shur decomposition of the matrix of interest and allow the user to focus on a reduced number of eigenvalues. We are using 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.

Note that the assumption that B does not depend on ω is not necessary for Eq. 29 to define a quadratic eigenvalue problem. Indeed, the matrix B being multiplied by ω in Eq. 29, the problem remains quadratic in ω as long as the impedance takes the following form:

1/Z = 1/Z 0 + Z 1 ω + Z 2 /ω, (34)
where Z 0 , Z 1 and Z 2 are complex valued constants. The eigenvalue problem Eq. 29 can be written as:

AP + ωBP + ω 2 CP = 0, ( 35
)
where the matrices A, B and C are defined as:

A ij = A ij + ∂Ω Z iZ 2 ρ 0 c 0 φ i φ j dσ, B ij = ∂Ω Z i ρ 0 c 0 Z 0 φ i φ j dσ, (36)
C ij = C ij + ∂Ω Z iZ 1 ρ 0 c 0 φ i φ j dσ,
and do not depend on the pulsation. The linearization process leading to Eq. 33 can then be applied to A, B and C respectively. Eventually, frequency dependent boundary impedance of the form Eq. 34 can be accounted for by a linear eigenvalue problem of size 2 × N .

IV. Accounting for the active flame effect Solving Eq. 25 as in the previous section means finding the eigenmodes of the burner, taking into account the presence of the flame through the mean temperature field but neglecting its unsteady effects. The boundary conditions are also accounted for and this approximation can provide relevant information about the shape and real frequency of the first modes of the combustor. However, since there is no coupling between the acoustics and the flame, there is no hope to discriminate between stable and unstable modes, which is the ultimate objective when addressing thermo-acoustic problems. The objective of this section is to discuss how the effect of the flame/acoustic coupling can be accounted for. In this case, the appropriate equation is Eq. 17 (or Eq. 19) instead of Eq. 25.

A. Asymptotic method

In this approach, [START_REF] Benoit | Numerical assessment of thermo-acoustic instabilities in gas turbines[END_REF] the flame is considered as a perturbation of the situation without combustion. Under this assumption, eigenfrequencies and eigenmodes (ω, p) of Eq. 17 are sought as first order expansions of the form :

ω = ω 0 + νω 1 + O(ν 2 ), p = p0 + ν p1 + O(ν 2 ), (37)
where (ω 0 , p0 ) are solutions of the homogeneous eigenvalue problem Eq. 25 and ν is a small parameter of the problem and O(ν 2 ) stands for terms of order ν 2 or higher. The ν parameter should be a measurement of the amplitude of the flame response so that the asymptotic analysis naturally retrieves the passive flame results when the flame/acoustic coupling tends to zero. For example in the case of the local n -τ model Eq. 16, it is natural to define ν as:

ν = Ω |n u (x)| dx, (38)
As suggested by global acoustic energy budgets, [START_REF] Poinsot | Theoretical and numerical combustion[END_REF] the eigenmode stability is determined by the whole spatial distribution of the flame and the boundary losses. This means that the asymptotic expansion should not be applied on the local equation 17 but on a global balance equation over the domain. Multiplying Eq. 17 by the complex pressure amplitude p and integrating over space, one obtains:

Ω p∇ • 1 ρ 0 ∇p dΩ + Ω ω 2 γp 0 p2 dx = Ω (γ -1) γp 0 q tot ρ 0 (x ref )U bulk n u (x)e iωτu(x )p∇p(x ref ) • n ref dx. (39)
Injecting Eq. 37 into Eq. 39 and keeping only first order terms in ν gives the following equation:

Ω p0 ∇ • 1 ρ 0 ∇ν p1 + ω 2 0 γp 0 ν p1 dx = -2ν Ω p2 0 ω 0 ω 1 dx + Ω (γ -1) γp 0 q tot ρ 0 (x ref )U bulk n u (x)e iω 0 τu(x) p0 ∇p 0 (x ref ) • n ref dx. (40)
The LHS term can be simplified by using a reduction order method [START_REF] Bender | Advanced mathematical methods for scientists and engineers[END_REF] in which p1 = p0 F 1 , F 1 being a function of x. Thanks to this relation, the LHS term of Eq. 40 becomes:

Ω p0 ∇. 1 ρ 0 ∇ν p1 + ω 2 0 γp 0 ν p1 dx = ν ∂Ω 1 ρ 0 p2 0 ∇F 1 • n BC dσ, (41)
which is obviously null on ∂Ω D since p0 = 0. Moreover, since the eigenmodes with flame (ω, p) and without flame (ω 0 , p0 ) satisfy the same boundary conditions, one can show that ∇F 1 • n BC = 0 on ∂Ω N and that the following relation is valid to first order in ν on ∂Ω Z ,

∇F 1 • n BC = iω 1 c 0 Z(ω 0 ) 1 - 1 Z(ω 0 ) ∂Z ∂ω (ω 0 ) . ( 42
)
Introducing this relation in the RHS term of Eq. 41, an expression for the perturbation νω 1 can be obtained:

νω 1 = q tot ρ 0 (x ref )U bulk Ω (γ -1)n u (x)e iω 0 τu(x) p0 ∇p 0 (x ref ) • n ref dx 2ω 0 Ω p2 0 dx + ∂Ω Z ic 0 p2 0 Z(ω 0 ) 1 -1 Z(ω 0 ) ∂Z ∂ω (ω 0 ) dσ . ( 43
)
In the case where the denominator of Eq. 43 is not null, this equation provides a simple way to check whether an eigenmode without combustion (ω 0 , p0 ) is made stable (ℑ(ω 0 + νω 1 ) < 0) or unstable (ℑ(ω 0 + νω 1 ) > 0) by the coupling with the unsteady flame. It generalizes the linear procedure proposed by Culick [START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems-An overview[END_REF] to cases where the boundary impedance is finite (|Z| = ∞ and Z = 0). Of course it is only valid for cases where the amplitude of the flame response is small, viz. ν << 1. Considering the limiting case of small frequency oscillations, it can be shown that this assumption is most likely not valid for practical cases where the burnt-to-fresh gas temperature ratio is of order 3 or larger. [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF] It is thus desirable to develop a more general method which would remain valid when the amplitude of the flame response is large. Deriving such method is the objective of the next section.

B. Iterative method

In order to handle cases where the amplitude of the flame response is not small, the finite element formalism applied previously to Eq. 25 is now extended to Eq. 17 and associated boundary conditions Eqs. 22-24. The LHS of Eq. 17 and boundary terms on ∂Ω Z will give rise to the LHS of Eq. 29 or Eq. 35 under the assumption Eq. 34. Using again the finite element formalism, the RHS term of Eq. 17 can be discretized to give rise to DP where D is a complex valued square matrix of size N whose generic element is

D kj = Ω γ -1 γp 0 q tot ρ 0 (x ref )U bulk n u (x)e iωτu(x) φ k ∇φ j (x ref ) • n ref dx. (44)
Note that D is a sparse, complex valued matrix which is neither symmetric nor hermitian.

Eventually, the discretized thermo-acoustic problem consists in the following eigenvalue problem of size N :

AP + ωBP + ω 2 CP = D(ω)P. (45)
In general, D depends non-linearly on the pulsation ω and no polynomial approximation can be given for this matrix. Indeed, typical models for the unsteady heat release involve a time delay (between a velocity and/or pressure perturbation and the response of the flame) which converts into an exponential term in ω in the Fourier space, thus in D. As a consequence, the eigenvalue problem given by Eq. 45 can hardly be re-written as a linear eigenvalue problem of larger dimension as done in the previous section and the classical linear algebra algorithm cannot be used to solve the problem efficiently. An option is then to use an iterative method to solve Eq. 45, the k th iteration consisting in solving the quadratic eigenvalue problem in ω k defined as:

(A -D(ω k-1 )) P + ω k BP + ω 2 k CP = 0. ( 46
)
A natural initialization is to set D(ω 0 ) = 0 so that the computation of the modes without 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 pulsation and associated mode. Finally, the algorithm to solve Eq. 45 is formally given by: Algorithm I:

1. Find the first few eigenfrequencies ω 1 0 , ω 2 0 , ω 3 0 , ... without flame coupling by solving Eq. 29, 2. Select one of these frequencies and call it ω 0 , 3. Set D(ω 0 ) = 0 and k = 1, 4. Solve Eq. 46,

Assess the error η

= |ω k -ω k-1 |/|ω 0 |,

iterate on k until η is small enough

To study the convergence of this algorithm, it is appropriate to consider the following eigenvalue problem P ω :

(P ω ) (A -D(ω)) P + ΩBP + Ω 2 CP = 0, ( 47
)
where ω is a complex number and Ω and P are the eigenvalue and eigenvector respectively. Consider also the operator F defined as :

F : C -→ C ω -→ F(ω), (48)
where F(ω) is the eigenvalue of P ω which minimizes the quantity |F(ω) -ω|. It is obvious from the definitions of the problem P ω and operator F that any complex number ω satisfying F(ω) = ω is an eigenfrequency of the non-linear eigenvalue problem 45. As a consequence, if algorithm I converges, it provides a solution to the thermo-acoustic problem. Besides, thanks to the fixed point theorem, algorithm I converges as soon as F is a contracting operator in a neighborhood of ω 0 . The validity of this condition can hardly be shown in the general case since F strongly depends on the flame response through the D term in the problem P ω which stems from the discretization of the coupling term. Indeed, the behavior and regularity of the local interaction index n u (x) and time delay τ u (x) with respect to both x and ω are far from understood. The only way to gain information about this term seems to post-process experimental [START_REF] Paschereit | Measurement of transfer matrices and source terms of premixed flames[END_REF] or numerical data. [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] In any case, obtaining general results about the contracting properties of the operator F from physical arguments is out of reach of the current understanding of the thermo-acoustic instabilities. However, in all the practical cases considered so far, [START_REF] Benoit | Prédiction des instabilités thermo-acoustiques dans les turbines à gaz[END_REF][START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF][START_REF] Martin | Analysis of acoustic energy and modes in a turbulent swirled combustor[END_REF] the above algorithm proved to be stable and to converge in less than 4 -5 iterations, suggesting good contracting properties for F.

V. Numerical results

The methodologies presented in the previous sections have been implemented in a 3D finite element based in-house acoustic solver called AVSP. Computations without active flame effects in 3D geometries and fully-reflecting boundary conditions were presented in Selle et al. [START_REF] Selle | Joint use of Compressible Large-Eddy Simulation and Helmoltz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF] where the ability of the approach to describe turning modes in a swirled combustor was also established. Computations where the acoustic flame coupling is accounted for in 3D geometries were discussed in Martin et al., [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] showing the potential of the method to discriminate between stable and unstable thermo-acoustic modes in a swirled staged com-bustor. These previous studies were illustrations of the joint use of acoustic and LES solvers to predict thermo-acoustic instabilities and understand their connections with hydrodynamic instabilities and acoustic tones. The present paper focuses more on the numerical strategy on which the thermo-acoustic solver is based and the objective of this section is to illustrate the potential of the method. The selected test cases include:

• an isothermal 2D domain with various boundary conditions,

• an anisothermal 2D domain with an active flame and simple boundary conditions,

• an anisothermal 3D industrial configuration with various boundary conditions.

• an anisothermal 3D domain with an active flame and complex valued boundary conditions,

The first two test cases allow a precise comparison of the numerical results with available analytical solutions in cases where 1D and 2D modes are subjected to partially reflecting boundary conditions and where 1D modes are coupled to an infinitely thin active flame. The third one illustrates the 3D capabilities of the formulation and the importance of accounting for complex valued boundary impedances. The fourth one is an illustration of acoustic calculations with complex valued impedance and 3D active flame.

A. Isothermal academic calculations

The computational domain is a rectangular cavity of length L = 0.5 m and height h = 0.1 m where the speed of sound is c 0 = 450 m/s (Fig. 1). The unstructured mesh contains 11367 triangular cells and 5813 nodes and is uniform over the domain. Two sets of boundary conditions are considered in the following sub-sections.

Longitudinal 1D modes

Longitudinal modes are considered first. A homogeneous Neumann condition (|Z| = +∞) is imposed at the bottom and top boundaries. The same condition is used at the left boundary while a complex impedance (Z = a + ib) is imposed at the right side (see Fig. 1). In this case, the problem admits a set of 1D solutions (longitudinal modes). Frequencies: Seeking for 1D harmonic solutions, viz. p(x, t) = ℜ [p(x)exp(-iωt)], the problem can be formulated as follows : and the eigenfrequencies are equal to

∂ 2 p ∂x 2 + ω 2 c 2 0 p = 0 with ∂ p ∂x = 0 for x = 0 and ∂ p ∂x - iω c 0 Z p = 0 for x = L, (49)
f m = m c 0 2L + c 0 2πL arctan( -i Z ), m ∈ N. (50)
In the case where Z is a purely reactive impedance, viz. Z = ib with b ∈ R, Eq. 50 shows that the eigenfrequencies are real, consistently with the fact that the eigenmodes are marginally stable in this case because there is no acoustic flux at the boundaries. On the other hand, when Z is a purely resistive impedance, viz. Z = a with a ∈ R, the eigenfrequencies are complex valued and equal to:

f = m c 0 2L -i c 0 4πL ln a + 1 a -1 , m ∈ N. (51)
If (a + 1)/(a -1) > 0, viz. |a| > 1, the real part of f is a multiple of the half-wave mode frequency c 0 /2L while in the case where (a + 1)/(a -1) < 0, viz. |a| < 1, ℜ(f ) belongs to the quarter-wave modes family. In any case, the real part of the frequency does not depend on a and its imaginary part does not depend on the mode index m. The numerical results are in full good agreement with these analytical findings as illustrated in Figs. 2 and3. Orthogonality of the modes: In general, the eigenmodes of the acoustic problem with finite complex valued impedance are not orthogonal (see appendix B for more details). This is illustrated in table 2 for the first four eigenmodes corresponding to Z = 0 and Z = i. Recall that the inner product of two eigenmodes pn (x) and pm (x) is defined as:

< pn , pm >= Ω pn (x)p * m (x)dx, ( 52
)
and is zero if and only if pn (x) and pm (x) are orthogonal. It follows from Appendix B that < pn , pm >= 0 when Z = 0 as soon as m = n. The computed products are indeed very small as shown in table 2. In the case Z = i, the inner products are not expected to be zero and the following expression can be derived:

| < pm , pn > | = 1 2(m + n) -1) × 2(4n -1) (4n -1)π -2 × 2(4m -1) (4m -1)π -2 . ( 53
)
Table 2 shows that the agreement with computed inner products is again very good. Note that large values (e.g. < p1 , p2 >≃ 0.15) are obtained. Assuming that the acoustic modes form an orthogonal functional basis is thus not justified and can lead to large errors in the case of finite valued boundary impedances.

Case < p1 , p2 > < p1 , p3 > < p1 , p4 > < p2 , p3 > < p2 , p4 > < p3 , p4 > Z = 0 2 × 10 -5 7 × 10 -4 2 × 10 -5 4 × 10 -5 10 -3 6 × 10 -5 Z = 0 -Exact 0 0 0 0 0 0 Z = i 0.

2D modes

This section proposes a comparison between numerical and analytical results in the case of 2D modes. An homogeneous Neumann condition (|Z| = +∞) is now imposed at the bottom, left and right boundaries while a complex impedance (Z = a + ib) is imposed at the top edge of the computational domain (see Fig. 4). Seeking for 2D harmonic solutions, viz. p(x, y, t) = ℜ [p(x, y)exp(-iωt)], the problem can be formulated as follows : The relevant dispersion relation is then :

∂ 2 p ∂x 2 + ∂ 2 p ∂y 2 +
e 2jkyh (k y - k Z ) -(k y + k Z ) = 0, k y = k 2 - mπ L 2 , m ∈ N. ( 55
)
In the case of a purely reactive impedance, viz. Z = ib with b ∈ R, there is no acoustic flux at the boundaries and consistently, the wave number k = ω/c 0 given by Eq. 55 is real (see Fig. 5). On the other, when Z is a purely resistive impedance, viz. Z = a with a ∈ R, the acoustic flux at the top boundary depends on the mode index and both the real and imaginary parts of the eigenfrequencies depend on a (see Fig. 6). In any case, the comparison between the analytical (Eq. 55) and numerical eigenfrequencies shows very good agreement.

B. Active flame calculation in an academic test case

This section presents a calculation where the influence of the flame considered as an acoustic element is accounted for. The configuration is 1D and consists of a duct of length L = 0.5 m and constant cross section where the fresh gas are separated from the hot gas by an infinitely thin flame located at the middle of the duct (Fig. 7). Modeling the unsteady effects of the c =694.36 m/s T =1200 K flame thanks to the classical n -τ model [START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. Part I[END_REF][START_REF] Crocco | Aspects of combustion instability in liquid propellant rocket motors. Part II[END_REF] (Eq. 9) and assuming that the hot-to-fresh gas temperature ratio is T 2 /T 1 = 4 (which means that the speed of sound in the hot gas c 2 is twice as large as in the fresh gas c 1 while the density ratio is ρ 2 /ρ 1 = 0.25), the dispersion relation takes a rather simple form provided Dirichlet or Neumann boundary conditions are used at x = 0 and x = L. [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF] For example, in the case where the acoustic velocity is zero at the inlet (x = 0) and the pressure is fixed at the outlet (x = L), the eigenfrequencies of the problem are solutions of the following dispersion equation:

cos L 4c 1 ω cos 2 L 4c 1 ω - 1 4 Γ(1 + ne iωτ ) -1 Γ(1 + ne iωτ ) + 1 - 3 4 = 0, ( 56
)
where the dimensionless coefficient Γ is the acoustic impedance ratio

Γ = ρ 2 c 2 ρ 1 c 1 . ( 57
)
Note that a more general solution of this problem is given in Kaufmann et al. [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF] where a jump in the duct cross section is allowed at the flame position. The solutions of Eq. 56 are also discussed in details in this reference, especially in the limit case n -→ 0. In the more general case where the interaction index is not small, Eq. 56 can be numerically inverted for any values of the unsteady flame model parameters n and τ . In this section, these exact eigenfrequencies are compared to the numerical results obtained by the asymptotic method and the iterative algorithm described in sections IV-A and IV-B respectively. Although analytically tractable, this test case is very demanding for any finite element based method. Indeed, discrepancies between the numerical and analytical results can be expected for two main reasons: [START_REF] Benoit | Numerical assessment of thermo-acoustic instabilities in gas turbines[END_REF] • in the simulations, the flame front is not infinitely thin as in the analytical model; it is at least as thick as the typical cell size in the interface region. This is a classical difficulty when Dirac delta functions are discretized,

• the reference point needed in the flame model Eq. 9 must belong to the fresh gas side to be consistent with the philosophy of the n -τ model. Since the flame thickness is not zero, the fresh gas region does not extend to the middle of the duct and the reference point cannot be strictly located at x = L/2 = 0.25 m.

Two 2D unstructured meshes have been considered, one coarse mesh (M1) and one fine mesh (M2). To keep the expected errors as small as possible, both meshes are refined near the middle of the duct, as shown in Fig. 8 for mesh M1. Using the finer mesh M2 should reduce the analytical/numerical discrepancies because a) the classical numerical errors stemming from the spatial gradient approximations are reduced and b) the flame thickness can be made smaller and the reference position can be chosen closer to its theoretical value L/2. The main characteristics of meshes M1 and M2 are gathered in table 3 where the minimal flame thickness and maximal reference position allowed by each mesh are also given. Regarding the flame characteristics, two different values of the interaction index have been considered, viz. n = 0.01 and n = 5, corresponding to a small and large amplitude of the flame response respectively. Regarding the time delay, a value leading to both stable and unstable modes among the first few eigenfrequencies has been selected, namely τ = 10 -4 s. For each value of n, a field of local interaction index n u (x) has been constructed so that the very surprising since the 'small' parameter of the problem, Eq. 38, is indeed larger than unity and the developments given in section IV-A are not justified. A better agreement is obtained when the iterative algorithm of section IV-B is used instead. The error increases with frequency, which is classical for finite element based methods, and is larger for the imaginary part of the eigenfrequencies. For example, the error on the real part of the third mode is close to 3.3 % while it is more than 22 % for the imaginary part. In order to investigate the sources of the error, the finer mesh M2 was used instead of M1. Relative errors on the frequency f of the third mode are given in Fig. 10 for different calculations based on the iterative algorithm. R1 and R2 corresponds to calculations using respectively M1 and M2 with the corresponding flame front description and reference point definition (see table 3). R21 corresponds to a simulation where M2 was used together with the flame front and reference point definition allowed by M1. Thus comparison between R1 and R21 gives indications about the numerical errors associated with the finite element approximation while comparing R21 and R2 gives indications about the effect of not considering an infinitely thin flame. In any case, the error on the real part ℜ(f ) of the frequency is very small and virtually not affected by the level of approximation. The error on the imaginary part ℑ(f ) remains significant in the case R21, although smaller than in the case R1, indicating that the numerical errors are significant when M1 is used. Using the fine mesh and associated physical parameters description, run R2, the error is again reduced significantly, indicating that the finite flame thickness effects are not negligible. Note that the algorithm requires only 1 iteration to provide a good approximation of ℜ(f ) with less than 1 % error while 2 iterations are required to reach the same accuracy on ℑ(f ).

C. Steady flame calculations in an annular test rig

The 3D acoustic solver is now used for characterizing the acoustic modes of an industrial demonstrator developed by the SNECMA MOTEURS company (SAFRAN Group). More precisely, an annular combustor consisting of a plenum connected to a combustion chamber with 18 swirled burners is considered. The geometry also contains 216 dilution holes evenly dispatched around the annular combustion chamber. The 3D unstructured mesh contains roughly 3.8 × 10 5 nodes and 10 6 tetraedra (see Fig. 11). The fields of density and speed of sound are inputs of the acoustic solver and they were obtained by performing a Reynolds-Averaged Navier-Stokes calculation over a 360/18 ≃ 20 degrees sector, the result of which being duplicated to produce the requested 360 degrees fields. Figure 12 illustrates the 3D nature of the data provided to the acoustic solver and the necessity to compute them by using an appropriate 3D CFD tool. In the same way, accounting for the acoustic/flame coupling would require complete fields of interaction index and time delay as discussed in section II-C. This would essentially require performing a LES of a 20 degrees sector of the annular combustor and post-processing the results in order to optimize n u (x) and τ (x). [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] Such simulation has not been performed yet on the geometry of interest and only passive flame acoustic calculations will be presented in this section. Note also that the adiabatic coefficient γ was considered constant in the following calculations. Two runs have been performed in order to illustrate the importance of accounting for complex boundary impedance when dealing with annular combustors.

• Run A: An homogeneous Neumann condition is used for the solid walls and inlets (entrance of the plenum and the dilution holes) while the classical zero pressure fluctuations condition is imposed at the outlet of the combustion chamber,

• Run B: The same conditions as in Run A are used for the solid boundaries and the inlets but for the outlet of the combustion chamber a frequency dependent impedance is now prescribed in order to account for the sonic surface that is present downstream the combustion chamber when the combustor is connected to a turbine. The frequency dependent impedance was obtained from classical acoustic theory 51 for nozzles assuming that the first compressor blades row is choked. Real and imaginary part of the corresponding impedance are shown in figure 13. are fully reflecting in this case), the imaginary part of the corresponding eigenfrequencies is always zero. For Run B, the outlet impedance mimics the presence of a chocked nozzle downstream the combustion chamber and part of the acoustic energy is leaving the computational domain through this boundary. Accordingly, all the modes are damped in Run B.

In the case where the classical outlet condition is used (Run A), all the modes are essentially longitudinal in the combustion chamber. This result is in contradiction with the fact that turning modes are often observed in annular combustors. [START_REF] Evesque | Spinning and azimuthally standing acoustic modes in annular combustors[END_REF][START_REF] Krebs | Thermoacoustic analysis of annular combustor[END_REF][START_REF] Schuermans | Non-linear Combustion Instabilities in Annular Gas-Turbine Combustors[END_REF] It is actually an artifact of the classical outlet condition which imposes zero pressure fluctuations and prevents modes with phase differences within the outlet plane from existing. This is confirmed by the results of Run B where three families of modes were found :

• orthoradial modes where the acoustic pressure field in the plenum and combustion chamber is virtually axi-symmetric (not exactly axi-symmetric because the geometry is truly 3D because of the injectors and the dilution holes) -see Fig. 14, modes a and b,

• longitudinal modes where the acoustic pressure field depends mostly on the axial distance along the 'symmetry' axis -see Fig. 14, mode c,

• coupled modes where the pressure field in the plenum and combustion chamber depends both on the axial distance and the azimuthal angle -see Fig. 14, mode d, These modes structures are not obtained if the outlet impedance is set to zero, as it is usual for outlets in acoustic simulations (Run A). Moreover, the lowest frequency mode, close to 350 Hz, is not captured in the Z = 0 case. This illustrates the importance of prescribing proper acoustic boundary conditions for predicting the eigenmodes of annular combustors.

D. 3D Active flame calculations

Three dimensional acoustic calculations based on the AVSP code and including complex impedance and active flame were discussed recently by Martin et al. [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] These authors con- sidered a swirled staged combustor which proved to be unstable both experimentally and numerically by LES. The acoustic calculations showed that :

• when no coupling with the flame is considered and Neumann/Dirichlet boundary conditions are used at the inlet/outlet, three acoustic modes (1L, 2L, 3L) are found below 1000 Hz, with frequencies close to 312, 431 and 841 Hz respectively,

• when complex valued impedance is prescribed at the outlet, all the acoustic modes are damped because of the acoustic losses at the boundary. The decay rates of the first three modes are close to -48, -32 and -4.5 rad/s respectively,

• when a complex valued impedance is used and the flame is made active, the real frequencies of the first three modes are only slightly shifted. Moreover, the 1L mode is heavily damped with a decay rate of -588 rad/s, 3L becomes virtually marginal and 2L becomes unstable with a growing rate close to 600 rad/s,

• these results compare favorably with the LES calculations performed under the same conditions and showing that the turbulent flame oscillates at approximately 380 Hz, the shape of the LES pressure fluctuations along the combustor axis being very close to the shape of the 2L unstable acoustic mode. Moreover, the 2L mode growing rates from the acoustic solver and the LES code are in fair agreement.

The study of Martin et al. [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] demonstrated the ability of the method presented in section IV-B to detect the existence of an unstable mode in the case of an unstable configuration. It is also important that the method predicts stability in the case of a stable configuration, even if the flame-acoustic coupling is accounted for thanks to the local n -τ model of section II-C. As a complement to the results [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] summarized above, a stable configuration is thus considered in the remaining of this section in order to better illustrate the the acoustic calculation procedure with 3D active flame.

The burner considered is a modified version of a single Siemens hybrid burner operated at 17 Mpa. Fully premixed methane / air is injected through two coaxial swirlers (diagonal and axial) at an equivalence ratio of 0.42 for the axial and 0.53 for the diagonal swirler. The flow rates injected in the combustion chamber are respectively 1.7 kg/s for the axial and 17.7 kg/s for the diagonal swirler. More details about the geometry of the burner can be found in previous references dedicated to the structure of the corresponding turbulent steady flame, [START_REF] Selle | Compressible Large-Eddy Simulation of turbulent combustion in complex geometry on unstructured meshes[END_REF] the effect of piloting on the flame, [START_REF] Sengissen | Large eddy simulation of piloting effects on turbulent swirling flames[END_REF] the relationship between hydrodynamic and acoustic fluctuations. [START_REF] Selle | Joint use of Compressible Large-Eddy Simulation and Helmoltz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF] The same burner was also used by Giauque et al. [START_REF] Giauque | System identification of a large-scale swirled partially premixed combustor using les and measurements[END_REF] to assess flame transfer functions from LES, including local interaction index and time delay maps. Note also that in the following the Siemens burner is not mounted on a cylindrical or square combustion chamber as in the previous studies; [START_REF] Selle | Compressible Large-Eddy Simulation of turbulent combustion in complex geometry on unstructured meshes[END_REF][START_REF] Selle | Joint use of Compressible Large-Eddy Simulation and Helmoltz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF][START_REF] Sengissen | Large eddy simulation of piloting effects on turbulent swirling flames[END_REF] a geometry more representative of the actual combustion chamber is used instead. To the author's knowledge, no instability were found experimentally under the considered operating conditions. The computational domain is shown in Fig. 15. The tetrahedral mesh is refined in the injection region in order to better capture the mean temperature gradient and flame response. As shown in Fig. 15, only a 360/24 = 15 degrees sector of the full configuration is considered. Note however that no major issues are expected in performing the acoustic calculation of the 360 degrees geometry (such calculation has been presented in section V-C for a passive flame case). The only technical difficulty to perform a 360 degrees active flame calculation is in defining one reference position per sector, which does change neither the mathematical nature of the underlying eigenvalue problem nor the related numerical issues. Although no information about the azimuthal modes can be obtained from 15 degrees calculations, such computations are presented below to illustrate the method. As in section V-C, the fields of density and speed of sound are inputs of the acoustic solver and they were obtained by performing a Reynolds-Averaged Navier-Stokes calculation over the computational domain. The flame transfer function was assessed by introducing acoustic perturbations at the diagonal swirler entrance as described in Giauque et al. [START_REF] Giauque | System identification of a large-scale swirled partially premixed combustor using les and measurements[END_REF] The difference with this reference however is that a broad band excitation was used instead of an harmonic signal. The Wiener-Hopf 42 equation was then used in order to obtained the local flame transfer function in the range 100 -600 Hz by post-processing the LES results. Typical maps of interaction index and time delays are given in Fig. 16 while the frequency dependence of these parameters is depicted in Fig. 17. Significant variations of the overall amplitude of the flame response are obtained when the frequency evolves from 100 to 600 Hz.

The frequencies associated with the first five acoustic modes are gathered in table 6. Results from both passive and active flame simulations are included. In both case, Neumann and axi-periodic conditions were used for the liner and lateral boundaries respectively. Complex valued impedances were prescribed at the combustion chamber inlet and outlet. In absence of experimental measurements, values consistent with the partially reflecting boundary conditions used in the LES performed to assess the flame response (shown in Figs. 16 and 17) were chosen, viz. 1/Z = -1 -7540i/ω for the outlet 55 and 1/Z = -1.67 for the inlet. With these values, the boundaries have a low-pass filter behavior, 55 the cut-off frequency being close to 600 Hz. Note also that the boundary impedances are of the form Eq. 34 so that all the passive flame modes have been obtained by solving the quadratic eigenvalue problem Eq. 35 only once. For each mode, the iterative procedure of section IV-B has then been used to obtain the active flame modes. In general, two iterations proved enough to reach a fair convergence of the complex frequency except for the first mode for which a third iteration has been necessary. Differently from the case investigated in Martin et al., [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF] the flame-acoustic coupling in the present configuration does not modify the acoustic modes substantially. The real frequency shift is 2 % or less for all modes except for the first one for which the frequency shift is more significant. All the passive flame modes are damped because of the acoustic radiation through the inlet and outlet boundaries. Note that the acoustic losses are less significant for the first radial mode (the 5th mode in 6). The active flame calculations show that the flame acoustic coupling has either a stabilizing effect (modes 1 and 3) or a destabilizing effect (modes 2, 4 and 5). Even if this effect is sometimes significant (see mode 4), it is always compensated by the damping effect of the boundaries so that the growing rate keeps a negative value for each mode, consistently with the overall stable behavior of the test rig.

VI. Conclusion

A finite element based methodology has been developed in order to determine the thermoacoustic modes of industrial systems like annular combustors. The method allows to handle frequency dependent complex valued impedance at the boundaries as well as 3D, distributed unsteady flames. Two ways of accounting for the flame/acoustic coupling have been developed and proved to give similar results in the limit case of small flame response amplitude. Several academic cases have been considered in order to validate the procedure and its ability to handle complex geometries has been illustrated by considering an annular combustor with 18 swirled injectors. A 15 degrees sector of another annular combustor has also be considered as an example of acoustic calculation including acoustic/flame coupling. Coupled to RANS/LES codes to provide the required input data, the numerical approach described in this paper can be used to predict the thermoacoustic modes of industrial systems and associated stability. consistent with the zero Mach number assumption can be obtained afterward.

For example, the linearization of the continuity equation, Eq. 1, gives:

∂ρ 1 ∂t + u 1 • ∇ρ 0 + ρ 0 ∇ • u 1 + u 0 • ∇ρ 1 + ρ 1 ∇ • u 0 = 0, (59)
Note that the first three terms do not depend on the mean velocity field and were already present in Eq. The two other terms of Eq. 59, namely u 0 • ∇ρ 1 and ρ 1 ∇ • u 0 , involve the mean velocity field and can be neglected if they are small compared to the time derivative, viz. small compared to ǫρ 0 c 0 /L a . Their order of magnitude being ǫρ 0 M c 0 /L a and ǫρ 0 M c 0 /L f respectively, it follows that the zero mean flow assumption is justified as long as the Mach number is such that M ≪ 1 and M ≪ L f /L a .

The same methodology can be applied to Eqs. 2 and 4. Linearizing Eq. 2 without assuming that the mean flow is at rest, one obtains Eq. 6 enriched by the terms ρ 0 u 1 • ∇u 0 , ρ 0 u 0 • ∇u 1 and ρ 1 u 0 • ∇u 0 at the LHS. By comparison with the time derivative term in Eq. 6, these terms are negligible if M ≪ L f /L a , M ≪ 1 and M ≪ L f /L a respectively. Regarding Eq. 4, the only new term when the mean flow is not at rest is u 0 •∇s 1 which can be neglected compared to ∂s 1 /∂t when M ≪ 1. Eventually, Table 7 gathers all the constraints that must be matched on the mean Mach number for the zero mean flow assumption to hold. For low frequency modes of typical combustors, the flame brush thickness is smaller than Equation Constraint mass M ≪ 1 and M ≪ L f /L a momentum M ≪ L f /L a , M ≪ 1 and M ≪ L f /L a entropy M ≪ 1

Table 7. Summary of the constraints on the mean Mach number implied by the zero mean flow assumption.

the acoustic wavelength and M ≪ L f /L a is most likely the most restrictive constraint on M .

By definition, the eigenvectors p1 (x) and p2 (x) are orthogonal if their 1/γ-weighted inner product is zero, viz. < p1 , p2 > γ = 0. Note that this definition of the orthogonality of the eigenvectors reduces to the more classical one, viz. < p1 , p2 >= 0, when γ is constant over the flow domain (see Eq. 52).

Taking the complex conjugate of Eq. 15 applied to (p 2 (x), ω 2 ) one obtains p *

2 γ = - p 0 ω * 2 2 ∇ • 1 ρ 0 ∇p * 2 -i γ -1 γω * 2 q * , (61)
and the < p1 , p2 > γ inner product is equal to :

< p1 , p2 > γ = - Ω p 0 ω * 2 2 p1 ∇ • 1 ρ 0 ∇p * 2 dx -i Ω γ -1 γω * 2 p1 q * dx. ( 62
)
The second integral strongly depends on the unsteady heat release. In the case of a complex turbulent flame and given the current knowledge and understanding of flame acoustic coupling, there is no reason to believe that this term should be zero or that it can be compensated by the first integral which does not depend on combustion. Thus in the case of an unsteady flame, the thermoacoustic modes are most likely non-orthogonal.

Let us now consider the more interesting passive flame case for which the orthogonality of the modes is sometimes assumed. [START_REF] Culick | Combustion instabilities in liquid-fueled propulsion systems-An overview[END_REF] Integrating the first integral of Eq. 62 by parts leads to < p1 , p2 > γ = p 0 ω * 2

2 Ω

1 ρ 0 ∇p 1 • ∇p * 2 dx - p 0 ω * 2 2 ∂Ω 1 ρ 0 p1 ∇p * 2 • n BC dσ, (63)
Integrating again by parts and using Eq. 15 applied to (p 1 , ω 1 ) lead to:

1 - ω 2 1 ω * 2 2 < p1 , p2 > γ = p 0 ω * 2 2 ∂Ω Z 1 ρ 0 (p * 2 ∇p 1 -p1 ∇p * 2 ) • n BC dσ, (64)
Note that the boundary integral is taken over ∂Ω Z only because neither ∂Ω D nor ∂Ω N contribute given that either p1 = p * 2 = 0 or ∇p 1 • n BC = ∇p * 2 • n BC = 0 on these boundaries. At this point, it is useful to consider the two distinct cases: Case 1: ∂Ω Z is empty: In this case, only Neumann or Dirichlet conditions are prescribed and the angular frequencies are real valued (ω * 2 = ω 2 ) since there is no acoustic radiation through the boundaries. Eq. 64 reduces to

1 - ω 2 1 ω 2 2 < p1 , p2 > γ = 0, (65)
which shows that p1 and p2 are orthogonal as soon as ω 1 = ω 2 . Since one can always find orthogonal eigenvectors sharing the same eigenfrequency by using the Gram-Schmidt orthogonalisation procedure, the acoustic problem admits a set of orthogonal eigenvectors when ∂Ω Z is empty.

Case 2: ∂Ω Z is not empty: Excluding the particular case where ω 1 = ω * 2 , the eigenvectors p1 and p2 are orthogonal only if the boundary integral in Eq. 64 turns out to be zero. In other words, the eigenvectors of the acoustic problem with finite impedance are non orthogonal in the general case. A sufficient condition for zeroing the RHS of Eq. 64 is obviously:

p * 2 ∇p 1 • n BC = p1 ∇p * 2 • n BC , on ∂Ω Z , (66)
Making use of Eq. 24, this constraint can be re-written as

ω 1 Z(ω 1 ) = - ω * 2 Z * (ω 2 ) , (67)
The following analytic form for Z(ω) ensures that Eq. 67 is satisfied for any values of ω 1 and ω 2 :

Z(ω) = iR 0 ω, R 0 ∈ R, (68)
Thus a sufficient condition for the acoustic problem to admit orthogonal eigenvectors is that the boundary impedance, if finite, is given as in Eq. 68. Note that a particular case is the Dirichlet condition Z = 0.

1 k

 1 Square sparse matrices of size N which do not depend on ω A, B, C, D Square sparse matrices of size N which do not depend on ω B, R Square sparse matrices of size N which depend on ω n ref Reference unit vector n BC Outward unit vector normal to the boundary P Complex valued column vector of size N . u Velocity vector x Position vector a, b Real and imaginary parts of the reduced impedance Z. c Speed of sound C p Heat capacity per unit mass at fixed pressure C v Heat capacity per unit mass at fixed volume f Frequency i Imaginary complex number √ -Wave number ω/c 0 . k y Wave number associated to the y direction, k y = k 2 -(mπ/4L) 2 , m ∈ N. L, h Length and height of the 2D computational domain. L a Characteristic length scale of the acoustic perturbations L f Characteristic length scale of the turbulent flame M Mach number m Mode index. N Size of the non-linear eigen value problem, N = N v -N d n Interaction index N d Number of vertices on ∂Ω D N v Number of vertices in the unstructured mesh describing Ω n u Field of interaction index p Static pressure q Heat release per unit volume q tot Volume averaged heat release. Scaling purpose r Heat capacities difference C p -C v s Entropy per mass unit S ref Area of the burner mouth at the reference position t Time U bulk Bulk velocity of the burner. Scaling purpose x f Infinitely thin premixed flame position Z Complex impedance Z 0 , Z 1 , Z 2 Complex valued constants S v Set of vertices of the mesh which do not belong to ∂Ω D S Ω Set of all the tetrahedral elements in the mesh Subscripts 0 Time averaged variable or passive flame mode (asymptotic method) 1 Fluctuating quantity or first order term in the asymptotic method k, j Indices BC Boundary value ref Reference for the flame model Conventions φ Complex amplitude of the fluctuating quantity φ Lp Linear mapping from C to C Lu Linear mapping from C 3 to C ℑ(z) Imaginary part of the complex valued quantity z ∇f Spatial gradient of f ℜ(z) Real part of the complex valued quantity z Symbols δ f Effective 1D premixed flame thickness. ǫ Small parameter of order p 1 /p 0 γ Adiabatic coefficient C p /C v Ω Computational domain ω Angular frequency ∂Ω Boundary of the computational domain ∂Ω D Subset of the boundary where p = 0 ∂Ω N Subset of the boundary where ∇p • n BC = 0 ∂Ω Z Subset of the boundary where the reduced impedance Z is imposed φ j Classical linear by parts form function for all node v j belonging to S v ρ Mixture density τ Time delay τ u Field of time delay

Figure 1 .

 1 Figure 1. Schematic of the computational domain for computing the longitudinal modes of the isothermal 2D domain.

Figure 2 .

 2 Figure 2. Influence of purely reactive impedance Z = ib on the 1D modes eigenfrequencies. Left: real part of the frequency. Right: Imaginary part. • numerical results, -analytical solution. Abscissa represent b.

Figure 3 .

 3 Figure 3. Influence of purely resistive impedance Z = a on the 1D modes eigenfrequencies. Left: real part of the frequency. Right: Imaginary part. • numerical results, -analytical solution. Abscissa represent a.

Figure 4 .

 4 Figure 4. Schematic of the computational domain for computing the 2D modes of the isothermal 2D domain.

  p = 0 for y = h.

Figure 5 .

 5 Figure 5. Influence of purely reactive impedance Z = ib on the 2D modes eigenfrequencies. Left: real part of the frequency. Right: Imaginary part. • numerical results, -analytical solution. Abscissa represent b.

Figure 6 .

 6 Figure 6. Influence of purely resistive impedance Z = a on the 2D modes eigenfrequencies. Left: real part of the frequency. Right: Imaginary part. • numerical results, -analytical solution. Abscissa represent a.

Figure 7 .

 7 Figure 7. 1D model problem for the acoustic/flame coupling

  2498 ≤ x ≤ 0.2502 0.2495

Figure 8 .

 8 Figure 8. 2D unstructured mesh M1 for the 1D active flame calculations.

Figure 10 .

 10 Figure 10. Relative errors (%) on the frequency of the 3 rd mode from runs R1 ( ), R21 ( ) and R2 ( ) versus the iteration number. Left: Real part, Right: Imaginary part of the frequency. The corresponding exact values are ℜ(f ) = 1227.3 and ℑ(f ) = 41.6. Iteration 0 corresponds to passive flame calculations.

Figure 11 .

 11 Figure 11. 3D unstructured mesh of the annular combustor.

Figure 12 .

 12 Figure 12. Temperature field in the combustor. Scale is from 450 (black) to 2500 K (white).

Figure 13 .

 13 Figure 13. Frequency dependent reduced boundary impedance imposed at the combustion chamber exit for Run B. Left: Real part, Right: Imaginary part

a:Figure 14 .

 14 Figure 14. Modulus of the pressure amplitude for the first four modes of the annular combustor -Run B. Scale is from 0 (white) to 1 (black).

Figure 15 .

 15 Figure 15. 3D unstructured mesh for the acoustic analysis of the Siemens combustor.

Figure 16 .

 16 Figure 16. Typical maps of the interaction index (left column) and time delay (right column) at 300 Hz (top row) and 500 Hz (bottom row). The time delay is only plotted in regions where the amplitude of the flame response (interaction index) is significant.

Figure 17 .

 17 Figure 17. Volume averaged interaction index (left) and time delay (right) versus frequency.

  Frequency and decay rate of the first 5 eigenmodes of the Siemens combustor for both the passive and active flame calculations. The shape of the passive flame modes are shown in the right column. Modulus of the pressure amplitude is shown with scale from 0 (white) to 1 (black). Active flame modes have very similar shape.

Table 2

 2

		15	0.10	0.08	0.076	0.056	0.051
	Z = i -Exact	0.151	0.106	0.081	0.076	0.062	0.052

. Inner products from the first four modes for Z = 0 and Z = i. The eigenmodes are normalized by imposing < pm , pm >= 1. Indices correspond to the mode number in Eq. 50.

Table 3 .

 3 Characteristics of meshes M1 and M2 and associated best physical values.

Table 4 .

 4 Frequency and decay rate of the first 4 eigenmodes of the 3D annular combustor.

	Run A

Table 5 .

 5 Frequency and decay rate of the first 4 eigenmodes of the 3D annular combustor. Run B

Table 6 .

 6

Acknowledgments

The authors acknowledge SNECMA for its support and the French national computing center CINES for giving access to parallel computers.

models given in Eqs. 9 and 11 correspond to the same flame response amplitude. A simple way to proceed is to assume that n u (x) is constant by parts, equal to η in the flame region and 0 anywhere else. Then, integrating Eq. 11 over space and equating the result to Eq. 9 gives the following connection between the classical and local interaction index:

where δ f is the thickness of the effective flame region in the computation and x f is the mean position of the flame, viz.

When the amplitude of the flame response is small, figure 9 (left) shows that both the asymptotic method and the iterative algorithm provide results in good agreement with the analytical solution, even if the coarser mesh M1 is used. With the selected time delay, the flame acoustic coupling makes the first and fourth modes stable while the second one is not modified (this mode admits a velocity node at the middle of the duct which zeroes the coupling term [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF] ) and the third mode is amplified by the presence of the flame. The passive flame calculation leads to marginal modes (consistently with the purely reflecting boundary conditions) in good agreement with the theoretical real eigenfreqencies. The stabilizing/destabilizing effect of the unsteady flame is well captured by both numerical methods, showing that the linearization process proposed in section IV-A is valid in this case. When the amplitude of the flame response is large (n = 5), figure 9 (right) shows that only the third mode is amplified by the flame/acoustic coupling, the first and the forth being damped and the second one being left unchanged. The passive flame calculation still produces marginal modes but the amplitude of the flame response being large, the associated real eigenfreqencies do not match the theoretical values very well. Large differences between the asymptotic method and the analytical solution also appear in this case. This is not

Appendix A: Validity of the zero Mach number assumption

To quantify the validity of the zero Mach number assumption, we first note that the time and length scales for any fluctuation (density, velocity or entropy) are L a /c 0 and L a respectively where L a is the acoustic wavelength of the perturbation and c 0 the speed of sound. For example, the ∂ρ 1 /∂t in Eq. 5 is of order ǫρ 0 c 0 /L a . Besides, we note that the steady part of the fields of density, velocity or entropy can evolve significantly over the flame brush thickness L f . Considering this quantity as the relevant length scale for assessing gradients of time averaged fields, the order of magnitude of the terms appearing in the linear counterpart of Eqs. 1, 2 and 4 can be assessed. Of course, the linearization is performed without assuming u 0 ≃ 0 so that constraints on the mean Mach number M = √ u 0 • u 0 /c 0 Appendix B: Orthogonality of the eigenmodes Consider p1 (x) and p2 (x) two eigenmodes of Eq. 15 associated to angular frequencies ω 1 and ω 2 . The objective of this appendix is to precise the conditions under which these two modes are orthogonal. We first extend the classical inner product in order to account for the spatial variations of γ, defining the following 1/γ-weighted inner product: