Camilo F Silva
email: kamilosilva@gmail.com.

Ignacio Duran
email: duran@cerfacs.fr.

Franck Nicoud
email: franck.nicoud@univ-montp2.fr

Stéphane Moreau
email: stephane.moreau@usherbrooke.ca.

Ignacio Durán

Boundary conditions for the computation of thermo-acoustic modes in combustion chambers

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. Introduction

The prediction of combustion instabilities in combustion chambers embedded in complex systems, such as aeronautical engines, is extremely difficult because different phenomena must be taken into account at once: the strong interaction between the flame and the most representative turbulent length scales of the reacting and burnt gases, the two-way interaction between the flame and the radiated acoustic waves it produces, the possible coupling between entropy waves caused by the hot spots generated by the flame and the acoustic waves produced by the non-homogeneities in the mean flow, and evidently the role played by all boundary conditions (inlets, outlets and walls).

Different strategies have been developed over the years [START_REF] Rayleigh | The explanation of certain acoustic phenomena[END_REF][START_REF] Poinsot | Interactions between acoustics and combustion[END_REF][START_REF]Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms and modeling[END_REF][START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF][START_REF] Rao | Use of finite element methods in frequency domain aeroacoustics[END_REF] to understand the physics of all these interactions and to create methodologies to control or avoid them. In recent years, reacting and compressible Large Eddy-Simulations (LES) have shown their capability to study the dynamics of turbulent flames [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] 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][START_REF] Roux | Studies of mean and unsteady flow in a swirled combustor using experiments, acoustic analysis and large-eddy simulations[END_REF][START_REF] Pitsch | Large Eddy Simulation of turbulent combustion[END_REF][START_REF] Schmitt | Large Eddy Simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high pressure burner[END_REF]. Due to its intrinsic nature (resolution of the unsteady 3D Navier-Stokes equations) LES is able to predict the interactions between the flame, the turbulent flow and the acoustic modes of the chamber, and for some particular configurations to distinguish between stable and unstable combustion systems [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF]. Yet, LES presently remains very CPU demanding and its use for parametric studies of aeronautical engines design is still excluded.

For this reason, simplifications on the geometry of the combustor as well as on the physics of the reactive turbulent flow may lead to interesting and useful strategies to study combustion instabilities. Following this philosophy, a wide spread methodology has been developed over the years which is based on the so called low-order models. A given system is characterized as a network of individual one dimensional/ two dimensional axisymmetric acoustic elements in which only perturbations of first order are kept [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]. Since acoustics is considered to be the only cause of these perturbations, all viscous terms in the Navier Stokes equations are neglected. The acoustic elements are then described by 1D Linearized Euler Equations (LEE) and they are associated with each other by relations of conservation of mass, momentum and energy. These methods aim to study the low frequency unstable modes that are originated by the coupling of plane acoustic waves with compact heat release sources. For complex geometries, as actual aeronautical combustion chambers, this 1D modeling strategy might remain too simplistic: the full three dimensional geometry might play an important role not only on the values of the predicted resonance frequencies but also on the growth rate of the acoustic modes, as shown for example in [START_REF] Silva | Combining a helmholtz solver with the flame describing function to assess combustion instability in a premixed swirled combustor[END_REF]). Instead of solving the LEE equations for a set of interconnected acoustic elements, it might be more precise to account for the actual 3D geometry. The development of LEE solvers for complex geometries associated to combustion chambers is still in a nascent state. In such an approach the stability not only of acoustic modes, but also vortical and entropy modes is studied [START_REF] Gikadi | Effects of the mean flow field on the thermo-acoustic stability of aero-engine combustion chambers[END_REF]. Nevertheless, the mean flow velocity field in typical combustion chambers is very low, and assuming that convective terms do not significantly contribute in the analysis of combustion instabilities is reasonable. Therefore, under the low Mach approximation, the LEE equations can be recast into a wave equation that if expressed in the frequency domain results in the so called Helmhotz Equation. The respective numerical tools, commonly known as Helmholtz solvers, have proved their ability to predict combustion instabilities [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF][START_REF] Selle | Joint use of compressible Large-Eddy Simulation and Helmholtz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF][START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF] in real combustion chambers.

The ability of the strategies above described (LES, low-order model techniques and Helmholtz solvers) to correctly estimate combustion instabilities depends strongly on the acoustic description at the system boundaries. The main objective when describing the acoustic behaviour of such boundaries is to be able to predict accurately the transmission and reflection of acoustic waves caused by impinging acoustic waves at these frontiers. The general idea is to solve the quasi-1D LEE in a domain that represents the regions either upstream (fan + compressor stages) or downstream (turbine stages + nozzle) of the combustion chamber. Several works have made important contributions in the study of acoustics of quasi 1D elements. The first significant study is due to Tsien [START_REF] Tsien | The transfer functions of rocket nozzles[END_REF]. In this work, the propagation of acoustic waves through a quasi 1D system is studied. The mean flow is considered to be linearly distributed along the axis. Subsequently, Marble and Candel [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] proposed analytical expressions that estimate the transmission and reflection of acoustic waves due to either acoustic or entropy waves that are injected in a quasi 1D system. The main assumption, in addition to consider plane acoustic waves, is to consider the system as compact, i.e. as a system whose length is considered small compared with a characteristic acoustic wavelength. These relations have been demonstrated to be very useful in the low frequency limit. In the work of Stow et al. [START_REF] Stow | Reflection of circumferential modes in a choked nozzle[END_REF], this theory is extended to second-order through an asymptotic analysis. A second-order correction is added through an effective length which is function of the mean flow spatial distribution. Consequently, the transmission and reflection coefficients are well predicted within a low frequency band. In the study of Moase et al. [START_REF] Moase | The forced response of choked nozzles and supersonic diffusers[END_REF], transmission and reflection coefficients are derived for a quasi-1D system whose mean flow is modelled to be piece-wise linear. As a result, an acoustic characterization of any quasi 1D geometry can be performed for any frequency band of interest as long as the mean flow can be described by piece-wise linear functions. Recently, a remarkable work has been carried out by Duran and Moreau [START_REF] Duran | Solution of the quasi-one-dimensional linearized euler equations using flow invariants and the magnus expansion[END_REF] in which the quasi 1D LEE are solved analytically through an asymptotic analysis based on the Magnus expansion of the quasi 1D LEE. In addition to analytical methods, numerical approaches to solve the quasi 1D LEE system can be also used to study the acoustic properties of a given quasi 1D system. Corresponding numerical schemes have demonstrated to be useful when evaluating the acoustic impedance of isentropic nozzle and diffusers to subsequently study combustion instability in complex 3D combustion chambers, as in the work of Lamarque and Poinsot [START_REF] Lamarque | Boundary conditions fo acoustic eigenmode computations in gas turbine combustion chambers[END_REF].

In the present study, a numerical approach is used to solve the quasi 1D LEE in the case where neither the total enthalpy nor the total entropy are constant over the entire domain. In the first section, the quasi 1D LEE system is presented. Compared to the system of equations exposed by [START_REF] Nicoud | About the zero Mach number assumption in the calculation of thermoacoustic instabilities[END_REF], additional terms are added accounting for the compression or expansion stages in the system. Subsequently, analytical solutions are derived to validate the quasi 1D LEE numerical tool in two specific cases: isentropic compact nozzles [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] and a one dimensional compressor. Finally the reflection coefficient at the inlet air circuit of an helicopter combustion chamber is computed and used to evaluate the first two acoustic modes of a 3D combustion chamber.

II. The quasi 1D Linearized Euler Equations

A quasi 1D Euler system of equations in the conservative form is presented (Nicoud and Wieczorek [START_REF] Nicoud | About the zero Mach number assumption in the calculation of thermoacoustic instabilities[END_REF]). This set of equations states that the flow is strictly one-dimensional and therefore only changes over the stream direction x are of relevance. It implies that only plane acoutic waves can be retrieved. The quasi 1D Euler system of equations reads:

∂ ∂t (ρA) + ∂ ∂x (ρuA) = 0 (1)
∂ ∂t (ρuA) + ∂ ∂x (ρu 2 A) + A ∂p ∂x -AF = 0 (2) ∂ ∂t (ρh t A) + ∂ ∂x (ρuh t A) -A ∂p ∂t -A Q -AW = 0 (3)
where ρ stands for the density of the fluid, u the velocity in the x direction, p the thermodynamic pressure and h t the total enthalpy. In addition the variables A, Q, F and W stand for the crosssection area, an external energy source, the volume-specific force added by the compressor or added to the turbine, and the volume-specific work done by the compressor or supplied to the turbine, respectively. In the following the conservation of energy (3), which is function of the total enthalpy h t , is replaced by the entropy (s) equation:

ρT ∂s ∂t + u ∂s ∂x -Q = 0 (4)
All variables in Eqs. ( 1), ( 2) and ( 4) are decomposed in mean and fluctuating quantities (φ(x, t) = φ(x) + φ ′ (x, t)) and only first order terms are kept. Subsequently, harmonic perturbations are considered (φ ′ = φe -jωt ). The quasi 1D Linearized Euler Equations LEE for a calorically perfect gas read:

