Muriel Amielh

Fabien Anselmet
email: anselmet@irphe.univ-mrs.fr

Yuchao Jiang

Ulf R Kristiansen

Pierre-Olivier O Mattei

Daniel Mazzoni

Cédric Pinhede

C Pinhède

Aeroacoustic source analysis in a corrugated flow pipe using low-frequency mitigation

Keywords: Corrugated pipe, acoustic mitigation, turbulent flow

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

Introduction

Due to their flexibility and local stiffness, corrugated pipes are used in many engineering applications. It is known that if the pipes carry a flow of dry gas, whistling sounds associated with the pipes' length resonances may develop. This phenomenon is used positively in musical toys like the Hummer [START_REF] Serafin | The voice of the dragon: A physical model of a rotating corrugated tube[END_REF][START_REF] Nakiboglu | Aeroacoustics of swinging corrugated tube: Voice of the Dragon[END_REF], but represents a problem in engineering applications as high sound levels and associated structural vibrations may develop. A well-known example is the "singing riser" phenomenon observed on some natural gas installations in the North Sea. Scientific study of the phenomenon dates back almost 100 years. Burstyn [START_REF] Burstyn | Eine neue pfeife (a new pipe)[END_REF] discussed briefly the phenomenon at the 1921 meeting of the German Society for Technical Physics. Around the same time, Cermak [START_REF] Cermak | Über die tonbildung bei metallschläuchen mit eingedrücktem spiralgang (on the sound generation in flexible metal hoses with spiraling grooves)[END_REF] studied the tones produced in short (0.5m) tubes having a narrow bore. He established that the frequencies excited are a little lower than the ones expected in smooth pipes, and that this could not be an effect of the end conditions. He also found that the frequencies heard correspond approximately to the "bump" frequencies based on flow velocity and corrugation pitch. More recent investigations show that the frequency effect is due to a lowering in the wave speed caused by the acoustic stiffness of the cavities ( [START_REF]Lecture Notes on the Mathematics of Acoustics[END_REF], see also section 3.2), and that the phenomenon occurs at a constant Strouhal number if the length parameter is taken as the cavity width plus the radius of curvature of the upstream cavity edge [START_REF] Nakiboglu | Aeroacoustic power generated by a compact axisymmetric cavity: Prediction of self-sustained oscillation and influence of the depth[END_REF]. Examples of experimental studies giving information on Strouhal and Reynolds numbers, necessary flow velocities, discussions on source positions, etc, for different test pipes, are found in references [7 -11]. The works cited above are performed on relatively short narrow bore pipes. A study directly related to the singing riser problem is the one by Belfroid et al. [START_REF] Belfroid | Flow induced pulsations due to flexible risers[END_REF].

Theoretical work on the subject includes the work by Debut et al. [START_REF] Debut | Flow-acoustic interaction in corrugated pipes: Time domain simulation of experimental phenomena[END_REF] who modelled the acoustic sources as a distribution of van der Pol oscillators along the pipe, and Goyder [START_REF] Goyder | On the modelling of noise generation in corrugated pipes[END_REF] who studied the distribution of sources and sinks along an acoustic wavelength in a corrugated pipe. More direct numerical studies have also appeared in recent years. Popescu et al. [START_REF] Popescu | Flow-induced acoustics in corrugated pipes[END_REF] used a turbulent model in Fluent to study the aero acoustic interaction in a short pipe. In parts of the pipe they found a strong interaction between cavity vortices and the shear layer separating the cavity fluid from the free stream. Nakibouglou et al. [START_REF] Nakiboglu | Aeroacoustics of swinging corrugated tube: Voice of the Dragon[END_REF][START_REF] Nakiboglu | Aeroacoustic power generated by a compact axisymmetric cavity: Prediction of self-sustained oscillation and influence of the depth[END_REF] modelled the flow in a pipe section including cavities using an incompressible Fluent version with a superposed oscillation corresponding to the frequency of a pipe mode. Their studies quantified the amount of sound power generated by the interaction of vortices shed at the cavities leading edges and the acoustic field.

Related to this, a study by Kopiev et al. [16] should be noted. In their study, a broad band signal is generated by a loudspeaker in a plenum chamber at one end of a corrugated pipe. They demonstrated that the acoustic signal would increase or decrease when propagating through the pipe depending on frequency and flow velocity.

Flow above single cavities with an associated sound production is a well-established field of study [START_REF] Howe | Acoustics of fluid-structure interactions[END_REF]. Of special interest to the present study is the finding that in cases where the pipe's entry section is sharp enough for the flow to separate, a small cavity a short distance downstream may give rise to strong longitudinal pipe modes [START_REF] Mazzoni | Aero-acoustics of a flow pipe having a single small cavity[END_REF][START_REF] Kristiansen | Aeroacoustic investigation of a flow pipe with a small cavity using the lattice Boltzmann method[END_REF]. The distance between the pipe's entry section and the single cavity's downstream corner is thought to be the relevant distance for Strouhal number estimations. As the corrugated test pipe used in the present study also has a sharp edged, partly flanged, entry section, it is expected that this will represent a possible source region.

The interested reader will also benefit from the recent review by Rajavel and Prasad [START_REF] Rajavel | Acoustics of corrugated pipes: A review[END_REF]. This paper summarizes the main available theories used to describe the singing phenomenon in corrugated pipes. It also reviews various experimental and computational works carried out for this problem and highlights remaining open questions such as the missing of the fundamental tone of corrugated pipes and the lack of prediction tools for the whistling effective acoustic sound pressure level.

The objective of the present research is to study both the acoustic and aerodynamic fields within a corrugated plane channel under singing conditions using both synchronised acoustic pressure and velocity measurements and particle image velocimetry (PIV).

We want to investigate more fully a previously published finding, mainly that adding a strong signal at a frequency very much lower than the whistling frequency of the channel, under some conditions, will quench the whistling [START_REF] Kristiansen | Experimental study of the influence of low frequency flow modulation on the whistling behavior of a corrugated pipe[END_REF]. In this previous study, an intermittent residual whistling noise remained. This whistling appeared and disappeared within a cycle of the low frequency signal. In the present study, we want to use this residual whistling in order to analyse the aeroacoustic source.

A model of the flow in a corrugated plane channel, based on a lattice Boltzmann technique, suitable for compressible flow modelling at low Mach numbers is faced to experimental results.

The paper is organized in the following way: section 2 describes the experimental arrangement and procedures, section 3 is devoted to the identification of sound pressure field characteristics inside the channel, section 4 reports the corrugated channel flow main characteristics, section 5 presents the analysis of the experimental results, with particular attention to the acoustic mitigation, section 6 presents the Lattice Boltzmann simulations, while section 7 presents discussion of the Lattice Boltzmann results. Finally, the final discussion and conclusion are presented in section 8.