∂ ū ∂x + ū ∂ ∂x + ū A ∂A ∂x ρ + ∂ ρ ∂x + ρ ∂ ∂x + ρ A ∂A ∂x û -jω ρ = 0 (5) 1 ρ ∂c 2 ∂x + ū ρ ∂ ū ∂x + c2 ρ ∂ ∂x ρ + ∂ ū ∂x + ū ∂ ∂x û + (γ -1) T 1 p ∂ p ∂x + ∂ ∂x ŝ - F ρ -jω û = 0 (6) γr Q ρp ρ + ∂s ∂x û + ū ∂ ∂x + (γ -1) Q p ŝ - r p Q -jωŝ = 0 (7)
where the linearized state equation and entropy expression (Eq. ( 8)) have been considered, i.e.

p p - ρ ρ -T T = 0 and ŝ c p = p γ p - ρ ρ . (8)
This system of equations can be expressed as an algebraic linear system in the form

A         ρ û ŝ        =         0 F ρ r p Q        (9)
where

A =         ∂ ū ∂x + ū ∂ ∂x + ū A ∂A ∂x -jω ∂ ρ ∂x + ρ ∂ ∂x + ρ A ∂A ∂x 0 1 ρ ∂c 2 ∂x + ū ρ ∂ ū ∂x + c2 ρ ∂ ∂x ∂ ū ∂x + ū ∂ ∂x -jω (γ -1) T 1 p ∂ p ∂x + ∂ ∂x γr Q ρ p ∂ s ∂x ū ∂ ∂x + (γ -1) Q p -jω         (10)
Note that this system of linearized equations models propagation of plane acoustic waves through a calorically perfect medium. These expressions are, nevertheless, very suitable when modeling the upstream or downstream region of typical aeronautical combustion chambers. Combustion instabilitiy, a phenomenon highly dependent on boundary conditions, takes place generally at low frequencies, i.e. at frequencies where only plane acoustic waves are present. Moreover, under normal operating conditions, reactions only takes place within the confinement of a combustion chamber.

Accordingly, assuming a constant heat capacity c p either upstream or downstream of the combustion chamber is fairly correct. Note however that the assumption can be easily relaxed as far as Eq. ( 9) is solved numerically.

Within the computational domain, the Quasi-1D LEE (Eqs. 5, 6 and 7) are discretized by finite differences (FD) with a second order centered scheme. At the inlet, a first order FD downwind scheme is applied, whereas at the outlet a first order FD upwind scheme is used. A staggered grid arrangement (velocity fluctuations stored at the cell edges, density and entropy fluctuations stored at the cell centers as shown in Fig. 1.) has been used in order to avoid the pressure field to be contaminated by the classical odd-even decoupling phenomenon [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Wieczorek | Numerical study of Mach number effects on combustion instability[END_REF]. For the (i + 1)/2 cell, equations of continuity, momentum and entropy are discretized as follows:

i i +1 i -2 i +2 i -1 1 2 n -1 n û ρ, ŝ (i + 1)/2 -1 (i + 1)/2 n/2 1 cell cell cell cell ∆x
• continuity equation (node i + 1) ρi+1 ∂ ūi+1 ∂x + ūi+1 (ρ i+3 -ρi-1 ) 4∆x + ρi+1 (û i+2 -ûi ) 2∆x + (û i+2 + ûi ) 2 ∂ ρi+1 ∂x + (û i+2 + ûi ) 2 ρi+1 A i+1 ∂A i+1 ∂x + ρi+1 ūi+1 A i+1 ∂A i+1 ∂x -jω ρi+1 = 0 (11)
• momentum equation (node i)

(ρ i+1 + ρi-1 ) 2ρ i ∂c 2 i ∂x + (ρ i+1 + ρi-1 )ū i 2ρ i ∂ ūi ∂x + c2 i ρi (ρ i+1 -ρi-1 ) 2∆x + ûi ∂ ūi ∂x + ūi,1 (û i+2 -ûi-2 ) 4∆x + (γ -1) Ti (ŝ i+1 + ŝi-1 ) 2p i ∂ pi ∂x + (γ -1) Ti (ŝ i+1 -ŝi-1 ) 2∆x -jω ûi = Fi ρi (12)
• entropy equation(node i + 1)

γr qi+1 ρi+1 pi+1 ρi+1 + (û i+2 + ûi ) 2
∂s i+1 ∂x + ūi+1 (ŝ i+3 -ŝi-1 ) 4∆x + (γ -1) qi+1 pi+1 ŝi+1 - r pi+1 qi+1 ŝi+1 -jωŝ i+1 = Qi+1 (13)
The resulting discretized matrix A for the previous system of equations is a sparse rectangular matrix of size m × m size, in which m = 3n/2 and n/2 is the number of cells of the system. A Jacobi preconditioner is used to improve the conditioning of matrix A and the linear system is solved by inverting the preconditioned matrix through the LU factorization. The resulting numerical tool is termed 'SNozzle'.

III. Definition of waves

In the frequency domain, acoustic and entropy waves ŵ are defined by an amplitude B and a phase φ, so that ŵ = B e jφ . They read [START_REF] Moase | The forced response of choked nozzles and supersonic diffusers[END_REF][START_REF] Leyko | Comparison of direct and indirect combustion noise mechanisms in a model combustor[END_REF]

ŵ+ = B + e jωx/c(1+ M) = p γ p + û c (14)
ŵ-= B -e -jωx/c(1-M) = p γ p - û c ( 15
) ŵS = B S e jωx/c M = ŝ c p (16)
where ŵ+ and ŵ-represent the non-dimensional acoustic waves propagating respectively in the direction of the flow at speed (ū + c), and against the flow at speed (ūc). ŵS stands for the entropy wave that is convected at the mean flow velocity (ū). Inversely, the variables p and û are expressed in terms of acoustic waves as

p γ p = 1 2 ŵ+ + ŵ- and û c = 1 2 ŵ+ -ŵ- (17)
The waves ŵ, that can appear in a 1D system, are shown in Fig. [START_REF] Poinsot | Interactions between acoustics and combustion[END_REF]. Indices 1 and 2 identify the upstream and downstream state respectively. The reflection R and transmission T coefficients are now defined as shown in table 1. The coefficient R AA is the reflection coefficient that measures the amplitude of the outgoing upstream traveling acoustic wave ŵ-1 with respect to an incoming downstream traveling wave ŵ+ 1 when no incoming entropy waves are assumed present ( ŵS 1 = 0). Similarly R SA relates the amplitude of an upstream outgoing traveling acoustic wave ŵ-1 with respect to an incoming downstream traveling entropy wave ŵS 1 . In this case, no incoming acoustic waves are present ( ŵ+ 1 = 0). The transmitted waves, defined by the tranmission coefficient T , are determined in a similar way. Being these acoustic parameters defined for a downstream boundary condition, they are used to characterize systems such as nozzles or turbines. The definition of these parameters should be inverted if an upstream boundary condition is taken as the reference, which is the case for intakes and compressors.

RAA → ŵ- 1 / ŵ+ 1 → with ŵS 1 = 0 RSA → ŵ- 1 / ŵS 1 → with ŵ+ 1 = 0 TAA → ŵ+ 2 / ŵ+ 1 → with ŵS 1 = 0 TSA → ŵ+ 2 / ŵS 1 → with ŵ+ 1 = 0
In the next two sections attention is paid to two of these cases. First, nozzles are treated as isentropic systems in which no enthalpy expansion is present, i.e. no turbines stages are considered.

In this manner, analytical solutions found in the literature for the acoustic response of compact nozzles to incoming acoustic or entropy waves [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] can be compared to the solution of the quasi 1D LEE in the zero frequency limit (ω ≈ 0). Second, the acoustic response of an ideal compact compressor embedded in a constant section duct is evaluated. This simplified case permits to compare the results obtained by the quasi 1D LEE to the solutions given by analytical relations which are indeed a generalization of analytical models such as the one exposed in [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF][START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF].

IV. Response of isentropic nozzles

As already stated in the introductory part of this work, determining appropriate acoustic boundary conditions for combustion chambers is crucial in order to correctly predict combustion instabilities. These acoustic boundary conditions depend, among others, on the mean flow either upstream (compressor stages) or downstream (turbine stages). As a first approximation, it is useful to consider the flow in these regions as isentropic ( Q, Q = 0) and isenthalpic ( F, F = 0). This means that the combustion chamber is only surrounded by perfect diffusers and nozzles.

A. Compact Nozzles

The continuity and energy relations in the framework of the 1D linearized Euler equations for compact systems (see Eqs. ( 62) and (64) in appendix) reduce to:

ρ1 ρ1 + û1 ū1 = ρ2 ρ2 + û2 ū2 (18)
Tt1 Tt1 = Tt2 Tt2 (19)
in which it has been assumed that the fluid has a constant specific heat c p and therefore ĥ = c p T .

Expressing the instantaneous total temperature as

T t = T (1 + [(γ -1)/2]M 2 ), in which M = u/c
represents the instantaneous Mach number, the following relation results

Tt

Tt = T T + (γ -1) M M 1 + γ-1 2 M2 = T T + (γ -1) 1 + γ-1 2 M2 M2 û ū -M2 ĉ c , (20)
where ĉ/c = T /2 T . Replacing Eq. ( 20) in Eq. ( 19) and considering also the relations exposed in Eq. ( 8), the following expression yields:

1 1 + γ-1 2 M1 2 M1 û1 c1 + ρ1 ρ1 + γ (γ -1) ŝ1 c p = 1 1 + γ-1 2 M2 2 M2 û2 c2 + ρ2 ρ2 + γ (γ -1) ŝ2 c p (21)
Finally, substituting the definition all variables into a wave description (Eq. ( 17)) at x = 0 into Eqs. ( 18) and ( 21) results in:

1 2 (1 + M1 ) M1 ŵ+ 1 - 1 2 (1 -M1 ) M1 ŵ- 1 + ŵS 1 = 1 2 (1 + M2 ) M2 ŵ+ 2 - 1 2 (1 -M2 ) M2 ŵ- 2 + ŵS 2 (22) 1 2
(1 + M1 )

ζ 1 ŵ+ 1 + 1 2 (1 -M1 ) ζ 1 ŵ- 1 + 1 ζ 1 (γ -1) ŵS 1 = 1 2 (1 + M2 ) ζ 2 ŵ+ 2 + 1 2 (1 -M2 ) ζ 2 ŵ- 2 + 1 ζ 2 (γ -1) ŵS 2 ( 23
)
in which the notation ζ = 1 + [(γ -1)/2] M2 has been introduced. Now the analytical relations concerning the acoustic response of compact nozzles [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF][START_REF] Leyko | Comparison of direct and indirect combustion noise mechanisms in a model combustor[END_REF] are

recalled. An unchoked nozzle, which is characterized by a subsonic flow everywhere, needs to be acoustically constrained both at the inlet and at the outlet. In order to acoustically characterize the nozzle independently from boundary conditions, the outlet is fixed as totally non-reflecting ( ŵ-2 = 0).

Note that this procedure also simplifies the final mathematical expressions. If the acoustic wave ŵ-

2
is of relevance, it can be considered in the equations through a reflection coefficient

R = ŵ- 2 / ŵ+ 2 .
Particularly, an expression for R AA for an arbitrary value of R = ŵ-2 / ŵ+ 2 can be found in the work of Lamarque and Poinsot [START_REF] Lamarque | Boundary conditions fo acoustic eigenmode computations in gas turbine combustion chambers[END_REF]. Since neither heat nor work is added to the system, entropy waves remain constant along the nozzle ( ŵS 1 = ŵS 2 ). Eqs. ( 22) and ( 23) are combined in order to obtain the following expressions:

R AA = M1 -M2 M1 + M2 M1 + 1 M1 -1 1 -1 2 (γ -1) M1 M2 1 + 1 2 (γ -1) M1 M2 with ŵS 1 , ŵ- 2 = 0 (24)
R SA = - M2 -M1 1 -M1 M1 1 + 1 2 (γ -1) M1 M2 with ŵ+ 1 , ŵ- 2 = 0 (25)
T AA = 2 M2 1 + M2 1 + M1 M1 + M2 1 + 1 2 (γ -1) M2 2 1 + 1 2 (γ -1) M1 M2 with ŵS 1 , ŵ- 2 = 0 (26)
T SA = M2 -M1 1 + M2 M2 1 + 1 2 (γ -1) M1 M2 with ŵ+ 1 , ŵ- 2 = 0 (27)
In a choked nozzle, when assuming the outlet pressure low enough, the flow becomes supersonic downstream the throat. In this case only acoustic boundary conditions at the inlet of the system should be imposed. Imposing acoustic boundary conditions at the outlet leads to an ill-posed problem: setting a given reflecting condition (for instance ( ŵ-2 = 0)) is physically not correct. In order to find the different wave coefficients shown in table 1, the problem is divided into two. First, the subsonic region is solved. Accordingly, the throat of the nozzle is established as the outlet with M2 = 1. Second, the entire domain is solved by accounting for the results obtained in the first step.

The final coefficients are then obtained yielding:

R AA = 1 -1 2 (γ -1) M1 1 + 1 2 (γ -1) M1
with ŵS 1 = 0 (28)

R SA = -M1 1 + 1 2 (γ -1) M1
with ŵ+ 1 = 0 (29)

T AA = 1 + 1 2 (γ -1) M2 1 + 1 2 (γ -1) M1
with ŵS 1 = 0 (30)

T SA = 1 2 ( M2 -M1 ) 1 + 1 2 (γ -1) M1 with ŵ+ 1 = 0 (31)
Equations ( 24) to ( 31) are classical equations for compact nozzles, first encountered in the work of Marble and Candel [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF]. It is interesting to observe that the coefficients R AA and R SA do not depend on the Mach number M2 at the supersonic outlet. Indeed, as expected, they only depend on M1 and on the flow conditions at the the throat, where M = 1.

B. Non compact nozzles

The compact assumption is very useful to study the acoustic response of a nozzle when the acoustic wave length λ is much larger than the nozzle length. Nevertheless, when acoustic and entropy disturbances fluctuate at large frequencies the jump conditions showed previously become unsuited since they do not account for phase changes across the nozzle. In order to study the influence of 'non-compactness', different methods are found in the literature [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF][START_REF] Stow | Reflection of circumferential modes in a choked nozzle[END_REF][START_REF] Moase | The forced response of choked nozzles and supersonic diffusers[END_REF][START_REF] Candel | Analytical Studies of Some Acoustic Problems of Jet Engines[END_REF]. Marble and Candel [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF] stated that the main influence of the nozzle finiteness is to alter the phase between the two wave trains. They showed it analytically by considering that the mean flow velocity changes linearly through the nozzle. Nevertheless, this assumption is too restrictive for real systems. In this section, the response of non compact nozzles is computed by solving numerically the Linearized Euler Equations under the quasi-1D assumption (see section II).

Two different nozzle configurations are studied. They are shown in Figure 3 and the mean flow and geometrical parameters are given in table 2.

A 1 A c A 2 M1 M2 A 1 A 2 M1 M2 (a) (b)
Fig. In order to test the reliability of the SNozzle tool (see section II), the results given by this solver are compared to those obtained by solving the complete set of Euler Equations using a numerical solver for the unsteady fluid flow equations. Besides, they are compared to the analytical results obtained under the compact assumption. More precisely, the following approaches are compared in what follows:

• Analytics: The acoustic parameters R AA , R SA , T AA , T SA are computed through the relations proposed by Marble & Candel [START_REF] Marble | Acoustic disturbances from gas nonuniformities convected through a nozzle[END_REF]. Equations ( 24) to [START_REF] Lamarque | Boundary conditions fo acoustic eigenmode computations in gas turbine combustion chambers[END_REF], and ( 28) to (31) are then used.

• Quasi-1D LEE (SNozzle): The inlet Mach number M1 , temperature T1 and pressure p1 are given as inputs to SNozzle as well as the geometry (x coordinate and variable section area A). The outlet is taken as totally non reflecting ( ŵ-2 = 0) for the non-choked case and is simply non constrained for the choked case. At the inlet either ŵ+ 1 or ŵS are imposed. Two thousand points are used for both geometries. Note that it was checked that grid convergence is reached for this level of spatial resolution up to a reduced frequency Ω = 1. Results were validated against the analytical quasi-1D solver proposed in the work of Duran and Moreau [START_REF] Duran | Solution of the quasi-one-dimensional linearized euler equations using flow invariants and the magnus expansion[END_REF].

• The 3D Euler Equations (AVBP): The LES solver AVBP [START_REF] Avbp | [END_REF] is used, in which the viscous terms of the Navier Stokes equations have been neglected. This numerical tool, developed at CERFACS, is second order accurate in space and time. The Euler equations are spatially discretized under the finite volume formulation with a cell-vertex approach. Care has to be taken specially with the use of the inlet/outlet boundary conditions, since both the mean flow and either the entropy or acoustic wave must be imposed [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF]. Special attention should be paid on how these waves are introduced into the problem since, on the one hand, a totally non-reflecting inlet/outlet are sought (the imposed and reflected waves, ŵ+ 1 and ŵ-1 , should be completely evacuated) and, on the other hand, the mean flow must remain constant. In AVBP, a NSCBC approach [START_REF] Poinsot | Boundary conditions for direct simulation of compressible viscous flows[END_REF] is used for the treatment of waves crossing the boundaries.

C. Results

The value of both transmission and reflection coefficients are shown in Fig. [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF]. They are plotted against the non-dimensional frequency Ω = f l n /c 1 where f is the frequency, l n the length of the nozzle and c1 the speed of sound at the inlet. It is observed that analytical results are well suited when considering small frequencies (Ω → 0). Both numerical tools converge to Marble & Candel expressions when λ ≫ l n for all coefficients.

At small frequencies, it is observed that the main influence of a choked nozzle is the increase of noise ( ŵ+ 2 ր) due to entropy waves, as already stated in [START_REF] Leyko | Comparison of direct and indirect combustion noise mechanisms in a model combustor[END_REF]. This can be seen in the value of T SA ( ŵ+ 2 / ŵS 1 ) which is as twice as large when the nozzle is choked. The reflected wave ŵ-1 also increases, but to a lesser extent. On the other hand, the acoustic response of the nozzle caused by an incoming 14 When the reduced frequency increases, the value of the acoustic coefficients starts to move away smoothly from the analytical results. On the one hand, when taking into account the unchoked case, it can be observed that results coming from both numerical solvers can be considered similar up to Ω ≈ 0.2. For larger frequencies, the 3D Euler solver curves representing R AA , T SA and R SA start to decrease faster, the largest differences appearing for R AA . On the other hand both numerical solvers give similar results when considering the choked case for the values of Ω under consideration. It is interesting to notice that AVBP, a numerical tool used in this study to solve the 3D Euler equations with second order of accuracy in space and time, is able to retrieve quite well the acoustic response of a chocked nozzle whereas underestimates the acoustic reflection R AA occurring in an unchoked nozzle. This is probably due to the numerical scheme used in the 3D Euler solver, which appears to be dissipative for small length waves, which is the case for the downstream travelling waves ŵS 1 and ŵ+ 1 in an unchoked nozzle for high values of Ω. Recall that the mesh is kept unchanged when Ω increases, so that numerical errors increase.

! !"# !"$ !"% !"& '! !# '! !' '! ! '! ' ! Unchoked Case M2 =0.

V. Response of an ideal compressor

In the previous sections, focus was made on the acoustic response of ducts with varying sections causing a change in the mean flow: the nozzle case. The purpose of this section is to focus on the acoustic response of a system when the changes in the mean flow are caused by enthalpy jumps.

An ideal compact compressor is considered here. The first step when modeling a compressor is to consider it as an element that creates a difference in both the total pressure and the kinetic energy in the flow, i.e., an element that exerts a work on the flow by changing its total enthalpy. As represented in Fig. 5, a jump of total enthalpy is considered in the middle of the constant section duct. This jump comes from the characteristic curve of the turbo-machinery.

w -

2 w + 2 M1 M2 Compact compressor -→ Imposed → Atmosphere - → R in l 1 1 2

Fig. 5 Scheme of the compact compressor configuration

In this configuration an upstream traveling acoustic wave ŵ-2 is imposed at the outlet while an acoustic reflection R in is imposed at the inlet. This analytical/numerical setup is a simplified representation of the inlet air circuit and compressor of an aeronautical engine. The acoustic waves, produced in the combustion chamber, travel upstream through the inlet air circuit and the compact compressor until reaching the atmosphere which is considered as a partial reflecting acoustic condition. A downstream acoustic traveling wave is subsequently reflected and propagates all the way through the inlet air circuit and compact compressor until arriving at the inlet of the combustion chamber. The acoustic reflection coefficient R out = ŵ+ 2 / ŵ-2 is then estimated. As done before, this study is carried out both analytically and numerically (SNozzle).

It should be noted that, although not explicited in this paper, a similar methodology can be applied for the acoustic response of a system composed by several turbines stages and a nozzle.

Yet, some simplifications in the acoustic model can be made when evaluating the impedance at the combustion chamber outlet. Indeed, the flow is generally choked at the first high pressure turbine stage [START_REF] Gordon | Aerothermodynamics of aircraft engine components[END_REF]. In such a case only the upstream flow influences the acoustic impedance at the outlet of the combustion chamber. Moreover, a simple expression for the acoustic impedance of such a boundary condition can be retrieved analytically in the low frequency limit [START_REF] Lamarque | Boundary conditions fo acoustic eigenmode computations in gas turbine combustion chambers[END_REF]. A stronger simplification and still satisfactory model can be made as well if acoustics of this choked outlet is seen as a wall: an acoustic boundary where perturbations of velocity are not allowed.

A. Modeling of the compact compressor element

The compressor may be seen as an element that exerts a change in the total enthalpy of the flow.

In a two dimensional model, this jump of total enthalpy through the compressor can be properly modeled by the conservation of rothalpy I as follows:

I = h t,1 -U 1 v θ,1 = h t,2 -U 2 v θ,2 or h t,1 -h t,2 = U 1 v θ,1 -U 2 v θ,2 (32)
where U is the angular velocity of the compressor and v θ is the projection of the velocity vector v on the plane of U . In a one dimensional model, however, tangential velocities v θ , which are normal to the flow stream u, do not count among the unknowns and cannot be used to model the rotor effect. The modeling strategy carried out in this work relies on the characteristic curves, i.e. the compressor map, which relates the total pressure ratio π c to the mass flow ṁ [START_REF] Leonard | A quasi-one-dimensional CFD model for multistage turbomachines[END_REF] as shown in Fig. 6. O p e r a t i n g l i n e The total enthalpy ratio is then given by:

h t,2 h t,1 = T t,2 T t,1 = π T where π T = π (γ-1)/ηpγ c ( 33
)
where π T and η p stand for the total temperature ratio and the polytropic efficiency respectively.

The total pressure ratio π c is defined by a second order polynomial π c = a ṁ2 + b ṁ + c and clearly depends on the performance characteristics of the compressor. The logarithmic derivative of T t can be then computed as:

η p γ γ -1 Tt,2 Tt,2 - Tt,1 Tt,1 = πc πc where πc πc = 2a m2 + b m a m2 + b m + c m m = K m m and m = ρūA (34)
Adding a compact one dimensional compressor to the system is therefore equivalent to add an additional term to the left hand side of the energy equation ( 19), as follows :

Tt,1 Tt,1 + γ -1 η p γ K ρ1 ρ1 + û1 ū1 = Tt,2 Tt,2 (35)
in which the constant K characterizes the compressor at a given operating point. As function of the Mach number M and the frequency ω of the acoustic perturbations, the continuity (Eq. ( 18))

and energy (Eq. ( 35)) equations read respectively:

(1 + M1 ) M1 B + 1 e jωx/c1(1+ M1) - (1 -M1 ) M1 B - 1 e -jωx/c1(1-M1) = (1 + M2 ) M2 B + 2 e jω(x-l1)/c2(1+ M2) - (1 -M2 ) M2 B - 2 e -jω(x-l1)/c2(1-M2) (36)
(1 + M1 )

ζ 1 + K(1 + M1 ) η p γ M1 B + 1 e jωx/c1(1+ M1) + (1 -M1 ) ζ 1 - K(1 -M1 ) η p γ M1 B - 1 e -jωx/c1(1-M1) = (1 + M2 ) ζ 2 B + 2 e jω(x-l1)/c2(1+ M2) + (1 -M2 ) ζ 2 B - 2 e -jω(x-l1)/c2(1-M2) (37)
The inlet acoustic condition (at x = 0) is given by the inlet reflection coefficient as

R in = B + 1 /B - 1 .
Combining Eqs. ( 36) and [START_REF] Leonard | A quasi-one-dimensional CFD model for multistage turbomachines[END_REF] with the inlet acoustic condition and making the incoming perturbation to be ŵ-2 = 1 leads to a 4 × 4 system as follows:

C             B + 1 B - 1 B + 2 B - 2             =             0 0 0 1             (38)
where the matrix C, for a compact compressor placed at x = l 1 , is given by

C =             1 -R in 0 0 (1+ M1) M1 e jωl1/c1(1+ M1) -(1-M1) M1 e -jωl1/c1(1-M1) -(1+ M2) M2 (1-M2) M2 (1+ M1) ζ1 + K(1+ M1) ηpγ M1 e jωl1/c1(1+ M1) (1-M1) ζ1 -K(1-M1) ηpγ M1 e -jωl1/c1(1-M1) -(1+ M2) ζ2 -(1-M2) ζ2 0 0 0 1             (39)

B. Inlet acoustic condition

The acoustic flux f ac , normalized by the mean value of the energy flux, is defined as

f ac mh t = ṁ′ m T ′ t Tt (40)
Replacing now the definition of f ac by the left handside of Eqs. ( 22) and ( 23) and neglecting entropy waves ( ŵS = 0) leads to a definition of the acoustic flux written in the form:

f inlet ,ac m1 h1,t = (1 + M1 ) 2 M1 ŵ+ 1 - (1 -M1 ) 2 M1 ŵ- 1 (1 + M1 ) 2ζ 1 ŵ+ 1 + (1 -M1 ) 2ζ 1 ŵ- 1 (41)
In order to focus on the effect of the acoustic energy generated inside the system on the outlet reflection coefficient R out , it is mandatory to assure that no acoustic flux is being imposed at the inlet of the configuration. For that purpose, one may assume either ṁ′ = 0 or Ṫt ′ = 0. For open apertures (not a wall), the second alternative is better suited. Satisfying T ′ 1,t / T1,t = 0 leads to a ratio between the downstream and upstream acoustic waves given by:

R in = w + w -= B + /B -e jφ = (1 -M1 ) (1 + M1 ) e jπ (42)
This reflection coefficient acts as a totally reflecting condition (R in = 1e jπ ) for zero Mach number flows and tends to zero for inlet flows reaching a choked condition.

C. Compressor modeling in Snozzle. The compressor is modeled in Snozzle as an isentropic compact element, which is placed between the nodes i -1 and i as illustrated in Fig. 7. Between these two nodes, the fluctuating mass and energy are related by:

ρi-1 ρi-1 + ûi-1 ūi-1 = ρi ρi + ûi ūi (43)
1 ζ i-1 Mi-1 ûi-1 ci-1 + ρi-1 ρi-1 + γ -1 η p γ K ρi-1 ρi-1 + ûi-1 ūi-1 = 1 ζ i Mi ûi ci + ρi ρi (44)
where ζ = 1 + 1 2 (γ -1) M2 and K, defined in Eq. [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF], is a constant related to the characteristics of the compressor and the mass flux at the operating point. The fluctuating velocity at the compressor inlet ûi,2 and the fluctuating density at the compressor outlet ρi+1,1 are expressed as function of neighboring parameters:

ûups = ûi-1 = Γ • ûi + Λ • ρi-1 (45) ρdown = ρi = Υ • ûi + Ω • ρi-1 (46)
After combining Eqs. ( 43) and (44) with Eqs. ( 45) and ( 46) , it can be found that the constants Γ, Λ, Υ and Ω are defined by:

Γ = ūi-1 c2 i-1 ζ i-1 η p γ(ū 2 i -c2 i ) ūi c2 i Θ (
Λ = ūi-1 c2 i-1 [η p γ(ζ i-1 -ζ i ) -ζ i ζ i-1 K(γ -1)] ρi-1 Θ (48) Υ = ρi ū2 i c2 i-1 ζ i-1 η p γ -ū2 i-1 c2 i ζ i η p γ -c2 i c2 i-1 ζ i ζ i-1 K(γ -1) c2 i ūi Θ (49) Ω = ρi ζ i η p γ(ū 2 i-1 -c2 i-1 ) -c2 i-1 K(γ -1)(ζ 2 i -ζ i ζ i-1 ) ρi-1 Θ (50)
where

Θ = ū2 i-1 ζ i η p γ + c2 i-1 ζ i ζ i-1 K(γ -1) -c2 i-1 ζ i-1 η p γ (51)
These parameters are injected into the discretized equations of continuity Eq. ( 11) and momentum Eq. ( 12) at nodes i -1 and i. In addition to validate the numerical tool Snozzle by comparing its results to the analytical ones, this section aims at studying the effect of a compact compressor on the outlet reflection coefficient defined as R out = ŵ+ 2 / ŵ-2 of a constant section intake. Five cases are evaluated as shown in table 3 and Fig. 8. No πc > 0 is considered here as it corresponds to unstable operating points. For the cases shown in table 3, the downstream mean values of temperature, pressure and Mach number are T2 = 600 K, p2 = 200 kPa and M2 = 0.1, respectively. 3). They are given by the intersection between the polynomial curves and the vertical line at ṁ = 0.1425 kg/s. Dashed lines corresponds to the cases for which the slope of πc is πc = -1.0 (A04) and πc = -1.5 (A05).

D. Results

Ω = f l n /c 2 .
As a first observation it is seen that numerical results coming from the Helmholtz solver match the analytical results very well for all cases. On the one hand, the modulus of R out seems to be a quantity that does not depend on the mean value of the pressure ratio at all: |R| out = 0.81 for all cases shown in Fig. 9. This result is actually expected since Eq. ( 35) tells us that the normalized fluctuation energy Tt / Tt must remain constant through the compact compressor as long as there is no pressure ratio fluctuations (π c = 0 or K = 0). Figure 9 also shows, on the other hand, that the argument of R out strongly depends on the value of πc . In a reference case (π c = 1), the argument of R out evolves linearly as expected for an homogeneous mean field. However, as soon as πc increases, the argument of R out moves away from this linear behavior and takes values that can be actually phase shifted up to π with respect to the reference case. It will be observed in the next section that small changes in the argument of R for an intake outlet can lead to very different estimations of the acoustic modes of a given combustion chamber.

In the low frequency limit, the downstream acoustic wave ŵ+ 2 can be analytically described by the following relation:

ŵ+ 2 = (1 -M2 )(Kζ 2 + γ M2 ) (1 + M2 )(Kζ 2 -γ M2 ) ŵ- 2 (52)
According to Eq. ( 52), values of ŵ+ 2 , and therefore R out , can be strongly amplified when K → γ M2 /ζ 2 . Compressor operating point is most of the time established for cases in which πc < 0 (see Fig. 6 ). For these cases the definition of the downstream acoustic wave ŵ+ 2 (Eq. ( 52)) does not present any singularity since K is always negative. Attention is only paid then to πc ≤ 0. Fig. 10 shows values of πc going from 0 until -1.5. On the one hand, it is observed from Fig. 10 that the values of R out , for the modulus, are always below those for which πc = 0. On the other hand, values for the argument of R out vary very little when decreasing the fluctuating quantity πc .

E. Model extension to a more general compressor

In the present work the whole compressor is modeled as a compact element that adds energy to the system by modifying its total pressure and kinetic energy. This can be extended to a non compressor by considering a sum of rotor/stator combinations. On the one hand, each rotor row could be modeled as a compact element that adds total enthalpy to the system. On the other hand, each stator row could be modeled as an isentropic element that adds a phase lag to the incoming acoustic waves. Finally, the sum of the contributions of each individual element would characterize the acoustics of an entire compressor stage. This a convenient approach in the framework of linear acoustics.

VI. Computation of the Reflection Coefficient on the Inlet air circuit of an Helicopter combustor chamber

The numerical tool described above can be used to compute the equivalent impedance of any air-intake duct with or without total enthalpy change at the compressor location. It is well known that the resonant modes of any combustion chamber are strongly dependent on the inlet and outlet acoustic boundary conditions used [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF]. On the one hand, the flow at the outlet of the combustion chamber (the high pressure distributor) is generally choked [START_REF] Gordon | Aerothermodynamics of aircraft engine components[END_REF]. In this case, and for low frequencies, it is known that the reflection coefficient can be considered close to that of a rigid wall. On the other hand, the same cannot be stated at the inlet, as acoustic properties of the combustion chamber inlet are fully dependent on the upstream configuration of the combustor. Wrong predictions of the resonant frequencies can be caused by simplified treatment of the inlet impedance.

Aeronautical engines are complex systems in which the combustion chamber is present just in a small region after the compressor stages. From the point of view of the combustor, the upstream region is then composed by several compression stages in which acoustics depends on enthalpy jumps (compressors), the intake inlet acoustic condition (usually related to the atmosphere) and the geometry of the entire intake system. In the remaining of this paper, the numerical tool discussed in sections II-V is used to assess the acoustic reflection coefficient at the inlet of the combustion chamber, i.e outlet of the air-intake duct R out . The intake belongs to an helicopter combustor and is shown in Fig. 11. Once this acoustic boundary condition is evaluated, it can be imposed in a Helmholtz solver computation or in a Large Eddy Simulation [START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF][START_REF] Wolf | Azimuthal instabilities in annular combustion chambers[END_REF] Low Pressure High Pressure Compressor Intake's inlet R out ?

x 0

x LP x HP Fig. 11 Air-intake duct configuration.

A. Mean parameters of the Intake Duct

SNozzle requires the mean temperature and pressure, the mean value of the Mach number at the outlet and the section area at several points of the curvilinear axis of the air-intake duct. The position and the total pressure jump of the compressors is also needed in order to characterize the mean flow. It is assummed, when describing the mean flow along the intake duct, that the compressors are acoustically compact. It means that the curvilinear length of the compressors (see Figs. 11 and14) is small compared to the acoustic wavelengths of interest. For this reason, the evolution of the mean flow along the compressors is left aside and, instead, a sudden jump in the total mean flow quantities is considered at the middle of the curvilinear axis of each compressor (dashed lines in Fig. 11).

In the present helicopter system, the compression of gases is not performed in an isentropic manner: instead of one, the polytropic coefficient has a value around η p = 0.88 for both compressor stages. This value of η p will be accounted for in one of the three cases considered for the evaluation of the mean field. This is done in order to illustrate the effect of small changes in the mean field on the estimation of the reflection coefficient. It should be pointed out that although the compressors are not anymore considered isentropic for the mean field, they are still taken as isentropic for the fluctuating field since no entropy waves are assumed to be generated (and therefore η p = 1 in Eq. ( 44)). Figure 12 shows the mean field of total pressure and total temperature corresponding to the helicopter intake considered for three different cases: a) constant total pressure along the intake duct, b) isentropic mean flow with isentropic total pressure jumps due to compressors and c) isentropic mean flow with non-isentropic total pressure jumps due to actual compressors with