Experimental arrangement

The experimental facility under study is made of a low speed wind tunnel whose fan rotation speed ensures a flow velocity up to 25 m/s. The tunnel is connected to a corrugated rectangular channel of length L = 2 m, prolonged by a loudspeaker itself ended by an open resonator of 17 m length in order to obtain high level acoustic oscillations at low frequency.

The rectangular corrugated channel is made of Plexiglas plates allowing direct inspection of the aerodynamic field by a PIV technique. The corrugated channel has lateral width D = 20mm and height B = 100mm. The lateral plates are machined in order to make a wall with a rectangular corrugation of pitch length Pt=20mm that is repeated all along the channel. The corrugation occupies the full height of the lateral walls of the channel. As shown in Figure 1, the corrugations are rectangular with a depth H=10mm, and a width w=10mm; the upstream corner of the corrugation is rounded with a curve with curvature radius r = 2.5mm. The top and bottom walls of the channel are perfectly smooth. The corrugations occupy about 83% of the wall surface of the channel. The entry of the corrugated channel is connected to a fan, as indicated in Figure 2, using a rectangular coupling box (see tick  in Figure 2). This box ensures, not only air tightness of the system corrugated channel/fan, but also that the channel has acoustically open flanged ends. The exit end of the corrugated channel is connected to a loudspeaker coupling box (see ticks and  in Figure 2) extended by the 17 m resonator whose resonance frequency is around 10 Hz. This resonator ensures that there is a strong coupling between the corrugated channel and the loudspeaker around its resonance frequency.

The loudspeaker is made of two push-pull sub-woofers SONIC-10MK2 (bandwidth f  [START_REF] Burstyn | Eine neue pfeife (a new pipe)[END_REF]300] Hz), mounted face to face and supplied with the same in-phase signals from a power amplifier ECLER MPA 4-150R. Velocity measurements are performed by hot-wire anemometry for temporal/spectral analysis and by PIV for spatial analysis. Two single hot-wires (DANTEC 55P11, 1.2mm length, 5m diameter) are used to characterise the flow velocity: the first one is located 10mm upstream of the vein entry (see tick  in Figure 2 and Figure 3) and the other one is downstream, inside the vein, 19mm just upstream of its exit section (see tick  in Figure 2). Hot-wire anemometry is controlled by the DANTEC Streamline system using the Streamware software. When not synchronously acquired with microphone signal by the Metravib NetdB system (24 bits), the hot-wire analogous signals are acquired by a 12 bit acquisition data card with a frequency of 25 kHz. For each hot-wire, 2 20 samples are acquired.

B

Figure 2. Photography of the experimental arrangement. Exit part of the corrugated channel and loudspeaker mitigation system; the fan is located on the left of the photography.

Hot-wire measurements are used to obtain the relation between the fan motor rotation and the flow velocity inside the vein (Figure 4). The upstream hot-wire is located in the vena contracta that feeds the corrugated channel, so that the velocity detected by this probe is lower than that measured by the downstream hot-wire. PIV is used to characterize the velocity field inside a cavity. The investigated cavity is the 35 th one, located at 700mm from the vein entry, ie at approximately one third of the vein length. Since the acoustic sound pressure inside the channel shows a succession of maxima and minima, (see Figure 6), this particular location was chosen because it is in the vicinity of a minimum of the sound pressure corresponding to a maximum of the acoustic velocity for the

Characterisation of the sound pressure field inside the pipe

Modal identification

Both ends of the corrugated pipe can be considered as open. Acoustic boundary conditions can then be described as Dirichlet conditions including a classical correction length to account for sound radiation [START_REF] Kristiansen | Experimental study of the influence of low frequency flow modulation on the whistling behavior of a corrugated pipe[END_REF][START_REF] Pierce | Acoustics. An Introduction to its Physical Principle and Applications[END_REF]. Such pipes present a discrete set of modes. Each mode is associated to a frequency f n given by

eff n c n, n 1 , 2 , . . . f 2( L 2 L )   
. c eff is the effective sound celerity in the corrugated pipe, it is slightly lower than the sound celerity inside the fluid due to the corrugations [START_REF] Elliott | Corrugated pipe flows[END_REF], and L is the correction length.

To our knowledge, no theory gives a sufficiently precise value of c eff for the present configuration. An estimate is proposed here through experimental investigation.

The first step consists in measuring the acoustic pressure at the 25 different positions described above. This is achieved by maintaining the fan rotation of 37 Hz that gives a constant flow speed of about 19.5 m/s inside the corrugated channel, allowing it to "whistle" at 592 Hz. Each record has a duration of 60 s. The 25 L p j rms pressure levels are estimated using the usual formula (1)

i=N 2 pj i j i=1 1 L = 10 log (p (t )) N  , (1)
where p j (t i ) is the acoustic pressure at instant t i , at point x j (x 0 = 195.1cm, x 1 = 194.1cm, …,

x 24 = 171.1cm), where the origin of abscissa x = 0 is the entry end of the channel. N = 2 18 ensures a nearly perfect estimate of the rms pressure level (the difference between the computation made with N = 2 18 and with N = 2 19 is about 0.05 dB).

An example of the spectrum of the recorded signal at x = 175.1cm is given in Figure 5 A comparison is presented in Figure 6 between the measured pressure and the identified one. This correction length is the simplest way to account for the effect of sound radiation. At low frequency for flanged circular pipe of radius R, the correction length is given by the well known result obtained by Rayleigh [START_REF] Strutt | Lord Rayleigh, The Theory of Sound[END_REF] :

L= 0.82 R.

(3) In the rectangular case, such a simple result is not known [START_REF] Boij | Flow effects on the acoustic end correction of a sudden in-duct area expansion[END_REF][START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF]. Nevertheless, it is possible to obtain accurate estimate of correction length for a rectangular channel by using results reported by [START_REF] Boij | Flow effects on the acoustic end correction of a sudden in-duct area expansion[END_REF][START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF]. For the presently considered rectangular geometry, with an aspect ratio B/D=5, from the numerical results presented in Figure 6 of [START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF], one obtains:

L 0.35 B/2   (4)
from which the correction length is L ≈ 17.5mm, which is in reasonably good agreement with the value estimated from relation [START_REF] Nakiboglu | Aeroacoustics of swinging corrugated tube: Voice of the Dragon[END_REF] if one considers that the present channel is connected to the loudspeaker coupling box and not to an infinite open space as assumed in [START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF]. This validates the experimental identification of the mode 8. From this, one can deduce the experimental effective sound celerity.

Effective sound celerity