η p < 1.
It should be noted that only by modeling the total enthalpy jumps through the compressors it is possible to estimate a reliable mean field, notably the Mach number that in fact determines the acoustic wavelengths. Most of existing quasi 1D -LEE solvers, however, do not account for any change in total enthalpy [START_REF] Moase | The forced response of choked nozzles and supersonic diffusers[END_REF][START_REF] Lamarque | Boundary conditions fo acoustic eigenmode computations in gas turbine combustion chambers[END_REF]. An illustration of the Mach number (normalized by the Mach number at the outlet of the intake) for the three cases already mentioned is shown in Fig. 13. As expected the introduction of the two compressor rows induces a strong variation of Mach number.

Accounting for the polytropic efficiency only reduces this parameter slightly mainly in the low pressure compressor.

Normalized curvilinear axis (m)

Air-Intake duct without total pressure jump Air-Intake duct with isentropic total pressure jump Air-Intake duct with non-isentropic total pressure jump Fig. 13 Mach number normalized by the corresponding values at the outlet of the air-intake duct.

B. Results

As described in section V, the reflection coefficient at the outlet of the intake duct R out is estimated by evaluating the ratio between the reflected acoustic wave ŵ+ 2 and the imposed acoustic wave ŵ-2 as illustrated in Fig. 14. The acoustic reflection at the inlet of the air-intake R in is partially reflecting and is associated to Eq. ( 42), as done in section V for the ideal 1D intake. The two compressors present in the inlet-air duct are assumed acoustically compact and the associated total pressure ratio π c , which is direclty related to the mass flow by the characteristic curves (see Fig. 6), is given by Eq. [START_REF] Kaufmann | Flow forcing techniques for numerical simulation of combustion instabilities[END_REF]. From section V it was observed, for an ideal 1D intake, that no significant changes in the estimation of R out (mainly phase) are obtained when either considering or not the fluctuating contribution πc . Therefore, in this study only the mean contribution of the total pressure ratio πc of each compressor is accounted for. Figure 15 shows both the estimated modulus and argument of R out for the above three cases studied. These two parameters are function of the frequency that ranges from 0 to 600 Hz. This upper value is taken assuming that a compressor remains acoustically compact if its length is ten times smaller than the acoustic wave length. The module of R out has a constant value of 0.9 when no total pressure ratio is considered whereas it smoothly oscillates between 0.9 and 0.98 when accounting for compact compressors. In spite of this, it can be considered that there is not any significative difference between the three cases studied for values of the module of R out . On the contrary big changes are found, as expected from the results of the previous section, when observing the values of the argument of R out : the phase difference between ŵ+ 2 and ŵ-2 increases faster as the frequency grows for the two cases accounting for total pressure changes. This fact may influence the resonance frequency of the combustion chamber dramatically. The influence of considering a polytropic coefficient when computing the mean field is finally observed. It is seen in Fig. 15 that almost no differences exist for both modulus and argument of R out if the mean field of the temperature is modified by applying this value of η pol .