It is well known that for acoustically rigid wall corrugated channels, the wave propagates along the channel at a slower, compared to a smooth one (i.e. non corrugated channel), apparent celerity c eff [START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF][START_REF] Elliott | Corrugated pipe flows[END_REF]. This effective sound speed can easily be deduced from experiments. The measurements show that the dominant frequency is close to 592 Hz. Since it corresponds to the eighth mode, the measured effective sound speed m eff c is computed by

m eff 2(L+2 L) c = 592 299 m/s 8   . (5)
Since the measurements are made at T = 293 °K, the true sound wave celerity is given by c 0 ≈ 20T 1/2 = 342 m/s. This apparent contradiction is based on the fact that inside the channel, it is difficult to define a clear propagating path for the wave since the apparent length of a corrugated channel is much longer than the corresponding smooth one. The simplest way to account for this is to use the model proposed by [START_REF] Elliott | Corrugated pipe flows[END_REF]. It considers that the corrugations act as if the walls of the channel were covered with spring-like, purely imaginary impedance with no energy loss/gain. This gives a value of c eff given, when the Mach number influence is neglected, by [26, p. 219, eq. 50-53]

eff 0 c p 1 c = c V Pw 1+ SP tS , (6)
where P = 2(B+D) is the perimeter of the tube, S= D.B the section of the tube, V c = 2w.H.B is the volume cavity, S p = 2w.B is the area presented by the cavity to the channel and Pt is the pitch length. With the experimental geometric data, one obtains an apparent velocity c eff ≈ 270 m/s that is clearly underestimated compared to the measured one, m eff c given by equation ( 5). This approximation can be improved by considering that this model supposes that the corrugations occupy the whole section of the channel, while in our experiments the corrugations occupy only 83 % of the inner surface of the channel since its top and bottom parts are smooth. A second improvement is obtained by noting that the spring-like effect of the cavity depends strongly on its size and, as the upstream edge of the cavity is rounded with a curvature radius r, then, the area presented by the cavity to the channel must be corrected to S p = 2 (w+r).B. Then, one obtains

2 eff 0 1 c= c w.H w ( B D ) 1+ 0.83 Pt w r D.B   , (7)
which gives an effective sound celerity c eff ≈ 288 m/s close to that deduced from experiments,

m eff c .
This formula can be used to estimate the effective sound celerity in rectangular partially corrugated channels.

Flow description

Since the whistling of corrugated pipes is induced by flow inside the pipe, let us describe its characteristics.

The flow velocities measured at the vein exit are varied in the range 0-25m/s according to the rotation frequency of the fan motor (Figure 4). For the previously described arrangement, the whistling begins to be detected on the upstream hot-wire signal at 16 m/s, corresponding to a fan rotation of around rot= 30 Hz. The upstream hot-wire captures the pure tone of the oscillations while the turbulence that develops all along the corrugated channel is superimposed to this sinusoidal oscillation on the downstream hot-wire signal (Figure 7). The whistling is associated to both a strong amplification of the velocity root-mean square measured by the upstream hot-wire and an increase of the sound pressure level from 50 up to 110 dB (Figure 8). The response of the system in presence of a flow is checked to be in agreement with the predicted eigen-frequencies. The whistling is detected on both velocity and pressure signals for which spectral analyses associate the strong values observed in figure 8 to the frequencies associated to 6 th to 9 th modes as shown in Figure 9. In this Figure, the Strouhal number based on these main peak frequencies f, the characteristic dimension (w+r) and the mean velocity U measured by the downstream hot-wire is shown to be constant around a mean value of St=0.36. This value is in agreement with the peak-whistling Strouhal numbers registered in corrugated pipes in the range 0.32-0.5 [START_REF] Nakiboglu | Whistling behavior of periodic systems: corrugated pipes and multiple side branch system[END_REF]. As indicated by Lattice Boltzmann simulations of the flow over a single cavity [START_REF] Kristiansen | Aeroacoustic investigation of a flow pipe with a small cavity using the lattice Boltzmann method[END_REF], the noise emission occurring when a corrugated pipe is submitted to an internal flow is associated to the formation of vortices. These vortex structures are convected by the flow to the downstream edge of the cavity and contribute to give a non-null time averaged acoustic power. PIV measurements are used to investigate the flow inside a cavity and detect possible structures. The instantaneous PIV fields show the presence of a vortex that occupies the whole cavity, but no advected structures are highlighted, probably due to the low data acquisition rate of PIV. The statistical analysis indicates that the mean vortex is almost centred in the cavity (Figure 10).

St v = f R vortex /V max .
The profiles of the mean vertical velocity of the vortex inside the cavity are extracted for each flow condition and are presented in Figure 11. The maximal velocity V max detected on each profile can be related to the modal frequency by the building of a new Strouhal number, 12) when R vortex = 4mm is the mean radius of the vortex deduced from the location of the extrema of V in Figure 11. Note that the slope l of the straight line that describes the relation between f and V max is homogeneous to a length and is estimated to l =10mm.

St v = f R vortex /V max .This Strouhal number is almost constant around St v =0.42 (Figure

Analysis of the experimental results

How an added low-frequency single tone sound reveals the origin of the riser song

As described in a previous study [START_REF] Kristiansen | Experimental study of the influence of low frequency flow modulation on the whistling behavior of a corrugated pipe[END_REF], a low-frequency acoustic field (around 10 Hz) leads to a dramatic reduction of the corrugated channel whistling. During the mitigation, as shown below in Figure 18, the very-low frequency carrier modulates the whistling in such a way that it appears only when the remaining low-frequency modulated velocity inside the channel is larger than a certain level.

In the velocity range 20-25m/s measured by the downstream hot-wire, four modes are simultaneously detected by the spectral analysis of the upstream hot-wire signal (see Figure 9). Above 17m/s, only one strong mode is highlighted that is to be related to the sudden increase of sound pressure level as shown in Figure 8. This pure tone is checked to be completely extinguishable by a low frequency signal emitted at 14 Hz by the two loudspeakers as shown on the downstream velocity spectral analysis, and without any obvious modification of the turbulence at the channel exit (Figure 13). As it can be inferred from Figure 13, the injected low frequency acoustic fluctuation level is rather high. In the work reported here, no particular optimization on both frequency and level of the signal imposed to loudspeakers has been performed. Upon analysing this phenomenon, we have observed, as shown below for the velocity signal in Figure 14, very short-time "wave packets" originating at the upstream entry of the channel when the flow speed reaches a sufficient level. These short-time wave packets will be shown to propagate toward the downstream end of the channel at the effective sound celerity, reflect back, propagate upstream and amplify the acoustic field inside the corrugated channel by successive reflections. This observation confirms that a sound source lies close to the pipe's entry section. It may also allow mitigation by controlling the low-frequency acoustic field to limit the amplification through wave superposition of the various reflected wave packets.

Two experimental mitigations at 14 Hz and 6 Hz are considered below. In both cases, the rotation speed of the fan is fixed at 37 Hz, allowing an upstream velocity (at 10mm in front of the channel entry) of 9.5 m/s and inside the channel, near its downstream end, of 19.5 m/s. For the first one, the frequency of the acoustic mitigation is 14 Hz and the remaining whistling is significantly reduced. In the second experiment, the mitigation frequency is set to 6 Hz with a remaining whistling level three times higher than that of the 14 Hz case. The two signals analysed and compared are the fluid velocity signal delivered by the entry hot wire probe and the acoustic sound pressure signal given by the downstream microphone probe. The results given in Figure 6 show that the acoustic pressure taken at x 0 = 1951mm is close to the maximum value and the minimum sound pressure close to the exit (the pressure acts similarly at the entry) ensures that the velocity fluctuation measured by the hot wire is maximum. In both measurements, the fluctuations created by the waves are then measured with a high signal to noise ratio. The analysis of signal data set is based on the comparison of the fluctuating parts of the flow velocity and acoustic pressure: continuous and low-frequency carrier contributions are removed from both signals. The remaining fluctuations are time frequency analysed and the time-shift correlation between velocity and acoustic pressure is computed.

In both cases (14 Hz and 6 Hz), the correlation has a maximum for a propagation time of 0.007 s very close to the time necessary for a wave to propagate from the entry to the exit at 299 m/s. This is an obvious proof that the noise source is at the entry of the channel. The small fluctuations created at the entry are amplified by successive reflections toward both corrugated channel ends leading to a limited or strong acoustic level depending on the mitigation frequency.

Acoustic mitigation at 14 Hz

The best mitigation results are obtained for an acoustic frequency close to 14 Hz.

In Figures 14 and15, as observed in our previous experiment [START_REF] Kristiansen | Experimental study of the influence of low frequency flow modulation on the whistling behavior of a corrugated pipe[END_REF], it can be seen that some remaining high frequency fluctuations are superposed over the velocity and pressure carriers. The study of these fluctuations is the core of our analysis. The continuous and carrier signals for velocity V c (t) and pressure P c (t) are described as a three term (higher order terms are negligible) Fourier expansion : The time-frequency representation (using a Gabor wavelet of frequency 6 [START_REF] Wolfram | Wolfram Mathematica[END_REF]) of the very beginning of both rectified signals is given in Figure 16. This shows, not only an obvious similarity between the two signals, but also a small time shift of the acoustic pressure compared to the velocity. A rough estimate gives this shift to be close to 0.007 s. Since 2/299 ≈ 0.0067 s, one can infer that it is precisely the time necessary for an acoustic perturbation to travel from the entry to the exit. To confirm this, one computes the Pearson correlation coefficient ( [START_REF] Wolfram | Wolfram Mathematica[END_REF]) between SP c (t) and SV c (t+) as a function of time delay . This correlation has been computed using 2 17 samples on both signals. As can be seen in Figure 17, the maximum value of correlation r, of about |r|=0.45, is reached at time  = 0.00689 s. This flying time is precisely that estimated by considering that the acoustic sources are located close to the entry of the channel. It is worth noting that the correlation fluctuates strongly at 669 Hz that is at the frequency of the remaining whistling.

V c (t) = A 0 +A 1 cos( f t+ 1 ))+A 2 cos(2 f t+ 2 )), with  f = 214

Acoustic mitigation at 6 Hz

Here also, the continuous and carrier signals for velocity V c (t) and pressure P c (t) are described as a three term (higher order terms are also negligible) Fourier

expansion: V c (t) = A 0 +A 1 cos( f t + 1 ))+A 2 cos(2 f t + 2 )),
with The frequency representation (using again a Gabor wavelet of frequency 6) of the very beginning of both rectified signals is given in Figure 20. This shows, not only a short shift of the acoustic pressure compared to the velocity similar to that observed for the 14 Hz mitigation, but also an interesting fine time structure of the four wave packets shown here (the same time structure can be observed over the complete record). As previously, each "big wave packet" of velocity and pressure is separated one from the others by again about 0.007 s, again corresponding to the time necessary for a wave to travel from the entry to the exit. Each big packet is made of a set of shorter packets each separated from the previous and following ones by a time of about 0.014 s (i.e. twice the flying time to the channel exit). Each of it has an intensity that increases until it reaches a maximum and decreases after that maximum. The explanation is that a small wave packet is emitted at the entry by the flow-cavity interaction. It propagates toward the exit end, at an apparent velocity of about 300 m/s (taking 0.007 s to reach the exit). It is reflected and back propagated to the entry. It is reflected again and increased in intensity by coupling with the entry source, and so on, until the entry source vanishes because of the flow mitigation. After that, the wave packet continues to reflect at both ends of the corrugated channel, but loses part of its energy by sound radiation and thus decreases in amplitude. For 14 Hz, this phenomenon is less evident due to the shorter time period of the mitigation. flying time is very close to that estimated if the acoustic sources are located close to the entry of the channel. The correlation is as twice bigger than previously. This is due to whistling that has much more time to be set on because the acoustic/flow mitigation is almost twice slower than previously and gives time to the acoustic field to take place inside the channel. It is worth noting that the correlation fluctuates strongly again at 669 Hz, that is at the frequency of the remaining whistling.

 f = 26 rad/s, A 0 = 9.34 m/s, A 1 = 2.
The acoustic mitigation could be optimized according to that fact. The frequency should be chosen in order to limit the wave reflection at the ends of the channel.

Lattice Boltzmann simulation