VII. Helmholtz Solver Computation

AVSP [START_REF] Martin | LES and acoustic analysis of combustion instability in a staged turbulent swirled combustor[END_REF][START_REF] Selle | Joint use of compressible Large-Eddy Simulation and Helmholtz solvers for the analysis of rotating modes in an industrial swirled burner[END_REF][START_REF] Nicoud | Acoustic modes in combustors with complex impedances and multidimensional active flames[END_REF][START_REF] Benoit | Numerical assessment of thermo-acoustic instabilities in gas turbines[END_REF] is the Helmholtz solver used here to compute the different acoustic modes of the aero-engine of interest. Most of the time the outlet of the combustion system is choked, which corresponds, for low frequencies, to an acoustic condition similar as u ′ = 0. The acoustic condition at the inlet is most of the time unknown. Usually, this boundary condition is also approximated by u ′ = 0 even though it is actually not correct. The purpose of this section is then to show how the first two acoustic modes of the chamber under study vary according to the imposed inlet acoustic condition R out [ †] .

Figure 16 shows the values of the resonance frequency as function of the modulus and argument of R out . First, it is assumed an inlet reflection coefficient in which arg(R out ) = 0 and |R out | varies from 0.1 to 10. This large range is studied to gain an overview of the associated effects on the combustion chamber acoustics. Clearly, large values of |R out | would mean that external acoustic energy enters the system. Second, |R out | is maintained equal to 0.9 and arg(R out ) is varied from -π to π. It is seen that a significant variation is produced on the resonant frequency of the first two acoustic modes when the phase of R out is changed at the inlet of the combustor. This is notably the case for the first mode since this frequency can vary up to 400 when modifying the phase from -π to π. The same shift of phase can modify the resonant frequency of the second mode until 100

Hz. These results are obtained maintaining |R out | = 0.9. Similar results are obtained for a fixed modulus |R out | = 1.0, and therefore it can be stated that results shown in Fig. 16 are representative

[ †] In this section Rout is associated with the inlet of the combustion chamber, i.e. the outlet of the air-intake duct (see Fig. 11 and Fig. 18) . of the range of |R out | of interest (Fig. 15).

Varying the modulus of the reflection coefficient R out does not imply a significative change in any of the first two acoustic modes of the combustor. These results are expected, since only the phase of R out can notably alter the wavelength of the acoustic fluctuations inside the system, and, therefore, the corresponding resonance frequency. The modulus, on the contrary, does not alter considerably this wave length but, instead, acts on the amount of acoustic energy subtracted from the system. This should be clear if the modulus and argument of R out are thought as the acoustic resistance (real part of the acoustic impedance) and acoustic reactance (imaginary part of the acoustic impedance) of the acoustic boundary, respectively.

Results shown in Fig. 16 are obtained for a long range of values of R out by changing either the associated module or phase. Now, the reflection coefficient R out obtained by the procedure of section VI is used to compute the acoustic modes of the combustion chamber. The effect of compressors at the air-intake duct can then be assessed.

With this in mind, and in order to get an overall view of the possible resonant frequencies of the

VIII. Conclusions

A numerical tool has been developed to compute the acoustic field of quasi-1D configurations accounting for total enthalpy jumps. This kind of situation is in fact present in the inlet-air circuit of actual aeronautical engines. Modeling this total pressure jump is not an easy task in the framework of a quasi 1D system: the concept of rothalpy cannot be developed. A model is then considered in which the total pressure ratio π c of each compressor is related to the characteristics curves of this element at a given operating point. In this way, it is possible to express the mean and fluctuating values of π c as a function of the mass flow and, as a consequence, to formulate two conservation equations through this interface. Analytical solutions are subsequently developed for compact compressors immersed in constant section ducts. These solutions are compared to the results obtained by the Helmholtz solver for validation purposes. Finally, the outlet reflection coefficient R out of an air-inlet circuit of an actual combustion chamber is computed. Significative differences on the estimate of R out results when comparing computations in which total pressure jumps are either considered or not. It has been also observed that the argument of R out , on the contrary to the modulus of R out , plays an important role on the value of the eigenfrequencies of the combustion chamber studied.

Appendix

In order to validate SNozzle, analytical solutions are derived for compact regions where the mean enthalpy or entropy may change. The starting point for this derivation are the 1D Euler equations for mass, momentum and total enthalpy. The contribution of the RHS terms on the above equations goes to zero for a compact region (x 1 → x 2 ) as long as only bounded fluctuations exist within this region, as also assumed in the work of Dowling [START_REF] Dowling | The calculation of thermoacoustic oscillations[END_REF]. Equations ( 59 where the indices 1 and 2 denote the regions upstream or downstream of any infinitely thin compact element; Q, F and Ŵ being the total heat, force and work associated to this element.

∂ρ

Fig. 1

 1 Fig. 1 Staggered grid

Fig. 2

 2 Fig. 2 Configurations under the compact assumption. (a) subsonic-subsonic. (b) subsonicsupersonic.

8 'Fig. 4

 84 Fig. 4 Reflection and Transmission coefficients for both unchoked and choked cases

Fig. 6

 6 Fig. 6 Map of a typical compressor

Fig. 7

 7 Fig. 7 Compact compressor in Snozzle

Table 3

 3 Five different configurations as function of πc and πc. For this cases the polytropic efficiency ηp is considered equal to one. Values of πc and πc are taken from typical values of actual aeronautical compressors.

Figures 9 and 10 Fig. 8

 108 Figures 9 and 10 show the different values of the reflection coefficient R out = ŵ+ 2 / ŵ-2 for the cases shown in table3. This reflection coefficient is plotted against the non-dimensional frequency

Fig. 9

 9 Fig. 9 Reflection coefficient at intake's outlet. The fluctuating value π ′ c is set to zero while the mean value πc varies from 1.0 until 5.0. Symbols and straight lines represent numerical and analytical results respectively.

Fig. 10

 10 Fig. 10 Reflection coefficient at intake's outlet. The mean value of the pressure ratio remains constant at πc = 2.0 while the fluctuating value πc varies from 0.0 until -1.5. Symbols and straight lines represent numerical and analytical results respectively.

Fig. 12

 12 Fig.[START_REF] Schmitt | Large Eddy Simulation and experimental study of heat transfer, nitric oxide emissions and combustion instability in a swirled turbulent high pressure burner[END_REF] Total pressure and total temperature profiles normalized by the corresponding values at the outlet of the air-intake duct.

Fig. 14

 14 Fig. 14 Model of air-intake duct.

Fig. 15 Reflection

 15 Fig. 15 Reflection Coefficient at the outlet of the helicopter intake duct.

! 9 Fig. 16

 916 Fig.[START_REF] Evesque | Low-order acoustic modelling for annular combustors: Validation and inclusion of modal coupling[END_REF] Dependance on the inlet reflection coefficient (Rout of the Air-intake) of the combustion chamber resonance frequencies of the first two acoustic modes -Helmholtz solver results.

Fig. 17

 17 Fig.[START_REF] Evesque | Spinning and azimuthally standing acoustic modes in annular combustors[END_REF] Reflection coefficient Rout at the outlet of the helicopter intake duct (grey lines).Resonance frequency of the combustion chamber when varying the argument of the reflection coefficient at its inlet (red lines). Intersection of red and grey lines (marked with crosses) represent expected resonance frequencies of the combustion chamber when coupled with the air intake under study.

  . The points of intersection between the red and grey lines, marked with a cross in Fig.17, may represent the resonance frequency of the combustion chamber when coupled with the air intake under study. The frequency associated with the second acoustic mode is practically not altered if R out corresponds to an intake either with or without compressors: in both cases it has a value around 500 Hz. In contrast, the frequency corresponding to the first acoustic mode is clearly influenced by the way R out is evaluated. Whereas it has values of ≈ 40 Hz if compressors are included in the computation of R out , it presents values of ≈ 80 Hz if compressors are not included in the computation. In order to verify this graphical analysis, computations are carried out with the Helmholtz solver. Results show that, on the one hand, when compressors are considered to estimate R out , the resonant frequencies are 51.5 Hz and 507.4 Hz for the first and second acoustic mode, respectively. The second acoustic mode is illustrated in Fig.18. On the other hand, when compressors are not accounted for in the air-intake, the combustion chamber resonant frequencies are 86.8 Hz and 506.6 Hz, for the first and second acoustic mode respectivelty. It is then observed that the value of R out used to estimate the first acoustic mode of the combustion chamber is not trivial. The approach described in the present study might be useful in order to reach an accurate prediction of this frequency of oscillation.

Fig. 18

 18 Fig. 18 Second acoustic mode of the Aeroengine combustor (507.4 Hz). Rout for an air-intake duct with compressors.

  notations have been introduced: Qdx = Q, Fdx = F and Wdx = W. After linearizing and considering harmonic oscillations (φ ′ = φe -jωt ), Eqs. (56), (57) and (58) become

Table 1

 1 Wave coefficients definition

Table 2

 2 3 Choked and Unchoked configurations studied Mean flow inlet conditions and geometrical parameters

		Unchoked Choked
		Case	Case
	A1/Ac	5.651	5.873
	A2/Ac	1.000	1.032
	M2	0.8	1.2
	M1	0.1
	T1	1300K
	p1	800 kPa

  Integrating over the 1D domain between positions x 1 and x 2 results in

					∂t	+	∂ ∂x	(ρu) = 0	(53)
				∂ρu ∂t	+	∂ ∂x	(ρu 2 ) = -	∂p ∂x	+ F	(54)
	∂ρh t ∂t	+	∂ ∂x	(ρuh t ) -	∂p ∂t	= Q + W	(55)
	x2 x1	∂ρ ∂t	dx +	x2 x1	∂ ∂x	(ρu)dx = 0
						x2 x1	∂ρ ∂t	dx + (ρu)

  ), (60) and (61) are then simplified to

									ρ1 ρ1	+	û1 ū1	=	ρ2 ρ2	+	û2 ū2	(62)
	ρ1	ū2 1	ρ1 ρ1	+ 2	û1 ū1		+ p1	p1 p1	+ F = ρ2	ū2 2	ρ2 ρ2	+ 2	û2 ū2	+ p2	p2 p2	(63)
	ρ1 ū1 ht1	ρ1 ρ1	+	û1 ū1	+	ĥt1 ht1	+ Q + Ŵ = ρ2 ū2 ht2	ρ2 ρ2	+	û2 ū2	+	ĥt2 ht2	(64)