This section presents a numerical study based on the use of the Lattice Boltzmann method (LBM) to predict the noise emission from the corrugated channel flow. The Lattice Boltzmann method differs from traditional patterns by numerically solving an equation based on the physical statistics. Today, although the Lattice Boltzmann method is well established in the domain of Computational Fluid Mechanics, there are few studies based on its acoustic performance. Although it is limited to weakly compressible flows, this method is unsteady and compressible; thereby it can be used to simulate the acoustic waves generated by a turbulent flow. The present approach is a pragmatic one, as it consists in using the Fourier analysis on the results issuing from a Lattice Boltzmann computation to extract acoustic components from the flow data. This Lattice Boltzmann computation is performed by using the Palabos open source code (http://www.palabos.org).

Geometry of the computed problem

The geometry for the Lattice Boltzmann numerical simulation concerns a two-dimensional plane channel of 210 mm long corrugated on both inner sides with 10 square cavities of dimensions 10x10 mm 2 and 11 square corrugations of size 10x10 mm 2 (see Figure 22).

The Lattice Boltzmann Method

The Lattice Boltzmann method is implemented on the domain presented in Figure 22. At the entrance of the domain (x = 0), there is a non-homogeneous Dirichlet boundary condition for the fluid velocity. The bounce-back approach is applied for this boundary condition. At the outlet (x = 420 mm) a Neumann boundary condition is applied for the fluid velocity. Accurate initial conditions are crucial in the Lattice Boltzmann simulation. Therefore, a typical approach is to use an equilibrium distribution function of the fluid velocity. In the present numerical simulation, the initial flow in the channel is given by the classical 1/7 th power law turbulent mean velocity profile [START_REF] De Chant | The venerable 1/7 th power law turbulent velocity profile: a classical non-linear boundary value problem solution and its relationship to stochastic processes[END_REF]. In the upstream tranquilizer chamber the flow streamlines are contracted to adapt the flow to the change of section. This contraction is performed by using a shape-function S y . The 1/7 th power law turbulent mean velocity profile is multiplied by this shape-function whose role is to "curve" the streamlines to adapt smoothly the flow to the change of section. This shape function acts from the middle of the first tranquilizer chamber (x ≈ 50 mm) where the section height is D i = 50 mm to the entrance of the channel (x = 105 mm) where the section height is D = 20 mm. This shape function is given for the flow lowest streamline (y = -D i /2 in the entrance domain):

 i 4 4 ii D/2 HD S:y s i n 1 ; w h e r e : 0 22 2                    (8)
where H i = 15 mm is the upward step at the entrance of the channel. A scale factor is applied in order to match  = 0 with roughly the middle of the upstream tranquilizer domain  ensures a smooth initialization for the flow at the change of section and reduces the compressibility errors on this initial flow. This wariness is essential for providing accurate initial conditions.

The collision model for the numerical simulation is the Smagorinsky eddy-viscosity model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. The basic experiment[END_REF] which assumes that the turbulence small scales induce a viscosity correction proportional to the norm of the strain-rate tensor      

2 ts x,t C S x,t    , where   ij ij Sx , t SS 
is the magnitude of the strain-rate tensor,  is the bandwidth of the LES filter and C s is a nondimensional coefficient called the Smagorinsky constant: It is an adjusted parameter ensuring the numerical stability of the model by increasing the viscosity as the rate-of-strain of the flow increases. The value of this Smagorinsky coefficient ranges from 0.1 to 0.2 [START_REF] Meyers | On the model coefficients for the standard and the variational multi-scale Smagorinsky model[END_REF]. For our numerical simulation, the value of the Smagorinsky coefficient that seems to be most convenient in terms of flow description is C s = 0.12.

General data for the computed LB model

The Palabos code works with specific Lattice Boltzmann discrete units [START_REF] Latt | Choice of units in lattice Boltzmann simulations[END_REF]. As a first step, the physical system must be converted to a dimensionless one (Table 1 in Appendix 1). As a second step, this dimensionless problem is converted to discrete units. The dimensions of the data are in this last step specific to the Lattice Boltzmann model. A characteristic length 0 L and a characteristic time delay 0 T are defined as follows: 0t L P  is the length of one corrugation plus one cavity; this is the pitch herein. The characteristic time delay is

LB LB 0 eff c uV c  .
The flow is computed during LB T =120 in the discrete Lattice Boltzmann units. The data of the computed problem are listed in Appendix 1 (Table 1). This computation results in 0.7 s of flow for five days of computation using one processor (10 cores) of a computer with eight processors Intel Xeon E7-8870 running at 2.4 GHz. The sampling rate for the data output is 37300Hz .

Lattice Boltzmann results and discussion

The results of this Lattice Boltzmann simulation consist in a fluid density ρ and a fluid velocity vector V at each point of the mesh grid (Table 1, Appendix 1), at each time step. All these results are in discrete Lattice Boltzmann units. In the present section, three observation domains are defined. Ω 1 is a discrete set of 100 measurement points M i = (x i ,0) (i=1, …, 100) uniformly distributed along the axis of the channel from the entrance (x = 0) to the outlet of the domain (x = 420 mm). Ω 1 represents an antenna of microphones partitioned on the axis, it is used to identify mode-shapes of the acoustic pressure inside the channel. Ω 2 is a rectangular domain corresponding to the first cavity located at the inlet of the corrugated part. Ω 3 is a domain covering two cavities located in the middle of the corrugated part. The frequency peaks plotted in Figure 23 do not correspond to the experimental frequency peaks presented in Figure 9. At the velocity used for the computation (ie 9 presents an experimental maximum sound pressure emission at the sixth mode of the channel ( 6 / 0 = f f ) meanwhile Figure 23 shows a maximum peak at the first mode of the domain 0 f /f= 1 . 1 . There are too many differences in the geometry (length of the domains) and in the boundary conditions between the experiment and numerical simulation to conclude on these mismatches. No matter these mismatches, the numerical example proposed in the present paper must be considered more as a way to analyze the physical phenomena involved in the whistling than as an exact representation of the experiment. The sound celerity c inside the corrugated part is calculated from the measurement of the wavelengths λ that appear in Figures 25 and 26 by using the well-known formula c = f λ which gives: One of the major problems encountered in this numerical simulation comes from the fact that there is no acoustic absorbing boundary condition proposed in the present Palabos version. So, overcoming the boundary conditions of the corrugated part's outer domains is impossible. Moreover, these boundary conditions are described in terms of fluid mechanics and not for an acoustic behavior. As a result, the eigen-modes observed in Figures 24,[START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF] S masks the much lower acoustic one. This behavior also appears on the experiment: in Figure 8, the downstream hot-wire measurements (white triangles) cannot detect the whistling phenomenon because it is covered by the turbulent noise. As expected, the cross-spectra S  and Vx S are significantly reduced in the regions of the tranquilization chambers (Figures 24, 25 and 26): the thin and deep cavities which cover the walls of these regions seem to actually play their role of dampers.

1 0 V = 16ms  ) Figure
The next observations concern images of the vorticity field calculated in the domains Ω 2 (first cavity) and Ω 3 (middle of the corrugated part). In Ω 2 domain, Figures 27a andb illustrate the vena contracta phenomenon that occurs at the inlet of the corrugated part where the flow separates from the wall. At this flow separation, vortex structures are created in the region very close to the inlet. These structures grow along the separation bubble and detach at a somewhat regular frequency in the region close to the upstream rounded cavity edge before being carried downstream by the flow. Another noticeable phenomenon is the two contrarotating vortices inside the cavity, and the expulsion of fluid from the upstream rounded cavity edge that feeds the region of the separation bubble. . Similarly, the convection velocity at the cavity mouth can be estimated experimentally from velocity V max in the cavity ( §4) and mean velocity U on the corrugated part axis ( §2). In the region of cavity mouth, the velocity is identified as the inside vortex maximal velocity V max , so that in the shear layer U c =(V max +U)/2. The results are presented in Figure 29. It is shown that the convection velocity deduced from the LBM simulation is quite in a good agreement with its experimental estimation for U=16m/s. Moreover, the convection velocity is linearly related to the velocity on the axis. According to numerical results, vortices are regularly emitted from the rounded edge and cross the cavity mouth over a distance equal to the pitch length t P . Consequently an emission frequency can be estimated by ct U/ P and it is experimentally obtained to be equal to the whistling frequency f (see Figure 29).

Final discussion and conclusion

The present experimental facility and the associated numerical tools were shown to be well suited for investigating flow acoustic situations. A series of combined acoustic sound pressure and fluid velocity measurements were conducted under different configurations of mitigation by a very low frequency acoustic wave. The aim of these experiments was to identify clearly the acoustic sound field that takes place inside the corrugated channel under the flow excitation. A low frequency acoustic mitigation induces a flow speed modulation around the mean velocity. This modulation dramatically changes the nature of the sound fields created by the corrugated channel, from intense tone to much quieter periodic whistling. Under this mitigation, a residual fluctuation was measured on both the acoustic field inside the channel and the flow velocity at its entry. By eliminating the low frequency forcing term on both signals and by correlating them, it was shown that the sound source is not only located near the entry end of the corrugated channel but also that the residual acoustic field is due to the superposition of waves reflected at both ends of the channel.

The strong coupling that exists between the acoustics and the flow was also investigated by numerical simulation using the Lattice Boltzmann Method. The results of this numerical simulation were analyzed by a Welch method which permits to identify acoustic eigen-modes of the channel. Limitations arise when analyzing the acoustic behavior because the boundary conditions involved in the lattice Boltzmann schemes are mainly designed for the description of the flow, their acoustic properties are not considered. As a result, the tranquilizer zones located at the inlet and at the outlet of the corrugated part participate strongly to the acoustic eigen-modes of the computed domain. In the present numerical simulation, a solution consisting in a multitude of thin and deep cavities was proposed to create two soundabsorbing regions at the ends of the domain. This method appeared to be of limited efficiency because the mode-shapes present a Neumann boundary condition at the inlet of the domain in place of the expected Dirichlet boundary condition at the inlet of the corrugated part.

Recently, some definitions of absorbing layers for nonreflecting acoustic boundary conditions appeared in the literature [START_REF] Xu | Analysis of the absorbing layers for the weakly-compressible lattice Boltzmann methods[END_REF]. An outflow sponge layer is also proposed in the version 1.4 of Palabos, but only for three-dimensional problems, and the present model is a two-dimensional one. The incorporation of such a non-reflective sponge layer in future work should improve greatly the accuracy of the simulation, thus allowing quantitative acoustic comparison with experimental data.

Acknowledgments.

The work presented herein was supported by Institut Carnot STAR through the CREATUCE Project. The contribution of the mechanics staff of IRPHE in designing and machining the corrugated channel is acknowledged.

Figure 1 .

 1 Figure 1. Geometry of the corrugated channel. a) Schematic view, b) Photography of details of the face to face vertical corrugated plates.