Georg T Becker

Sherman C Lo

David S De Lorenzo

Per K Enge

Christof Paar

Secure Location Verification A Security Analysis of GPS Signal Authentication

The use of location based services has increased significantly over the last few years. However, location information is only sparsely used as a security mechanism. One of the reasons for this is the lack of location verification techniques with global coverage. Recently, a new method for authenticating signals from Global Navigation Satellite Systems(GNSS) such as GPS or Galileo has been proposed. In this paper, we analyze the security of this signal authentication mechanism and show how it can be used to establish a secure location verification service with global coverage. This new security service can be used to increase the security of various different applications, even if they are not directly connected to navigation or positioning.

Introduction

The use of location based services has grown tremendously in recent years. One of the reasons is that the GNSS receivers have become very small and cheap. Many mobile phones, for example, already have GPS capability. There has also been research in the area of location based access control, mainly for wireless networks. However, there has been a lack of technologies that can provide location verification services with global coverage. But location has great potential to be used as a cryptographic primitive. In many cases, a communication partner can be identified primarily by its location. In other cases the location can be used as an additional authentication parameter to increase the strength of an authentication protocol. In [START_REF] Denning | Location-based authentication: grounding cyberspace for better security[END_REF] the idea of using GPS to establish a location verification service was first introduced. However, this approach relied on the fact that the verifier as well as the user both have a trusted GPS-authentication device. Furthermore, this approach is vulnerable to spoofing, especially since the selective availability was turned off. Recently, Lo et al. introduced a new way to authenticate GPS signals in [START_REF] Lo | Signal authentication -a secure civil gnss for today[END_REF] called SAGA. In this paper we analyze the security of this new signal authentication mechanism and evaluate its usability. The advantage of this mechanism is that it does not only increase the security of location self-verification, but it can also be used for location verification. Previous approaches of location verification either relied on trusted devices or on bidirectional systems with only local coverage such as distance bounding [START_REF] Chandran | Position based cryptography[END_REF]. In SAGA on the other hand, the location and time of the reception of any GPS signal can be securely determined. This enables a location verification service that does not need any trusted devices and that does not need to have a bidirectional communication with any location service provider. Therefore, SAGA can be used as a building block to set up different secure location based services, ranging from traditional applications such as secure positioning and secure tracking to new applications such as location based access control. The main difference between this paper and [START_REF] Lo | Signal authentication -a secure civil gnss for today[END_REF] is that in [START_REF] Lo | Signal authentication -a secure civil gnss for today[END_REF] SAGA is described from a technical and navigational point of view. It includes a proof of concept implementation and complexity estimations but does not explain any possible attacks on the system nor states any security assumptions. In this paper we will look at SAGA from the computer security point of view. The second chapter is aimed to provide enough information so that non-navigation experts can understand the functionality of SAGA and our security analysis. In the third chapter we then come to the security analysis of SAGA by explaining the possible attacks on the system. From this threat analysis we then derive the security assumptions under which the system is secure. The conclusions from this security analysis are drawn in the last section. As location verification is not commonly defined in the literature yet, we start this paper with a definition of location verification and location self-verification: In secure location verification from A to B, B can be sure that entity A was at location L A at time t. In secure location self-verification for A, A can be sure that A was at location L A at time t.

Secure Authentication for GNSS Applications (SAGA)

In this section, we give an overview how Secure Authentication for GNSS Applications (SAGA) works. A more technical description of SAGA with test results can be found in [START_REF] Lo | Signal authentication -a secure civil gnss for today[END_REF]. In this paper we explain SAGA using GPS, but as other GNSS systems such as GALILEO or COMPASS work similarly, they can be used with SAGA as well.

GPS background

In GPS, the position is determine by measuring the arrival times of signals from different GPS satellites. With these arrival times, pseudo ranges between the receiver and the satellites are determined. Trilateration is used to calculate the position of the receiver out of these pseudo-ranges. Signals from four different satellites are needed to solve the trilateration equations, three to determine the 3dimensional position and one to determine the accurate time. In the current GPS constellation (typically 24-32 satellites), all satellites transmit signals on at least two frequencies: L1 and L2 (at 1575.42 MHz and 1227.60 MHz, respectively). The GPS satellites transmit a civilian signal on the L1 frequency. This C/Acode sequence (for Coarse Acquisition) is publicly available free-of-charge to any user worldwide. GPS satellites also transmit a secret military signal on both L1 and L2. This P-code (for Precision) sequence is encrypted to deny access to unauthorized users, becoming the P(Y)-code. As all satellites transmit on the same frequency, code-division multiple access is used to ensure that the satellites do not interfere with each other. Figure 1 illustrates the GPS signal structure. The 50 Hz data signal gets added (XOR) with the 1.023 MHz C/A code sequence. This C/A code is a pseudo-random code sequence with a length of 1023 bits (also called chips). Each satellite has its own unique pseudo-random code sequence. The resulting combination of the 50Hz data message with the 1.023 MHz C/A code is transmitted using binary phase modulation on a sinusoidal 1575.42 MHz carrier frequency. By adding the 50 Hz data signal with the 1.023 MHz C/A code, the signal is spread over a wider bandwidth. In this way, the signal can be recovered, although the signal is transmitted roughly 20 DB below the thermal noise floor. To recover the signal, a code replica of the C/A code is correlated with the received signal. A correlation peak shows the presence of the C/A code in the received signal. A positive correlation peak indicates that the currently transmitted data bit is a '0', a negative correlation peak indicates a '1'. The military P(Y) code works in nearly the same way. Instead of the 1.023 MHz C/A code, a 10.23 MHz P code is used. This has the effect that the code is hidden deeper in the noise. To restrict access to this signal, the P code gets encrypted to the P(Y) code by adding a secret and very long pseudo-random sequence W to the P code. Both the length of the code and its hidden nature make it very hard to discover the P(Y) code. The C/A and P(Y) code are both transmitted on L1 with a frequency of 1575.42 MHz. The P(Y) code is shifted by 90 degree (π/2) in phase compared to the C/A code,(also called phase quadrature) so that the C/A code has its minimum and maximum when the P(Y) code is zero and vice versa.

Functionality of SAGA

We now describe how the military P(Y) signal can be used to securely authenticate the GPS signal, even without the knowledge of the secret W code. The main idea is to exploit the fact that the P(Y) code sequence received at location 1 is identical to the sequence received at location 2, except for the difference in satellite-to-receiver signal travel time (and some frequency differences due to the receiver clock and Doppler effect). The idea is to cross-correlate the samples taken at location 1 and location 2 with each other. This will result in a correlation peak when the phases of the P(Y) codes of the two samples are aligned with each other. The presence of this peak indicates the verifier that the same code is hidden in both samples. Of course, the C/A code in both samples can create a correlation peak as well. But keep in mind that the C/A and P(Y) code are orthogonal to each other. Hence, the verifier knows exactly where the P(Y) code should be located in the noise. The verifier takes sample points when the sine term of the C/A carrier goes to zero, so that the carrier of the P(Y) signal goes to one or minus one. 4 Hence, only the P(Y) code creates a correlation peak. The C/A code is known for each satellite. This enables the verifier to match each correlation peak to the different satellites, as the P(Y) code is in phase with the C/A code of the same satellite.

The verifier can precisely measure the time offsets between the satellite signals of the two samples to determine the arrival times of the different signals.

If the verifier knows the arrival times of at least four satellites, the verifier can determine the exact location and time at which the sample was taken using the same trilateration methods that are used for the normal location determination in GPS. Hence, this system can be used to provide location verification from A to B. To provide location verification, the verifier B needs a signal S B that is valid and has not been spoofed. When A sends a signal S A to B, B can use the reference signal S B and the signal S A to verify the location and time at which S A was received.

The accuracy of this methods lies within a magnitude of the normal GPS positioning determination. Hence, the verifier can determine the position where the signal was received with an accuracy in the meter range.

Threat Model Against Location Verification Using SAGA

In this section, the possible attacks on location verification using SAGA are described.

Signal-synthesis attacks using the secret code: The hidden signals are generated using a secret code. In the case of GPS, this code is the military P(Y) code. An attacker who possesses this code can generate valid signals for every position he wants. Therefore, the system is only secure as long as the code is kept secret. If the attacker does not know the secret code used to generate the signals, he can try to guess the code. However, if the code is long enough and pseudo-randomly generated, it is computational infeasible for an attacker to successfully guess the secret code (which is true for the P(Y) code). It might be possible to use high-gain steerable antennas to raise the P(Y) code above noise. However, for every satellite one expensive and big (more than 10 meters) highgain steerable dish antenna is needed. Note that if the P(Y) code is revealed, the anti-spoofing capability and the restricted use of the military GPS signals would be broken. Hence, the same security assumption is needed for the military GPS anti-spoofing mechanism.

Signal-synthesis attack without the secret code: An attacker can try generating and transmitting his own navigation signals. Such an attack is called signal-synthesis attack. However, the attacker does not know the secret P(Y) code and therefore the attacker's signals will not match with the verifier's signals. Hence, the attacker would need to attack the verifier's reference signal as well. An attacker can try a signal-synthesis attack by inserting a hidden signal h into the verifier's reference signal. To be successful, the attacker needs centimeter knowledge of the verifier's antenna. The attacker sends a hidden signal to the verifier's antenna that is buried deep in the noise so that it does not interfere with the GPS signals. As the C/A code is not changed at all, and the hidden signals are buried well below the thermal noise floor, a verifier cannot detect the existence of these hidden signals. The attacker can now create a signal that the verifier falsely accepts as valid. To do this, the attacker first generates the C/A codes for the different satellites like they are expected at the wanted spoofing position. With the knowledge of the location of the verifier's antenna, the attacker can determine the travel time of the hidden signal to the verifier's antenna and therefore the offset between the hidden signal h and the C/A codes of each satellite in the verifier's data sample. Using the same offsets, the attacker aligns a copy of the hidden signal h with the C/A code of each satellite in his data sample. When the verifier correlates the data sample of the attacker with his data sample, the hidden signal h in the attackers's data sample and in the verifier's data sample correlates. The verifier cannot distinguish a correlation peak that is generated by the hidden signal h from a correlation peak that is generated by the P(Y) code and will therefore falsely accept the signal as valid. Without knowing the P(Y) code, the verifier will not be able to detect this attack. Signal observation techniques of the C/A code will be useless, as the original GPS signals are kept untouched by this attack. Using directional antennas to get the reference signal and shielding the antenna can make the attack much more complicated, as it would be more difficult for the attacker to insert the hidden signal into the verifier's signal. The insertion of the hidden signal can also be mitigated by collecting data samples from antennas at closely related locations (e.g. 3-5 meters). The attacker would need to align a hidden signal h i for each used antenna i. Furthermore, each antenna would receive the hidden signals h i with a different phase. Using cross-correlation techniques the presence of these signals is detectable.The verifier can also use reference signals from different places to increase the security. By using signals from different locations a web-of-trust can be build. This would significantly increase the attack complexity as the attacker would need to spoof each of these locations. Note that signal-synthesis attacks become very complicated in case an attacker needs to attack a receiver over-the-air, e.g. when he attacks location self-verification. In this case the attacker needs to somehow get rid of the original C/A code and P(Y) code (if cross-correlation is used to detect spoofing), which can be very complicated without physical access to the user's receiver.

Delay attack: In a delay attack, the attacker delays the incoming signals for the same amount of time. If all signals are delayed for the same amount of time, this has no impact on the position computation. However, this results in a clock offset at the receiver. B can still validate at which time A was at the location L A . But the clock of A and B are not synchronized. Therefore, delay attacks are very powerful against time synchronization, but do not have a direct impact on secure location verification. But if B's clock is not synchronized to a standard time reference such as UTC (from US Naval Observatory, GPS, etc.), e.g. because B is being attacked by a delay attack as well, B might falsely accept an old signal as fresh. Hence, as a requirement B's clock must be securely synchronized if B needs to decide whether the signal is fresh or not.

Selective-delay attack: In a selective delay attack, the attacker delays each satellite signal for a different amount of time so that a false position is calculated. This is a very powerful attack against navigation systems.[4] However, to be able to delay each signal for a different amount of time it must be possible to separate the signals from each other. But this is very difficult for the P(Y) signals as they are hidden in the noise. It might be possible to separate the signals by using high gain directional antennas for each satellite. Using a directional antenna pointing at one satellite, the C/A and P(Y) code of the target satellite are stronger than the signals from the other satellites.(But might be still below the thermal noise floor) If you combine signals from two directional antennas that target different satellites, a verifier might be able to detect the signals of the two satellites, while the signals of the other satellites might be too weak. Using this method the C/A and P(Y) signal from one satellite can be separated from the signals from the other satellites. But note that this attack needs at least four very good high-gain directional antennas and quite some knowledge in signal processing. Furthermore, this attack needs to be done in real-time, as the verifier can precisely determine the freshness of the signal. Whether this attack is successful depends on the verifier's ability to detect the signals from the not-targeted satellites. This strongly depends on the attacker's as well as the verifier's antennas and the verifier's effort to find these signals.

Relaying attack (wormhole attack): This is the most powerful attack against location verification with SAGA. In a relaying attack, the attacker relays the signals S v received at location L v to the attacker's location L A . As S v is a valid signal, B will falsely validate A's position as L v . This kind of attack is the biggest problem for all passive location verification techniques, as these techniques only verify the location of the received signal, and not of the receiver. So in passive location verification services that use GNSS techniques, it will only be possible to proof that an entity has access to a receiver (signal) at the claimed location, but not that he is actually there. As the exact reception time (less than a millisecond) of the signal is known these relay attacks can be made more difficult by setting up sharp bounds of the freshness of the signal. Note, that the accuracy of SAGA lies within the low meter range, hence, an attacker can only collect valid signals if he is within a few meters from the valid location.

Security Assumptions for Location Verification

We will now summarize the needed security assumptions in order to provide location verification services with SAGA.

1. B can be sure that the signals he has received are valid and no other signal than the P(Y) code is hidden in the noise. 2. An attacker does not have a signal from the claimed location L A for the claimed time period. 3. It is impossible to separate the signals from the different satellites from each other, so that they cannot be delayed for different amounts of time. 4. The attacker does not possess the secret code needed to generate the hidden signals.

Additional security assumption for location self-verification:

To prevent delay attacks, A either needs to be securely synchronized with GPS time or A needs to be sure that B's signal S B is fresh.

Conclusion

The new mechanism to authenticate GPS signals is very promising to enable secure location verification services. As GNSS signals cover great areas, only about 6 reference stations can provide reference signals that enable location verification with global coverage. However, looking at the security assumptions it is clear that careful consideration is needed for every application to decide whether or not secure location verification is possible with SAGA. The key assumption for SAGA is assumption number 2, that an attacker does not have a signal from a valid location. This assumption is not just limited to SAGA but is rather a general shortcoming of location verification: If a malicious user has a collaborator at the claimed location, the verifier cannot distinguish whether the received location signals are the user's or the collaborator's signal. Hence, he will not be able to know whether the user is at the claimed location or some collaborator. Therefore, SAGA should be used in applications where it is very unlikely that an attacker has access to a signal at a valid location at the claimed time. As an example application, a server with confidential information might restrict its access only to the company area and maybe the home of some employees that sometimes work from home. In this case, SAGA can be used for location based access control as an additional security mechanism for this server. Of course an attacker could try to circumvent this security mechanism by collecting a location signal at a valid location during the attack. However, this would significantly increase the complexity of the attack, especially as in many cases an attacker could be living far away from the target, e.g. in another country. Furthermore, the fact that the verifier will have meter knowledge of the attackers position, as well as the fact that he needs to be very close (in the meter range) to a specific location increases the chance that the attacker gets caught significantly. So location verification would not make the system unbreakable, but it could significantly increase the complexity of an attack. Hence, for many real-world systems, the proposed location verification techniques can significantly increase security. Furthermore, location verification can be a security tool that provides security in situation where traditional security mechanism such as passwords often fail. The main reason for this is that the security of location verification with SAGA does not rely on any secret information that can be lost and reused for later attacks. If a location signal is not fresh it is of no use for an attacker.

If SAGA needs to be resistant against very sophisticated attackers, the assumption that the P(Y) signals cannot be separated from each other might not be true, as an attacker could use very sophisticated high-gain steerable antennas. But in most cases, a possible adversary does not have access to such technology.

It should be further noted, that there currently does not exist any alternative to SAGA for using civil GNSS signals for location verification. Civil GNSS signals do not have any security mechanism so that spoofing can easily be done. [START_REF] Humphreys | Assessing the Spoofing Threat[END_REF] Hence, the security of SAGA far exceeds the security of normal civil GNSS services. Especially in applications, such as tracking, where location verification is the primary goal there is no alternative right now to SAGA with a comparable level of security when using GNSS. SAGA is also currently the most secure civil system for using GPS for location self-verification.

Fig. 1 .

 1 Fig. 1. The GPS signal structure on the L1 frequency.

In communication jargon, the verifier separates out the quadrature (sine) component from the in phase (cosine) component of the signal.