Figure 3 .

 3 Figure 3. Photography of the entry hot-wire.

  wire (fit : U =1.47+0.487 rot) upstream hot-wire (fit : U =1.18+0.236 rot) mean velocity U (m/s) rot (Hz)

Figure 4 .

 4 Figure 4. Relation between the fan motor rotation and the mean velocity at the entry and exit of the corrugated channel measured by hot-wires.

2 MP

 2 is the amplitude measured, n = 8 indicates that the mode under consideration is the 8-th channel longitudinal mode. One obtains L = 1 cm as the correction length of the channel.

Figure 5 .

 5 Figure 5. Spectrum of acoustic sound pressure at position x =175.1cm for a fan rotation of 37 Hz.

Figure 6 .

 6 Figure 6. Variation of the pressure along the corrugated channel (cm) for a fan rotation of 37 Hz, Black dots: measurement, red curve: identification.

Figure 7 .

 7 Figure 7. Response of the corrugated channel to an internal flow. Velocity time series from hot-wires with whistling.

Figure 8 .

 8 Figure 8. Response of the corrugated channel to an internal flow. Amplification of U rms measured by hot wire at the vein entry and sound pressure level measured by microphone inside the vein.

  /s) Figure 9. Response of the corrugated channel to an internal flow. Whistling frequencies (from hot-wire measurements) and associated Strouhal numbers, St=f (w+r)/U.

Figure 10 .

 10 Figure 10. Mean vortex in the cavity: Velocity vectors measured by PIV and location of the extracted V-profiles.

Figure 11 .

 11 Figure 11. Mean vortex in the cavity: extracted mean velocity profiles from PIV measurements.

Figure 12 .

 12 Figure 12. Mean vortex in the cavity: associated modes detected by hot-wire measurements and Strouhal number based on the vortex characteristics: St v = f R vortex /V max .

Figure 13 .

 13 Figure 13. Whistling mitigation (extinction) by a low-frequency acoustic field, downstream hot-wire velocity signal spectral analysis.

Figure 14 .

 14 Figure 14. Time records, in s, of velocity (left) SV o (t), in m/s, and sound pressure (right) SP o (t), in Pa. The black curves are the original data, the red ones represent the identified carrier wave.

  rad/s, A 0 = 9.28 m/s, A 1 = 3.46 m/s,  1 = -6.04 rad, A 2 = -0.042 m/s,  2 = 0.13 rad and P c (t) = B 0 +B 1 cos(( f t+ 1 )+B 2 cos(( f t+ 2 ), with B 0 = 1.4 Pa, B 1 = -1469 Pa,  1 = -4.7 rad, B 2 = 84.25 Pa,  2 = -3.51 rad. The amplitude and phase coefficients are identified using the standard Mathematica [28] procedure. It is worth noting that rectified signals SV c (t) = SV o (t) -V c (t) and SP c (t) = SP o (t) -P c (t) in Figure (15) are alike.

Figure 15 .

 15 Figure 15. Time records, in s, of velocity (left) SV c (t), in m/s, and sound pressure (right) SP c (t), in Pa, of the rectified data.

Figure 16 .

 16 Figure 16. Time-frequency representation of the signal using Gabor wavelet transform of velocity (left) SV c (t) and sound pressure (right) SP c (t).

Figure 17 .

 17 Figure 17. Pearson's correlation coefficient for about half a second delay (left) and 0.1 second delay (right). The diamond marks the maximum absolute value.

Figure 18 .

 18 Figure 18. Time records, in s, of velocity (left) Sv o (t), in m/s, and sound pressure (right) Sp o (t), in Pa. The black curves are the original data, the red ones represent the identified carrier wave.

Figure 19 .

 19 Figure 19. Time records, in s, of velocity (left) SV c (t), in m/s, and sound pressure (right) SP c (t), in Pa, of the rectified data.

Figure 20 .

 20 Figure 20. Time-frequency representation of the signal using Gabor wavelet transform of velocity (left) SV c (t) and sound pressure (right) SP c (t).

Figure 21 .

 21 Figure 21. Pearson's correlation coefficient for about one second delay (left) and 0.1 second delay (right). The diamond marks the maximum absolute value. The Pearson correlation between SP c (t) and SV c (t+) as a function of  is again computed using 2 17 samples on both signals. The result is presented in Figure 21. One obtains the maximum value of Pearson's correlation r of about |r|=0.82 at time  = 0.0062 s. Again, this

Figure 22 .

 22 Figure 22. Geometrical arrangement for the Lattice Boltzmann simulation.The upstream and downstream edges of the corrugations are rounded with a radius equal respectively to 1 mm and 3 mm. The distance between two facing corrugations is 20 mm. This computed channel is a two-dimensional representation of the experiment detailed in the previous section where the length of the corrugated part is reduced from 2 m to 210 mm in order to reduce the computational effort. In the upstream and downstream regions of this corrugated part, there are two tranquilizer chambers of 105 mm long and 50 mm high. The role of these tranquilizer chambers consists in introducing a flow step at the inlet and outlet sections of the channel to simulate a channel with acoustically open flanged ends. At the air inlet, this upward step reduces sharply the section of the inflow of 60%, it results in an acceleration of the fluid velocity and a separation of the flow from the walls: A separation bubble forms at the air inlet that reattaches to the wall a few centimeters downstream in the region of the first corrugations. The walls of the tranquilizer region are made of a multitude of very thin (≈ 0.1 mm) and deep (10 mm) cavities. These cavities are thin enough so that no flow structure can develop in their interior. Their low thickness promotes viscous effects: their role is to attenuate the acoustic modes in the tranquilizers. The x-axis is classically in the middle of the domain with its origin at the left end side.

-

  of the corrugated part (x = 105 mm). performs rising of the lowest streamline, equal to the value of the upward step H i over the distance 50 mm. The highest streamline i with respect to the channel axis. For the streamlines which are between these two extreme streamlines, the shape-function y S( ) , with ii -D /2 < y < D /2, is calculated by linear interpolation between the functions i . The x-component V x of the fluid velocity in the change of section is given by the flow conservation. The y-component of the velocity V y in this change of section is given by the derivation of the shape-function

  / V, this is the delay for the fluid particles to cross the pitch t P . The Mach number of the flow is0e f f M = V /c .Considering the sound celerity on the lattice LB c1 / 3  and the flow velocity on the lattice LB u , the Mach number is expressed in Lattice Boltzmann units as LB LB M u/ c  . The fluid velocity on the lattice LB u is thus calculated from these two definitions of the Mach number:

  Ω 2 and Ω 3 are two domains on which the flow behavior is investigated and compared to the PIV measurements.The first comparison between the simulation and the experiment focuses on the interior acoustic field. Let us define an arbitrary reference point M ref located in the middle of the corrugated part. A Welch analysis is performed on the antenna Ω 1 . The results of this Welch analysis are the cross-spectrum for the fluid density , f between the points M i = (x i , 0) and M ref , at frequency f. For the analysis of the results, the value of f issuing from this Welch analysis of the Lattice Boltzmann results is expressed in dimensionless units. The dimensioning factor is the first eigen-frequency f 0 of a corrugated domain of length T L with Neumann-Dirichlet boundary conditions (here 0e f fT f c/ ( 4 L )1 7 8 H z ). Figure23presents the sum of all the 100 spectra of the fluid density  …, 100) calculated on the antenna Ω 1 and plotted as a function of the dimensionless frequency 0 / f f . The frequency peaks plotted in Figure 23 correspond roughly to the first eigen-modes of a corrugated domain of length T L with Neumann-Dirichlet boundary conditions. The n-th eigen-frequency of this corrugated domain of length T L with Neumann-Dirichlet boundary conditions is equal to in agreement with the expected ones according to the spectral resolution. This confirms that, for the first eigen-frequencies, the computed channel behaves as a corrugated domain of length T L with Neumann-Dirichlet boundary conditions.

Figure 23 .

 23 Figure 23. Normalized cross-spectrum of the fluid density .

Figures 24 ,

 24 Figures 24, 25 and 26 present the fluid density spectrum

1 Figure 24 . 7 Figure 25 . 11 Figure 26 .

 1247251126 Figure 24. Normalized cross-spectra for the frequency

  calculated from the wavelength measurements in Figures 25 and 26 is very close to the experimental value eff c = 299 m/s . This result confirms that the mode-shapes observed on the density power spectrum in Figures 25 and 26 are indeed acoustic modes.

  and 26 are not eigen-modes of a channel with open ends, but eigen-modes of a channel with unknown acoustic boundary conditions. Nevertheless, these acoustic boundary conditions can be estimated visually from the shape of the function S  presented inFigures 24,[START_REF] Duan | A hybrid finite element approach to modeling sound radiation from circular and rectangular ducts[END_REF] and 26. These functions S  correspond to a Neumann acoustic boundary condition at the inlet of the domain (x = 0): the non-homogeneous bounce-back method that has been used at the inlet to create the incoming flow seems to induce at this section an acoustic Neumann boundary condition. Note that the bounce-back method is also used to model the walls of the corrugated part for which a Neumann boundary condition is also expected. The outflow condition is a zerogradient boundary condition for the fluid velocity vector V. The acoustic boundary condition induced by this outflow boundary condition cannot be identified clearly from the observation of the mode-shapes in Figures 24, 25 and 26. The velocity spectra Vx S plotted on the bottom of Figures 24, 25 and 26 are very noisy. It is likely that the large band turbulence contribution to the spectrum Vx

Figure 27 .

 27 Figure 27. Vorticity in the inlet section (domain  2 ). (a) at time 0.43124 s, (b) at time 0.43237 s, (see associated video vorticity_inlet.avi)

Figure 28 .Figure 29 .

 2829 Figure 28. Vorticity in domain   (a) at time 0.42785 s, (b) at time 0.42898 s, (c) at time 0.43011 s, (d) at time 0.43124 s. (see associated video vorticity_middle.avi)

  , showing that the main dominant pressure peak has a frequency of about 592 Hz with a sound pressure level L p = 23.7 dB (ref 1 Pa) while the sound pressure level of the whole signal is close to L p 20 = 27.4 dB (ref 1 Pa). Even if this location is close to the sound pressure minimum, this measurement shows that the signal to noise ratio is sufficient to ensure a good estimate of the whistling. The correction length L has been estimated by fitting the experimental results with the equation (2):

Table 1. Data of the computed problem.


Romain Garby
email: romain.garby@imft.fr

Laurent Selle
email: laurent.selle@imft.fr

Thierry Poinsot
email: thierry.poinsot@imft.fr

Combustion, flow and spray dynamics for aerospace propulsion Large-Eddy Simulation of combustion instabilities in a variable-length combustor Simulation aux grandes échelles des instabilités de combustion dans un brûleur de longueur variable

Keywords:

le calcul axi-symétrique reproduit qualitativement certaines caractéristiques de l'instabilité.

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

r é s u m é

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

Introduction

Most high-performance propulsion systems exhibit high-amplitude pressure oscillations at some stage in their design process. Unsteady heat release from combustion can couple with acoustic eigenmodes of the combustion chamber: this phenomenon is called combustion instability [START_REF] Clavin | Turbulent premixed flames and sound generation[END_REF][START_REF] Harrje | Liquid propellant rocket instability[END_REF][START_REF] Culick | Unsteady Motions in Combustion Chambers for Propulsion Systems[END_REF]. These instabilities can lead to important stresses and heat fluxes on the chamber walls and injectors, sometimes even resulting in the destruction of the engine [START_REF] Searby | A parametric acoustic instability in premixed flames[END_REF][START_REF] Candel | Flame interactions as a source of noise and combustion instabilities[END_REF]. Combustion instabilities are very difficult to predict because of the wide range of physical phenomena involved (acoustics, turbulence, chemical reactions, etc.) and they are often discovered at an advanced phase of the conception program.

A trial-and-error methodology with incremental modifications of the engine is usually necessary to recover stability but the associated cost is prohibitive. Therefore, the prediction of combustion instabilities through numerical simulation is a major issue. For gas turbine or aeronautic propellers, several studies have already been conducted [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] Lieuwen | Combustion Instabilities in Gas Turbine Engines: Operational Experience, Fundamental Mechanisms and Modeling[END_REF][START_REF] You | A generalized model of acoustic response of turbulent premixed flame and its application to gas-turbine combustion instability analysis[END_REF][START_REF] Pitsch | Large eddy simulation of turbulent combustion[END_REF][START_REF] Schildmacher | Unsteady flame and flow field interaction of a premixed model gas turbine burner[END_REF][START_REF] Noiray | Investigation of azimuthal staging concepts in annular gas turbines[END_REF] but for rocket engines the number of unsteady simulations is smaller [START_REF] Menon | Large-eddy simulations of combustion instability in an axisymmetric ramjet combustor[END_REF][START_REF] Venkateswaran | Computational fluid dynamic analysis of liquid rocket combustion instability[END_REF][START_REF] Oefelein | Modeling high-pressure mixing and combustion processes in liquid rocket engines[END_REF][START_REF] Oefelein | Mixing and combustion of cryogenic oxygen-hydrogen shear-coaxial jet flames at supercritical pressure[END_REF][START_REF] Zong | Cryogenic fluid jets and mixing layers in transcritical and supercritical environments[END_REF][START_REF] Schmitt | Large-Eddy Simulation of oxygen/methane flames under transcritical conditions[END_REF]. This could be explained by the extreme thermodynamic conditions in rocket engines and the wide range of length scales: the characteristic length of the injector is a fraction of a millimeter while the flame length can be several dozens of centimeters.

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

Experimental setup

Description

The burner studied in this paper was designed and operated at Purdue University (Indiana, USA) [START_REF] Yu | Experimental study of high-frequency combustion instability in a continuously variable resonance combustor[END_REF][START_REF] Sisco | Examination of mode shapes in an unstable model combustor[END_REF][START_REF] Yu | Spontaneous longitudinal combustion instability in a continuously-variable resonance combustor[END_REF]. The system is called Continuously Variable Resonance Chamber (CVRC) because the length of the oxidizer injector can be varied continuously, allowing for a detailed investigation of the coupling between the acoustics and combustion in the chamber. The CVRC, presented in Fig. 1, is an axi-symmetric burner with a 15-inch long (38.1 cm) combustion chamber closed by a choked nozzle. The reactants are fed into the chamber by a coaxial injector: oxidizer at the center and fuel at the periphery.

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

Operating point

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

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

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

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

Stability

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

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

Numerical method

Solver and models

The simulations are carried out with the AVBP code developed at CERFACS [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] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] that solves the compressible Navier-Stokes equations on unstructured grids. The numerical scheme is the second-order in space and time Lax-Wendroff formulation [START_REF] Lax | Systems of conservation laws[END_REF]. Turbulence is modeled using the classical Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: 1. The basic experiment[END_REF]. The flame/turbulence interaction is accounted for by the Dynamically Thickened Flame (DTF) model [START_REF] Légier | Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion[END_REF]. In this approach, the flame is artificially thickened to be resolved on the LES mesh and subgrid reaction rate is recovered through an efficiency function [START_REF] Charlette | A power-law wrinkling model for LES of premixed turbulent combustion: Part I -Non-dynamic formulation and initial tests[END_REF].

The CH 4 /O 2 kinetics are modeled using a reduced two-step mechanism (BFER [START_REF] Franzelli | A two-step chemical scheme for Large-Eddy Simulation of kerosene-air flames[END_REF]) that accounts for six species (CH 4 , O 2 , CO 2 , CO, H 2 O and N 2 ) and two reactions:

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

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

Computational domain and boundary conditions

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

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

Results and discussion

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

Mean flow field

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

Flame stabilization and dynamics

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

Analysis of the unstable acoustic mode

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

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

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

1 τ t+τ t V γ -1 γ p p ′ ω′ dv dt and F = 1 τ t+τ t Σ p ′ u ′ .n dσ dt ( 6
)
where V is the volume of the domain and Σ its boundary; ω′ is the fluctuation of heat release rate and n the outward normal of Σ .

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

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

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

Conclusion

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

Fig. 1 .

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

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

Fig. 2 .

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

Fig. 3 .

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

Fig. 4 .

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

Fig. 5 .

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

Fig. 6 .

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

Fig. 7 .

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

Fig. 8 .

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

Fig. 9 .

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

Fig. 10 .

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

Table 1

 1 Operating point.

	Experiment	LES

c ≈ 1904 K 2512 K Chamber static pressure P c 1.34 MPa 1.54 MPa

Table 2

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

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

< 10 -4

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

Table 3

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

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

Acknowledgement

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

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