Thanks Lukas

Anaïs, Karim, Xavier, Tim Romain Houssam

Marie Wafa

Matteo Fernando

Thibaud

Rudy Benjamin

Marie Christine Tristani

white Cette thèse a été réalisée grâce à l'aide, au soutien et la présence de nombreuses personnes.

Je remercie les professeurs

Azeddine Kourta et Christophe Airiau, mes directeurs de thèse, pour m'avoir accueilli au sein du groupe EMT2, à l'IMFT et pour m'avoir donné l'opportunité de participer à cette expérience internationale. Je leur suis reconnaissant pour avoir assuré la direction de mes travaux et pour m'avoir fait partager leur expérience dans la recherche avec enthousiasme, patience et motivation, ainsi que pour l'aide à la rédaction de la thèse, et des résumés en français en particulier. Je remercie à nouveau Christophe Airiau pour son encadrement pendant ces années ainsi que pour sa confiance et sa bonne humeur. Il a pris le temps de relire et corriger ce manuscrit. Je remercie Thierry Poinsot pour avoir accepté de présider mon jury de thèse. Je remercie Jean-Christophe Robinet et Christian Tenaud pour avoir évalué mes travaux de thèse en tant que rapporteurs, et Aloïs Senginssen pour avoir accepté de prendre part à mon jury. Merci à tous pour vos observations pendant la soutenance.

"Materialism and morality have an inverse relationship. When one increases, the other decreases".

Introduction

From the day one of the human life on earth, man started exploring surroundings to understand its activities to adapt his way of life just to live. With the time, his thinking has evolved to develop and modify things to handle the difficulties around him. He started to develop his abilities to protect himself from other humans, rain, snow, sun or fire. With the increase of population, he found ways to share water, food and shelter with his community. With the discovery of wheel and his skill of taming animals around him, he started exploring the land which was spread before him . His ability improved with time and led him to cross water bodies with wooden logs, then with boats and ships. He discovered fossil fuels and invented ways to use them to burn, to produce energy and to cater his needs from cooking fast in homes to move fast on vehicles on rails or on roads. After making a lot of trials (sometimes fatal) while exploring the space above him, he found a way to fly heavier body faster than birds. Now, the man with his fast paced community concentrates on safety and comfortable journey on road, rail, on or under water, in the space above him or the space outside his planet. Now he realises that he has more responsibilities on surroundings while making his life, journey safe and comfortable. He creates lots of institutions on various disciplines to observe, study and analyse the eco-cycle and to recommend the community about the perturbations in the ecosystem. He recently found that his unoptimised or careless usage of natural resources led to huge concern on the perturbations on the environment. Now he communicates with other communities to decrease the exploitation of the natural sources without compromising the quality of life. Governments on the world started initiating lots of projects for reducing the pollution in any form. Every project got a collection of experts from science, engineering and technology to perform research to find every possible solution to reduce the pollution. More concern and restrictions are laid before the (terrestrial, airborne and nautical) vehicle industry to reduce the carbon emission.

Members of the European Commission initiated numerous projects with environmental concern in their mind for the betterment of climate and the quality of life for the present and future generations to follow. The PhD work presented here is a part of a project whose aim is to reduce air and noise pollution which are generated from civilian aircrafts, vehicles on rails and roads.

Introduction

Figure 1.1: An aircraft with landing wells during take off. 1

Noise

The noise, defined as unwanted, excessive, uncomfortable sound, is a major problem in day to day life. Researchers have known for years that exposure to excessively-loud noise can cause changes in blood pressure as well as changes in sleep and digestive patterns -all signs of stress on the human body. The very word "noise"itself derives from the Latin word noxia, which means injury or hurt.

In 1996, the European Commission published the Green paper [START_REF]Future Noise Policy[END_REF] which showed that about 20% of the population in the European Union live in the so called 'grey areas'

where the noise exposure exceeds an equivalent noise level of 65 dB at daytime. The same document discuss about variety of topics related to noise pollution such as exposure of population to the noise level in their surroundings. European Union estimates the external cost of noise pollution vary between 0.2% and 20% of Gross Domestic Product.

It mentions about the noise pollution from the vehicular transport apart from industrial noise pollution. The European Commission considers living close to an airport to be a possible risk factor for coronary heart disease and stroke, as increased blood pressure from noise pollution can trigger more serious maladies and European Union terms the neighbourhood as unhealthy and unacceptable place to live.

The most common noise sources can be divided into aerodynamic and mechanical.

There are various noise -generating elements on aircrafts which includes the engines, the engine housing and airframe. Though sound produced by the engine is high. But during take off and landing, the sound generated by the airframe components: landing gears, high-lift devices, fuel vents are dominant. In general, noise control is an active or passive means of reducing sound emissions, often incentivized by personal comfort, environmental considerations or legal compliance. Effective noise control focuses on reducing the noise from these sources as near of the source as possible. Locating the source of sound and reducing the intensity or loudness of the sound can be done only after a series of research by moving the engineers from industries and research scholars from universities to work together.

tional time-averaged or statistical approaches and introducing time-dependent methods for aeronautical research at an early stage of career development. One common research topic was chosen by these four institutes to solve the unsteady flow over airflow fuel vents (i.e cavities of rectangular and cylindrical shapes) by different investigating methods such as experiments and/or numerical simulation.

In Institut de Mécanique des Fluides de Toulouse, usage of numerical simulation is proposed as the method of choice to identify the flow events that are related to the dynamics of coherent structures and are the main acoustic noise sources in the cavity flow.

Motivation and Objectives

Cavity flows are studied for practical purposes such as reduction of drag, energy consumption and unnecessary noise. Cavities represents the landing gears, fuel vents of air borne vehicles and in the terrestial vehicles they are found as windows of train coaches, space between wagons of trains, sun roof and windows of cars. These unavoidable cavities generates noise inside the cabin space giving high human discomfort and also disturbs the ecological system which involves humans, birds and animals who live near the roadways, railways and runways. Figure 1.2 shows a civilian aircraft with landing gears engaged.

Not only by the noise pollution, the ecological system is also affected by the air pollution. Carbon emission from the engine of the vehicles is increased due to more 1. Introduction consumption of fossil fuel with increase of drag on the vehicle with cavities.

Apart from noise and air pollution, the structure and components of the vehicles fail without warning due to the fatigue. The load or drag weakens the mechanical properties of the material. This increases the weight and volume of the vehicle in whole.

Cavity flows contains a wide range of physical phenomenon like unsteady shear layer, vortex shedding, recirculation eddies, instabilities and three dimensional effects. The main objective of this PhD work is to determine the sound sources at near and far field of the cavity with an incoming thick turbulent boundary layer and to investigate the turbulence of the cavity flow using numerical methods. An acoustic analogy is coupled with Large Eddy Simulation. Large Eddy Simulation is performed on a significant number of test cases to obtain hydrodynamic pressure values in the domain. Sound Pressure Level is obtained using acoustic analogy from the hydrodynamic pressure values determined by Large Eddy Simulation method.

Plan of the thesis

The contents of this thesis are organised as follows:

• Chapter 2 : This chapter is devoted to discuss about the literature related to cavity flows, turbulence, Direct Numerical Simulation and Large Eddy Simulation and about Aeroacoustics which includes acoustic analogy and the procedures that are followed to determine the sound pressure level of the noise generated by the cavity flow.

• Chapter 3 : This chapter starts with the description about the inflow condition, asymptotic modeling and contains sections to explain the modeling of turbulent boundary layer, mixing length model, zero pressure gradient boundary layer and adverse pressure gradient boundary layer.

• Chapter 4 : provides a general overview of main features of the AVBP solver.

The chapter begins with the description of the governing equations of Large Eddy

Simulation and the method used to discretisation of the governing equations. Description of boundary conditions is also included. The test cases, geometries, meshing, challenges while performing simulations are also discussed in this chapter. Results are summarised and analysis on boundary layer turbulence in and downstream of the cavity are included here. At the end of this chapter results obtained from acoustic analogy are presented and analysed.

• Chapter 5 : In this final chapter of the thesis, observations and conclusions are laid. Future extension of this PhD work and perspectives are listed at the end.

Cavity flow, turbulence and aeroacoustics

Plusieurs études numériques se sont intéressées à l'aéroacoustique de la cavité 2D ou 3D: simulation numérique directe (Gloerfelt et al [START_REF] Gloerfelt | Aerodynamic noise induced by laminar and turbulent boundary layers over rectangular cavities[END_REF]), modèles de Boltzmann sur réseau (Ricot et al [START_REF] Ricot | Numerical simulation of unsteady cavity flow using Lattice Boltzmann Method[END_REF]), simulations LES (Larchevêque et al [START_REF] Larchevêque | Largeeddy simulation of a compressible flow past a deep cavity[END_REF]).... Ainsi on a pu constater que la largeur de la cavité (dans la direction de l'envergure) modifie les oscillations de pression (effet 3D), que les niveaux sonores les plus importants sont obtenus pour les cavités les plus larges et que l'épaissseur de la couche limite amont pilote la croissance des oscillations.

Cavité-Oscillations de l'écoulement associées

Les oscillations induites par une cavité peuvent être classées en trois catégories (figure 2.2):

1. Fluide-élastique: elles apparaissent quand les modes propres de la surface de la cavité sont forcés (élasticité de la paroi)

2. Fluide-résonant: il existe une oscillation auto-entretenue a une longueur d'onde équivalente aux dimensions de la cavité. Il y a couplage entre les modes acoustiques de la cavité et la couche cisaillée au-dessus de la cavité.

3. Fluide-dynamique: elles sont liées à un mécanisme de feedback. Ce régime implique l' amplification des instabilités de la couche cisaillée provoquées par le retour de l'interaction de la couche cisaillée avec le coin ou la paroi aval.

En raison de

Ecoulement de cavité bidimensionnel et tridimensionnel

Il est généralement convenu de considérer l'écoulement de cavité comme essentiellement bidimensionnel. Cependant, on sait que les tourbillons longitudinaux au sein de la couche limite, des effets de bord sur la couche de cisaillement et dans la cavité, ainsi que des instabilités de type Taylor-Görtler dues à la forte courbure de la recirculation peuvent induire une tridimensionnalisation de l'écoulement. Des études expérimentales (Rockwell et Knisley [START_REF] Rockwell | Observations of the three-dimensional nature of unstable flow past a cavity[END_REF]) et numériques (Rizzetta and Visball [START_REF] Rizzetta | Large-eddy simulation of supersonic cavity flowfields including flow control[END_REF], Larchevêque et al [START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a three-dimensional open cavity at high Reynolds number[END_REF],

Chang et al [START_REF] Chang | Analysis of the flow and mass transfer processes for the incompressible flow past an open cavity with a laminar and a fully turbulent incoming boundary layer[END_REF]) ont analysé l'aspect tridimensionnel à l'intérieur de la cavité.

Cavity flow, turbulence and aeroacoustics

Ecoulement de cavité cylindrique

La cavité cylindrique a également été étudiée. L'expérience de Hiwada et al [START_REF] Hiwada | Some characteristics of flow pattern and heat transfer pas a circular cylindrical cavity[END_REF] montre l'évolution de l'écoulement en faisant varier le rapport d'aspect de 0,1-1 (voir tableau 2.1). Rona [START_REF] Rona | The acoustic resonance of rectangular and cylindrical cavities[END_REF] Une revue détaillée des calculs aéroacoustiques peut être trouvée dans Larsson [START_REF] Larsson | Computational Aero Acoustics for Vehicle Applications[END_REF] et Tam [START_REF] Tam | Computational Aeroacoustics: An overview of computational challenges and applications[END_REF].

2. Cavity flow, turbulence and aeroacoustics

Introduction

This chapter is devoted to discuss the following topics elaborately: Cavity flows, direct numerical simulation, Navier-Stokes equations, Turbulence, RANS and Aeroacoustics.

Earlier and recent studies on two-dimensional and three-dimensional cavities related to aspect ratio, Mach number and other parameters are reviewed along with the physical phenomenon occuring in the cavity flows under different conditions. Navier-Stokes equation for direct numerical simulation are given elaborately. Turbulence and its quantities are presented. Finally under computational aeroacoustics, Lighthill-Curle's analogy is discussed and the equation required to determine sound pressure level is also derived.

Cavity flows 2.2.1 Physical phenomenon, Resonance

Karbon [START_REF] Karbon | Simulation and design of automobile sunroof buffeting noise control[END_REF] observes that when a vehicle moving at 50km/h, with the sunroof open, the noise in the cabin space can reach more than sound pressure level of 98dB which will bring stress and fatigue to the travellers. Lid-driven cavity flow does not considers the interaction between the shear layer and the recirculating flow but just models the flow field inside the cavity [START_REF] Sinha | Laminar separating flow over backsteps and cavities Part II: Cavities[END_REF]. showing the incoming boundary layer at the leading edge of the cavity, shear layer over the cavity and the pressure perturbation from the trailing edge of the cavity due to the impingement of the shear layer on the corner of the downstream of the cavity. Eddy or eddies are created inside the cavity depending on various parameters which will be dicussed inside the chapter.

Karamcheti [START_REF] Karamcheti | Sound radiation from surface cutouts in high speed flow[END_REF] reported that there is a minimum cavity length needed for generation of cavity noise, depending on the Mach number of the flow and whether the approaching boundary layer is turbulent. If the cavity length is less than the minimum length, the flows will not oscillate. Sarohia [START_REF] Sarohia | Experimental investigation of oscillations in flows over shallow cavities[END_REF] stated that the parameters cavity depth D and initial momentum thickness θ 0 at the leading edge also are as important as the minimal cavity length L. Gloerfelt et al [START_REF] Gloerfelt | Aerodynamic noise induced by laminar and turbulent boundary layers over rectangular cavities[END_REF]

Cavity-related flow oscillations

The understanding of cavity-related flow oscillations was simplified by Rossiter and Naudascher. They divides them into three categories 1. Fluid-elastic oscillations: They occur when a cavity surface itself is forced into oscillation. In other words, this regime encompasses flows that are affected by the elastic boundaries of the cavity.

2. Fluid-resonant oscillations: These are caused when a self sustaining oscillation in the flow has a wavelength of the same order as one of the cavity dimensions. This regime couples the acoustic modes of the cavity and shear layer over the deeper cavities and for the cavities subject to high Mach number flow.

3. Fluid-dynamic oscillations: These are related to the cavity feedback resonance mechanism. This regime involves shear-layer instability amplification due to feedback from interaction of the shear layer. These interactions occurs for low-speed flow past shallow cavities.

Due to the self-sustaining nature of the feedback mechanism, acoustic pulses are generated periodically and a narrow band acoustic tone results. A semi-empirical formula to predict the frequency of this tone was predicted by Rossiter:

St n = f n L U = n -γ M + 1 κ n = 1, 2, ..., (2.1)
where St n is the Strouhal number corresponding to the n th mode frequency f n , and κ = 1 1. [START_REF] Kraichnan | Eddy viscosity in two and three dimensions[END_REF] and γ = 0.25 are empirical constants corresponding to the average convection speed of disturbances in the shear layer, and a phase delay.

Classification and main results

Open and closed cavities

Earlier, according to Sarohia [START_REF] Sarohia | Experimental investigation of oscillations in flows over shallow cavities[END_REF] shallow cavities have aspect ratios L D less than unity whereas deep cavities have L D ratios greater than unity. Rossiter [START_REF] Rossiter | Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF] defines the cutoff to be a ratio of 4.0. observations. Karamcheti noticed that, for a fixed freestream Mach number M ∞ and depth D, there exists a minimum cavity length L min below which no sound emission is noticed. For a fixed cavity, experimental results further showed a minimum Mach number below which no sound emission was noticed. For a given flow, the prerequisite of a minimum length L min for the onset of cavity oscillations strongly suggests that the mechanism of cavity oscillations depends upon the stability characteristics of the shear layer. Rockwell [START_REF] Rockwell | Prediction of oscillation frequencies for unstable flow past cavities[END_REF] and Rockwell and Naudascher [START_REF] Rockwell | Review-self-sustaining oscillations of flow past cavities[END_REF] clarified the significant parameters for this oscillation type as Re, δ 2 , θ 0 /L, L/W . Rockwell and Naudasher [START_REF] Rockwell | Review-self-sustaining oscillations of flow past cavities[END_REF] predicted the main oscillatory frequency for incompressible flow over two-dimensional cavities based on linear inviscid stability theory. The predictions agreed well with the [START_REF] Srinivasan | Navier-stokes calculations of transonic flow past cavities[END_REF]. Recently, Tracy and Plentovich [START_REF] Tracy | Measurements of fluctuating pressure in a rectangular cavity in transonic flow at high Reynolds numbers[END_REF] and Raman et al [START_REF] Raman | Tone noise and near field pressure produced by jet-cavity interaction[END_REF] have concluded that the disagreement found in the literature stems from the dependence of the cavity flow type on Mach number as well as L D . It was shown that the boundary layer thickness at the cavity lip is also an important parameter [START_REF] Ahuja | Effects of cavity dimensions, boundary layer and temperature on cavity noise with emphasis on benchmark data to validate computational aeroaoucstic codes[END_REF], [START_REF] Tracy | Measurements of fluctuating pressure in a rectangular cavity in transonic flow at high Reynolds numbers[END_REF]. Colonius [START_REF] Colonius | Computational aeroacoustics: progress on nonlinear problems of sound generation[END_REF] states that the momentum thickness θ 0 at the leading edge of the cavity plays a vital role in the selection of the modes and in governing the growth of the shear layer [START_REF] Rowley | On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[END_REF], [START_REF] Tam | Computational Aeroacoustics: An overview of computational challenges and applications[END_REF] that spans an open cavity [START_REF] Charwat | An investigation of separated flows. Part I: The pressure field[END_REF]. Gharib and Roskho [START_REF] Gharib | The effect of flow oscillations on cavity drag[END_REF] specified the threshold for self sustained oscillation and the wake mode. They also found L θ for lower limit for the cavity resonance to be approximately

Shear and wake mode

For the low Mach number flows, the cavity has been classified as shear mode or wake mode according to the shear-layer on the cavity. In a shear layer mode, the length of the cavity plays an important role. According to Rossiter [START_REF] Rossiter | Wind-tunnel experiments on the flow over rectangular cavities at subsonic and transonic speeds[END_REF], the oscillations generated are driven by the vortices from the shear layer. The wave length of this periodic oscillation is usually close to the cavity length or 1 N of the cavity length. The oscillation of the shear layer is confined within a narrow region near the straight line between the leading and trailing edge of the cavity. The recirculation flow inside the 2.3. Classification and main results cavity is usually relatively quiescent and the interaction between the shear layer and flow inside the cavity is weak. In a cavity with shear-layer mode, the shear layer spans the mouth of the cavity and stagnates at the downstream wall. Both fluid-resonant and fluid-dynamic regimes can be found in the cavity with shear-layer mode. When the shear layer oscillates in the shear layer mode, multiple discrete and high magnitude peaks will be present in the pressure spectra. These peaks are the cavity tones. There is usually one tone with higher magnitude than the rest of the spectrum as it so that it possesses most of the energy. This tone is referred to as the dominant tone or the fundamental frequency. Karamcheti [START_REF] Karamcheti | Sound radiation from surface cutouts in high speed flow[END_REF] discovered that the frequency of the dominant tone is inversely proportional to the cavity length. As the length of cavity becomes even longer, the fundamental frequency disappears and strong intermittencies will overcome the coherent oscillation. The feedback mechanism becomes ineffective at this point. This mode of oscillation is called a wake mode. This mode is identified by the stagnation of the flow prior to the downstream wall. Gharib and Roshko [START_REF] Gharib | The effect of flow oscillations on cavity drag[END_REF] noted the flow looked similar to a bluff-body wake and named the mode as wake mode. In the wake mode, self oscillations cease, the cavity flow becomes unstable on a large scale, and the drag increase with the presence of the cavity. The depth of the cavity becomes more important in this type of mode. Direct numerical simulations by Rowley et al [START_REF] Rowley | On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[END_REF] showed similar results for a two-dimensional rectangular cavity. In this mode, the vortex grows near the leading edge of the cavity until it fills the cavity, then it sheds downstream, collides onto the rear wall, and ejects out of the cavity. The region of the shear layer oscillation is much larger, up to the depth of the cavity. Three-dimensionality has been shown to play a role in suppressing the wake mode. Wake mode is less likely to appear in threedimensional flows and at higher Reynolds numbers, for example Rowley et al [START_REF] Rowley | On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[END_REF]. Large eddy simulations by Shieh and Morris [START_REF] Shieh | Parallel computational aeroacoustic simulation of turbulent subsonic cavity flow[END_REF] showed that two-dimensional cavities in wake mode return to shear-layer mode when three-dimensional disturbances are present in the incoming boundary layer. Suponitsky et al [START_REF] Suponitsky | On three-dimensionality and control of incompressible cavity flow[END_REF] showed that the development of a threedimensional flow field, generated by the introduction of the random in flow disturbance into a two-dimensional cavity oscillating in wake mode, yielded the transition to the shear-layer mode, regardless of the amplitude and shape of the inflow disturbance. [START_REF] Hiwada | Some characteristics of flow pattern and heat transfer pas a circular cylindrical cavity[END_REF] vided a transition between two-dimensional and three-dimensional flow. They reported three-dimensionality in the mean flow, and much lower (about 15dB) acoustic loads than the two-dimensional flow. The three-dimensional cavity flow have been studied

Two dimensional and three dimensional cavity flow

using Large Eddy Simulation approach by Rizzetta and Visbal [START_REF] Rizzetta | Large-eddy simulation of supersonic cavity flowfields including flow control[END_REF], Larchevêque et al [START_REF] Larchevêque | Large-eddy simulation of a compressible flow in a three-dimensional open cavity at high Reynolds number[END_REF] and Chang et al [START_REF] Chang | Analysis of the flow and mass transfer processes for the incompressible flow past an open cavity with a laminar and a fully turbulent incoming boundary layer[END_REF]. These studies have been mainly focused on the frequencies of oscillation and coherence of the Rossiter modes. The three-dimensional incompressible LES, coupled with the Lighthill-Curle acoustic analogy is used by Suponitsky et al [START_REF] Suponitsky | On three-dimensionality and control of incompressible cavity flow[END_REF], to investigate the oscillation mechanism and sound source of a two-dimensional cavity with a length-to-depth ratio of L D = 4 and Reynolds number of Re D = 5000. At the inflow boundary a streamwise velocity profile is specified as a power law of

u u ∞ = y δ 1 7
Simulations without and with inflow disturbance are carried out. More evidence of three-dimensional structures in cavity flows have been presented in the work of Faure et al [START_REF] Faure | Visualizations of the flow inside an open cavity at medium range Reynolds numbers[END_REF]. They investigated experimentally the interaction between a laminar boundary layer and an open cavity L D = 0.5 -2 for medium range Reynolds numbers. In their work, they relate the three-dimensional structures to the primary vortex inside the cavity.

High Mach number cylindrical cavity flow

Cylindrial cavity flows are more complex than the rectangular cavities. Hiwada et al [START_REF] Hiwada | Some characteristics of flow pattern and heat transfer pas a circular cylindrical cavity[END_REF] performed experiments on cylindrical cavities with ratio cavity depth / cavity diameter 0.1 to 1.0 and the observations can be found in the table 2.1. Dybenko et al [START_REF] Dybenko | Turbulent flow over circular cylindrical cavities with varying depth to diameter ratio[END_REF]observed that the symmetric flow is related to the occurence of an acoustic feedback mechanism.

Rona [START_REF] Rona | The acoustic resonance of rectangular and cylindrical cavities[END_REF] developed an analytical model to investigate oscillations in circular cavities and he predicts the asymmetric modes being oriented in one or the other direction. Recent experiments on cylindrical cavities are done by Marsden et al [START_REF] Marsden | Investigation of flow features and acoustic radiation of a round cavity[END_REF] [START_REF] Marsden | Investigation of flow features and acoustic radiation of a round cavity[END_REF].

ically to approach experimental results. Preliminary numerical results are presented in [START_REF] Marsden | Investigation of flow features and acoustic radiation of a round cavity[END_REF] for the flow configuration of 90 m/s. An Euler numerical method was used by Grottadaurea and Rona [57]

Direct Numerical Simulation

It becomes important to discuss about Direct Numerical Simulation (DNS) as few test cases were perfomed with this numerical approach. Computing power in the recent because the unsteady eddy motions of turbulence appear over a wide range.

Navier Stokes Equations

Conservative form

The Navier stokes equations for the direct numerical simulation which are presented in this section are followed in the solver AVBP. More details about AVBP is given in chapter 4. Conservation equations describing the evolution of a compressible flow with chemical reactions of thermodynamically active scalars reads,

∂ρu i ∂t + ∂ ∂x j (ρu i u j ) = - ∂ ∂x j [P δ ij -τ ij ] (2.2) ∂ρE ∂t + ∂ ∂x j (ρEu j ) = - ∂ ∂x j [u i (P δ ij -τ ij ) + q j ] + ωT + Q r (2.3) ∂ρ k ∂t + ∂ ∂x j (ρ k u j ) = - ∂ ∂x j [J j,k ] + ωk (2.4)
It should be noted that index k is reserved to refer to the k th species and will not follow the summation rule unless other specified or implied by the sign.

In Eqs 2.4-2.4 respectively corresponding to the conservation laws for momentum, total energy and species Y , the following symbols denote respecticely ρ, u i , E, ρ k , density, the velocity vector, the total energy per unit mass and ρ k = ρY k for k = 1 to N (N is the total number of species Y ). The source term in the total energy equation 2.4, is decomposed for convenience into a chemical source term and a radiative source term such that: S = ωT + Q r . Corresponding source terms in the species transport equations 2.4 are ωk . It is usual to decompose the flux tensor into an inviscid and a viscous component.

They are respectively noted for the three conservation equations:

Inviscid terms

The inviscid terms from above equations are grouped in matrix form as

    ρu i u j + P δ ij (ρE + P δ ij ) u j ρ k u j     (2.5)
where the hydrostatic pressure P is given by the equation of state for a perfect gas (Eq. 2.13).

Navier Stokes Equations

Viscous terms

Similarly, the components of the viscous flux tensor take the form:

    -τ ij -(u i τ ij ) + q j J j,k     (2.6)
where J k is the diffusive flux of species k and is presented in section 2.5.4 (equation 2.24).

The stress tensor τ ij is explained in section 2.5.5 (equation 2.26). The section 2.5.6 is devoted to the heat flux vector q j (equation 2.28).

Thermodynamical variables

The standard reference state used is P 0 = 1 bar and T 0 = 0K. In AVBP solver (for details see chapter 4), the sensible mass enthalpies (h s,k ) and entropies (s k ) for each species are tabulated for 51 values of the temperature (T i with i = 1...51) ranging from 0K to 5000K with a step of 100K. Therefore these variables can be evaluated by:

h s,k (T i ) = T i T 0 =0K C p,k dT = h m s,k (T i ) -h m s,k (T 0 ) W k
, and (2.7)

s k (T i ) = s m k (T i ) -s m k (T 0 ) W k , with i = 1...51 (2.8)
The superscript m corresponds to molar values. The tabulated values for h s,k , (T i ) and s k (T i ) can be found in the JANAF tables [START_REF] Stull | JANAF thermochemical tables[END_REF]. W k is molecular weight of the species Y k With this assumption, the sensible energy for each species can be reconstructed using the following expression :

e s,k (T i ) = T i T 0 =0K C v,k dT = h s,k (T i ) -r k T i i = 1...
e s,k (T ) = e s,k (T i ) + (T -T i ) e s,k (T i+1 ) -e s,k (T i ) T i+1 -T i (2.10)

Cavity flow, turbulence and aeroacoustics

The sensible energy and enthalpy of the mixture may then be expressed as:

ρe s = N k=1 ρ k e s,k = ρ N k=1 Y k e s,k (2.11
)

ρh s = N k=1 ρ k h s,k = ρ N k=1 Y k h s,k (2.12
)

The equation of state

The equation of state for an ideal gas mixture is given as:

P = ρ r T (2.13)
where r is the gas constant of the mixture dependant on time and space: r = R W where W is the mean molecular weight of the mixture:

1 W = N k=1 Y k W k (2.14)
The gas constant r and the heat capacities of the gas mixture depend on the local gas composition as:

r = R W = N k=1 Y k W k R = N k=1 Y k r k (2.15) C p = N k=1 Y k C p,k (2.16
)

C v = N k=1 Y k C v,k (2.17)
where R = 8.3143 J/mol.K is the universal gas constant. The adiabatic exponent for the mixture is given by γ = C p C v . Thus, the gas constant, the heat capacities and the adiabatic exponent are no longer constant. Indeed, they depend on the local gas composition as expressed by the local mass fractions Y k (x, t):

r = r(x, t), C p = C p (x, t), C v = C v (x, t), and γ = γ(x, t) (2.18)
The temperature is deduced from the the sensible energy, using equations 2.10 and 2.11.

finally boundary conditions make use of the speed of sound of the mixture a ∞ which is given by:

a 2 ∞ = γrT (2.19)
2.5. Navier Stokes Equations

Conversation of Mass: Species diffusion flux

In multi-species flows the total mass conservation implies that:

N k=1 Y k V k i = 0 (2.20)
where V k i are the components in directions (i = 1, 2, 3) of the diffusion velocity of species k. They are often expressed as a function of the species gradients using the Hirschfelder Curtis approximation:

X k V K i = -D k ∂X k ∂x i , (2.21)
where X k is the molar fraction of species k :

X k = Y k W W k .
In terms of mass fraction, the approximation 2.21 may be expressed as:

Y k V k i = -D k W k W ∂X k ∂x i , (2.22)
Summing equation 2.22 over all k's shows that the approximation 2.22 does not necessarily comply with equation 2.20 that expresses mass conservation. In order to achieve this, a correction diffusion velocity V c is added to the convection velocity to ensure global mass conservation (see [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]) as:

V c i = N k=1 D k W k W ∂X k ∂x i (2.23)
and computing the diffusive species flux for each species k as:

J i,k = -ρ D k W k W ∂X k ∂x i -Y k V c i (2.24)
Here, D k are the diffusion coefficients for each species k in the mixture (see 2.5.7); J i,k is computed. Using equation 2.24 to determine the diffusive species flux implicitly verifies equation 2.20.

Viscous stress tensor

The stress tensor τ ij is computed and is given by the following relations: where µ is the shear viscosity (see 2.5.7).

τ ij = 2µ S ij - 1 3 δ ij S

Heat flux vector

For multi-species flows, an additional heat flux term appears in the diffusive heat flux.

This term is due to heat transport by species diffusion. The total heat flux vector then writes:

q i = -λ ∂T ∂x i Heat conduction -ρ N k=1 D k W k W ∂X k ∂x i -Y k V c i h s,k Heat flux through species diffusion = -λ ∂T ∂x i + N k=1 J i,k h s,k (2.28)
where λ is the heat conduction coefficient of the mixture (see 2.5.7). The second term is added to the classical heat flux vector.

Transport coefficients

In CFD codes for multi-species flows the molecular viscosity µ is often assumed to be independent of the gas composition and close to that of air so that the classical Sutherland law can be used. Same assumption is proposed for the multi-gas AVBP(see chapter 4), yielding:

µ = c 1 T 3/2 T + c 2 T ref + c 2 T 3/2 ref (2.29)
where c 1 and c 2 must be determined so as to fit the real viscosity of the mixture. For air at T ref = 273K, c 1 = 1.71 × 10 -5 kg/m.s and c 2 = 110.4K (see [START_REF] Stull | JANAF thermochemical tables[END_REF]). A second law is available, called Power law:

µ = c 1 T T ref b (2.30)
with b typically ranging between 0.5 and 1.0. For example b = 0.76 for air.

The heat conduction coefficient of the gas mixture can then be computed by intro-2.6. Turbulence ducing the molecular Prandtl number of the mixture as:

λ = µC p P r (2.31)
with P r supposed as constant in time and space.

The computation of the species diffusion coefficients D k is a specific issue. These coefficients should be expressed as a function of the binary coefficients D ij obtained from kinetic theory (Hirschfelder et al. [START_REF] Hirschfelder | Molecular Theory of Gases and Liquids[END_REF]). The mixture diffusion coefficient for species k, D k , is computed as (Bird et al. [START_REF] Bird | Transport Phenomena[END_REF]):

D k = 1 -Y k N j =k X j D jk (2.32)
The D ij are complex functions of collision integrals and thermodynamic variables. For a DNS code using complex chemistry, using Eq 2.32 makes sense. However in most cases, DNS uses a simplified chemical scheme and modeling diffusivity in a precise way is not needed so that this approach is much less attractive. Therefore a simplified approximation is used in AVBP for D k . The Schmidt numbers S c,k of the species are supposed to be constant so that the binary diffusion coefficient for each species is computed as:

D k = µ ρS c,k (2.33)
Note that the Schmidt number for each species k is assumed to be constant in time and space and is given as input parameter. P r and S c,k model the laminar (thermal and molecular) diffusion. Usual values of Schmidt and Prandtl numbers for premixed flames are those given by PREMIX in the burnt gas. The kinetics, radiation are not discussed here in the work.

Turbulence

Horace Lamb said:

"I am an old man now, and when I die and go to heaven there are two matters on which I hope for enlightenment. One is quantum electrodynamics, and the other is the turbulent motion of fluids. And about the former I am rather optimistic".

The earliest identification of turbulence as a prominent physical phenomenon had already taken place during the time of Leonardo da Vinci (circa 1500). But there seems to have been no significant progress in understanding until last part of 19 th century. The figure 2.4 is a rendition of one found in a sketch book of Leonardo da Vinci. Turbulence has been described as a random/irregular motion, both in time and space.

In 1877, Boussinesq [START_REF] Boussinesq | Essai sur la théorie des eaux courantes[END_REF] hypothesis that turbulent stresses are linearly proportional to mean strain rates is still the base of most turbulence models. Laminar, transitional and turbulent are the three regimes which can generally be observed in a flow field. The flows in the laminar regime are smooth, streamlined and the adjacent layers of fluid slide past each other in an orderly manner. The transition is due to the instability of the laminar state which makes the change from a laminar flow to a turbulent one. Most flows encountered in nature and in industrial applications are turbulent. Turbulence is characterised by disorganised motion over a range of length and time scales. The concept of an energy cascade wherein an equilibrium exists between drawing energy from mean flow gradients and dissipating that energy at the smallest scales allows one to estimate the range of scales for a given flow. However, due to its complexity, understanding of turbulent flows or turbulence is still incomplete. The understanding of turbulent behaviour in flowing fluids is one of the most intriguing, frustrating and important problems in all of classical physics.

Osborne Reynolds [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF] [START_REF] Chapman | Observations, theoretical ideas, and modeling of turbulent flows-past, present and future[END_REF], [START_REF] Mcdonough | Introductory lectures on Turbulence-Physics[END_REF].

nate over the diffusion forces. Turbulent flows dissipate energy. The kinetic energy of the fluid is converted into heat, at the smallest scales, due to viscous effects. Poincaré [START_REF] Poincaré | Les methodes nouvelles de la mécanique celeste[END_REF] found that relatively simple nonlinear dynamical systems were capabale of exhibiting chaotic random-in-appearance behaviour that was in fact, completely deterministic.

Lorenz [START_REF] Lorenz | Deterministic nonperiodic flow[END_REF] in 1963 was first to propose the connection between deterministic chaos and turbulence. Chapman and Tobak [START_REF] Chapman | Observations, theoretical ideas, and modeling of turbulent flows-past, present and future[END_REF] divide the century between Reynolds experiments in 1883 to the then present time 1984 into three overlapping "movements"that they term statistical, structural and deterministic. Figure 2.5 provides a sketch similar to the one presented in [START_REF] Chapman | Observations, theoretical ideas, and modeling of turbulent flows-past, present and future[END_REF]. McDonough [START_REF] Mcdonough | Introductory lectures on Turbulence-Physics[END_REF] has discussed more about statistical, structural and deterministic movement and included additional entries to the figure 2.5.

The first major result was obtained by Prandtl [115] in 1925 in the form of a prediction of the eddy viscosity (introduced by Boussinesq) that took the character of a "first-principles "physical result, and as such no doubt added significantly to the credence of the statistical approach. Prandtl's "mixing-length theory"was based on an analogy between turbulent eddies and molecules or atoms of a gas and purportedly utilized kinetic theory to determine the length and velocity (or time) scales needed to construct an eddy viscosity (analogous to the first-principles derivation of an analytical description of molecular viscosity obtained from the kinetic theory of gases). The smallest scales associated with turbulence are much larger than the molecular mean free path. Thus, turbulence is a continuum process. Other significant characteristic is that turbulence is three dimensional, an important role in setting up and maintaining the continuum of scales characteristic of a high Reynolds number turbulence being the vortex stretching.

In two space dimensions vortex stretching cannot occur. Turbulence is a diffusive process since it causes rapid mixing and increases the rates of mass, momentum, and heat transfer. It should be noted that it is not a property of the fluid, but it is a property of the flow. A turbulent flow field at high Reynolds number consists of vortices (eddies) of various sizes, from the largest to the smallest ones. Each eddy can be related to a scale of velocity, time and length. These initial large vortical scales will break up due to vortex stretching to develop smaller and smaller scale structures. Four main sets of scales in a turbulent flow (there may be more if other physical phenomena, e.g., heat transfer and/or combustion are important); these are:

1. the large scale, based on the problem domain geometry, 2. the integral scale, which is an O(1) fraction (often taken to be ∼ 0.2) of the large scale, 3. the Taylor microscale which is an intermediate scale, basically corresponding to Kolmogorov's inertial subrange, and 4. the Kolmogorov (or "dissipation") scale which is the smallest of turbulence scales.

Turbulence

Turbulence energy dissipation rate ǫ is given as

ǫ = 2ν S 2 (2.35)
where S is the strain rate tensor. The length and time scales are derived which are associated with the Taylor microscale. The definition for the Taylor microscale length provided in [START_REF] Lesieur | Turbulence in Fluids[END_REF]:

λ 2 = | u ′ | S (2.36)
| u ′ | is the square root of the turbulence kinetic energy k. The Taylor microscale length expressed as

λ = ν | u ′ | ǫ (2.37)
The Taylor microscale length is roughly consistent with the Kolmogorov inertial subrange scales. The smallest scales of turbulence were derived by Kolmogorov under the assumption that dissipation wouls be important at these scales. The two physical parameters needed to describe behavior are viscosity ν and dissipation rate ǫ of turbulence kinetic energy. The length scale given by Kolmogorov is

η k = ν 3 ǫ 1/4 (2.38)
and Kolmogorov time scale is

τ k = ν ǫ 1/2 (2.39)
It is of interest to compare some of these various scales. We observe the length scales and Reynolds numbers can be related as follows. First, we can compare the Kolmogorov length scale η k with the integral scale length ℓ using equation 2.38 with length scale Unsteady-RANS [START_REF] Singer | Development of computational aeroacoustics tools for airframe noise calculations[END_REF] and Detached Eddy Simulation [START_REF] Hedges | Detached-eddy simulation over a simplified landing gear[END_REF] are other time-dependent numerical predictions used in complex geometries from industries. Unclosed terms arise when introducing operators on the set of compressible Navier-Stokes equations and these terms from these manipulations and models need to be supplied for the problem to be solved. In RANS the operation consists of a temporal or ensemble average over a set of realisations of the studied flow (Pope [START_REF] Pope | Turbulent Flows[END_REF] and Chassing [START_REF] Chassaing | Turbulence en mécanique des fluides, analyse du phénomène en vue de sa modélisation à l'usage de l'ingénieur[END_REF]). The unclosed terms represents the physics taking place over the entire range of frequencies present in the ensemble of realisations under consideration. Figure 2.6 illustrates the range of wavenumber modeled or/and resolved by DNS, RNS and LES. Large eddy simulation (LES) is dicussed more in detail in the chapter 4.

ℓ = | u ′ | 3 /ǫ solved for ǫ to write η k ∼ ν 3 | u ′ | 3 /ℓ ; η k ℓ ∼ ν 3 | u ′ | 3 ℓ 3 ∼ Re -3/4 ℓ , ℓ η k ∼ Re 3/

Aeroacoustics

Acoustics is the science of sound that includes its production, propagation, and its effects. The acoustic theory from fluid mechanics focuses on the mathematical description of sound waves. The flows are governed by a nonlinear system of equations. This is the fact responsible for the complexity of fluid dynamics research and consequently for flow acoustics. Flow-generated sound is a one of the problems in many engineering

I = p ′ 2 ρ ∞ a ∞ (2.44)
The presence of cavities in vehicles changes the drag and heat transfer and may cause intense periodic oscillations, which in turn may lead to severe buffeting of aerodynamic structure and generation of sound [START_REF] Sarohia | Experimental investigation of oscillations in flows over shallow cavities[END_REF]. In aeroacoustics, turbulence is the principal source of broadband noise. With the recent increase in the performance of computers to perform numerical simulation of sound, Computational aero acoustics has become very popular. 2.9 Computational Aeroacoustics 2.9.1 Generalities Computational Aeroacoustics combines the classical approaches of flow field computation with acoustics. Computational methods for flow-generated sound can be divided into two kinds: direct computation and indirect, or hybrid computation. The direct approach computes the sound together with its fluid dynamic source field by solving the governing equations without modelling. The direct approach is very expensive than DNS. Because it solves all the scales from Kolmogorov microscale to distance covered by the sound waves in the computational domain. Other than very accurate numerical schemes, this approach needs high quality grids with less than 1% stretching. In the hybrid approach, the computation of flow is decoupled from the computation of sound, which can be performed during a post-processing stage based on aeroacoustic analogy.

The far-field sound is obtained by integral or numerical solutions of acoustic analogy equations using computed source field data. Figure 2.7 shows the main computational approaches which may be used when evaluating the sound field generated by turbulent flows. From the figure 2.8, the flow region is dominated by hydrodynamic phenomena.

The pressure fluctuations which are present in this region is due to turbulence or large structures. Consider a source region of characteristic length scale L source containing individual sources (eddies) of size ℓ ed . The hydrodynamic pressure fluctuations dominate, since the energy of the acoustic field is 1% of the total energy [START_REF] Ribner | The generation of sound by turbulent jets[END_REF]. The far field is a region where the turbulence is less and the mean flow field is typically homogeneous.

The far field and the source region is separated by a distance d. The only phenomena in this region is acoustic wave propagation. In the integral forms of acoustic analogies, the use of leading-order terms in an acoustic far-field expansion (with respect to λ ac d , where λ ac is the acoustic wavelength) leads to much simpler evaluations of sound. For small amplitudes and low Mach numbers M , far field can be described by a linear ho-

Acoustic analogy

The notion of "analogy"refers to the idea of representing a complex fluid mechanical process that acts as an acoustic source by an acoustically equivalent source term. The first major step in the development of acoustics was done by Sir James Lighthill [START_REF] Lighthill | On sound generated aerodynamically: I. General Theory[END_REF], [START_REF] Lighthill | On sound generated aerodynamically: II. Turbulence as a source of sound[END_REF] who published in 1952 his "acoustic analogy". This represents one of the first theories on aerodynamic noise generation for describing the radiation of the sound field generated by a turbulent flow. Hybrid method is where flow field is resolved using CFD solver and the flow field is employed as found in Lighthill [START_REF] Lighthill | On sound generated aerodynamically: I. General Theory[END_REF] where an analogy between the propagation of sound in an unsteady unbounded flow to that in an uniform medium at rest, generated

∂ 2 ρ ∂t 2 -a 2 ∞ ∂ 2 ρ ∂x 2 i = ∂ 2 T ij ∂x i ∂x j (2.47)
is non homogeneous wave equation and called as Lighthill equation. where analogy. Its use is justified at low Mach number flow where source-propagation ambiguities diminish and additional approximations can make it analytically more tractable. The Lighthill's analogy does not include the effect of solid boundaries in the flow, thus it considers only aerodynamically generated sound without solid body interaction. The formulation was extended by Curle [START_REF] Curle | The influence of solid boundaries upon aerodynamic sound[END_REF] and Ffows Williams and Hawkins [START_REF] Williams | Sound generation by turbulence and surfaces in arbitrary motion[END_REF] to take into account the generation and the scattering mechanisms when solid bodies are present.

T ij = ρu i u j -τ ij + (p -a 2 ∞ ρ)δ ij is
The solution to Lighthill's equation was given by Curle which is

ρ(x, t) -ρ 0 = 1 4πa 2 ∞ ∂ 2 ∂x i ∂x j V T ij r dV (y) Volume contribution - 1 4πa 2 ∞ ∂ ∂x i S n j r (pδ ij -τ ij )dS(y) Surface contribution (2.48)
x is the observer position, y is the source position and r =| x-y | is the distance between them. τ = t -r a ∞ is the retarded time, which is the time of the emission of a signal that reaches the observer location at time t. The displacement between the observer and the source can be expressed as r = (t -τ )a ∞ . If the observer in equation 2.48 is located in region where the flow is isentropic, the density fluctuation at this location can be written as

ρ(x, t) -ρ 0 = p(x, t) -p a 2 ∞ (2.49)
Using the following derivatives

∂f (τ ) ∂x i = ∂f ∂τ ∂τ ∂x i = - 1 a ∞ ∂r ∂x i ∂f ∂τ (2.50)
and

∂r ∂x i = ∂ (x j -y j ) 2 ∂x i = (x i -y i ) (x j -y j ) 2 = x i -y i r = l i (2.51)
equation 2.48 can be written as p(x, t) -p 0 = 1 4π

∂ ∂x i V -l i Ṫij a ∞ r + T ij r 2 dV (y) - 1 4π S -l i n j ṗδ ij -τij a ∞ r + pδ ij -τ ij r 2 dS(y) = 1 4π V l i l j Tij a 2 ∞ r + 2 Ṫij a ∞ r 2 + 2 T ij r 3 - ∂l i ∂x i Ṫij a ∞ r + T ij r 2 dV (y) + 1 4π S l i n j ṗδ ij -τij a ∞ r + pδ ij -τ ij r 2 dS(y) (2.52)
The derivative ∂l j ∂x i is expanded as

∂l j ∂x i = ∂ ∂x i x j -y j r = δ ij -l i l j r (2.53)
Inserting this expansion into equation 2.52

p(x, t) -p 0 = 1 4π V l i l j a 2 ∞ r Tij + 3l i l j -δ ij a ∞ r 2 Ṫij + 3l i l j -δ ij r 3 T ij dV (y) + 1 4π S l i n j ṗδ ij -τij a ∞ r + pδ ij -τ ij r 2 dS(y) (2.54)
The above derivation is followed from the work of Larsson et al [START_REF] Larsson | Aeroacoustic investigation of an open cavity at low Mach number[END_REF] and they identified the surface pressure dipole as the dominating terms for an open cavity. In the presence of walls, the sound radiation by turbulence is enhanced. The compact bodies radiate a dipole sound field associated with the force theory exert on the flow as a reaction to the hydrodynamic force of the flow applied on them. Sharp edges are particularly efficient radiators. In low Mach number flows, the main source of sound generation is due to the interaction of the flow with the cavity walls. The vortices shed from the cavity leading edge create pressure fluctuation when they impinge onto the cavity vertical wall. These surface pressure fluctuations make these surface integral contribution to far field noise dominant with respect to that of the volume integral. Larsson et al [START_REF] Larsson | Aeroacoustic investigation of an open cavity at low Mach number[END_REF] investigated in their numerical study this assertion by evaluating all the terms in Curle's acoustic analogy applied to a cavity flow and concluded that the volume integral contribution is indeed negligible. Curle's dimensional analysis [START_REF] Curle | The influence of solid boundaries upon aerodynamic sound[END_REF] also reports that the dipole sources along the wall becomes increasingly important at low Mach number over quadrupole sources. For performing numerical simulation, it is better to retain the spatial derivatives inside the integral. If the dipole terms are the main contributors to the radiated noise 2. Cavity flow, turbulence and aeroacoustics and neglecting the viscous term in the equation 2.54

ρ(x, t) -ρ ∞ = 1 4πa 2 ∞ S l i n j ṗδ ij a ∞ r + pδ ij r 2 dS(y) (2.55)
The pressure fields obtained from the simulation are only available in a two dimensional plane. Therefore the equation 2.55 is integrated in the out-of-plane direction from -∞ to +∞ yielding

p(x, t) -p 0 = 1 4π L l i n j π ṗδ ij a ∞ + 2 pδ ij r 2 dL(y) (2.56)
The flow is two-dimensional and uniform in the spanwise direction. The surface integral becomes a line integral along the cavity walls. Another two-dimensional form of the Curle's equation will be used, where equation 2.55 is integrated in the z-direction from -w to +w, where w is half the cavity spanwise extension, yielding

p(x, t) -p 0 = 1 4π L l i n j 2 arctan w r ṗδ ij a ∞ + 2w pδ ij r 2 dL(y) (2.57)
A general overview and details about other computational techniques used in CAA can be found in works of Larsson [START_REF] Larsson | Computational Aero Acoustics for Vehicle Applications[END_REF], Tam [START_REF] Tam | Computational Aeroacoustics: An overview of computational challenges and applications[END_REF] and Large-Eddy simulation for acoustics [START_REF] Wagner | Large-Eddy Simulation for Acoustics[END_REF].

Conclusion

Literature related to cavity flows is huge. Topics concerned to the present work is reviewed. Governing equations for the Direct numerical equations are given. Turbulence and RANS are also discussed. Lighthill-Curle's equation in two-dimensional form is derived to compute the sound pressure levels. The pressure field has to be determined by the numerical simulation. The pressure fields are fed as input to the equation (derived from the acoustic analogy) to measure the sound generated in the cavities. The last mentioned part will be dealt in the final Chapter 5.

Chapter 3

Inflow conditions and asymptotic modelling Contents

La section 3.6 présente la validation et l'analyse de l'approche asymptotique pour une couche limite turbulente d'équilibre en présence d'un gradient de pression adverse.

Il apparaît clairement que le modèle de Michel ne permet pas de bien prendre en compte les effets de gradient de pression, et qu'il faut modifier le paramètre n, jusqu'à n = 24 (approximativement) pour obtenir un bon accord entre l'approche asymptotique et les résultats de DNS proposés par Skote [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF], en particulier sur la contrainte turbulente.

Deux gradients de pression faibles et modérés ont été comparés. On observe que la viscosité turbulente sans dimension dépend fortement du gradient de pression dans le cas de l'approche asymptotique, alors que la dépendance est plus faible dans le cas des simulations numériques directes. Un post-traitement et une analyse des données DNS

disponibles met en évidence la difficulté de calculer avec précision le rapport u e /uτ , le nombre de Reynolds Re τ ou l'épaisseur de la couche limite dans le cas d'un gradient longitudinal de pression non nul. Des bosses sont observées en particulier sur le profil de la viscosité turbulence au voisinage de la fin de la couche limite, introduisant une erreur importante sur le calcul du Re τ . La sensibilité du modèle asymptotique déficitaire à la valeur de Re τ semble, en plus, plus importante dans le cas du gradient de pression adverse. Un autre modèle de longueur de mélange avec plus de paramètres pourrait certainement s'avérer nécessaire, en particulier lorsqu'on se rapproche du décollement.

En conclusion l'approche asymptotique, pour une couche limite d'équilibre est capable de fournir des profils de vitesse longitudinale d'excellente qualité pour servir de conditions d'entrée dans des simulations, pour peu que le modèle de longueur de mélange représente bien la physique. Ainsi la viscosité turbulente doit être calibrée vis-à-vis

du gradient de pression de l'écoulement. Par contre, du fait de ses limites, l'approche asymptotique ne peut fournir de profils de vitesse ou de contraintes turbulentes très convenables à faibles nombres de Reynolds.

Introduction

This Chapter is devoted to the determination of the inflow velocity profiles which are required for the test cases carried out with LES simulations where a turbulent boundary layer interacts with a cavity. Then, a brief introduction on boundary layers is given the analytical approach to provide velocity profiles, based on the defect law and the wake correction of the turbulent logarithmic law. Then, some turbulent boundary layer profiles are produced using Successive Complementary Expansion Method(SCEM). An alternate blending function is discussed under the mixing length model section and the validations with experimental are given. Zero pressure gradient and adverse pressure gradient cases were simulated using asymptotic approach and validated against the Direct Numerical Simulation data of Skote [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF] and experiments of Klebanoff [START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF],Townsend [START_REF] Townsend | The structure of the turbulent boundary layer[END_REF].

Boundary Layer

Laminar boundary layer

The presence of a wall has a dominant effect on the processes that produce turbulence.

The external flow is determined by the displacement of streamlines about the body and in which viscosity is negligible (potential flow) and the pressure field is developed.

But boundary layers are thin regions in the flow where viscous forces are important.

Although the name boundary layer originally referred to the layer of fluid next to the wall. The essential ideas are that the layer is thin in the direction across the streamlines and that viscous stresses are important only within the layer and the velocity satisfies the no-slip condition at the wall.

" A very satisfactory explanation of the physical process in the boundary layer between a fluid and a solid body could be obtained by the hypothesis of an adhesion of the fluid to the walls, that is, by the hypothesis of a zero relative velocity between fluid and wall. If the viscosity was very small and the fluid path along the wall not too long, the fluids velocity ought to resume its normal value at a very short distance from the wall. In the thin transition layer, however, the sharp changes of velocity, even with small coefficient of friction, produce marked results ". Prandtl postulated that the strain rate very near to the surface would become as large as necessary to compensate for the vanishing effect of viscosity, so that at least one viscous term remained. This very thin region near the wall became known as Prandtl's boundary layer, and the length scale characterising the necessary gradient in velocity became known as the boundary layer thickness [START_REF] George | Lectures in turbulence for the 21 st century[END_REF]. The boundary layer thickness δ(x) defined as the y value at which

u (x, y) | y=δ = 0.99u e (x) (3.1)
u e is the velocity outside the boundary layer, where fluid can be considered as inviscid.

In case of flat plate boundary layer with zero incidence, u e is constant and equals to the upstream velocity u ∞ . The thickness depends on small velocity differences. In other words, it is the layer where viscous effect continue to be important. As the layer is thin, the derivatives across the flow direction might be expected to be larger than derivatives in the flow direction. More reliable ways to characterise the thickness of boundary layer are displacement thickness δ * , momentum thickness θ and shape factor H. The flow near the surface is retarded, so that the streamlines must be displaced outwards to satisfy continuity. To reduce the total mass flow rate of a frictionless fluid by the same amount, Lele [START_REF] Colonius | Computational aeroacoustics: progress on nonlinear problems of sound generation[END_REF] that the value of momentum thickness θ at the cavity leading edge plays a vital 3.2. Boundary Layer role in the selection of the modes and in governing the growth of the shear layer(see Colonius & Lele [22], Rowley et al [START_REF] Rowley | On self-sustained oscillations in two-dimensional compressible flow over rectangular cavities[END_REF] and Tam [START_REF] Tam | Computational Aeroacoustics: An overview of computational challenges and applications[END_REF]) that spans an open cavity(see Charwat et al [START_REF] Charwat | An investigation of separated flows. Part I: The pressure field[END_REF]).

Turbulent boundary layer

Most of the flow around any body are turbulent in nature. For example turbulent boundary layer flow occurs on a high speed train, where the gap between the coaches build the cavity and the boundary layer developing along the train may have a size comparable to the cavity depth. In aeroplanes, these turbulent boundary layer flow occur during taking off, flying at high velocity, landing and taxing. Many researchers worked on the turbulent flows and turbulent boundary layers . Turbulent flows over (rough) walls have been studied by Hagen [START_REF] Hagen | Über den Einfluss der Temperatur auf die Bewegung des Wassers in Röhren[END_REF] in 1854 and Darcy [START_REF] Darcy | Recherches expérimentales relatives au mouvement de leau dans les tuyaux[END_REF] in (1857), who were concerned with pressure losses in water conduits. Study and analysis on the turbulent boundary layers were started while performing measurements in wind-tunnel experiments. Experiments performed by Schultz-Grunow [START_REF] Schultz-Grunow | Neues reibungswiederstandsgesetz für glatte platten[END_REF], Ludwieg and Tillman [START_REF] Ludwieg | Investigations of the wall shearing stress in turbulent boundary layers[END_REF], Klebanoff [START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF] and Smith & Walker [START_REF] Smith | Skin-friction measurements in incompressible flow[END_REF] were noteworthy. The first Direct Numerical simulation of a turbulent boundary layer was performed by Spalart [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF]. Skote et al [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF] obtained turbulent boundary layers at different pressure gradients. The overall structure of turbulent boundary layers can be found in textbooks for instance by Townsend [START_REF] Townsend | The Structure of Turbulent Shear Flows[END_REF]. In this work, mean thick turbulent boundary layer profiles are produced to impose on inlet of the computational domain to simulate cavity flows at different velocities and at different Reynolds number.

Power law

The algebraic law for a flat plate turbulent boundary layer under zero pressure gradient known as power law is given here. Because this approach has been initially used to generate inflow conditions in the simulation of cavity flows in this work. Then in the simulations at the later part, mean turbulent boundary layer profile were imposed on the computational domain of the cavity.

Consider an incompressible flow over a smooth flat plate (zero pressure gradient).

Simpler, but less accurate, relations between δ, δ * , θ and H can be obtained if one uses the power-law assumption for the velocity distribution in which one assume

u u ∞ = y δ 1 n (3.4)
Here the exponent n is about 7 in a constant pressure boundary-layer, increasing slowly with the Reynolds number. Using (3.4) and the definitions of δ * , θ and H, one can show 3. Inflow conditions and asymptotic modelling that

δ * δ = 1 1 + n θ δ = n (1 + n)(2 + n) (3.5) H = 2 + n n
Other formulas obtained from power-law assumptions, given by Schlichting [START_REF] Schlichting | Boundary layer theory[END_REF] are:

δ x = 0.37 U ∞ x ν ∞ -1/5 (3.6) θ x = 0.036 U ∞ x ν ∞ -1/5 (3.7)
Those equations are valid for Reynolds number Re x = U ∞ x ν ∞ , between 5 × 10 5 and 10 7 . The dimensionless skin friction co-efficient is

C f = τ w 1 2 ρu 2 ∞ (3.8)
At higher Reynolds numbers the boundary layer thickness can be calculated more accurately by the following empirical formula given by Granville δ x = 0.0598 logRe x -3.170 (3.9) This equation was obtained on the assumption that the boundary layer is turbulent from the leading edge onwards.

Analytical method

The boundary layer is described by a two-layer structure (see Mellor [START_REF] Mellor | The large Reynolds number, asymptotic theory of turbulent boundary layers[END_REF] and Yajnik [START_REF] Yajnik | Asymptotic theory of turbulent shear flows[END_REF]). The overall description of a turbulent boundary layer is dependent on two separate inner and outer length scales:

1. The outer length scale is commonly taken as the thickness of the boundary layer δ in outer layer where convective transport terms are important 2. an inner layer whose thickness is of order ν u τ , where u τ is friction velocity and is

u τ = τ w ρ
where τ w is the wall shear stress.

Analytical method

In between these layers, there is an overlap layer where both the convective transport and the viscous term are negligible. This is the logarithmic overlap region. The wall layer is further divided into a viscous sublayer where visous shear stress dominates and turbulent stresses are unimportant and into a buffer layer where both stresses have to be taken into account. These layers are well represented in the figure 3.2.

Normally in the boundary layer, viscosity effect is dominant below y + ≈ 5. The most active part of the flow lies between 10 y + 100 which is called the buffer region. The buffer layer is difficult to analyse theoretically since both viscous and turbulent stresses are important. For example, from the DNS of unsteady channel flow from Jiménez and P. Moin [START_REF] Jiménez | The minimal flow unit in near wall turbulence[END_REF] observed that in moderate Reynolds number flows, this buffer region generates most of the turbulent energy as it contains the nonlinear self-sustaining cycle.

Since the above said layers have different length scales, the whole turbulent boundary layer can never be self similar. The wall layers alone are self similar. The outer layers of so called equilibrium boundary layers can be considered approximately self similar.

To describe the mean velocity profile in a turbulent boundary layer, similarity solutions are sought in the inner and the outer regions. In the inner region, the mean stream wise velocity u scales with the wall friction velocity u τ and with the viscous length scale

l = ν u τ , so that u + = u u τ = f yu τ ν (3.10)
In outer region, the velocity profile is described by the velocity defect law

u e -u u τ = f y δ (3.11)
In eqs. 3.10 and 3.11, u + is the normalised stream wise velocity, u e is the free-stream velocity, ν is the kinematic viscosity, y is the wall-normal distance and δ is the boundary layer thickness, which is taken as the wall-normal distance at which u = u e . The outerlayer velocity distribution depends also on the external pressure gradient. Based on the existence of an overlap region between the inner and the outer regions, Coles [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF] proposed the following additive law of the wall and law of the wake in non-dimensional form:

u + = 1 κ ln y + + B + Π κ f (η) f (η) = 1 -cos (πη) (3.12) y + = yu τ ν , η = y δ
where y + is the non-dimensional wall-normal distance (also called inner variable), η is the non-dimensional wall-normal distance (also called outer variable), κ the von Kármán constant, B the logarithmic law constant and Π is the wake parameter. The wake parameter Π represents the effect on the outer layer dynamics. Coles [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF] determined 3. Inflow conditions and asymptotic modelling the wake parameter as

Π = κ 2 u + e - 1 κ ln Re τ -B (3.13) u + e = u e u τ , Re τ = δu τ ν
where u + e is the normalised free-stream velocity and Re τ is the boundary layer Reynolds number which defines the scale separation between the outer and inner lengths. For a given Reynolds number, this Π parameter cancels and then the changes the sign (becomes negative) for

Re τ = exp κ u + e -B
It characterises the deviation of log law profile at η → 1. At distances from the wall of the order of boundary layer thickness, the size of the structures is limited by δ, which becomes the relevant length scale. Let

f (η) = A 1 η 2 + A 2 η 3 (3.14)
be a cubic polynomial approximation to f (η) in eq. 3.12. Substituting the boundary conditions u| y=δ = u e and ∂u ∂y y=δ = 0 (3.15) in eq. 3.12, with f (η) from eq. 3.14, gives

A 1 = 6 1 + 1 6Π and A 2 = -4 1 + 1 4Π ,
with Π defined by eq. 3.13. The law of the wake of eq. 3.12 then becomes

u + = Log-law of the wall 1 κ ln y + + B + 1 κ (η) 2 (1 -η) Pure wall flow +2 Π κ (η) 2 (3 -2η)
Pure wake component

. ( 3

Successive Complementary Expansion Method

In this section, an approach is mainly developed for the boundary layer, but many extensions can be found for other flow such as channel flow.

Successive Complementary Expansion Method

According to Cousteix & Mauss [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF], Successive Complementary Expansion Method (SCEM) discusses about "singular perturbation problems" with a small parameter ǫ,

where when ǫ → 0, the solution does not tend uniformly towards the corresponding reduced problem obtained for ǫ = 0. It is necessary to observe that the non-uniformity occurs in a domain whose dimension is smaller than the initial domain. The principle of SCEM is to find an "uniformly valid approximation" which is uniformly valid in the whole flow field with an improved approximation near the walls. This improved approximation can be attained by adding a correction which takes into account the effects of viscosity. The successive complemetary expansion method consists here in seeking contiguous asymptotic matches between the inner and the outer regions of an incompressible turbulent boundary layer. This approach has been initially introduced by Schlichting [START_REF] Schlichting | Boundary layer theory[END_REF], Clauser [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF], Mellor & Gibson [START_REF] Mellor | Equilibrium turbulent boundary layers[END_REF] and Bradshaw [START_REF] Bradshaw | The turbulence structure of equilibrium boundary layers[END_REF].

Mixing length model

ρu ∂u ∂x + ρv ∂v ∂y = - ∂P ∂x + ∂ ∂y µ ∂u ∂y -ρ < u ′ v ′ >
where the pressure gradient is given by

∂P ∂x = dP dx = -ρ e u e u e x
because in a boundary layer flow, the pressure gradient across the flow is zero i.e ∂P ∂y = 0.

With zero incidence of the flat plate, the streamwise pressure gradient is zero as well and the pressure is constant.

Across the boundary layer, the local shear stress is given by where u ′ and v ′ are the time-dependent fluctuations of the streamwise and flow-normal velocity components and are unknown. To avoid having to resolve these unknowns, the Reynolds shear stress is evaluated using Prandtl's mixing length model [115] ℓ, with the Van Driest [START_REF] Van Driest | On dimensional analysis and the presentation of data in fluid-flow problems[END_REF] near-wall damping correction F. This gives

τ = τ turb + τ lam = -ρu ′ v ′
τ t = -ρu ′ v ′ = ρ F2 ℓ 2 ∂u ∂y ∂u ∂y (3.18)
F = 1 -exp - y + 26 (3.19)
In the inner region, ℓ = κy is linear, while in the outer region, ℓ/δ → 0.085 as y → δ.

These two asymptotic behaviour can be merged analytically into a single distribution for the mixing length ℓ across the whole boundary layer by the using a "blending" function.

Michel et al. [START_REF] Michel | Application d'un schéma de longueur de mélange à l'étude des couches limites turbulentes d'équilibre[END_REF] used a blending function which is

ℓ(η) = δ c ℓ tanh κη c ℓ (3.20)
with c ℓ = 0.085 and κ = 0.41. In [START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF] Airiau propose an alternative blending function which improves the prediction of the turbulent shear stress profile at the interface between the inner and the outer layer, at low Reynolds numbers Re τ . This is

ℓ(η) = δ κη 1 + κη c ℓ n 1 n (3.21)
For 2.6 < n < 2.7, the ℓ(η) profile from equation 3.21 almost matches that from equation 3.20.

Inner region velocity profile

Inner region is given by y + ∈ [0 + , 50 -100]. Normalising the local shear stress τ by τ w = ρu 2 τ and assuming a monotonic velocity profile, from eq. 3.18,

1 = ∂u + ∂y + + ℓ +2 F2 ∂u + ∂y + 2 (3.22)
where ℓ + = ℓu τ ν . In the viscous layer, the zone close to the wall with y + < 1, the velocity is really small and ℓ + is linear with respect to y + . The turbulent shear stress is negligible comparing to the viscous laminar shear stress. Then

1 = ∂u + ∂y + ⇒ u + = y +
This approximation falls in the range 10 < y + < 40. F (y + ) → 1 when y + > 60-80 and the viscous term is neglected:

1 = κ 2 y + 2 ∂u + ∂y + ⇒ ∂u + ∂y + = 1 κy +
then the velocity satisfies a logarithmic law :

u + = 1 κ log y + + C (3.23)
This region where the log-law is true is called the logarithmic region. Integration of equation (3.22)(for 40 < y + < 100 -1000) gives (see Cousteix [START_REF] Cousteix | Turbulence et couche limite[END_REF], Schlichting [START_REF] Schlichting | Boundary layer theory[END_REF])

C ≈ 5.25.
Recent calculation from Cousteix [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF] produces the value 5.28 (calculations peformed in this work produce 5.28). Equation 3.22 is a quadratic in ∂u + ∂y + so the analytical solution is given by :

∂u + ∂y + = 2 1 + 1 + 4 ℓ + (y + ) F (y + ) 2 (3.24)
Integrating equation 3.24 with respect to y + with the boundary condition u + (x, 0) = 0 gives the inner layer tangential velocity profile that asymptotes to the log-law of the wall in equation 3.16 for y + → ∞.

Outer region velocity profile

In the outer region, the Reynolds stress component is dominant over the laminar shear stress where viscous stress is negligible, so τ ≃ τ t . From eq. 3.18, with the van Driest damping constant F → 1 at y + ≥ 100. The shear stress is :

τ = τ t = ρℓ 2 ∂u ∂y 2 (3.25)
In an equilibrium turbulent boundary layer, the similarity solution for the outer layer can be expressed in terms of the velocity defect F ′ (η) = u + e -u + and the shear stress is obtained from the integration of the streamwise momentum equation :

τ + = τ τ w = 1 - F F 1 + 1 F 1 + 2β ηF ′ (3.26)
where

F (η) = η 0 F ′ (ξ) dη F 1 = F (1) β = - δ u τ du e dx (3.27)
The shear stress, from equation 3.25 is expressed as a function of the derivative of the defect law F :

τ τ w = ℓ δ 2 F ′′2
where F ′′ = dF ′ dη . Substituting for τ τ w in eq. 3.26, the similarity solution for the outer region becomes ℓ δ

2 F ′′2 = 1 - F F 1 + 1 F 1 + 2β ηF ′ (3.28)
The parameter β represents the pressure gradient. Clauser had defined the factor β c

β c = δ * u τ dP dx (3.29)
and is related to β as

β c = β δ * δ u e u τ (3.30)
For zero pressure gradient flow, β = 0.

To determine the most important boundary layer quantities, it is necessary to calculate the values of F 1 , F 2 and G as

F 1 = 1 0 F ′ dη, F 2 = 1 0 F ′2 dη, G = F 2 F 1
For β = 0, and with Michel's mixing length model, F 1 = 3.15 and G = 6.13.

In the neighbourhood of η = 0, it is easy to demonstrate that F ′ (η) becomes loga-

Successive Complementary Expansion Method

rithmic from the equation (3.28)

ℓ δ 2 F ′′2 = 1, L(η) = ℓ δ = κη F ′ (η) = - 1 κ log η + D v (β)
For a zero pressure gradient boundary layer flow (β = 0), Cousteix & Mauss [START_REF] Cousteix | Asymptotic Analysis and Boundary Layers[END_REF] give

D v = 1.76.

Asymptotic matching of the inner and outer profiles

A matching condition is sought for the velocity profiles of the inner and outer regions, solutions of equations 3.24 and 3.28. This is obtained from standard asymptotic analysis (Cousteix & Mauss [25]) by considering y + → ∞ in equation 3.24 and η → 0 in equation 3.28. That gives respectively

u + = 1 κ ln y + + C (3.31) u + e -u + = - 1 κ ln η + D v (3.32)
Adding eq. 3.31 to eq. 3.32 gives

u + e = 1 κ ln Re τ + C + D v (3.33)
Equation 3.33 can be re-casted as function of the wall skin friction coefficient

C f = 2 τ w (ρu 2 e ) = ρu 2 τ 1 2 ρu 2 e = 2γ 2 (3.34)
that is imposed with same value in the inner and outer regions and provides the matching criterion for the two profiles

2 C f = 1 κ ln Re τ + C + D v (3.35)

Boundary layer quantities

It is possible to calculate analytically the displacement thickness, the momentum thickness and the shape factor of the boundary layer as soon as the velocity profile is known.

γ = C f 2 = u τ u e = 1 u + e

Inflow conditions and asymptotic modelling

The displacement thickness δ * is given by

δ * = δ 0 1 - u u e dy =⇒ δ * δ = γF 1 with γ = C f 2 = u τ u e = 1 u + e .
The momentum thickness θ is determined by

θ = δ 0 u u e 1 - u u e dy =⇒ θ δ = γF 1 -γ 2 F 2
The shape factor H of the boundary layer is

H = δ * θ = 1 1 -γG

Turbulent shear stress and turbulent viscosity

To compare results with experimental or numerical data, the turbulent shear stress and the turbulent viscosity are converted to non-dimensional form, in the inner and outer region.

Turbulent shear stress

The turbulent shear stress values are calculated using mixing length model.

τ t = -ρ < u ′ v ′ >= ρ F 2 y + ℓ + ∂u ∂y 2
In the inner region, the non-dimensional shear stress is

τ t τ w = F y + L y + R τ R τ ∂u + ∂y + 2
Usually, for the defect zone, the damping function is not taken into account F , but here F is retained. Because F = 1 only for y + ≤ 100 :

τ t τ w = F y + L(η)F ′′ (η) 2 (3.36)

Derivative ∂u ∂y

In the internal layer, the velocity derivative in inner variable y + is given by:

∂u ∂y = u τ l + ∂u + ∂y + = u 2 τ ν ∂u + ∂y + and l + = ν u τ (3.37)

Successive Complementary Expansion Method

For the external layer, in outer variable η:

∂u ∂y = - u τ δ F ′′ (η) = u 2 τ ν ∂u + ∂y + and ∂u + ∂y + = - 1 Re τ F ′′ (η) (3.38)
Turbulent dynamic viscosity ν t By the definition, the turbulent dynamic viscosity ν t is given in dimensional form as

ν t = F 2 ℓ 2 ∂u ∂y (3.39) in non-dimensional form νt νt = ν t u τ δ = F 2 ℓ δ ℓ u τ ∂u ∂y (3.40)
In the internal layer, the above expression is written in variable

y + as νt = F 2 ℓ δ ℓ u τ u 2 τ ν ∂u + ∂y + = Re τ F 2 (ℓ + (y + )) 2 ∂u + ∂y + = Re τ F 2 (κy + ) 2 ∂u + ∂y + (3.41)
For the external layer, with the variable η, it yields

νt = F 2 L 2 (η) F ′′ (η) = Re τ F 2 L 2 ∂u + ∂y + (3.42)
The relationship between the non-dimensional turbulent viscosity νt and the nondimensional turbulent stress τ t τ w is deduced:

• In the inner region, in variable y + :

τ + t = νt Re τ ∂u + ∂y + (3.43)
• In the outer region, in variable η :

τ + t = νt F ′′ (η) (3.44)

Numerical implementation

Expliciting the outer region velocity profile poses several challenges. Equation 3.28 is non-linear and is ill-defined in at the upper boundary layer limit, at η → 1, where F ′′ → 0, and at the lower boundary layer limit, at η → 0, where ℓ/δ → 0 and F ′′ → ∞.

To solve the problem, auxiliary approximate solutions are imposed on the floor of the laminar sub-layer and at the edge of the boundary layer, as shown in figure 3.4 so that the edges of the inner and of the outer regions are modelled analytically while the overlap region is resolved numerically.

Inflow conditions and asymptotic modelling

η = ǫ 0 η = 1 η = 1 -ǫ 1
numerical integration of eq. Let

f (η) = F (η) F (1)
. On the floor of the laminar sub-layer, imposing η = 0 and ℓ = κy, as in section 3.4.1, eq. 3.28 becomes

κ η F 1 f ′′ (η) 2 = 1 -f (η) + (1 + 2βF 1 ) ηf ′ (η) (3.45)
with the boundary condition f (0) = 0. Let β = 2βF 1 . In a zero pressure gradient boundary layer, β = 0 by eq. 3.27, for which eq. 3.45 has the explicit solution

f (η) = η 2 4α 2 - η ln η α + C η f ′ (η) = η 2α 2 - 1 + log η α + C f ′′ (η) = 1 2α 2 - 1 αη with α = F 1 κ. The integration constant C is determined by evaluating f ′ (η) at η = ǫ 0
on the floor of the laminar sub-layer. In a non-zero pressure gradient boundary layer, βηf ′ → 0 as η → 0, so the zero pressure gradient profile is used on the floor of the laminar sub-layer.

At the edge of the boundary layer, close to η = 1, eq. 3.28 becomes

ℓ 1 F 1 f ′′ (η) 2 = 1 -f (η) + 1 + β ηf ′ (η) (3.46)
with the boundary conditions f (1) = 1, f ′ (1) = 0, f ′′ (1) = 0 and ℓ 1 is evaluated from eq. 3.21 at η = 1. Cousteix [START_REF] Cousteix | Outer boundary layer self-similar solution[END_REF] proposed the solution for eq. 3.46 :

f (η) = 1 - (1 -η) 3 3 (3.47) f ′ (η) = (1 -η) 2 f ′′ (η) = -2 + 2η
For β = 0, the analytical solution has the attractive property of being independent from F 1 and ℓ 1 . The same solution is used in case of pressure gradient flow (β = 0), as βηf ′ (η) = 0 by the boundary condition f ′ (1) = 0 in eq. 3.46.

3.5 Zero pressure gradient boundary layer

Comparison of velocity profiles

The analytical and numerical methods for predicting a boundary layer mean turbulent velocity profile are tested against a range of streamwise velocity reference data (experiments and numerical simulations) from zero pressure gradient boundary layers of Spalart [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF], Erm and Joubert [START_REF] Erm | Low-Reynolds-number turbulent boundary layers[END_REF], De Graaff and Eaton [START_REF] De Graaff | Reynolds number scaling of the flatplate turbulent boundary layer[END_REF] and Österlund [START_REF] Österlund | Experimental studies of zero pressure-gradient turbulent boundary layer flow[END_REF], over the range Re θ ∈ [300, 31000]. Table 3.1 lists the values of u + e , Re τ and Π at each Re θ of the reference velocity records. The values of u + e and Re τ (column 1 to 3) are the ones reported in publications [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF][START_REF] Erm | Low-Reynolds-number turbulent boundary layers[END_REF][START_REF] De Graaff | Reynolds number scaling of the flatplate turbulent boundary layer[END_REF][START_REF] Österlund | Experimental studies of zero pressure-gradient turbulent boundary layer flow[END_REF] while Π (column 4) has been obtained by fitting eq. 3.16 using the least squares fit. The normalised mean streamwise velocity u + is plotted against the normalised wall-normal distance y + in figure 3.5 for different Reynolds numbers. The symbols used in figure 3.5 are measured values of [START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF][START_REF] Erm | Low-Reynolds-number turbulent boundary layers[END_REF][START_REF] De Graaff | Reynolds number scaling of the flatplate turbulent boundary layer[END_REF][START_REF] Österlund | Experimental studies of zero pressure-gradient turbulent boundary layer flow[END_REF] at different Re θ , labelled as in table 3.1. The continuous lines show the fitted analytical profiles for the outer layer. For clarity, an incremental shift of u + = 2.5 is applied to all curves. The three "0"labels on the vertical axis of figure 3 [START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF].

Re θ = 5200, and Re θ = 31000 respectively. The quality of the predictions is quantified by evaluating the mean square percentage error ǫ for each profile

ǫ = 1 N N i=1 u + e -u + ref u + ref 2 (3.48)
where u + a is the predicted value and u + ref is the corresponding reference (experimental, numrical) value for a given y + i in a discretized velocity profile of N points. The ǫ obtained at different Re θ with u + a evaluated from equation 3.16 is reported in table 3.1 (column 5). The maximum ǫ is 2.05% at Re θ = 31000. Such error enables the use of eq. 3.16 to predict the mean streamwise velocity of boundary layers in many common engineering applications, where an error margin of 5% is often acceptable. The reference data seem to be randomly distributed about the fitted curve with no underlying trend, suggesting that the curve fit has captured most of the u + dependence on δ, u e , u τ , and Re θ . velocity profile down to the wall. In the outer layer, the asymptotic method captures the Reynolds number dependent transition between the log-law and the constant freestream velocity for most of the curves. The free-stream velocity at Re θ = 22845, 12663 and 3654 appear to be under-predicted. In table 3.1, are found from column 7 to 9, the Reynolds number Re τ , the non-dimensional velocity u + e and the non-dimensional error ǫ num given from the asymptotic approach. The parameters Re θ , Re τ and u + e of C f are directly related for a given value of the pressure gradient coefficient β. Here the Reynolds numbers was choosen from the reference value, and Re τ and u + e were determined by an iterative Newton approach. The differences between the reference data and the present calculations are confirmed by the corresponding numerical mean square percentage error, ǫ num , which is computed by evaluating u + e in eq. 3.48. Specifically, the ǫ num at Re θ = 22845, 12663 and 3654 are higher than for some of the other Reynolds numbers, due to the difference in the normalised free-stream velocity between experiment and prediction.

Whereas, in general, the error from the numerical velocity profile is higher than that from the analytical profile, it is within the range for which the predictions can be used for engineering accurate predictions.

The difference between the normalised free-stream velocity from experiments and from the SCEM approach is further investigated in figure 3.7, where the outer layer portion of the predicted velocity profile for Re θ = 22845 is re-plotted on a larger scale.

The continuous black line is the numerical prediction obtained by matching the experimental value of Re θ in the matched complementary expansion, the red dash-dot line is obtained by matching the experimental value of Re τ , while the dashed blue line shows the predicted profile with a matched normalised free-stream velocity u + e . Matching the experimental Reynolds numbers seems to give similar profiles irrespective of whether the target Reynolds number is defined with respect to the momentum thickness, Re θ , or the friction velocity, Re τ . Fitting the outer profile by imposing the normalised free-stream velocity u + e seems to over-predict the boundary layer thickness, leading to a coarser fit with experiment compared to the numerical predictions obtained by matching the profile Reynolds number.

Validation of the new mixing length model with experiments

The optimised value of the n parameter in the new mixing length model (eq. 3.21) has been determined to fit with the experimental measurements of the non-dimensional length ℓ reported in Klebanoff [START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF]. shown here) and it was found to have very little effect on the predicted normalised ν t , which is also the trend in experiment [START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF][START_REF] Townsend | The structure of the turbulent boundary layer[END_REF].

In Rona et al [START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF], no attempt have been made to predict the time-averaged velocity profiles of a boundary layers at Re τ < 300. A small explanation is required. In the asymptotic approach, with the skin friction value, u + e , is obtained by matching the outer layer velocity profile to the inner layer velocity profile in the logarithmic layer. When Re τ < 140, an overlap region in the form of a logarithmic layer is no longer present, which prevents the method from evaluating u + e . Here the matched complementary expansion method in its present formulation has reached its Re τ applicability limit. To illustrate this upper limit in Re τ , the Figure 3.9 shows the velocity profile created using asymptotic approach for Re τ = 900 in inner variable y + . In the u + profile, the log law region is shed on light in the interval y + ∈ [START_REF] Gharib | The effect of flow oscillations on cavity drag[END_REF]200]. The inner region velocity profile is obtained by integrating the equation 3.24 and the velocity profile of the outer region is determined The non-dimensional shear stress τ t τ w and the derivative du + dy + (see 3.4.6) are plotted against inner variable y + in the figure 3.10. The shear stress (continous line) which is obtained from the asymptotic method shows a discontinuity near y + ∼ 100, in the overlapping log-law region. We obviously observe that the continuity of the velocity and of the shear stress are fixed. The derivative du + dy + decreases smoothly with increase in y + and the mixing length continuously grows from zero at the wall to a constant value at the edge of the boundary layer which implies that the shear shear is maximum at a given distance from the wall (quite close to the wall), decreases away from wall and goes to zero in the external flow (i.e outside the boundary layer). This discontinuity on the non-dimensional shear stress results from the complex product of the decreasing function du + dy + and of the increasing function, the mixing length ℓ. This discontinuity neither exist in the reality (experiment) nor in the Direct Numerical Simulation (see later).

Comparison with Direct Numerical Simulation

For turbulent flat plate boundary, numerous accurate direct numerical simulations are not available, especially in the case of equilibrium boundary layer. For this comparison, numerical data (shear stress and profile) which is referred here as Skote's data is available on-line is taken as reference. Analysis and curves can be found in the PhD of Skote [START_REF] Skote | Studies of turbulent boundary layer flow through direct numerical simulation[END_REF] and in [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF]. Three cases are considered in equilibrium turbulent state: a turbulent boundary layer with the zero pressure gradient (ZPG) flow and two cases with small and moderate adverse pressure gradient (APG1 and APG2). A summary of the Skote data [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF] is given in the table 3.2 as described in the reference publication. The nondimensional pressure gradient is given by the Clauser parameter β c (eq. 3.29). In [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF] a different way is presented to evaluate the pressure gradient parameter and the equilibrium state is discussed as well. The β c value given here can be considered as mean value over a given range of Reynolds number Re θ .

From the numerical data, the main important quantities which characterise a turbulent boundary layer have been calculated. For instance, the following parameters F 1 , G and R θ Re τ have been determined from the following formula :

F 1 = 1 0 u + e -u + dη, G = u + e 1 - 1 H , R θ Re τ = u + e θ δ
Re tau and u + e have been read from the files (column 2 of the table 3 It should be noted that the DNS always produces higher shape factor H and parameter G values than in the asymptotic case. On the contrary, the non-dimensional external velocity u + e and consequently the Reynolds number Re θ are over-estimated with the asymptotic approach. The difference can have different reasons, from the difficulty to evaluate the exact value of the pressure gradient in DNS's data to the assumption made as the equilibrium state of the boundary layer. The difficulty of determining a right value of the Reynolds number form DNS's data is also discussed in a next section. [START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF] with Re τ = 1540, and experiment of Townsend [START_REF] Townsend | The structure of the turbulent boundary layer[END_REF] with Re τ = 2775. The figure shows the influence of the Re τ number on the turbulent viscosity. This influence seems to be higher with the asymptotic approach, if DNS results are considered as the reference. The discontinuity of slope in the asymptotic curves comes from the discontinuity observed in the turbulent shear layer.

For a given Re τ number and in the region close to the wall, the turbulent viscosity determined from asymptotic approach fits very well with that of DNS data, indicating a really good evaluation of the skin friction. In the region η ∈ [0.4, 0.8], all the turbulent viscosity curves from experiment, DNS and asymptotic approach fit together. But in the interval η ∈ [0.1, 0.4], turbulent viscosity from asymptotic approach is underestimated and a discontinuity appears in this region. The other region η ∈ [0.8, 1] experimental data and asymptotic approach curves are in good agreement except the curves from pressure gradient and dramatically it increases the drag with decreasing lift in turbulent flows. Investigation of adverse pressure gradient boundary layer and control of separation are the two main topics in aerodynamics. Here, focus is laid on small or moderate streamwise pressure gradient, before separation. In the present asymptotic approach, the wall shear stress is used as the reference quantity (or as parameter) which was the final output of the problem through the skin friction coefficient.

Experimental work of Clauser [START_REF] Clauser | Turbulent boundary layers in adverse pressure gradients[END_REF] and the work of Rotta [START_REF] Rotta | Turbulent boundary layers in incompressible flow[END_REF] demonstrated that equilibrium boundary layers in both zero and adverse pressure gradients could exist at least approximately for a certain distance along a smooth wall. Other notable experiments are of Herring & Norbury [START_REF] Herring | Experiments on equilibrium turbulent boundary layers in favourable pressure gradients[END_REF] in favorable pressure gradients and Bradshaw [START_REF] Bradshaw | The turbulence structure of equilibrium boundary layers[END_REF] in adverse pressure gradients.

Townsend [START_REF] Townsend | Equilibrium layers and wall turbulence[END_REF] tried to set out the necessary conditions for the existence of an equi- The two cases are the constant-pressure and adverse-pressure gradient flows that have investigated experimentally. For the adverse pressure gradient case from figure 3.16(b), the maximum u + value is found higher (u + ≈ 60) than the case with constant pressure (u + ≈ 30) which is observed in the figure 3.16(a).

Comparison with DNS

To validate the proposed approach and to test the new blending function with the parameter n in the mixing length model, the obtained results are compared with Skote's DNS results. As for the zero pressure gradient testcase, the tables 3. A good agreement is obtained from asymtotic approach and Skote's DNS is found. A main difference was observed in the outer part of the boundary layer, since asymptotic 3.6. Adverse pressure gradient boundary layer approach predicts slight over estimation of the non-dimensional external velocity u + e . The influence of the n parameter exists on the velocity profile even if it is really weak.

The agreement is improved by increasing the value of n.

In the zero pressure gradient case, the turbulent shear stress is the relevant quantity for comparison or validation. The shear stresses of APG1 and APG2 cases (which are calculated by asymptotic approach) are plotted in the figures 3. The moderate pressure gradient APG2 case is a little bit different (see figure 3. [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF]).

An inappreciable agreement was not observed. The asymptotic approach seems to follow the DNS in parallel, the discontinuity in the slope of the turbulent shear stress is less significant (the sign of the slope is the same) and just one global maximum is produced by the asymptotic analysis. However the same trend is observed in the both plots, a change in the slope occurs when y + ≈ 90 and the peak of the shear stress is predicted at the same location. As for the previous case, increasing the parameter n from 4 to 24 improve the agreement, especially in the overlapping region.

One of the conclusions could be that the parameter n is related to the pressure gradient. 1. The distance from the wall where the velocity defect u + d is smaller than a given small parameter ǫ (0.001 for instance).

Eddy viscosity

2. The distance from the wall where the non-dimensional turbulent shear stress τ + is smaller than a given small parameter ǫ (0.001 for instance).

3. The distance from the wall where the non-dimensional turbulent viscosity is smaller than a given small parameter ǫ (0.001 for instance).

In a 'perfect' turbulent boundary layer, these definitions should provide the same value. As explained earlier, the last definition corresponds to a limit of 0 over 0 in the edge of boundary layer which must be equal to 0. Theoretically it can be imposed, numerically the round off errors, or any other numerical approximation can not insure a convergence of the turbulent viscosity to zero at the edge of the boundary layer.

In this study, the guess of Re τ from the DNS data has followed the two first definitions. This led to the results provided here. Figure 3.22 represents the following function which is proportional to the non-dimensional turbulent viscosity : is observed in the pressure gradient case. Finally, from Skote et al [START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF], it can be concluded that the bumps on the turbulent viscosity could be an effect of a shorter height of computational domain and the boundary condition at this location.

Re τ νt = τ + du + dy + = Re τ ν t u τ δ 3
DNS data can provide a valuable element of comparison but can not be used as the only source for validation.

Conclusion

Analytical and asymptotic methods were followed to obtain the time averaged velocity profiles of a turbulent boundary layer and validation against reference data have been carried out. The analytical method given in this chapter is an extension to wake law from Coles [START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF] that matches both the free stream velocity and the velocity gradient at Finally, it seems that under the equilibrium assumption, the asymptotic approach is able to provide accurate velocity and turbulent stress profiles which can be imposed as inlet of the computational domain where DNS, RANS or LES is performed. For the present work, it is required for the cavity flow simulations.

These works have been published in the four proceedings [START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF], [START_REF] Gandhi | Analysis of turbulence and noise emission in open rectangular cavity and separated flows[END_REF], [44], [START_REF] Rona | On the generation of the mean velocity profile for turbulent boundary layers with pressure gradient under equilibrium conditions[END_REF] Les méthodes numériques sont basées sur des schémas de type Lax-Wendroff [START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF] ou de type éléments finis faible-dissipation Taylor-Galerkin(Donea [START_REF] Donea | Taylor-Galerkin method for convective transport problems[END_REF], Donea et al [START_REF] Donea | An analysis of time discretization in the nite element solution of hyperbolic problems[END_REF],

Quartapelle & Selmin [START_REF] Quartapelle | High-order Taylor-Galerkin methods for non-linear multidimensional problems[END_REF], Colin & Rudgyard [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF]). Un modèle de

L'approche des résidus du volume de contrôle

Pour la description de l'approche on considère les équations de Navier-Stokes laminaires sous forme conservatives ∂U ∂t

+ ∇ • F = 0 (4.1)
Les termes d'espace sont approximés à chaque volume de contrôle, pour obtenir le résidu

R Ω j = 1 V Ω j ∂Ω j F • ndS (4.
R Ω j = 1 N d V Ω j i∈Ω j F i .d S i (4.3)
Le volume V Ω j est défini par:

V Ω j = 1 N 2 d i∈Ω j x i .d S i (4.4) ∇ • x = N d .
Une fois le résidu calculé on obtient la forme semi-discrète:

dU k dt = - 1 V k j|k∈Ω j D k Ω j V Ω j R Ω j (4.
D k Ω j = 1 n n I + C δt Ω j V Ω j A Ω j • d S k (4.6)
Calcul des gradients

∂U ∂x C ≈ 1 V C ∂Ω C U • n∂S (4.7)
qui donne l'approximation suivante:

∇U Ω j = 1 V Ω j i∈Ω j U i d S i (4.8)
L'approximation du gradient au noeud est obtenue en réalisant la moyenne des gradients des cellules:

∇U k = 1 V Ω k j|k∈Ω j V j ∇U Ω j (4.9)
Calcul du pas de temps

U n+1 = U n + ∆t ∂U ∂t n + 1 2 ∆t 2 ∂ 2 U ∂t 2 n + O(∆t 3 ) (4.11)
Soit la forme conservative des équations de Navier-Stokes laminaires:

∂U ∂t + ∇ • F = 0 (4.12)
On a :

∂U ∂t = -∇ • F (4.13)
et:

∂ 2 U ∂t 2 = ∂ ∂t (-∇ • F) = -∇ • ∂F ∂t = -∇ • A ∂U ∂t = ∇ • [A (∇ • F)] (4.14)
si A = ∂F ∂U la matrice Jacobienne. La solution au pas de temps n + 1 est donnée par:

U n+1 = U n -∆t ∇ • F - 1 2 ∆t ∇ • [A (∇ • F)] -O ∆t 2 (4.15)
Le schéma à deux pas de Taylor-Galerkin (TTGC)

Il est difficile de développer des schémas d'ordres très élevés (en espace) pour des maillages structurés en volumes finis. La formulation aux noeuds peut être étendue à l'approche éléments finis où les ordres élevés sont possibles. Les schémas de Taylor-Galerkin (TG) développées initialement par Donea [START_REF] Donea | Taylor-Galerkin method for convective transport problems[END_REF][START_REF] Donea | An analysis of time discretization in the nite element solution of hyperbolic problems[END_REF] combinent le dveloppement de Taylor en temps et la discrétisation en espace de Galerkin. Colin et Rudgyard [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF] ont développé le schéma (TTGC) de troisième ordre en espace et en temps.

Ũn = U n + α∆t ∂U ∂t n + β∆t 2 ∂ 2 U ∂t 2 n (4.16
)

U n+1 = U n + ∆t ∂ Ũ ∂t n + γ∆t 2 ∂ 2 Ũ ∂t 2 n ( 4
Ũn = U n -α∆t∇ • F n + β∆t 2 ∇ • [A (∇ • F n )] (4.19) U n+1 = U n -∆t∇ • F n + γ∆t 2 ∇ • [A (∇ • F n )] (4.20)
Multipliant ces équations par une série de fonctions tests linéaires φ i ("redskin tent"functions) et intégrant le résidu sur le domaine du calcul Ω, nous obtenons cette formulation faible:

Ω Rn φ i dV = -αL i (U n ) -β∆tLL i (U n ) (4.21) Ω R n+1 φ i dV = -L i Ũn -γ∆tLL i (U n ) (4.22) avec Rn = Ũn -U n ∆t , R n+1 = U n+1 -U n ∆t (4.23) et L i (U) = Ω ∇ • F (U n ) φ i dV (4.24) LL i (U) = Ω A (∇ • F (U n )) ∇φ i dV LL 0 i (U n ) - ∂Ω φ i A (∇ • F (U n )) dS BT i (U n ) (4.25)
Le terme LL i peut être séparé en faisant une intégration par partie en supposant qua la normale à la surface dS est externe. La première contribution LL 0 i (U n ) est intégrable sur tout le domaine alors que la seconde, BT i (U n ), est non nulle uniquement sur les limites du domaine. La méthode de Galerkin est ensuite appliquée à la divergence du flux et aux résidus. Ainsi, elles peuvent être exprimées comme une somme de fonctions-forme linéaires ( identiques aux fonctions-test utilisées pour dériver la formulation faible), 4. Numerical simulation and LES models donnant:

R n = k R n k φ k (4.26) ∇ • F = k F k ∇φ k (4.27)
où F k le flux discrétisé en chaque point du domaine. Avec le choix des fonctions forme, les résidus sont exprimés ainsi:

Ω Rn φ i dV = k Ω φ i φ k dV Rn k = k M ik Rn k (4.28)
notant M ik comme les composantes de ce qui appelée la matrice de masse qui, dans AVBP, est inversée localement par la méthode de Jacobi itérative.

Dans la discrétisation, les contributions des intégrales dans l'équation 4.24 et 4.25 permettent d'avoir i comme seulement provenant des cellules adjacentes.

L i (U n ) = j|i∈Ω j L i (U n ) Ω j (4.29) LL i (U n ) = j|i∈Ω j LL i (U n ) Ω j (4.30)
En utilisant les équations 4.26,

L i (U n ) Ω j et LL i (U n ) Ω j on a: L i (U n ) Ω j = k|k∈Ω j F n k Ω j φ i ∇φ k dV (4.31) LL i (U n ) Ω j = A n Ω j k|k∈Ω j F n k Ω j ∇φ • ∇φ k dV -A n Ω j k|∈∂Ω j ∩∂Ω F n k ∂Ω j ∩∂Ω φ i ∇φ k dS (4.32)
Pour les éléments triangulaires et tétrahédriques le gradient de la fonction forme est constant sur chaque élément et l'intégrale de φ i prend une forme simple (see Colin & Rudgyard [19]).

∇φ k = - S k n d V Ω j (4.33) Ω k φ k dV = V Ω j n v(Ω j ) ∀k ∈ Ω j ( 4
L i (U n ) Ω j = k|k∈Ω j F n k ∇φ k Ω j φ i dV = (∇ • F n ) Ω j Ω j φ i dV = R n Ω j V Ω j n v(Ω j ) (4.35)
De même:

LL 0 i (U n ) Ω j = A n Ω j k|k∈Ω j F n k ∇φ k • Ω j ∇φ i dV = - 1 n d A n Ω j R n Ω j • S i|Ω j (4.36)
Pour plus d'informations sur les schémas dans AVBP vous pouvez consulter le chapitre 5 de la thèse de Lamarque [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF] .

Termes de diffusion

Dans les équations de Navier-Stokes d'espèces ou de modèles, on a des termes de diffusion qui ont la forme générale suivante:

∂u ∂t = ∇. (ν∇u) (4.37)
AVBP utilise deux differentes discrétisations du terme de diffusion: un opérateur 4△ pour la diffusion laminaire et un opérateur 2△ pour la diffusion turbulente. La discrétisation est aussi en volumes fines ou éléments finis. Pour plus de détails voir le manuel de AVBP [10]. Pour les volumes finis, on a:

u n+1 i -u n i ∆t = 1 V i ∇. (ν∇u) | i (4.38)
Avec les éléments finis, la matrice de masse est appliquée et on a:

u n+1 i -u n i ∆t = M -1 ∇. (ν∇u) | i (4.39)
Viscosité Artificielle

ζ J Ω j = max k∈Ω j ζ J k (4.40)
Pour les scalaires (la pression P ) :

ζ J k = | ∆ k 1 -∆ k 2 | | ∆ k 1 | + | ∆ k 2 | + | P k | (4.41) où ∆ k 1 et ∆ k 2 sont: ∆ k 1 = P Ω j -P k ∆ k 2 = ∇P k • x Ω j -x k (4.42)

Le capteur de Colin

Il est défini par:

ζ C Ω j = 1 2 1 + tanh Ψ -Ψ 0 δ - 1 2 1 + tanh -Ψ 0 δ (4.43) avec Ψ = max k∈Ω j 0, ∆ k | ∆ k | +ǫ 1 P k ζ J k (4.44) ∆ k = | ∆ k 1 -∆ k 2 | -ǫ k max | ∆ k 1 |, | ∆ k 2 | (4.45) ǫ k = ǫ 2 1 -ǫ 2 max | ∆ k 1 |, | ∆ k 2 | | ∆ k 1 | + | ∆ k 2 | +P k (4.46)
Dans AVBP les valeurs utilisées sont:

Ψ 0 = 2 × 10 -2 δ = 1 × 10 -2 ǫ 1 = 1 × 10 -2 ǫ 2 = 0.95 ǫ 3 = 0.5 (4.47)
Pour la viscosité artificielle utilisée il y a deux opérateurs: un de second ordre qui opère comme une viscosité "classique". Il lisse les gradients, et introduit de la dissipation.

Le quatrième ordre s'utilise comme un bi-Laplacien et assure le contrôle des hautes fréquences (voir les expressions dans partie qui suit en anglais).

Simulation de Grandes Echelles (LES)

Dans cette méthodes seules les grandes échelles énergétiques sont calculées et les effets de petites échelles sont modélisés. Elle permet de faire les calculs à des plus grands nombres de Reynolds que la Simulation de Grandes Echelles avec des coût plus faibles.

C'est une approche immédiate (voir Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF]) en comparaison à la simulation statistique calssique(RANS). Elle permet le calcul des écoulements instationnaires et peut donc être utilisée pour l'aéroacoustique, l'aéroélasticité ou le contrôle des écoulements.

Historiquement cette méthode a été utilisée en météorologie ( Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF], Lilly [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF], Deardorff [START_REF] Deardorff | Three-dimensional numerical study of the height and mean structure of a heated planetary boundary layer[END_REF] Mason [START_REF] Mason | Large-Eddy Simulation: A critical review of the technique[END_REF]). Et ensuite elle a été appliquée à d'autres cas de plus en plus complexes ( Kraichnan [START_REF] Kraichnan | Eddy viscosity in two and three dimensions[END_REF], Chasnov [START_REF] Chasnov | Simulation of the Kolmogorov inertial subrange using an improved subgrid model[END_REF], Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF], Schumann [START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF],

Moin and Kim [START_REF] Moin | Numerical investigation of turbulent channel flow[END_REF], Piomelli [START_REF] Piomelli | High Reynolds number calculations using the dynamic subgridscale stress model[END_REF], Akselvoll and Moin [START_REF] Akselvoll | Large-eddy simulation of turbulent confined coannular jets[END_REF], Haworth and Jansen [START_REF] Haworth | Large-eddy simulation on unstructured deforming meshes: towards reciprocating IC engines[END_REF])

. Les équations pour la LES sont obtenues en appliquant un filtre aux équations de Navier-Stokes compressibles. On résoud donc des équations filtrées. Pour tenir compte des échelles sous-mailles, on introduit des modèles. La résolution permet ainsi de déterminer le détacheemnt des tourbillons ou l'acoustique par exemple (voir Poinsot

& Veynante [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]).

Equations résolues en LES

Dans le cadre de ce travail l'écoulement est non-réactif. On utilise dans les équations filtrées relatives à cet écoulement.

Filtrage

Pour séparer les grandes échelles des petites, un filtre passe-bas G △ , est appliqué aux équations du mouvement(voir Leonard [START_REF] Leonard | Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows[END_REF]). Il s'agit d'un produit de convolution entre toute fonction, f , avec la fonction filtre G △ :

f = f (x ′ )G △ (x -x ′ )dx ′ (4.48)
La quantité filtrée, f , représente les structures de grandes échelles alors que les structures de tailles plus petit que la taille du filtre, △, sont contenues dans l'écoulement résiduel,

f ′ : f ′ = f -f (4.49)
Pour l'écoulement à densité variable ρ, une moyenne pondérée par la masse f (Favre [START_REF] Favre | Statistical equations of turbulent gases[END_REF] pour éviter l'apparition d'autres termes inconnus.

f = ρf ρ (4.50)

Numerical simulation and LES models

Navier-Stokes filtrées sans réactif

Les équations s'écrivent comme suit:

∂ ρ ∂t + ∂ ∂x j (ρũ i ) = 0 (4.51) ∂ ρũ i ∂t + ∂ ∂x j (ρũ i ũj ) = - ∂ ∂x j P δ ij -τij -τ t ij (4.52) ∂ ρE ∂t + ∂ ∂x j (ρEu j ) = - ∂ ∂x j u i (P δ ij -τ ij ) +

Termes non-visqueux

Les trois parties du flux non-viqueux sont:

f I =          ρũ 2 + P ρũṽ ρũ w ρ Ẽ ũ + P u ρk ũ          , ḡI =          ρũṽ ρṽ 2 + P ρṽ w ρ Ẽṽ + P v ρk ṽ          , hI =          ρũ w ρṽ w ρ w2 + P ρ Ẽ w + P w ρk w          (4.60)

Termes visqueux filtrés

Les composantes du tenseur des flux visqueux ont la forme suivante: Pour le cas non-réactif

FV =          -τ xx -τ xy -τ xz -(uτ xx + vτ xy + wτ xz ) + q x J x,k          , (4.61) ḠV =          -τ xy -τ yy -τ yz -(uτ xy + vτ yy + wτ yz ) + q y J y,k          , (4.62)
HV =          -τ xz -τ yz -τ zz -(uτ xz + vτ yz + wτ zz ) + q z J z,k          (4.63) avec τ ij = 2µ S ij - 1 3 δ ij S ll (4.64) approximation : τ ij ≈ 2µ Sij - 1 3 δ ij Sij (4.65) avec : Sij = 1 2 ∂ ũj ∂x i + ∂ ũi ∂x j (4.66) µ ≈ µ( T ) ( 4
J i,k = -ρ D k W k W ∂X k ∂x i -Y k V c i (4.71) approximation : J i,k ≈ -ρ D k W k W ∂ X k ∂x i -Y k V i c (4.72) avec : V i c = N k=1 D k W k W ∂ X k ∂x i (4.73) D k ≈ µ ρSc k (4.74)
Le flux de chaleur filtré q i

q i = -λ ∂T ∂x i + N k=1 J i,k h s,k (4.75)
approximation :

q i ≈ -λ ∂ T ∂x i + N k=1 J i,k hs,k (4.76)
avec : λ ≈ µC p ( T ) P r (4.77)

Termes de sous-maille

Les flux de la turbulence de sous-maille sont:

F t =          -τ xx t -τ xy t -τ xz t q x t J x,k t          , G t =          -τ xy t -τ yy t -τ yz t q y t J y,k t          , H t =          -τ xz t -τ yz t -τ zz t q z t J z,k t          (4.78) avec τ ij t = -ρ ( u i u j -u i u j ) (4.79) approximation :τ t ij = 2ρν t S ij - 1 3 δ ij S ll (4.80) avec : S ij = 1 2 ∂ u i ∂x j + ∂ u j ∂x i - 1 3 ∂ u k ∂x k δ ij (4.81)
Dans l'équation 4.81, τ t ij est le tenseur de sous-maille, ν t est la viscosité sous-maille, et S ij est le tenseur de déformation résolu. La modélisation de ν t est expliqué dans la section 4

J t i,k = ρ u i Y k -u i Y k (4.82) modélisé comme :J ij t = -ρ D t k W k W ∂ Xk ∂x i -Ỹk Ṽ c
q i t = ρ u i E -u i E (4.86)
modélisé comme :

q i t = -λ ∂ T ∂x i + k J i,k t hs,k (4.87)
avec : λ t = µ t c P P r t (4.88)

Modèle de sous-maille

Le modèle de sous-maille (SGS)s'écrit:

τ ij t = -ρ ( u i u j -u i u j ) (4.89) = 2 ρ ν t S ij - 1 3 τ ll t δ ij (4.90)
Le modèle de Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] repose sur la viscosité turbulente :

ν t = (C S △) 2 2 S ij S ij (4.91)
avec △ du filtre (△ = 3 √ △x△y△z), C S est la constante du modèle égale à 0.18 mais est comprise entre 0.1 et 0.18 suivant le cas (Lilly [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF], Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF]). Le modèle dynamique de Germano repose sur la détermination dynamique de C S comme fonction d'espace et du temps.

ν t = (C S D △) 2 2 S ij S ij (4.92) où C 2 S D = 1 2 M ij M ij L ij L ij (4.93)
4. Numerical simulation and LES models et

M ij = ∆2 2 < S ij >< S ij > L ij =< u i >< u j > -< u i u j > (4.94)

WALE

Pour obtenir la viscosité turblente proche d'une paroi, l'amortissement de Van Driest est souvent utilisé [START_REF] Van Driest | On the turbulent flow near a wall[END_REF]. Une autre voie est WALE (Wall-Adapting Local Eddy-viscosity)

proposée par Nicoud & Ducros [START_REF] Ducros | Wall-adapting local eddy-viscosity models for simulations in complex geometries[END_REF][START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]:

s d ij = 1 2 g2 ij + g2 ji + 1 3 g2 kk δ ij (4.95) ν t = (C w △) 2 s d ij s d ij 3 2 S ij S ij 5 2 + s d ij s d ij 5 4
(4.96)

avec g2 ij = ∂u i ∂x k ∂u k ∂x j , Cw = 0.4929.

Conditions aux limites

Les conditions aux limites jouent un rôle important en simulation numérique. Et surtout pour l'aspect propagatif contenu dans les équations résolues par AVBP (voir Schöfeld &

Rudgyard [START_REF] Schöfeld | Steady and unsteady flows simulations using the hybrid flow solver avbp[END_REF] et Poinsot & Veynante [START_REF] Poinsot | Combustion numérique-cours de DEA[END_REF]). Dans AVBP appliquer les conditions aux limites est équivalent à trouver les résidus R aux limites. Ce dernier esst obtenu en utilisant le schéma d'intégration de Runge-Kutta.

U n+1 = U n -R∆t (4.97)

Comme le code est explicite seule la solution d'indice n est utilisée. Pour corriger le résidu aux limites deux méthodes sont utilisées: sans relation caractéristique ou avec (méthode (NSCBC)(voir Moureau et al [102] et Poinsot & Lele [113]). Dans ce dernier cas contrairement au premier aucune valeur n'est imposée seules les amplitudes des ondes sont spécifées. A chaque pas de temps, on utilise l'approche normale qui utilise les relations caractéristiques au flux normal des dérivées du résidu. La formulation des conditions aux

limites a fait l'objet d'études poussées ( Nicoud [103]).

Les relations caractéristiques

Le traitement des conditions aux limites dans AVBP est résumé sur la figure 4.4.

L'avancement explicite en temps dans AVBP donne U n+1 pred :

∂U = U n+1 pred -U n = -R P ∆t (4.98)
Le résidu total R P peut être divisé en deux parties :

∂U = -∆t R P BC + R P U (4.

The AVBP solver

The chapter is totally devoted to explain the numerical solver which is used to simulate the cavity flow. The chapter carries details about the numerical methods, artificial viscosity, derivation of governing equations for LES, models used in large eddy simulation (LES) and usage of characteristic boundary conditions in the simulation.

The AVBP (A Very Big Project) was historically motivated by the idea of building a modern software tool for Computational Fluid Dynamics (CFD) of high flexibility, efficiency and modularity. It was started at CERFACS in January 1993 as an initiative of Michael Rudgyard and Thilo Schönfeld. The aim was to create an unstructured solver capable of handling grids of any cell type. The use of these so-called hybrid grids is motivated by the efficiency of unstructured grid generation, the accuracy of the computational results (using regular structured elements) and the ease of mesh adaptation. The philosophy of building AVBP upon software libraries was adopted to best meet the modularity requirement.

AVBP is a parallel CFD code that solves the laminar and turbulent compressible Navier-Stokes equations in two and three space dimensions. Steady state or unsteady flows may be simulated. For the prediction of unsteady turbulence, various Large-Eddy Simulation (LES) subgrid scale models have been implemented. AVBP was initially conceived for primarily stationary external flows for aerodynamics applications. Since the mid-nineties the emphasis of applications is on the modeling of unsteady turbulent flows (with and without chemical reactions) for mainly internal flow configurations. These activities are partially related to the rising importance of the understanding of the flow structure and mechanisms leading to turbulence. The prediction of these unsteady turbulent flows is based on the LES approach which has emerged as a prospective technique for problems associated with time dependent phenomena and coherent eddy structures.

An Arrhenius law reduced chemistry model allows investigating combustion for complex configurations.

The handling of unstructured or hybrid grids is one key feature of AVBP. With the use of these hybrid grids, where a combination of several elements of different type is used in the framework of the same mesh, the advantages of structured and unstructured grid methodologies are combined in terms of gridding flexibility and solution accuracy.In order to handle arbitrary hybrid grids, the data structure of AVBP employs a cell-vertex finite-volume approximation. The basic numerical methods are based on a Lax-Wendroff [START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF] or a finite-Element type low-dissipation Taylor-Galerkin(Donea [START_REF] Donea | Taylor-Galerkin method for convective transport problems[END_REF] method of AVBP (Rudgyard [START_REF] Rudgyard | Cell vertex methods for compressible gas flows[END_REF][START_REF] Rudgyard | Cell vertex methods for steady inviscid flow[END_REF]), the discrete values of the conserved variables are stored at the cell vertices (or grid nodes), while conservation relations are applied to the grid (or primary) cells. The advantages of using such a discretisation are:

• The native capability of handling unstructured hybrid meshes.

• An easy and efficient parallelisation.

• Increased accuracy without an important additional cost, can be obtained by using the same spatial differential operators in a finite element framework (see Section 4.2.6).

In the cell-vertex method employed within AVBP both solution and coordinate vectors are stored at the nodes of the grid. However, most of the operations are done on the elements and often a transfer from the cell vertices (the nodes) to the cell centers is required. This collecting of the nodal information to temporary arrays that contain the information of the vertices for an element is done in a so-called data gather operation Nomenclature: In the rest of the section the following subscripts are used:

• i ∈ [1, N node ]
is the index used for the global node numbering and the nodal values.

• j ∈ [1, N cell ] is used for the cell numbering.

• k ∈ [1, n v (Ω j )]
is the local numbering of the vertices of a cell Ω j , with n v (Ω j ) the number of vertices of the cell Ω j .

• Ω j is the index used to design a value at the centre or associated with the j-th cell.

Numerical method

• R i is the global nodal residual.

• R Ω j is the global cell residual.

• R i|Ω j is the part of the residual of element j to be scattered to node i.

Weighted Cell Residual Approach

For the description of the weighted cell-residual approach the laminar Navier-Stokes equations are considered in their conservative formulation:

∂U ∂t + ∇ • F = 0 (4.102)
where U is the vector of conserved variables and F is the corresponding flux tensor. For convenience, the latter is divided into an inviscid and a viscous part,

F = F I (U) + F V U, ∇U
. The spatial terms of the equations are then approximated in each control volume Ω j to give the residual

R Ω j = 1 V Ω j ∂Ω j F • ndS (4.103)
where V Ω j is volume of the control volume and ∂Ω j denotes the boundary of Ω j with normal n.

This cell-vertex approximation is readily applicable to arbitrary cell types and is hence straight-forward to apply for hybrid grids. The residual 4.103 is first computed for each element by making use of a simple integration rule applied to the faces. For triangular faces, a straightforward mid-point rule is used, which is equivalent to the assumption that the individual components of the flux vary linearly on these faces. For quadrilateral faces, where the nodes may not be co-planar, in order to ensure that the integration is exact for arbitrary elements if the flux functions do indeed vary linearly, each face is divided into triangles and then integrated over the individual triangles. The flux value is then obtained from the average of four triangles (two divisions along the two diagonals). This so-called 'linear preservation property 'plays an important part in the algorithm for ensuring that accuracy is not lost on irregular meshes. Computationally, it is useful to write the discrete integration of equation 4.103 over an arbitrary cell as

R Ω j = 1 N d V Ω j i∈Ω j F i .d S i (4.104)
where F i is an approximation of F at the nodes, N d represents the number of space dimensions and {i ∈ Ω j } are the vertices of the cell. In this formulation the geometrical information has been factored into terms d S i that are associated with individual nodes 4. Numerical simulation and LES models of the cell but not faces; d S i is merely the average of the area-weighted normals for triangulated faces with a common node i, i ∈ Ω j . Note, that for consistency one has i∈Ω j d S i = 0. A linear preserving approximation of the divergence operator is obtained if the volume V Ω j is defined consistently as

V Ω j = 1 N 2 d i∈Ω j x i .d S i (4.105) since ∇ • x = N d .
Once the cell residuals are calculated, one may then define the semi-discrete scheme

dU k dt = - 1 V k j|k∈Ω j D k Ω j V Ω j R Ω j (4.106)
where D k Ω j is a distribution matrix that weights the cell residual from cell centre Ω j to node k (scatter operation), and V k is a control volume associated with each node.

Conservation is guaranteed if k∈Ω

j D k Ω j = I.
In the present context, equation 4.106 is solved to obtain the steady-state solution using explicit Euler or Runge-Kutta timestepping.

The family of schemes of interest makes use of the following definition of the distribution matrix:

D k Ω j = 1 n n I + C δt Ω j V Ω j A Ω j • d S k (4.107)
where n n is the number of nodes of Ω j , δt Ω j is the cell 'time-step'and A is the Jacobian of the flux tensor. The simplest 'central difference'scheme is obtained by choosing C = 0 and is neutrally stable when combined with Runge-Kutta time-stepping. A Lax-Wendroff type scheme may also be formulated in which case C is chosen to be a constant that depends on the number of space dimensions and the type of cells used-it may be shown that this takes the simple form C = n 2 v 2 N d . If one replaces the cell 'timestep'δt Ω j by a matrix ΦΩ j with suitable properties, one may also obtain an SUPG-like scheme (for Streamwise Upwind Petrov-Galerkin) from Brooks & Hugues [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] which has slightly better convergence and shock-capturing behaviour, however, at some extra computational cost.

Computation of gradients

In order to recover the nodal values of the gradients ∇U a cell approximation ∇U Ω j is first calculated and then distributed to the nodes. The cell-based gradient is defined in a manner similar to the divergence equation 4.104 so as to be transparent to linear solution variations:

∂U ∂x C ≈ 1 V C ∂Ω C U • n∂S (4.108)
4.2. Numerical method which leads to the approximation:

∇U Ω j = 1 V Ω j i∈Ω j U i d S i (4.109)
A nodal approximation of the gradient is then obtained using of a volume-weighted average of the cell-based gradients:

∇U k = 1 V Ω k j|k∈Ω j V j ∇U Ω j (4.110)

Computation of time step

Temporal discretisation is explicit for all numerical schemes in AVBP. The practical implementation of this kind of approach is relatively straightforward and the computational cost per iteration is small. The main drawback of explicit codes is that the time step ∆t is limited for stability reasons:

∆t < CFL min (∆x) max | u | +a ∞ (4.111)
where u is the propagation speed of a perturbation in the flow, a ∞ is the sound speed, ∆x is the mesh size and CFL is the Courant-Friedrichs-Lewy number. The CFL value required for stability changes slightly depending on the scheme adopted. In AVBP, the CFL value is fixed to 0.7.

The Lax-Wendroff scheme

The main convective schemes are a finite volume Lax-Wendroff type scheme (LW) from

Lax & Wendroff [START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF] and a finite element two-step Taylor-Galerkin scheme (TTGC) from Colin & Rudgyard [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF]. These two schemes are respectively 2 nd and 3 rd order in time and space. The diffusive scheme is a typical 2 nd order compact scheme. Element types handled by AVBP are triangles and quadrangles in 2D and tetrahedrons, prisms, pyramids and hexahedrons in 3D. The time integration is fully explicit to maximise accuracy.

The form of the distribution matrix D i|Ω j (see equation 4.107) determines the different numerical schemes available in AVBP. In the following D i|Ω j is derived for the Lax-Wendroff scheme [START_REF] Lax | Systems of conservation laws[END_REF][START_REF] Lax | Difference schemes for hyperbolic equations with high order of accuracy[END_REF]. The Lax-Wendroff scheme (second order accurate in space and time) is based on a Taylor expansion in time of the solution U .

U n+1 = U n + ∆t ∂U ∂t n + 1 2 ∆t 2 ∂ 2 U ∂t 2 n + O(∆t 3 ) (4.112)

Numerical simulation and LES models

Considering the conservative formulation of laminar Navier-Stokes equation

∂U ∂t + ∇ • F = 0 (4.113)
the first temporal derivative can be expressed as:

∂U ∂t = -∇ • F (4.114)
In a similar manner, the second derivative can be recast as:

∂ 2 U ∂t 2 = ∂ ∂t (-∇ • F) = -∇ • ∂F ∂t = -∇ • A ∂U ∂t = ∇ • [A (∇ • F)] (4.115)
assuming that temporal and spatial derivatives can be exchanges and defining A = ∂F ∂U as Jacobian matrix. Hence, substituting equations 4.114, 4.115 into equation 4.112, the solution a time n + 1 can be written as:

U n+1 = U n -∆t ∇ • F - 1 2 ∆t ∇ • [A (∇ • F)] -O ∆t 2 (4.116)
In discrete form, remembering the basic principle of the cell-vertex approach, the nodal residual R i is obtained by summing the contributions of all the surrounding elements.

The value is then scaled by the nodal volume V i :

R i = 1 V i j|i∈Ω j R i|Ω j (4.117)
The residual contribution to node i of element j can be written as:

R i|Ω j = R Ω j V Ω j n v(Ω j ) -LW i|Ω j (4.118)
The first term in equation 4.118 is the cell residual computed as in equation 4.104. It is weighted by the volume of the cell divided by the number of vertices of the element.

The LW i|Ω j term is computed on the dual cell C i taking advantage of the Green-Gauss theorem:

LW i|Ω j = 1 2 ∆t Ω j T C i ∇ • [A (∇ • F )] dV = 1 2 ∆t ∂C i ∇ • [A (∇ • F )] dS (4.119)
This term is then discretised to give:

LW i|Ω j ≃ 1 2 ∆t [A (∇ • F )] Ω j • S i|Ω j n d (4.120)
where S i|Ω j is the normal associated with node i and cell j it is computed according to

R i|Ω j = I - ∆t 2 n d n v(Ω j ) V Ω j A Ω j • S i|Ω j R Ω j V Ω j n v(Ω j ) (4.121)
Recalling now equation 4.118 the distribution matrix takes the form:

D i|Ω j = 1 n v(Ω j ) I - ∆t 2n d n v(Ω j ) V Ω j A Ω j • S i|Ω j (4.122)

The TTGC numerical scheme

It is nearly impossible to develop schemes of higher order (in space) on unstructured meshes in a finite volume context. The cell-vertex formulation can be extended to a finite element approach, where higher order schemes are possible. Taylor-Galerkin (TG) schemes were first given by Donea [START_REF] Donea | Taylor-Galerkin method for convective transport problems[END_REF][START_REF] Donea | An analysis of time discretization in the nite element solution of hyperbolic problems[END_REF] coupling a Taylor expansion in time and a Galerkin discretisation in space. Colin and Rudgyard [START_REF] Colin | Development of high-order Taylor-Galerkin schemes for LES[END_REF] developed a two-step Taylor-Galerkin scheme (TTGC) that is third-order in space and time.

Ũn = U n + α∆t ∂U ∂t n + β∆t 2 ∂ 2 U ∂t 2 n (4.123) U n+1 = U n + ∆t ∂ Ũ ∂t n + γ∆t 2 ∂ 2 Ũ ∂t 2 n (4.124) α = 1 2 -γ and β = 1 6 (4.125)
first and second temporal derivatives can be replaced as done for the Lax-Wendroff scheme (see equation 4.114 and 4.115) giving:

Ũn = U n -α∆t∇ • F n + β∆t 2 ∇ • [A (∇ • F n )] (4.126) U n+1 = U n -∆t∇ • F n + γ∆t 2 ∇ • [A (∇ • F n )] (4.127)
Multiplying these equations by a set of linear test functions φ i ("redskin tent"functions) and integrating them over the computational domain Ω, leads to the following weak formulation:

Ω Rn φ i dV = -αL i (U n ) -β∆tLL i (U n ) (4.128) Ω R n+1 φ i dV = -L i Ũn -γ∆tLL i (U n ) (4.129)
4. Numerical simulation and LES models with

Rn = Ũn -U n ∆t , R n+1 = U n+1 -U n ∆t (4.130)
and

L i (U) = Ω ∇ • F (U n ) φ i dV (4.131) LL i (U) = Ω A (∇ • F (U n )) ∇φ i dV LL 0 i (U n ) - ∂Ω φ i A (∇ • F (U n )) dS BT i (U n ) (4.132)
The LL i term can be split by performing an integration by parts assuming the surface normal dS external. The first contribution LL 0 i (U n ) is integrated over all the computational domain while the second one, BT i (U n ), is non zero only at boundaries. It should be noticed that the LL i term involves second spatial derivatives (like the LW i term, see for example equation 4.120), that are not expected when dealing with convection problems. The Galerkin method is then applied to the flux divergence and to residuals.

Hence, they can be expressed as a sum of linear shape-functions (same functions as the test-functions used to derive the weak formulation), leading to:

R n = k R n k φ k (4.133) ∇ • F = k F k ∇φ k (4.134)
where F k is the discrete flux at each point of computational domain. With this choice of shape functions, the residuals are recast as:

Ω Rn φ i dV = k Ω φ i φ k dV Rn k = k M ik Rn k (4.135)
denoting M ik as the components of the so-called mass matrix which, in AVBP, is inverted locally by an iterative Jacobi method.

In the spatial discretisation, the contributions of integrals in equation 4.131 and 4.132

to node i come only from the adjacent cells.

L i (U n ) = j|i∈Ω j L i (U n ) Ω j (4.136) LL i (U n ) = j|i∈Ω j LL i (U n ) Ω j (4.137)

Numerical method

Taking advantage of equations 4.133, L i (U n ) Ω j and LL i (U n ) Ω j can be written as:

L i (U n ) Ω j = k|k∈Ω j F n k Ω j φ i ∇φ k dV (4.138) LL i (U n ) Ω j = A n Ω j k|k∈Ω j F n k Ω j ∇φ • ∇φ k dV -A n Ω j k|∈∂Ω j ∩∂Ω F n k ∂Ω j ∩∂Ω φ i ∇φ k dS (4.139)
For triangular and tetrahedron elements the gradient of the shape function is constant1 over each element and the integral of φ i takes a simple form(see Colin & Rudgyard [19]).

∇φ k = - S k n d V Ω j (4.140) Ω k φ k dV = V Ω j n v(Ω j ) ∀k ∈ Ω j (4.141)
Substituting relations 4.140 and 4.141 in equation 4.138 yields:

L i (U n ) Ω j = k|k∈Ω j F n k ∇φ k Ω j φ i dV = (∇ • F n ) Ω j Ω j φ i dV = R n Ω j V Ω j n v(Ω j ) (4.142)
Applying the same procedure to the first term of equation 4.139 leads to:

LL 0 i (U n ) Ω j = A n Ω j k|k∈Ω j F n k ∇φ k • Ω j ∇φ i dV = - 1 n d A n Ω j R n Ω j • S i|Ω j (4.143)
These two operators are therefore equivalent to the ones encountered in the cell-vertex finite volume discretisation (see equation 4.121). The scaling for the nodal volume does not appear explicitly in this derivation but it is taken into account in the mass matrix.

A more complete study of numerical schemes available in AVBP can be found from the Chapter 5 of Lamarque [START_REF] Lamarque | Schémas numériques et conditions limites pour la simulation aux grandes échelles de la combustion diphasique dans les foyers d'hélicoptère[END_REF] thesis.

It has a computational cost of approximately 2.5 times Lax-Wendroff (which is slightly less than the three-step Runge-Kutta). Achieving higher order in space is par-4. Numerical simulation and LES models ticularly useful for three-dimensional, unsteady simulations since it provides a much better accuracy on meshes already used for second-order simulations.

Diffusion terms

The

V i : u n+1 i -u n i ∆t = 1 V i ∇. (ν∇u) | i (4.145)
If it is associated to a FE scheme TTGC, the mass matrix is applied to the operator:

u n+1 i -u n i ∆t = M -1 ∇. (ν∇u) | i (4.146)
It is important to remark that a FV or FE convection schemes can be associated with any FV or FE diffusion scheme.

Artificial Viscosity

The numerical discretisation methods in AVBP are spatially centered. These types of schemes are known to be naturally subject to small-scale oscillations in the vicinity of [START_REF] Jameson | Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes[END_REF] and the "Colin-sensor"ζ C Ω j from Colin [START_REF] Colin | Simulations aux grandes échelles de la combustion turbulente prémélangé dans les statoréacteurs[END_REF].

Jameson cell sensor

For every cell Ω j , the Jameson cell-sensor ζ J Ω j is the maximum over all cell vertices of the Jameson vertex-sensor ζ k :

ζ J Ω j = max k∈Ω j ζ J k (4.147)
The Jameson vertex-sensor for scalar quantity, for pressure P is :

ζ J k = | ∆ k 1 -∆ k 2 | | ∆ k 1 | + | ∆ k 2 | + | P k | (4.148)
where the ∆ k 1 and ∆ k 2 functions are defined as:

∆ k 1 = P Ω j -P k ∆ k 2 = ∇P k • x Ω j -x k (4.149)
where a k subscript denotes cell-vertex values while Ω j is the subscript for cell-averaged values. ∇P k is the gradient of P at node k as computed in AVBP. ∆ k 1 measures the variation of p inside the cell Ω j (using only quantities defined on this cell). ∆ k 2 is an estimation of the same variation but on a wider stencil (using all the neighbouring cell of the node k).

The Colin sensor

The Jameson sensor is smooth and was initially derived for steady-state computations.

But for most unsteady turbulent computations it is however necessary to have a sharper sensor, which is very small when the flow is sufficiently resolved, and which is nearly 4. Numerical simulation and LES models maximum when a certain level of non-linearities occurs. The exact definition of the Colin sensor is:

ζ C Ω j = 1 2 1 + tanh Ψ -Ψ 0 δ - 1 2 1 + tanh -Ψ 0 δ (4.150) with Ψ = max k∈Ω j 0, ∆ k | ∆ k | +ǫ 1 P k ζ J k (4.151) ∆ k = | ∆ k 1 -∆ k 2 | -ǫ k max | ∆ k 1 |, | ∆ k 2 | (4.152)
ǫ k = ǫ 2 1 -ǫ 2 max | ∆ k 1 |, | ∆ k 2 | | ∆ k 1 | + | ∆ k 2 | +P k (4.153)
The numerical values used in AVBP are

Ψ 0 = 2 × 10 -2 δ = 1 × 10 -2 ǫ 1 = 1 × 10 -2 ǫ 2 = 0.95 ǫ 3 = 0.5 (4.

154)

The operators

There are two artificial viscosity operators in AVBP: a 2 nd order operator and a 4 th order operator. The 2 nd order operator acts like a "classical"viscosity. It smoothes gradients, and introduces artificial dissipation. It is thus associated to a sensor which determines where it must be applied. Doing this, the numerical scheme keeps its order of convergence in the zones where the sensor is inactive, while ensuring stability and robustness in the critical regions. It smooths any physical gradient. The 4 th order operator is a less common operator. It acts as a bi-Laplacian and is used to control spurious highfrequency wiggles. The user-defined parameters smu2 for 2 nd order operator and smu4 for 4 th order operator are selected before the start of simulation.

The 2 nd order operator A cell contribution of the 2 nd order artificial viscosity is first computed on each vertex of the cell Ω j :

R k∈Ω j = - 1 N v V Ω j ∆t Ω j smu2 ζ Ω j w Ω j -w k (4.155)
The nodal residual is then found by adding the surrounding cells contributions:

dw k = j R k∈Ω j (4.156)
For example, on a 1D uniform mesh, of mesh size ∆x, and for ζ Ω j = ζ:

dw k = - smu2 2 ζ ∆x ∆t (w k-1 -2w k + w k+1 ) (4.157)
4.2. Numerical method which can be interpreted as:

dw k = -ν AV (∆ k,∆x w) dx (4.158)
with:

ν AV = smu2 ζ∆x 2 2∆t = smu2 ζ∆x | u + c | 2 CF L and ∆ F D k,∆x w = w k-1 -2w k + w k+1 ∆x 2 (4.159)
where ∆ F D k,∆x is exactly the classical FD Laplacian operator evaluated at k and of size ∆x. This shows that ν AV can be seen as an artificial viscosity (it has the same units as a physical viscosity), which is controlled by the user-defined non-dimensional parameter smu2.

The 4 th order operator

The technique used for the 4 th order operator is identical to the technique of the 2 nd order operator. A cell contribution is first computed on each vertex:

R k∈Ω j = - 1 N v V Ω j ∆t Ω j smu4 ∇w Ω j • x Ω j -x k -w Ω j -w k (4.160)
The nodal value is then found by adding every surrounding cells contributions:

dw k = j R k∈Ω j (4.161)
For example, on a 1D uniform mesh, of mesh size ∆x, this yields:

R k∈Ω lef t = smu4 2 ∆x ∆t 1 2 w k -w k-2 2 ∆x + w k+1 -w k-1 2 ∆x • -∆x 2 - smu4 2
∆x ∆t

w k-1 + w k 2 -w k (4.162) R k∈Ω right = smu4 2 ∆x ∆t 1 2 w k+1 -w k-1 2 ∆x + w k+2 -w k 2 ∆x • -∆x 2 - smu4 2
∆x ∆t

w k + w k+1 2 -w k (4.163)
Adding these two contributions delivers:

dw k = smu4 ∆x 16∆t (w k-2 -4w k-1 + 6w k -4w k+1 + w k+2 ) (4.164)
4. Numerical simulation and LES models which can be interpreted:

dw k = κ AV ∆∆ F D k,∆x w dx (4.165) with κ AV = smu4 ∆x 4 16 ∆t = smu4 ∆x 3 | u + c | 16 CFL and ∆∆ F D k,∆x w = w k-2 -4w k-1 + 6w k -4w k+1 + w k+2 ∆x 4 (4.166)
where ∆∆ F D k,∆x is exactly the classical FD bi-Laplacian operator evaluated at k and of size ∆x. This shows that κ AV can be seen as an artificial 4 th order hyper-viscosity, which is controlled by the user-dened non-dimensional parameter smu4.

Artificial viscosity model

In "Colin"model, three sensors are used in conjunction. The first one is based on total energy, the second one is based on species densities, and the last one is the maximum of the two previous.

ζ COL E = ζ C Ω j (ρE) , ζ COL Y = max k=1,neqs ζ C Ω j (ρ k ) , and ζ COL max = max ζ COL E , ζ COL Y (4.167)
The sensor used in the 2 nd order operator is also the maximum sensor ζ COL max . The 2 nd and 4 th order operators are then applied on each species. This model is particularly dedicated to LES of reactive flows. The lack of 4 th order artificial viscosity on momentum allows to keep many small scale structures, while damping the wiggles on energy and species.

Schönfeld-Lartigue-Kaufmann (SLK model) is an improvement of the "Colin"model.

It is noticed that applying no 4 th order artificial viscosity at all on momentum can yet lead to non-negligible wiggles in some cases (often depending on the quality of the mesh) and this had to be avoided. This new "SLK"model is very similar to the "Colin"model, except that instead of setting the modified 4 th order coefficient to zero for momentum equation, it is set to 10% of the value used for the other equations. This model is very well suited for computations on poor quality meshes that exhibit velocity wiggles with the standard "Colin"model. Colin and SLK model are used in the present work. Values of the user defined parameters smu2 and smu4 are given in next chapter along with the details of the test cases.

Large Eddy Simulation

The LES can be considered as a midpoint between the RANS approach in which all the turbulent scales are modeled and the DNS in which all the turbulent scales are com-4.4. Governing equations for LES puted. In a LES simulation, only the largest scales -the scales that contain most of the energy -are computed; the effect of the smallest scales are modeled. The smallest scales have a more predictable behaviour and should be easier to model. LES has been highly developed by the engineering computational fluid dynamics community since its inception in 1970. Large-eddy simulation (LES) resolves only the dynamically important flow scales and models the effects of smaller scales whereas DNS resolves all flow scales as mentioned in the section 2.4. Because of its high computational cost, usage of DNS is limited to simple flow configurations at low to moderate Reynolds numbers. Large Eddy Simulation(see Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF]) is nowadays recognized as immediate approach in comparisons to the more classical Reynolds Averaged Navier-Stokes (RANS) methodologies. LES gives access to the dominant unsteady motion so that it can be used to study aeroacoustics, fluid-structure interaction or the control of turbulence by an appropriate unsteady forcing. Much of the pioneering work on LES (e.g., Smagorinsky [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF],

Lilly [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF], Deardorff [START_REF] Deardorff | Three-dimensional numerical study of the height and mean structure of a heated planetary boundary layer[END_REF]) was motivated by meteorological applications, and atmospheric boundary layers remain a focus of LES activities (e.g., Mason [94]). The development and testing of LES methodologies have focused primarily on isotropic turbulence (e.g., Kraichnan [75], Chasnov [START_REF] Chasnov | Simulation of the Kolmogorov inertial subrange using an improved subgrid model[END_REF]), and on fully-developed turbulent channel flow (e.g., Deardorff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF], Schumann [START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF], Moin and Kim [START_REF] Moin | Numerical investigation of turbulent channel flow[END_REF], Piomelli [START_REF] Piomelli | High Reynolds number calculations using the dynamic subgridscale stress model[END_REF]). A primary goal of work in this area is to apply LES to flows in complex geometries that occur in engineering applications (e.g., Akselvoll and Moin [START_REF] Akselvoll | Large-eddy simulation of turbulent confined coannular jets[END_REF], Haworth and Jansen [START_REF] Haworth | Large-eddy simulation on unstructured deforming meshes: towards reciprocating IC engines[END_REF]). The derivation of the new governing equations is obtained by introducing operators to be applied to the set of compressible Navier-Stokes equations. Unclosed terms arise from these manipulations and models need to be supplied for the problem to be solved. In LES, the operator employed in the derivation is a spatially localized filter of given size △ to be applied to a single realisation of the studied flow. Resulting from this "spatial average"is a separation between the large (greater than the filter size) and small (smaller than the filter size) scales. The unclosed terms are in LES representative of the physics associated with the small structures (with high frequencies) present in the flow. Due to the filtering approach, LES allows a dynamic representation of the large scale motions whose contributions are critical in complex geometries. The LES predictions of complex turbulent flows are henceforth closer to the physics since large scale phenomena such as large vortex shedding and acoustic waves are embedded in the set of governing equations (see Poinsot & Veynante [112]). The basic idea of LES is to resolve (large) grid scales (GS), and to model (small) sub grid-scales(SGS).

Governing equations for LES

The test cases presented in this work will not be reactive and therefore this section discusses only the filtered equations solved by AVBP for a turbulent non-reacting flow.

Numerical simulation and LES models

With this intention, the filtering procedure is presented in Subsection 4.4.1. Subsection 4.4.2 describes the equations solved for LES of non-reacting flows. Then, the different terms of the flux tensor are presented. finally, different models of the subgrid stress tensor available in AVBP are described.

Filtering procedure

To separate large and small scales, a low-pass (in wavenumber) filter, G △ , is applied to the equations of motion (see Leonard [START_REF] Leonard | Energy Cascade in Large-Eddy Simulations of Turbulent Fluid Flows[END_REF]). Mathematically, it consists of a convolution of any quantity, f , with the filter function G △ :

f = f (x ′ )G △ (x -x ′ )dx ′ (4.168)
The resulting filtered quantity, f , represents the large-scale structures of the flow whereas all the structures of size smaller than the filter length, △, are contained in the residual field, f ′ :

f ′ = f -f (4.169)
To apply this filtering procedure to the instantaneous balance equations [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], the filter G △ (typical a box or a Gaussian filter [START_REF] Sagaut | Introduction à la simulation des grandes échelles, mathématiques et applications[END_REF]) must satisfy some properties which are: conservation of constants, linearity and commutation with temporal and spatial derivatives. The latter is satisfied only for homogeneous filters (i.e., grid meshes). For the sake of simplicity, this property is assumed hereafter.

For variable density ρ, a density-weighted filter f (Favre [START_REF] Favre | Statistical equations of turbulent gases[END_REF] averaging) is used, in order to avoid modeling of additional terms introduced by density fluctuations:

f = ρf ρ (4.170)

Filtering Navier-Stokes equations for non-reacting flows

The balance equations (mass, momentum, energy and species) for large-eddy simulations are obtained by filtering the instantaneous balance equations (see Poinsot &

Veynante [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]):

∂ ρ ∂t + ∂ ∂x j (ρũ i ) = 0 (4.171) ∂ ρũ i ∂t + ∂ ∂x j (ρũ i ũj ) = - ∂ ∂x j P δ ij -τij -τ t ij (4.172) ∂ ρE ∂t + ∂ ∂x j (ρEu j ) = - ∂ ∂x j u i (P δ ij -τ ij ) + qj + qt j + ωT + Qr (4.173) ∂ ρk Ỹk ∂t + ∂ ∂x j ρ Ỹk ũj = - ∂ ∂x j Jj,k + Jt j,k + ωk (4.174)
where ũi , Ẽi and Ỹk denote the filtered velocity vector, total energy per unit mass and species mass fractions, respectively. A repeated index implies summation over this index (Einstein's rule of summation). Note also that the index k is reserved for referring to the k th species and does not follow the summation rule (unless specifically mentioned).

Writing the vector of the filtered conservative variables as follows: Ū = (ρũ, ρṽ, ρ w, ρ Ẽ, ρ Ỹk )

, equations (4.173)-(4.174), can be expressed as:

∂ Ū ∂t + ∇ • F = s (4.175)
where s is the filtered source term and F is the flux tensor which can be divided in three The cut-off scale corresponds to the mesh size (implicit filtering). As usually done, we assume that the filter operator and the partial derivative commute.

Inviscid terms

The three spatial components of the inviscid flux tensor based on the filtered quantities are:

f I =          ρũ 2 + P ρũṽ ρũ w ρ Ẽ ũ + P u ρk ũ          , (4.180) ḡI =          ρũṽ ρṽ 2 + P ρṽ w ρ Ẽṽ + P v ρk ṽ          , ( 4

Filtered viscous terms

The components of the viscous flux tensor take the form:

FV =          -τ xx -τ xy -τ xz -(uτ xx + vτ xy + wτ xz ) + q x J x,k          , (4.183) ḠV =          -τ xy -τ yy -τ yz -(uτ xy + vτ yy + wτ yz ) + q y J y,k          , (4.184) HV =          -τ xz -τ yz -τ zz -(uτ xz + vτ yz + wτ zz ) + q z J z,k          (4.185)
filtering the balance equations leads to unclosed quantities which need to be modeled.

The filtered diffusion terms are (see Poinsot & Veynante [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF], Chapter 4):

Laminar filtered stress tensor τ ij For non-reacting flows: 4.197) approximation :

τ ij = 2µ S ij - 1 3 δ ij S ll (4.186) approximation : τ ij ≈ 2µ Sij - 1 3 δ ij Sij (4.
J i,k = -ρ D k W k W ∂X k ∂x i -Y k V c i (4.193) approximation : J i,k ≈ -ρ D k W k W ∂ X k ∂x i -Y k V i c (4.194) with : V i c = N k=1 D k W k W ∂ X k ∂x i (4.195) D k ≈ µ ρSc k (4.196) Filtered heat flux q i q i = -λ ∂T ∂x i + N k=1 J i,k h s,k (
q i ≈ -λ ∂ T ∂x i + N k=1 J i,k hs,k (4.198)
with : λ ≈ µC p ( T ) P r (4.199)

These forms assume that the spatial variations of molecular diffusion fluxes are negligible and can be modeled through simple gradient assumptions.

Subgrid-scale turbulent terms

The three components of the turbulent subgrid-scale flux take the form:

F t =          -τ xx t -τ xy t -τ xz t q x t J x,k t          , G t =          -τ xy t -τ yy t -τ yz t q y t J y,k t          , H t =          -τ xz t -τ yz t -τ zz t q z t J z,k t          (4.200)

Numerical simulation and LES models

As highlighted above, filtering the transport equations leads to a closure problem evidenced by the so called "subgrid-scale"(SGS) turbulent fluxes. For the system to be solved numerically, closures need to be supplied. Details on the closures are:

The Reynolds tensor τ ij t

τ ij t = -ρ ( u i u j -u i u j ) (4.201) modeled as : τ t ij = 2ρν t S ij - 1 3 δ ij S ll (4.202) with : S ij = 1 2 ∂ u i ∂x j + ∂ u j ∂x i - 1 3 ∂ u k ∂x k δ ij (4.203)
In equation 4.203, τ t ij is the SGS tensor, ν t is the SGS turbulent viscosity, and S ij is the resolved strain rate tensor. The modeling of ν t is explained in section 4.4.5

The subgrid scale diffusive species flux vector J t i,k

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

J ij t = -ρ D t k W k W ∂ Xk ∂x i -Ỹk Ṽ c,t i (4.205) with : Ṽ c,t i = N k=1 D t k W k W ∂ Xk ∂x i (4.206) D t k = ν t Sc t k (4.207)
The turbulent Schmidt number Sc t k = 1 is the same for all species and is fixed in the source code (like P r t ). Note also that having one turbulent Schmidt number for all the species does not imply, Ṽ c,t i = 0 because of the W k W term in equation 4.206.

The subgrid scale heat flux vector q i t

q i t = ρ u i E -u i E (4.208) modeled as :q i t = -λ ∂ T ∂x i + k J i,

Subgrid scale model

In the broader context of turbulence modelling including LES, Pope [START_REF] Pope | Turbulent Flows[END_REF] suggests five criteria: level of description, completeness, cost and ease of use, range of applicability

τ ij t = -ρ ( u i u j -u i u j ) (4.211) = 2 ρ ν t S ij - 1 3 τ ll t δ ij (4.212)
In the equation, τ ij t is the SGS tensor to be modelled, ν t is the SGS turbulent viscosity, u i is the Favre filtered velocity vector (compressible flows) and S ij is the resolved strain rate tensor. The models in AVBP differs only in the way the turbulent viscosity value ν t is calculated.

Smagorinsky's Model

The Smagorinsky's model [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF] was developed in the sixties by Smagorinsky. An eddy viscosity was supposed to take into account subgrid-scale dissipation through a Kolmogorov k -5/3 cascade. It was heavily tested for multiple flow configurations. This closure has the particularly of supplying the right amount of dissipation of kinetic energy in homogeneous isotropic turbulent flows. Smagorinsky's eddy viscosity is

ν t = (C S △) 2 2 S ij S ij (4.213)
where △ denotes the filter characteristic length(i.e △ = 3 √ △x△y△z), C S is the model constant set to 0.18 but can vary between 0.1 and 0.18 depending on the flow configuration. The constant C S may determined in isotropic turbulence which was performed by Lilly [START_REF] Lilly | The representation of small-scale turbulence in numerical simulation experiments[END_REF]. In the Smagorinsky model, the coefficient C S appears only in the product, hence decreasing △ is equivalent to decreasing C S . Locality is however lost and only 4. Numerical simulation and LES models global quantities are maintained in this model. The Smagorinsky model is known as being "too dissipative"in presence of a wall and transitioning flows are not suited for its use. More details can be found in Lesieur [START_REF] Lesieur | Large-Eddy Simulations of Turbulence[END_REF] and Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF].

Dynamic Smagorinsky's Model

The dynamic model was proposed by Germano et al [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF], with important modifications.

It is constructed to determine "dynamically"of the Smagorinsky model constant C S as a function of space and time. The extensions of this model were provided by Lilly [START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF] and Meneveau et al [START_REF] Meneveau | A Lagrangian dynamic subgrid-scale model of turbulence[END_REF]. The dynamic procedure has been most successful in remedying the standard Smagorinsky model's serious deficiencies in laminar flows, transitional flows and in the viscous near-wall region. The starting point for development of the dynamic model is a special case of the Germano identity [START_REF] Germano | Turbulence: The filtering approach[END_REF], relating the Leonard stress to a similar twice-filtered tensor and a second filtering of the usual SGS tensor. For example, details of construction of dynamic models were detailed by Sagaut [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF]. Here in this model, the eddy viscosity is given by

ν t = (C S D △) 2 2 S ij S ij (4.214)
where

C 2 S D = 1 2 M ij M ij L ij L ij (4.215)
the tensors from the previous expression are given by

M ij = ∆2 2 < S ij >< S ij > L ij =< u i >< u j > -< u i u j > (4.216)
L ij is called Germano's identity. The terms on the right hand side have to be modeled.

The "test"filter of characteristic length ∆ equal to the cubic root of the volume is introduced and defined by all the cells surrounding the cell of interest. Clipping and smoothing ensures non negative values for C S D . It can be shown that the dynamic model gives a zero subgrid-scale stress at the wall, where L ij vanishes, which is a great advantage with respect to the original Smagorinsky model; it also gives the proper asymptotic behaviour near the wall [START_REF] Lesieur | Large-Eddy Simulations of Turbulence[END_REF]. Simulations with this dynamic model were discussed by Lesieur and Métais [START_REF] Lesieur | New trends in Large-Eddy Simulations of turbulence[END_REF].

WALE Model

To obtain the right scaling for the turbulent viscosity when approaching a solid boundary, the Van Driest damping function is often used [START_REF] Van Driest | On the turbulent flow near a wall[END_REF]. A more elegant way is the WALE (Wall-Adapting Local Eddy-viscosity) proposed by Nicoud & Ducros [START_REF] Ducros | Wall-adapting local eddy-viscosity models for simulations in complex geometries[END_REF][START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF]. It is an eddy viscosity model based on the square of the velocity gradient tensor and accounts 4.5. Boundary conditions for the effects of both the strain and the rotation rate to obtain the local eddy-viscosity.

It recovers the proper y 3 near-wall scaling for the eddy-viscosity without requiring a dynamic procedure. They replace the characteristic filtered rate of strain by a term that detects strong rates of deformation and/or rotation and not shear as in Smagorinsky model.

s d ij = 1 2 g2 ij + g2 ji + 1 3 g2 kk δ ij (4.217) ν t = (C w △) 2 s d ij s d ij 3 2 S ij S ij 5 2 + s d ij s d ij 5 4 (4.218) with g2 ij = ∂u i ∂x k ∂u k ∂x j
where △ denotes the filter characteristic length (cubic-root of the cell volume), Cw = 0.4929 is the model constant was calibrated numerically on isotropic decaying turbulence and gij denotes the resolved velocity gradient. This expression for ν t allows for the right scaling of turbulent viscosity when approaching walls and also for the prediction of transition.

Boundary conditions

Boundary conditions plays an important role in any numerical tool, and especially here in AVBP because of acoustics present in the governing equations (see Schöfeld & Rudgyard [START_REF] Schöfeld | Steady and unsteady flows simulations using the hybrid flow solver avbp[END_REF] and Poinsot & Veynante [START_REF] Poinsot | Combustion numérique-cours de DEA[END_REF]). Large-eddy simulation requires the setting of boundary conditions to fully determine the system and obtain a mathematically wellposed problem. LES contains a large number of degrees of freedom. So it needs a precise space-time deterministic representation of the solution at the computational domain boundaries. Applying boundary conditions in AVBP is equivalent to finding the residual R on the boundaries. The residuals in AVBP are obtained using a Runge-Kutta multi-step time integration. The derivation for a single-step Runge-Kutta scheme is given below. In such a scheme the solution at time t + ∆t(U n+1 ) can be derived from the solution at time t using:

U n+1 = U n -R∆t (4.219)
Note that the code is explicit therefore only the solution indexed n is used. For each node in the boundaries, a residual R computed by the scheme is corrected using the target values from the boundary conditions. The residuals can be corrected using two methods:

Non characteristic method : This method imposes directly the target conservative variables using the residual. In most cases, this means simply replacing the boundary value predicted by the scheme by the target value. Characteristic method : The correction is applied through the use of a wave decomposition. This method is called Navier-Stokes characteristic boundary condition(NSCBC)(see Moureau et al [START_REF] Moureau | Numerical methods for unsteady compressible multi-component reacting flows on fixed and moving grids[END_REF] and Poinsot & Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]). Here, no value is imposed directly to variables such as velocities or densities. Only waves amplitudes are specified.

Characteristic conditions for Euler equations were first derived by Thompson [START_REF] Thompson | Time dependent boundary conditions for hyperbolic systems[END_REF][START_REF] Thompson | Time dependent boundary conditions for hyperbolic systems, II[END_REF] the extension to Navier-Stokes equations is due to Poinsot & Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF]. In other words, compressible flows are characterised by waves whose physics is to be respected in numerical simulations. Characteristic boundary conditions allows for the correct treatment of waves impacting a boundary of the computational domain. At run time, either the widely used full residual or the normal approach is selected. In the normal approach, the characteristic boundary conditions are applied to the normal flux derivatives of the residual, which is in general more accurate and implies a decomposition of variations in the conservative variables into a set of ingoing and outgoing waves. The correct formulation of the boundary conditions has been the subject of an extensive investigation(refer to Nicoud [START_REF] Nicoud | Defining wave amplitude in characteristic boundary conditions[END_REF]).

Building the characteristic boundary condition

The Navier-Stokes equations which is written for a multigas flow have been discussed

for DNS in section 2.4and LES in section 4.4. To be able to apply boundary conditions, following mathematical operations becomes essential:

Boundary conditions

• all equations must be written in a reference frame linked to shape of the boundary section. This section has a normal vector n and two tangential vectors t 1 and t 2 .

• all equations must be written in characteristic variables, i.e. in variables such as

∂W 1 = ∂u n + ∂P ρa ∞
, ∂W 2 = ∂u n ..... etc., because these are the only variables which are significant in terms of waves. There are many ways to perform this and appendix of AVBP [10] describes the routine followed in AVBP. Basically, transformations start from conservative variables ∂U to primitive variables ∂V then to primitive variables in n, t 1 , t 2 ∂V n and finally to characteristic variables ∂W (see figure 4.3) that can be seen as waves carrying information normally to the boundary.

The derivation starts directly from the equation for V n .

∂V n ∂t + N ∂V n ∂ n + T 1 ∂V n ∂ t 1 + T 2 ∂V n ∂ t 2 + S = 0 (4.220)
where V n = (u n , u t 1 , u t 1 , P, ρ 1 , ....ρ N ) T is the primitive variable vector (N = 1 for a single species gas), assuming that dX n + dY t 1 + dZ t 2 = d n + d t 1 + d t 2 to simplify the notation. N represents the normal Jacobian, T 1 and T 2 the two tangential Jacobians along t 1 and t 2 . The S term sums all the contributions related to diffusion terms and chemical reactions. The principle of characteristic boundary conditions is to diagonalise the normal Jacobian N to write convection equations for characteristic variables W :

∂W ∂t + D ∂W ∂ n = S W -T W (4.221)
where D is the diagonal matrix containing the propagation velocity (eigen values of N ) of the waves and S W -T W are the sum of all non-hyperbolic terms associated to the wave: reaction, diffusion and tangential terms.

D =             u. n + a ∞ 0 . . . 0 . . . 0 0 u. n -a ∞ 0 . . . u. n . . . 0 u. n 0 . . . . . . . . . 0 0 . . . 0 . . . u. n             (4.222)
Variations of characteristic variables ∂W can be obtained from variations of primitive variables and vice versa, using L and R matrices (see figure 4.3):

∂W = L∂V n , ∂V n = R∂W (4.223)
4. Numerical simulation and LES models Similar relations hold for the passage from variations of conservative variables ∂U to variations of characteristic variables ∂W using L U and R U matrices:

∂W = L U ∂U, ∂U = R U ∂W (4.224)
The expressions of the variations of characteristic variables in terms of primitive variables are recalled :

         ∂W 1 ∂W 2 ∂W 3 ∂W 4 ∂W 4+k          =             ∂u n + 1 ρa ∞ ∂P -∂u n + 1 ρa ∞ ∂P ∂u t 1 ∂u t 2 - Y k a 2 ∞ ∂P + ∂ρ k             (4.225)
The first two characteristic variations represent acoustic disturbances, the third and the fourth are related to variations in shear velocity and the ∂W 4+k correspond to chemical species k. The associated propagation velocities are:

         λ 1 λ 2 λ 3 λ 4 λ 4+k          =          u n + a ∞ u n -a ∞ u n u n u n          (4.226)
The so called entropy wave can be recast by summing all the species waves as:

∂W S = N k=1 ∂W 4+k = - 1 a 2 ∞ ∂P + ∂ρ (4.227)
Some important inverse relations useful for the imposition of the incoming waves from the outgoing ones are recalled:

∂u n = 1 2 ∂W 1 -∂W 2 (4.228) ∂u t 1 = ∂W 3 (4.229) ∂u t 2 = ∂W 4 (4.230) ∂u = 1 2 n x ∂W 1 -∂W 2 + t 1x ∂W 3 + t 2x ∂W 4 (4.231) ∂v = 1 2 n y ∂W 1 -∂W 2 + t 1y ∂W 3 + t 2y ∂W 4 (4.232) ∂w = 1 2 n x ∂W 1 -∂W 2 + t 1z ∂W 3 + t 2z ∂W 4 (4.233) 4.5. Boundary conditions ∂P = 1 2 ρa ∞ ∂W 1 + ∂W 2 (4.234) ∂ρ k = ρ k 2a ∞ ∂W 1 + ∂W 2 + ∂W 4+k (4.235) ∂ρ = ρ 2a ∞ ∂W 1 + ∂W 2 + ∂W S (4.236) ∂Y k = 1 ρ ∂W 4+k -Y k ∂W S (4.237) ∂r = 1 ρ k r k ∂W 4+k -r∂W S (4.238) ∂T = βT 2a ∞ ∂W 1 + ∂W 2 - j r j T ρr ∂W 4+j (4.239) ∂ρu = ρ (u + a ∞ n x ) 2a ∞ ∂W 1 + ρ (u -a ∞ n x ) 2a ∞ ∂W 2 + ρt 1x ∂W 3 + ρt 2x ∂W 4 + u∂W S (4.240) ∂ρv = ρ (u + a ∞ n y ) 2a ∞ ∂W 1 + ρ (u -a ∞ n y ) 2a ∞ ∂W 2 + ρt 1y ∂W 3 + ρt 2y ∂W 4 + v∂W S (4.241) ∂ρw = ρ (w + a ∞ n z ) 2a ∞ ∂W 1 + ρ (u -a ∞ n z ) 2a ∞ ∂W 2 + ρt 1z ∂W 3 + ρt 2z ∂W 4 + w∂W S (4.242)
It is important to understand :

• the dealing with variations of characteristic variables (obtained by variations of conserved or primitive variables). The choice of calculating these variatiosn has to be done on the ∂W variables. In the following sections the two approaches: the spatial form and the temporal form are discussed.

• ∂W described here correspond exactly to strength variables in the AVBP coding.

• relations 4.229 to 4.242 are used for all boundary condition formulations in AVBP to prescribe incoming waves. As explained in the following sections, only the the evaluation of the wave amplitudes from the predicted values of U depends on the formulation chosen.

It should be noted that strength(4) is always the acoustic incoming wave (i.e. ∂W 1 for inflows and ∂W 2 for outflows) since it follows the sign of the vector normal to the boundary which, in AVBP, points always inward.

The procedure of Boundary Condition treatment in AVBP is summarised in figure 4.4.

Numerical simulation and LES models

Predicted values

strength p = strength p = strength p = LU λ ∂U ∂n δt LU (R p -R p t ) LU R p R in,P BC = RU • strength(in) P Boundary Conditions INLET WAVE UVW T Y INLET RELAX UVW T Y OUTLET RELAX P R in,C BC = RU • strength(in) C U n+1 = U n -∆t(R P BC -R in,P BC + R in,C BC + R P U )
Prediction of the numerical scheme compute waves amplitude contribution of "wrong" incoming waves Specify incoming vs outgoing waves Same approach for all formulations

Redistribute corrected incoming waves to residuals

Apply corrected residuals to advance the solution in time The explicit time advancement scheme of AVBP leads to the predicted value U n+1 pred :

∂U = U n+1 pred -U n = -R P ∆t (4.243)
The total residual R P can be split into two parts :

∂U = -∆t R P BC + R P U (4.244)
R P BC is the residual part which will be modified by the BC treatment and R P U the part which will be left unchanged. The objective of the BC treatment is to construct the final value of U at time n + 1 :

U n+1 U n+1 = U n -∆t R P BC + R P U (4.245)
where R C BC is the part of the residual which has been corrected using R P U , U n , the type of BC and the target values. The correction is made in the following way:

R C BC = R P BC -R in,P BC + R in,C BC (4.246) 4.5. Boundary conditions
i.e. by substituting the contribution on the residuals of the predicted "wrong"incoming waves R in,P BC by their correct values given by the boundary condition R in,C BC . A fundamental issue is to choose the residual part to update R P BC . There are two main methods in AVBP to update R P BC linked to the spatial and temporal formulation described in the next sections. Other ways to choose the part of update do exist, using:

• the advection terms of the bi-characteristic equations (from Hirsch [START_REF] Hirsch | Finite volume method and conservative discretization with an introduction to finite element method[END_REF]).

• a Fourier decomposition of the solution at the boundary (from Giles [START_REF] Giles | Non-reflecting boundary conditions for euler equation calculations[END_REF]).

• viscous and reacting terms (see Sutherland [START_REF] Sutherland | Improved boundary conditions for viscous, reacting, compressible flows[END_REF]).

• a decomposition between the convective and the acoustic part to build the waves(from Prosser [START_REF] Prosser | Improved boundary conditions for the direct numerical simulation of turbulent subsonic flows. I. inviscid flows[END_REF][START_REF] Prosser | Toward improved boundary conditions for the DNS and LES of turbulent subsonic flows[END_REF] These topics are not be presented here. More details and comparisons between all these methods can be found in Nicoud & Poinsot [START_REF] Nicoud | Boundary conditions for compressible unsteady flows[END_REF] and in the AVBP handbook [10].

Spatial formulation

In the spatial formulation, which is the initial form of the Navier-Stokes characteristic boundary condition method from Poinsot & Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF], the ∂W are defined from spatial gradients:

∂W = strength = -λ ∂W ∂n ∆t (4.247)
where λ is a vector containing the eigenvalues of the normal Jacobian, i.e. the propagation speed of the waves. This means that the variations of characteristic variables in the spatial formulation are proportional to normal spatial gradients of variables. Following the development by Poinsot & Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] we can introduce the L notation:

L = λ ∂W ∂n (4.248)
More informations on Navier-Stokes characteristic boundary condition and on the equivalence with the ∂W notation can be found in the AVBP handbook [10]. To build the boundary condition, variations of characteristic variables ∂W must be obtained from residuals. The computations of the strength from the residuals R P is then performed using the normal residual approach. This corresponds to the Navier-Stokes characteristic boundary condition formulation from Poinsot & Lele [START_REF] Poinsot | Boundary conditions for direct simulations of compressible viscous flows[END_REF] in which spatial derivatives normal to the boundary are used to update R P BC . To do this, the residual R P must be split in two parts :

R P = R P n normal part + R P t + R P Dif f usion + R P Chemistry non normal part (4.249)

Numerical simulation and LES models

The Navier-Stokes characteristic boundary condition method assumes that only the normal part must be updated :

R P BC = R P n (4.250)
while the non normal part is unchanged:

R P U = R P t + R P Dif f usion + R P Chemistry (4.251)
Therefore, variations of conservative variables linked to the normal residual can be written as

∂U = -R P n ∆t (4.252)
The normal part of the residuals can be defined in the following way:

R P n = N U ∂U ∂n (4.253)
where

N U = A U n x + B U n y + C U n z z is the normal Jacobian in conservative variables.
With wave decomposition from Handbook of AVBP [10], N U is:

N U = R U DL U (4.254)
where, as usual, D is the eigenvalues diagonal matrix. The values of predicted strength are obtained by :

strength P = ∂W = -L U ∂U = -L U R U DL U ∂U ∂n ∆t = -L U λ i ∂U ∂n ∆t (4.255)
Characteristic variables variations are therefore calculated using spatial normal derivates of conserved variables. The boundary condition are applied to impose the ingoing waves strength(in) C and the solution is projected back to the residuals according to equation 4.246:

R C BC = R P BC -R in,P BC + R in,C BC (4.256)
where:

∆tR P BC = R U strength P ∆tR in,P BC = R U strength(in) P ∆tR in,C BC = R U strength(in) C
The final value for U n+1 is then :

U n+1 = U n -∆tR C BC -∆t R P t + R P Dif f usion + R P Chemistry (4.257)

Boundary conditions

Note that this method does not enforce strictly the value of U n on the boundary since the tangential, viscous and chemical terms are not accounted for when assessing the corrected value of the incoming waves.

Colin [START_REF] Colin | Simulations aux grandes échelles de la combustion turbulente prémélangé dans les statoréacteurs[END_REF] developed an alternative method for calculating the normal part of the residuals (iwave = 3). The idea is to subtract the transverse part of the residual from the total residual.

R P BC -R P n = R P -R P t (4.258)
This transverse residual R P t can be calculated on the boundary using the "complete"centered numerical scheme. On the contrary, gradients normal to the wall used for iwave = 1, are calculated with a "truncated"and less precise scheme, since, on the boundary, we have access only to cells inside the domain.

Temporal formulation

Computing spatial derivatives as in the spatial form can be difficult. An alternative solution is to use time variations to evaluate R P BC : in the temporal formulation originally introduced by Thompson [START_REF] Thompson | Time dependent boundary conditions for hyperbolic systems[END_REF], the ∂W are defined as

∂W = ∂W ∂t ∆t = strength (4.259)
Characteristic variables variations are then calculated as a temporal variation (not a temporal derivative) of primitive (or conserved) variables. The computation of the variations of characteristic variables strength from the residuals R P is then performed using the full residual approach. In this case the total residual R P is used for R P BC so that R P U = 0 in equation 4.244 The predicted variations in conservative variables are

∂U = -R P ∆t (4.260)
where R P is the actual residual calculated by AVBP before the application of boundary conditions. This means that only time changes are used to compute waves and there is no need for normal spatial gradients. Now predicted variations of characteristic variables can be computed from the variations of conservative variables using the left passage matrix L U .

strength P = ∂W = L U ∂U = -L U R P ∆t (4.261)
All waves going out of the domain are left unchanged in strength P while corrected incoming waves strength(in) C are computed using the relations not detailed here but found in the AVBP handbook. Having modified strength(in) P , the corrected R C BC is obtained, as for the spatial formulation, by : where:

R C BC = R P BC -R
∆tR P BC = R U strength P ∆tR in,P BC = R U strength(in) P ∆tR in,C BC = R U strength(in) C
and U n+1 can finally be obtained by equation 4.245:

U n+1 = U n -R C BC ∆t (4.263)

Boundary conditions

Because of the wide range of applications of AVBP, variety of boundary conditions are available in AVBP. Boundary conditions which are used this work are given here in the table 4.2. Derivation, implementation are discussed in this section.

No-Slip Conditions

The presence of solid wall inhibits growth of small scales and modifies the turbulence dynamics in several ways. When a fluid flow over a solid surface, the layer next to the surface may come attached to it. This is called the 'no-slip condition'. In most references concerned with fluid mechanics, the only boundary condition discussed is the no-slip condition. This condition is the analog of the constitutive relations and therefore only holds when at least one material is a Navier-Stokes fluid. The no-slip boundary condition demands that the velocity component tangential to the wall be the same as the tangential velocity of the wall. If the wall is at rest relative, then the no-slip condition demands the tangential flow velocity be identically zero at the surface. The no-slip conditions are normally ignored when the inviscid approximation is made.

WALL WAVE NOSLIP ADIAB

This no-slip adiabatic wall boundary condition imposes zero velocity on the wall through a characteristic treatment of the acoustic and the shear waves (strength(4), strength(

and strength(3)) such as the one used with INLET WAVE UVW T Y. other waves are left unchanged. It also imposes a zero heat flux and species mass flux through the wall. The Von Neumann conditions are applied in a weak way. The "weak"part of the boundary condition expresses that the temperature, pressure and the species gradients normal to the wall are zero. In the "characteristic"part we let the scheme predict every variable, except for the velocity which must be corrected. Hence in terms of wave:

strength(4) = strength(5) + 2(U t n -U n ) strength(2) = U t t1 -U t1 strength(3) = U t t2 -U t2 (4.264) (4.265)

WALL WAVE NOSLIP ISOT

This no-slip isothermal wall boundary condition imposes both zero velocity and prescribed temperature on the wall through a characteristic treatment of the acoustic, shear and species waves ( strength(4), strength(2), strength(3) and strength(5 + k)) such as the one used with INLET WAVE UVW T Y. Hence:

strength(4) = strength(5) + 2(U t n -U n ) strength(2) = U t t1 -U t1 strength(3) = U t t2 -U t2 (4.266) strength(5 + k) = Y k strength(1) strength(1) = ρβ a ∞ strength(5) + (U t n -U n ) - ρ(T t -T ) T t
Since species are not imposed on walls (target values = predicted values) species contributions are neglected.

Inlet

In AVBP, a variety of inlet boundary conditions (introducing an acoustic excitation, turbulent perturbation, supersonic) are available. But inlet conditions relevant to the simulation are discussed in detail here. The table 4.3 gives the correspondance between the notation used in derviation and in AVBP.

INLET WAVE UVW T Y

This inlet characteristic boundary condition allows to impose the velocity components, the static temperature and the mass fraction at an inlet in a strong way. The ingoing waves are computed from the knowledge of the outgoing waves and of the current state in such a way that velocity and temperature are properly imposed.

∂u = 0, ∂v = 0, ∂w = 0, ∂T = 0, and
∂Y k = 0 (4.267)
For an inlet, it is necessary to impose values for the incoming acoustic wave, the two shear waves (in 3D) and all the species waves. Since, as explained before, the ∂W notation is employed for both spatial and temporal formulations to prescribe incoming waves, the derivation of all boundary conditions will be made with this notation. Moreover, some hints on the actual coding are added. The incoming acoustic wave can be derived

∂W 1 = ∂W 2 + 2∂u (4.268)
Using equations 4.268, 4.238 and 4.239 can be rewritten to give the entropy wave

∂W S = - ρ∂T T + ρβ a ∞ ∂W 2 + ∂u - ρ r ∂r (4.269)
In AVBP, the variations of the characteristic variables are called strength. The values of these variables are the same as ∂W but the numbering is different. Note that the incoming acoustic wave is always strength(4) since it is related to the boundary normal which, in AVBP, is always internal. According to this change of variables and to relations 4.267, waves should be written as:

strength(4) = strength(5) strength(2) = 0 strength(3) = 0 (4.270) strength(5 + k) = 0 strength(1) = ρβ a ∞ strength (5) 4.5. Boundary conditions
Velocity, temperature and species mass fractions should not change. To avoid drifts of imposed values, the variations of these variables are added. For example ∂u is approximated by the difference between the target value of the velocity and its actual value.

The boundary condition is finally written:

strength(4) = strength(5) + 2 U t n -U n strength(2) = U t t1 -U t1 strength(3) = U t t2 -U t2 (4.271) strength(5 + k) = ρ Y t k -Y k strength(1) = ρβ a ∞ strength(5) + U t n -U n -ρ T t -T T t - k ρW Y t k -Y k W k
To recast the last term of equation 4.269 as done in the actual coding the following expression is used:

k r k ∂Y k = k rW W k ∂Y k = ∂r (4.272)

INLET RELAX UVW T Y

This characteristic boundary condition INLET RELAX UVW T Y imposed at the inlet of any domain with relaxing co-efficients for velocity components, temperature, and species.

The purely non reflecting condition should impose no incoming wave at all, which means :

strength(2) = 0 strength(3) = 0 strength(4) = 0 and strength(5 + k) = 0 (4.273)
However, to keep the mean inlet variables (U n , U t1 , U t2 , T and Y k ) under control (around the target values U t n , U t t1 , U t t2 , T t and Y t k ), the condition is written in a different way: by specifying the ingoing waves to relax every variable toward its target: Each relaxation coefficient, K i which is homogeneous to a frequency, allows the boundary condition to act as a high frequency filter, with a cut frequency of the order of K i .

strength(4) = 2K Un ∆t(U t n -U n ) strength(2) = K Ut ∆t(U t t1 -U t1 ) strength(3) = K Ut ∆t(U t t2 -U t2 ) (4.274) strength(5 + k) = ρK Y ∆t(Y t k -Y k ) strength(1) = -ρK T ∆t (T t -T ) T 4.
Therefore,it is possible to keep the mean inlet variables around their target values and at the same time let the high frequency waves leave the domain. Choosing values for K i requires a priori evaluation which is described by Selle et al. [START_REF] Selle | Actual impedence of non reflecting boundary conditions: Implications for computaion of resonators[END_REF] and few tests. It should be noted when the relax coefficient is increased, INLET RELAX UVW T Y will behave like INLET WAVE UVW T Y. For all relaxed inlet boundary conditions, the entropy wave is calculated retaining only the contribution of temperature. This is an approximation that can be accepted, since a "soft"boundary condition is dealt, allowing the fluctuation of boundary values.

Outlet

Outlet boundary conditions followed in the simulation are discussed in detail here. Similar to the table from section4.5.5, the table 4.4 gives the correspondance between the notation used in derviation and in AVBP.

OUTLET WAVE P

This characteristic outlet boundary condtion OUTLET WAVE P allows to impose the static pressure at an outlet in a strong way. This means that the ingoing wave is computed from the knowledge of the outgoing waves and of the current state in such a way that the pressure is properly imposed. This outlet condition imposes no pressure variation on the boundary condition:

∂P = 0 (4.275)
This purely reflecting condition should be recast into:

strength(4) = strength(5) (4.276)
4.5. Boundary conditions but to avoid pressure drifts, the condition is corrected :

strength(4) = strength(5) + 2 ρa ∞ (P t P ) (4.277)
This ensures that P remains exactly equal to target pressure P t .

OUTLET RELAX P

The static pressure is imposed at an outlet in a soft way in OUTLET RELAX P characteristic outlet boundary condition. The OUTLET WAVE P boundary condition are perfectly reflecting, and thus do not allow the acoustic energy to leave the computational domain, which may lead to an accumulation of energy in the domain and to a non-physical behaviour. OUTLET RELAX P boundary condition impose quantities in a partially nonreflecting way. The amount of reflection is controlled by the "relax"coefficient. A relax equal to zero leads to a perfectly non-reflecting boundary condition while large values lead to nearly reflecting boundary condition. The maximum value allowed for the relax is 1 ∆t . With this maximum value, the RELAX boundary conditon acts as the WAVE boundary condition. Above this value, the boundary condition is unstable (over-relaxation).

The first "relaxed"boundary condition is a subsonic outlet boundary where pressure is imposed but with a non-reflecting condition. This is achieved by imposing a relaxation on the strength(4) wave. A perfectly non-reflecting boundary condition would be strength(4) = 0. However this formulation can be shown to be ill-posed (it leads to pressure drift, as no information is provided from the outside, (more details in Poinsot& Veynante [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]).

strength(4) = 2 K P ∆t (P t -P n ) ρa ∞ (4.278)
The K P parameter is the so-called "Pressure Relax"coefficient. This coefficient has the same unit as a pulsation (s -1 ) . It allows the boundary condition to act as a high frequency filter, with a cut-off frequency of the order of K P . 1 K P is thus a rough estimation of the relaxation time (the time needed to move from P n to P t ). When K P = 0, this formula gives a perfectly non-reflecting boundary condition. The incoming wave is independent of the outgoing wave, which is the definition of a non-reflecting boundary condition. Choosing a value of K P is equivalent to choosing a reflecting coefficient R. Selle et al [START_REF] Selle | Actual impedence of non reflecting boundary conditions: Implications for computaion of resonators[END_REF] have derived the following relation between R and K P :

R = R K -1 1 -i 2ω K P (4.279)
Finally, relation 4.278 is imposed only if the local Mach number is less than unity.

OUTLET R RELAX P

The classical outlet boundary condition OUTLET R RELAX P can be easily extended to a partially reflecting characteristic boundary condition. By writing

strength(4) = 2 K P ∆t (P t -P n ) ρa ∞ -strength(5) • R K (4.280)
a constant reflection coefficient of magnitude R K is imposed in the frequency range where the classical relax outlet is non-reflecting. Then the reflection coefficient is:

R = -R K + R K -1 1 -i 2ω K P (4.281)

Conclusion

The spatial discretisation, numerical schemes Lax-Wendroff and TTGC were discussed along with the artificial viscosity in this chapter. Governing equations, filtering procedure, turbulence models were given in detail for the large eddy simulation. The characteristic boundary conditions handled in AVBP solver were slso explained elaborately. With the solver, few simulations are conducted to simulate the cavity flow and the results are analysed in the next chapter.

Chapter 5

Analysis of the cavity flows

Introduction

In low Mach number flows, the ratio between the acoustic and hydrodynamic length scales is of order 10 -100. It explains why usually the numerical tools or approaches are different to compute hydrodynamic field and aeroacoustic perturbations. In this chapter, the flow field is determined by the large eddy simulation at very low Mach number in the range of incompressible flow. The acoustic pressure field is provided by the Lighthill-Curle's acoustic analogy but without volume integral.

In this chapter are successively presented and analysed the hydrodynamic over a two-dimensional and three-dimensional cavity flow. Geometry, mesh, boundary conditions and mean flow field and turbulent quantities are discussed. It is followed by the aeroacoustic study and a conclusion.

Two-dimensional cavity

Geometry and mesh

All simulations were performed on the cavity of aspect ratio (L/D) 4. Through out the work, the length of the cavity is maintained as 0.04 m and depth of the cavity as 0.01 m.

The coordinates are non-dimensionalised by depth of the cavity D. to the downstream corner of the cavity is higher than in the upstream region of the cavity.

In the literature, it is predicted that the mesh at corner and the vertical walls of the downstream edge of the cavity has to handle high gradients during the impingement of shear layer or the energetic eddies. Lean mesh density is at the top of the computational domain. The lean mesh density reflects the absence of eddies and gradients in top extreme region of the domain. No buffer region have been added. to solve accurately high gradients, bad outflow boundary conditions, waves reflection, etc ...). They will not be discussed here though a lot of time has been spend to solve them step by step.

Numerical schemes and LES Model

As mentioned in the chapter 4 that the solver AVBP handles many species (N

Inlet condition

The details of the two-dimensional test cases are given the table 5.4. The dimensional variables which characterise the cavity flow according to Colonius [START_REF] Colonius | An overview of simulation, modeling and active control of flow/acoustic resonance in open cavities[END_REF] are the length of the cavity L, the depth of the cavity D, free stream velocity u ∞ , momentum thickness θ, velocity of sound in the medium a ∞ , and kinematic viscosity ν ∞ . The role played by the momentum thickness θ at the leading edge of the cavity in the selection of modes 5.6) signifies that the ingoing waves which are computed pointwise so that the relaxation tends to drive the velocity components and temperature towards the exact profiles given by the reference state. At the inlet, the ingoing waves are taken proportional to the difference between the actual state at the boundary nodes and the reference velocity and temperature. Three integer parameters allow to fix the way the waves are calculated in the solver AVBP, the way the reference state is defined and the type of relaxation that is performed. Four real parameters are used to prescribe the values of the relaxation coefficients: relax on Un, relax on Ut, relax on T, relax on Y.

Temporal approximation is followed. The variations of the conservative variables are assessed from their time derivative and the strength of the wave. More details are found in the section 4.5.5.

With the details from the section 4.5.6, the top and right hand side edges of the computational domain are treated as outlet and characteristic boundary condition OUTLET RELAX P is applied with relax on P, the relaxation parameter on the pressure (see equation 4.280).

An improved normal approximation as in [START_REF] Colin | Simulations aux grandes échelles de la combustion turbulente prémélangé dans les statoréacteurs[END_REF] is followed and suits well for the outlet boundaries. Finally the bottom edges of the domain are treated as the solid walls with no slip and adiabatic boundary condition WALL NOSLIP ADIAB (see the subsection 4.5.4).

Boundary layer flow part

Velocity profiles

The incoming boundary layer of the test cases behaves, in the region

Cavity results

Time traces

The three cases U40, U20 and U5.8 are imposed with thick boundary layer with low

Mach number (see table 5.4).

Wake mode

For test cases U20 and U40, wake mode is observed. Figure 5.12 shows the instantaneous vorticity fields ω z over a period for the test U40. A vortex is formed from the trailing edge and fills the cavity region is shown in figure 5.12(a). A low pressure zone is created at the downstream wall of the cavity. In figure 5.12(b), the vortex detaches and impinges on the downstream corner of the cavity. Due to the impingement, the is ruptured and moves out of the cavity, while another eddy enter the cavity from the leading edge of the cavity (see fig. 5.12(c)). The eddy which is broken at this point of time moves downstream of the cavity, while another new eddy grows to fill the cavity is shown in the figure 5.12(d). The flow above the cavity region is affected by the flow from the cavity. The free stream flow is periodically directed into the cavity.

Near the upstream vertical wall of the cavity in the figure 5.12(b) and in the middle of the cavity of the figure 5.12(c), interaction of two counter rotating eddies produce high gradients, and invoke numerical errors resulting in blowing up of the solution. The inclusion of artificial viscosity (see subsection 4.2.7) tends to smooth these gradients and introduces artificial dissipation. Different values of artificial viscosity values have been tried for this configuration and the solution is converged with the following artificial viscosity values: smu4= 0.05 for 4 th order operator and smu2= 0.2 for 2 nd order operator.

The flow was found to be highly unsteady and strongly influenced by the behaviour of the shear layer. Larsson et al [START_REF] Larsson | Aeroacoustic investigation of an open cavity at low Mach number[END_REF] observed wake mode at M = 0.15 in his twodimensional direct numerical simulations. Only a small portion of the computational domain near the cavity is shown.

Shear mode

The test case U5.8 (with equilibrium turbulent boundary layer) oscillates in shear mode.

The turbulent boundary layer which separates from the leading edge of cavity forms an oscillating shear layer. Figure 5.13 shows the instantaneous vorticity fields ω z over a period T . Figure 5.13(a) shows the shear layer stretching from the upstream of the cavity and is parallel to the bottom of the cavity. Over the right upper corner of the cavity, shear layer with a tongue like structure extends to the downstream of the cavity from the vortex near to the vertical wall of the cavity. Figure 5.13(b) describes the complex interaction between the shear layer and the vortex at the downstream wall of the cavity. The incoming shear layer extends until the middle of the cavity region and the lip of the shear layer swipes on the vertical wall at the trailing edge of the cavity.

The swiping action cuts the tongue like shear layer to travel downstream of the cavity.

C f = τ w 0.5ρu 2 ∞
The shear stress on the bottom wall depends on the interaction of vortices and wall surface. The shear stress values found on upstream side of the cavity bottom (between x D = 5 and 5.5) for U5.8 are lesser than the values from the cavity U40. In this region, the cavity U5.8 which is operating in the shear mode, the interaction between the vortices and the bottom wall is less. A peak is observed near the downstream region of the cavity (at x D = 8.5) indicating the interaction of primary vortex and the bottom wall.

In case of cavity U40 which falls in wake regime, the peaks are observed in the upstream region and the downstream region of the cavity. It indicates the level of interaction between the vortices and the bottom wall. At x D = 8.5, the skin friction value of U5.8 is higher than the U40. This interaction reflects on the shear stress values and therefore should increase the overall value of cavity drag. From the figure 5.16, the energy cascade for the cavity U20 can be observed. The slope of the cascade fits well with the theoretical prediction of -

Two-dimensional cavity

Also shown in the figure 5.17 of test case U40 is a -5 3 Kolmogorov slope plotted by the solid line. Two apparent peaks corresponding to the dominant oscillation frequency and its first harmonic can be observed. A less prominent peak at the frequency of the second harmonic is also indicated. mode of the shear layer mode, and it is a higher frequency than the one found with the other two testcases U20 and U40 which are operating in wake mode. A maximum value v ′ 2 ≈ 0.2 is found for the mean fluctuating normal velocity (see figure 5.22). The role played by the u ′ 2 quantities is slightly less than the v ′ 2 quantities in the shear mode. In the wake mode (U40), the contour values of u ′ 2 are in the same order of the contours of v ′ 2 where as in the shear mode (U5.8), the order of u ′ 2 is less than the v ′ 2 . Increasing velocity from 5.8 m/s to 40 m/s increases the order of magnitude of the turbulent intensity by two. This clearly states the role and the amount of velocity fluctuations due to vortices in the cavities operating in wake mode. In the shear layer mode, the stress profiles (see figure 5.29) of the test case U5.8 are qualitatively similar to profiles from Bertier et al [START_REF] Bertier | Large eddy simulation of a subsonic flow over a cavity on general unstructured grids[END_REF] and are typically of separated flows.

Turbulent fluctuations

Mean fluctuating velocity field

q v ′2 q v ′2 0
u ′ u ′ u ′ v ′ v ′ v ′ ×10 -3 ×10 -3
u ′ u ′ u ′ v ′ v ′ v ′
u ′ u ′ u ′ v ′ v ′ v ′ ×10 -4 ×10 -5 ×10 -6
u ′ u ′ u ′ v ′ v ′ v ′ ×10 -3 ×10 -3

Aeroacoustics

The interaction of the vortex with the trailing edge of the cavity generates pressure waves which are radiated into the far field. These pressure waves are identified as aerodynamic noise. To determine the sound pressure level using the acoustic analogy, an acoustic

Geometry and mesh

The aim of this section is to determine whether a wake mode is observed in threedimensional simulation with a power law profile imposed at the inlet of the domain. For the two-dimensional test cases U20, U40 presented in the earlier section, wake 5. Analysis of the cavity flows mode was observed. The presence of wake mode is due to the absence of 3D effects, which is considered as an artifact of two-dimensional cavity flow. The both modes were analysed from the means and instantaneous turbulent velocity and vorticity field. The region of high turbulence activity have been observed through the iso-contours of main turbulent turbulent fluctuations.

Conclusion

Comparison with turbulent intensity field from Haigermoser [START_REF] Haigermoser | A timeresolved particle image velocimetry investigation of a cavity flow with a thick incoming turbulent boundary layer[END_REF] has shown that region of high intensity were given different by the simulation, but levels have been predicted well.

A three-dimensional simulation was performed to verify the presence of wake mode as found in the two dimensional test case. In the three-dimensional cavity, wake mode was not observed for the test case with u ∞ = 20 m/s. Due to the absence of three-dimensional effects in the two-dimensional flow, presence the wake mode (in twodimensional cases) can be considered as an artifact.

Sound pressure levels for the test case u ∞ = 5.8m/s and u ∞ = 40m/s are calculated.

The maximum SPL values (which are observed at the downstream vertical wall) are validated. The mode of propagation and the directivity of the propagation for these low Mach numbers are not prominent. But this study again proves the influence of the boundary layer thickness and velocity at the inlet in the generation of sound and its intensity.

Conclusions

Résumé étendu en français

Conclusions

This PhD is concerned with a numerical study of the two-dimensional turbulent cavity flows regarding the fluid dynamic and the acoustic aspects. It is integrated into a Marie Curie Early Stage Training Actions called AeroTraNet.

The numerical approach is based on Large Eddy Simulations which are resolved using the AVBP code from CERFACS. The acoustic analysis based on the Lighthill-Curle analogy has been performed.

The flow over a cavity can be divided into three regions : An incoming turbulent boundary layer which develops and grows before reaching the cavity, the turbulent flow inside and above the cavity which can be distinguished into a wake mode and a shear layer mode, and an unsteady boundary layer flow downstream of the cavity which is with or without separation and convects the structures ejected from the cavity.

In our case, an incoming turbulent boundary layer has been initially defined with a power law mean velocity profile.

In the case of an equilibrium boundary layer, extension of the usual asymptotic approach based on the defect law in the outer turbulent region and the usual friction The turbulent intensity is related to the mode, which depends on the ratio of the momentum thickness of the boundary layer to the cavity depth and on the Mach number.

Increasing the value of the incoming velocity from 5.8 m/s to 40 m/s increases the order of magnitude of the turbulent intensity by two.

At 20 m/s, the two-dimensional case has been extended to three-dimensional to demonstrate that the wake mode is a pure artifact of the numerical simulation, since only a shear mode is observed in three-dimensional cases. In the three-dimensional case, initially the flow remains uniform in the span wise direction and oscillates in wake mode. After an initial period, the wake mode ruptures into a shear mode. The incoming turbulent boundary layer leads to a low frequency flapping of the shear layer resulting in an irregular vortex shedding from the leading edge of the cavity. An explanation is due to the turbulent structures of the incoming boundary layer which continue to exist inside the cavity shear layer resulting in instability mechanisms.

The aeroacoustic study has been done to analyse the acoustic sources and propagation. A large set of instantaneous pressure field was determined from large eddy simulation and used as an input to the Curle's analogy, to compute the acoustic pressure of an observer region. The Sound Pressure Level (SPL), is related to the upstream velocity and the oscillating mode of the shear layer. In all the cases a weak directivity in the upstream direction has been observed. Wake modes with a complex vorticity dynamic directly impacts locally SPL and directivity.

Suggestions for future work

With the increased computing resources (for e.g. IDRIS, Paris and CALMIP, Toulouse), Large eddy simulation and direct numerical simulation could be performed to simulate cavity flows with huge domains and well refined grids. Improved sub grid models such as WALE, Filtered Smagorinsky model and Dynamic Smagorinsky model could be followed to allow a better representation of local phenomena typical of complex turbulent flows.

For a wide range of Reynolds number, thick turbulent boundary layer from experiments (with turbulent quantities) can be added to the existing boundary layer so that they could be imposed at the inlet of the two-dimensional and three-dimensional domains.

Studies over cavities with different aspect ratios could be made with the thick turbulent boundary layer of varying Reynolds number to analyse the parameters which influence the instability mechanism and mode switching.

Study on the three-dimensional rectangular cavities could be extended to cylindrical cavities which show more complex flow inside and over the cavity region with asymmetric acoustic near field.

In calculating the sound pressure level for the low Mach number flow, the volume Auteur : Thangasivam GANDHI Calcul et analyse de l'interaction aéroacoustique dans un écoulement turbulent subsonique affleurant une cavité L'objectif de cette thèse est d'étudier numériquement l'aéroacoustique à faibles nombres de Mach (M < 0.3) pour un écoulement de couche limite turbulente épaisse affleurant une cavité, sur la base de simulations numériques à grandes échelles (LES). Un profil de vitesse en loi puissance et pour une couche limite d'équilibre ont servi comme conditions en entrée du domaine de calcul. La couche limite d'équilibre, sans et avec gradient de pression adverse, a été résolue par une approche asymptotique basée sur une formulation déficitaire avec un nouveau modèle de longueur de mélange. Ce dernier a été validé pour améliorer les comparaisons avec les expériences et les simulations numériques directes. Des simulations LES ont permis de regarder l'influence de l'épaisseur de la couche limite turbulente amont sur le mode d'oscillation d'une cavité L/D = 4. Un accord satisfaisant avec les expériences d'Haigermoser et l'émergence du mode de cisaillement a été obtenu pour la vitesse amont de 5.8m/s. Le mode était de type sillage pour les deux autres cas tests (20 et 40 m/s). Finalement, une simulation 3D a montré que le mode de sillage est un artefact du calcul 2D. En utilisant l'analogie de Lighthill-Curle et les champs de pression instationnaire issus de la simulation, nous avons déterminé les niveaux de pression sonore dans le champ proche et lointain. Conformément aux expériences d'Haigermoser, une faible directivité vers l'amont est trouvée. Le mode de sillage influence très fortement les niveaux de pression acoustique. Mots clefs : Couche limite turbulente, longueur de mélange, cavité, LES, aéroacoustique.

Numerical investigation of aeroacoustic interaction in the turbulent subsonic flow past an open cavity

The objective of this thesis is to study numerically the aeroacoustics of low Mach number (M < 0.3) flow with thick turbulent boundary layer past a cavity based on Large Eddy Simulation (LES). Velocity profiles from power law and equilibrium turbulent boundary layer were imposed as inlet conditions on the computational domain. The equilibrium turbulent boundary layer profiles (zero and adverse pressure gradient) have been generated using asymptotic approach with an improved mixing length model. A good agreement is observed between the computed boundary layer profiles and the profiles obtained from experiments and direct numerical simulations. LES results present the influence of the thickness of the incoming turbulent boundary layers on the mode of oscillation in the shallow cavity of L/D = 4. An agreement with the experiments of Haigermoser and the shear mode have been found for the upstream velocity 5.8 m/s. Wake mode was observed for the other two test cases at 20 and 40 m/s. A 3D cavity simulation is performed to show that the wake mode observed in the 2D calculations is an artifact. The hydrodynamic pressure field obtained from the 2D simulation is used as an input to the acoustic analogy (Lighthill-Curle's analogy), to compute the acoustic pressure field at the near and far-field of the cavities. Conforming the experiments of Haigermoser, a weak directivity of sound propagation was observed. Shear mode influences the sound pressure levels strongly. Key words : Turbulent boundary layer, mixing length, cavity flow, LES, aeroacoustics.

Institut de Mécanique des Fluides de Toulouse, Allée du Pr. Camille Soula, 31400, Toulouse.

Figure 2 .

 2 Figure 2.1(a) illustrates the length L, depth D and width W in an experimental setup with respect to the stream wise flow direction and the figure 2.1(b) carries details

= 1 ,

 1 performed Direct Numerical Simulation (DNS) on two-dimensional cavity of L D = 2 with thick laminar upstream boundary layer and three-dimensional Large Eddy Simulation (LES) for higher Reynolds number on cavities of laminar and turbulent regime. On a cavity of aspect ratio L D Ricot et al [122] used Lattice Boltzmann method for aeroacoustic computations of low subsonic M = 0.044 flows. Gloerfelt et al [52] performed two dimensional Direct Numerical Simulation and hybrid methods to evaluate the far-field noise with a relative thick laminar incoming boundary layer on a cavity of aspect ratio L Schematic diagram of cavity flow

Figure 2 . 1 := 1 .

 211 Figure 2.1: Cavity flow.

Figure 2 . 2 :

 22 Figure 2.2: Sketch of open and closed cavity flow at subsonic speed.[START_REF]Aerodynamics and aero-acoustics of rectangular planform cavities. part I: Time-averaged flow. Data Item 02008[END_REF]

Figure 2 .Figure 2 . 3 :

 223 Figure 2.3: Sketch of open and closed cavity flow at super sonic speed. [109]

L θ ≈ 80 .= 4 .

 804 When the ratio of the cavity length to the momentum thickness of the incoming boundary layer ( L θ ) is in the range 80 < L θ < 120, the self-sustained oscillations take place in the shear layer mode. When L θ exceeds 120, the drag abruptly increases due to the onset of the wake mode. Grace et al [56] performed measurements of both laminar and turbulent upstream boundary layers cases with low Mach number. They found no evidence of self-sustained oscillations in streamwise velocity data obtained using a hotwire or in wall pressure fluctuation data obtained using a microphone when an incoming boundary layer is turbulent. They examined mean and turbulent flow fields in a shallow cavity with aspect ratio L D The laminar cases with L θ = 130 and 190 and the turbulent cases with L θ = 78 and 86 were performed with corresponding Re θ were 2892, 3949 for laminar cases and 6318, 12627 for turbulent cases respectively. A cavity with a laminar incoming boundary layer of ratio L D = 4 at very low Mach number was studied by Ozsoy et al [107]. The results brought observation of Reynolds number sensitivity on the mean and turbulent flow velocities and on the vortex characteristics. In spite of the large values of L θ ranging from 114 to 160 no feedback mechanism involving regular flow self-sustained oscillations were observed.

  times have become powerful to perform Direct Numerical Simulation of the Navier-Stokes equations for turbulent flows. They are restricted to low Reynolds number and on simple geometries, The three dimensional unsteady Navier-Stokes equations also apply to turbulent flow when the values of the dependent variables are understood as instantaneous values. Direct numerical simulation (DNS) resolves all flow scales including the small dissipative scales (see figure 2.6 and acoustic propagation. The simulation domain must be sufficiently large to include all the sound sources of interest and at least part of the acoustic near field. However, the important computational cost related to the strong requirements in terms of mesh resolution and temporal discretisation, prevents the DNS approach from being used for industrial applications. The Navier-Stokes equations completely describe turbulent flows. Therefore a DNS of turbulence does not need any modelling of turbulence. Turbulent flows are intrinsically unsteady and involve various length scales. Therefore an accurate simulation must provide sufficient spatial and temporal resolution. An estimate of the necessary spatial resolution is possible when assuming that the total number of necessary grid points N must at least be equal to 2. Cavity flow, turbulence and aeroacoustics the ratio of the integral turbulent length scale to the Kolmogorov length scale: Even for the restricted cases, DNS becomes difficult and extremely expensive computing problem

2 .Figure 2 . 4 :

 224 Figure 2.4: Leonardo da Vinci sketch of turbulent flow.

Figure 2 . 5 :

 25 Figure 2.5: Movements in the study of turbulence, as described by Chapman and Tobak [12], [95].

2 .

 2 Cavity flow, turbulence and aeroacousticsE(κ) ∼ κ -5/3 ln κ ln κ DNS ln κ LES

Figure 2 . 6 :

 26 Figure 2.6: DNS, RANS and LES on Energy spectrum

Sound generated by fluid

  flows is an area of research which has received an increasing amount of attention during the last 15 years. At this point clear definitions of sound, sound wave and pressure fluctuation should be made. Sound is defined to be the pressure fluctuation. Sound wave is defines to be the part propagating as waves at the velocity of sound in a medium where as hydrodynamic pressure fluctuation is defined to be the pressure fluctuations associated with turbulence . The sound sources are generated by motion, either by the free fluid motion, either by a solid body-fluid interaction. It is possible to split the acoustic problem into two parts: fluid flow and acoustic problem.

Ludwig Prandtl -Figure 3 . 1 :

 Prandtl31 Figure 3.1: Boundary layer with details.

Figure 3 . 2 :ρ e u 2 e θ = ∞ 0 (

 320 Figure 3.2: Log law of wall. McDonough [95].

. 16 )

 16 Equation 3.16 is validated over a relatively wide range of momentum thickness based Reynolds number Re θ = u e θ ν in section 3.4.1. To evaluate eq. 3.16, Rona et al[START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF] take κ = 0.41 and B = 5.0, as proposed by Coles[START_REF] Coles | The law of the wake in the turbulent boundary layer[END_REF].

Figure 3 .

 3 Figure 3.3 illustrates the shear stress τ = τ total near the wall. The shear stress is summation of laminar shear stress (τ lam ) and turbulent shear stress (τ turb ). The laminar stress is more dominant in the region very close to the wall (viscous layer). The dominancy and the influence decreases in the region away from the wall. The turbulent shear stress increases with the increase in y and decreases outside the boundary layer. Fundamental equations for incompressible turbulent boundary layer are given here ∂u ∂x + ∂v ∂y = 0

Figure 3 . 3 :

 33 Figure 3.3: Shear stress near the wall.

Figure 3 . 4 :

 34 Figure 3.4: Boundary layer decks.

  .5 correspond to Re θ = 300, 3. Inflow conditions and asymptotic modelling u +

Figure 3 .Figure 3 . 6 :

 336 Figure 3.6 compares velocity profiles obtained using the Successive Complementary Expansion Method of section 3.4.7 with the same reference data of figure 3.5. n = 4 was used for the numerical prediction of the mixing length in eq. 3.21. The symbols used in figure 3.16 are measured values[START_REF] Spalart | Direct simulation of a turbulent boundary layer up to Re θ = 1410[END_REF][START_REF] Erm | Low-Reynolds-number turbulent boundary layers[END_REF][START_REF] De Graaff | Reynolds number scaling of the flatplate turbulent boundary layer[END_REF][START_REF] Österlund | Experimental studies of zero pressure-gradient turbulent boundary layer flow[END_REF] at different Re θ , labelled as in table 3.1. The continuous lines show the normalised numerical velocity profiles. For clarity, the same incremental shift of u + = 2.5 as in figure3.5 is applied to all curves.The origin of the ordinate of figure 3.6 refers to the Re θ = 300 profile. Figure3.6 shows that the Successive Complementary Expansion Method of section 3.4.7 produces a full

Figure 3 . 7 :

 37 Figure 3.7: Outer layer profile determined from asymptotic approach. Re θ = 22845. (+) experiment, (-) SCEM approach.

Figure 3 . 1

 31 Figure 3.8(a) compares the normalised mixing length distribution across a zeropressure gradient boundary layer with ℓ (η) /(δF 1 ) obtained from measurements at Re τ = 1540 by Klebanoff [73], reported in Hinze [65]. The ℓ distribution (Michel's model, eq. 3.20) is shown by the continuous line while the dashed line shows the distribution from equation 3.21 with n = 4. n = 2.7 would provide the same plot as Michel's model case. At this Reynolds number, the new formulation appears to be a good improvement in the predicted the mixing length. No effort has been made to further optimise n ∈ ℜ

Figure 3 . 8 :

 38 Figure 3.8: Turbulence model variables. (•) experiment[START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF] at Re τ = 1540, ( ) experiment[START_REF] Townsend | The structure of the turbulent boundary layer[END_REF] at Re τ = 2775, (--) asymptotic approach at Re τ = 1000 with F 1 = 3.1479 from eq. 3.20 (Michel's model), (-) asymptotic approach at Re τ = 1000 with F 1 = 3.1044 from eq. 3.21. (present model)

Figure 3 . 9 :

 39 Figure 3.9: Velocity profile Re τ = 900

Figure 3 . 10 :Figure 3 . 11 :

 310311 Figure 3.10: Non-dimensional turbulent stress τ t τ w (continous line) and non-dimensional velocity slope du + dy + (dashed line) vs y + at Re τ = 900

Figures 3. 12

 12 Figures 3.12(a) and 3.12(b) show the streamwise non-dimensional velocity profiles

  Figures 3.12(a) and 3.12(b) show the streamwise non-dimensional velocity profiles u + produced by DNS of Skote and asymptotic approach. A very good agreement could be observed for the both Reynolds number Re τ = 220 and Re τ = 270. As shown in tables, the u + e value is a little bit over-estimated in the asymptotic approach. Two plots in the figure 3.13 show the comparison of the non-dimensional turbulent shear stress for the test cases with Re τ = 220 (see figure 3.13(a)) and Re τ = 270 (see figure 3.13). In the two cases, the shear stress curves are normally smooth with the DNS while the curves from asymptotic approach follows the DNS data all except a certain range of y + ∈ 95 to 100, in the log-law intermediate region. The function which could be more synthetic is the non dimensonal turbulent viscosity νt = ν u τ δF 1 . The figure 3.14 gives the turbulent viscosity of the zero pressure gradient

Figure 3 . 12 :

 312 Figure 3.12: Comparison of velocity profiles from Skote's DNS and asymptotic approach for n = 4

Figure 3 . 13 :

 313 Figure 3.13: Comparison of shear stress of test cases Re τ = 220 and Re τ = 270 of Skote's DNS and aysmptotic approach n = 4

3. 6 .Figure 3 . 14 :

 6314 Figure 3.14: Non-dimensional turbulent viscosity, zero pressure gradient, comparison with Skote's DNS.

3. 6

 6 Figure 3.15: Effect of pressure gradient on boundary layer profiles[START_REF] White | Fluid Mechanics[END_REF]

Figure 3 . 16 :

 316 Figure 3.16: Velocity profile: (a) Constant Pressure and (b) Adverse pressure gradient [170].

3. 6 .Figure 3 . 17 :

 6317 Figure 3.17: Comparison of velocity profiles of the test cases (a) APG1 and (b) APG2 with asymptotic approaches with n = 4 and n = 24

Figure 3 .

 3 Figure 3.16 presents comparison of velocity profiles(non-equilibrium) from the computations performed by Wilcox [170] and from the experimental work of So & Mellor [151].

Figure 3 . 18 :Figure 3 . 19 :

 318319 Figure 3.18: Comparison of turbulent shear stresses of test case APG1 from asymptotic approach with n = 4 and n = 24 and from DNS data of Skote

  18 and 3.19 respectively for the values n = 4 and n = 24 and they are compared with the shear stress values obtained from Skote's DNS. In the APG1 case (see figure 3.18), with n = 4, the shear stress curve traces the turbulent shear stress of the DNS throughout the boundary layer except the range y + ∈ 90 -110, where the discontinuity occurs. Even for the case with n = 24, such discontinuity is observed though the magnitude of the τ + is less than for the case with n = 4. Let's remain that the necessity, in the asymptotic approach, to overlap the inner and outer velocity profile of the boundary layer gives rise to this discontinuity which delimitates exactly the point where the both inner and outer regions are joined. The two local peaks in the shear stress curves (for n = 4 or 24) are shown with the present approach. They represent the maximum shear stress values (from equations in the sub section 3.4.6) which correspond to the inner and outer region of the boundary layer. Increase in the value of n shifts the discontinuity to the right hand side and improves the profile of shear stress curve (see figure 3.18), without suppressing the two peaks.

Figure 3 .Figure 3 . 20 :

 3320 Figure 3.20 compares the turbulent viscosity curves ν t u τ δF 1 of APG1 testcase with n = 4 and n = 24, with the Skote's DNS results. Experimental values from Townsend [161]and Klebnoff[START_REF] Klebanoff | Characteristics of turbulence in a boundary layer with zero pressure gradient[END_REF] in zero pressure gradient testcase (and a much larger R τ value) are given as reference. The DNS results exhibit a local maximum for the shear stress at the end of the inner region. Then shear stress should decrease to zero at the edge of the boundary layer (η = 1). It is a similar problem like the ZPG testcases. The asymptotic

3. 6 . 4 Figure 3 . 21 :

 64321 Figure 3.21: Comparison of turbulent viscosity: moderate pressure gradient test case -APG2

Figure 3 . 22 : 24 DNSFigure 3 . 23 :

 32224323 Figure 3.22: DNS analysis :Re τ ν t u τ δ = τ + / du + dy + One point on 3 is designed for the DNS curve

3 .

 3 Inflow conditions and asymptotic modelling the boundary layer edge. The method is shown to predict the outer region of turbulent boundary layers rather well for zero streamwise pressure gradient test cases over the Reynolds number range 300 ≤ Re θ ≤ 31000, with a maximum mean square percentage error of 2.05%. A modification was proposed to the Successive Complementary Expansion Method presented in Cousteix & Mauss [25], with an improved blending function for the mixing length model. Comparison against experimental data shows that this blending function improves the prediction of the mixing length and of the eddy viscosity in the log-law and outer region of a zero pressure gradient boundary layer. The new model is validated against experimental and numerical reference velocity profile data over the Reynolds number range 300 ≤ Re θ ≤ 31000 under zero streamwise pressure gradient and found to achieve engineering accurate predictions. The new blending function introduces an additional adjustable parameter n in the model that can undergo a more extensive calibration over a wider experimental dataset to improve the predictions. Velocity profiles, shear stress profiles and turbulent viscosity profiles for zero pressure gradient and adverse pressure gradient boundary layer have been compared to DNS data from Skote's work. A good agreement is observed in the non-dimensional streamwise velocity and turbulent shear stress in zero pressure gradient (ZPG) case than in adverse pressure (APG) gradient boundary layer. The n parameter of the mixing model is approximately 4 and 24 respectively in ZPG and APG cases. The limit of the asymptotic model has been discussed as well as the sensitivity if the Reynolds number Re τ . The plot of the non-dimensional eddy viscosity determined from DNS results has shed a light on the problem of accuracy in the Direct Numerical Simulation, at least for the adverse pressure gradient case.

4. 2

 2 Numerical method 4.2.1 The cell-vertex discretisation AVBP numerical schemes are based on the cell-vertex method which naturally ensures a high compactness. The flow solver used for the discretisation of the governing equations is based on the finite volume (FV) method (Hirsch [66]). There are three common techniques for implementing FV methods: the so-called cell-centered, vertex-centered and cell-vertex approaches. In the first two ones, the discrete values of the solution are stored at the centre of the control volume (grid cells for the cell-centered formulation and median dual cells for the vertex-centered one, see fig 4.1) and neighbouring values are averaged across the control volume boundaries in order to calculate the fluxes. In the alternative cell-vertex technique, used as underlying numerical discretisation 4. Numerical simulation and LES models Ω j

Figure 4 . 2 :

 42 Figure 4.2: Cell-vertex principle:(a) gather and (b) scatter operation. García [46].

(figure 4

 4 .2(a)). At this stage each cell has locally its information available at the vertices and for example can calculate the cell gradient. The cell quantity is then distributed back to the global nodes through an inverse so-called scatter operation (figure 4.2(b)).

  Navier-Stokes equations, species and model equations include diffusion terms which have the general form: ∂u ∂t = ∇. (ν∇u) (4.144) The diffusion term on the right-hand side of equation 4.144 can be discretised in many ways. The AVBP code uses two different discretisation of the diffusion term: a 4△ operator for the diffusion by laminar diffusivity and a 2△ operator for the diffusion by turbulent diffusivity. As the AVBP code is a finite volume/finite element solver, It was chosen to consider only finite volume (FV) or finite element (FE) discretisations of the diffusion operator. More details related to the operators are found in the handbook of AVBP [10]. The left-hand side (LHS) of 4.144 can be discretised in a FV or FE manner. If a FV convection scheme like the Lax-Wendroff(LW) scheme is used, the LHS is discretised by the mass lumped matrix. The RHS operator is then simply divided by

  I , ḡI , hI T (4.177) Viscous terms : FV = f V , ḡV , hV T (4.178) Turbulent subgridscale terms : Ft = f t , ḡt , ht T (4.179)

4. 4 .

 4 Governing equations for LES and accuracy. One needs a subgrid model to model the turbulent scales which cannot be resolved by the grid and the discretisation. LES models are derived on the theoretical ground that the LES filter is spatially and temporally invariant. Variations in the filter size due to non-uniform meshes or moving meshes are not directly accounted for in the LES models. Change of cell topology is only accounted for through the use of the local cell volume, that is △ = 3 √ V cell . The application of any spatial filtering operation to the Navier Stokes equations here done implicitly through the numerical approximation being tied to the cell size △ leads to the LES equations The filtered compressible Navier-Stokes equations exhibit subgrid scale (SGS) tensors and vectors describing the interaction between the non-resolved and resolved motions. The influence of the SGS on the resolved motion is taken into account in AVBP by a SGS model based on the introduction of a turbulent viscosity ν t . Such an approach assumes the effect of the SGS field on the resolved field to be purely dissipative. The previous hypothesis is essentially valid within the cascade theory of turbulence. The notion of turbulent viscosity can therefore be introduced and yields a general model for the SGS which reads

Figure 4 . 3 :

 43 Figure 4.3: Relation between different set of variables and intermediate matrices involved in the wave decomposition process.

Figure 4 . 4 :

 44 Figure 4.4: Scheme of the global procedure for characteristic boundary conditions from Staffelbach [155].

Figure 5 . 1 4 Table 5 . 1 :Figure 5 . 1 :Figure 5 . 2 :

 514515152 Figure 5.1: Schematic diagram of the computational domain

Figure 5 .

 5 Figure 5.2 shows the density of mesh resolution near the walls and in the cavity region. The boxed region which is highlighted at the upper left corner represents the high mesh density. It is the region where the shear layer and other important mechanisms begin for hydrodynamics and aeroacoustics. The mesh density at the region immediate

  Large eddy simulations were performed on the test cases with larger domain (mentioned in figure5.1) with coarser mesh and refined mesh. The final computational domain on which the simulations are performed has 202 × 121 cells in the cavity region and 824 × 244 cells in the upper part of the domain. It corresponds to 225 500 grid points.Many numerical problems have been encountered (numerical instabilities, difficulty

Figure 5 . 3 :

 53 Figure 5.3: Mean stream wise velocity profile for the turbulent boundary layer in external units, u ∞ = 20m/s

5 Figure 5 . 4 : 5 Figure 5 . 5 :Figure 5 . 6 Figure 5 . 6 :Figure 5 . 7 : 5 . 5 .Figure 5 . 8 :

 5545555656575558 figure 5.4 are similar to those given from the asymptotic approach, but not equal since the equilibrium boundary layer assumes a constant boundary thickness in zero pressure gradient flow. A more accurate analysis show higher non dimensional external velocity u + e in the simulated flow than in the asymptotic solution. As a simple conclusion, the comparison with asymptotic approach have demonstrated that in this section the boundary layer is no longer in equilibrium. The discrepancy increases with the inlet mean velocity.

Figure 5 . 1 ( 1 (Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :

 51159510511 Figure 5.9: Time traces of the velocity components u and v versus non-dimensional time T * for the test case U20.

4 Figure 5 . 12 :

 4512 Figure 5.12: Instantaneous vorticity fields ω z for wake mode (test case U40) at four different times (a-d) corresponding to approximately a quarter of a period of oscillations. Only a small portion of the computational domain near the cavity is shown.

4 Figure 5 . 13 :- 8 Figure 5 . 14 :

 45138514 Figure 5.13: Instantaneous vorticity fields ω z for shear mode (test case U5.8) at four different times (a-d) corresponding to approximately a quarter of a period of oscillations. Only a small portion of the computational domain near the cavity is shown.

3 Figure 5 . 15 :Figure 5 .

 35155 Figure 5.15: Comparison of skin friction on the bottom of the cavities U5.8 and U40

Figure 5 . 16 :

 516 Figure 5.16: Energy spectrum of velocity component in x-direction at x D = 8, y D = 0 for the test case U20.The sampling time of 125 ≤ T * ≤ 600

From figure 5 . 18 ,Figure 5 . 17 :Figure 5 . 18 :

 518517518 Figure 5.17: Energy spectrum of velocity component u at x D = 8, y D = 0 for the test case U40. The sampling time of 100 ≤ T * ≤ 800

Figures 5 .Figure 5 . 19 :Figure 5 . 20 :

 5519520 Figures 5.21 and 5.22 portray the mean fluctuating velocity fields u ′ 2 and u ′ 2respectively which are normalised with respect to the square of the free stream velocity u 2 ∞ for the test case U5.8. The contours which are extended from the trailing edge to the middle of cavity region and parallel to the bottom of the cavity describe the shear layer. The contours over the cavity and the downstream of the cavity represent the fluid

Figure 5 . 21 :Figure 5 . 22 :

 521522 Figure 5.21: Mean fluctuating velocity u ′ 2 field for the test case U5.8.

Figure 5 . 23 :

 523 Figure 5.23: Mean fluctuating velocities from experiments of Haigermoser et al[START_REF] Haigermoser | A timeresolved particle image velocimetry investigation of a cavity flow with a thick incoming turbulent boundary layer[END_REF]

Figures 5 .

 5 Figures 5.23(a) and 5.23(b) from Haigermoser et al [61] represent respectively the variance u ′ 2 and v ′ 2 of the longitudinal and of normal to the wall fluctuation velocities, normalised with respect to the square of the free stream velocity u 2 ∞ . Though the mean fluctuating velocity contours ( u ′ 2 and v ′ 2 ) from Haigermoser et al [61] do not match, the order of fluctuating velocity values fit very well with the values found in this test case. It should be noted that in this test case U5.8, no turbulent quantities were included in the incoming turbulent boundary layer. The mean fluctuating velocity u ′ 2 contours from figure 5.24 for the test case U40 show the maximum value of the fluctuating component near the bottom of the wall.In the wake mode, the vortex which fills the cavity brushes the bottom of the cavity before impinging on the vertical downstream of the wall. This phenomenon explains the presence of higher value of fluctuating component near the bottom of the cavity. The oscillation of the shear layer in the cavity region disturbs the fluid flow above the cavity

Figure 5 . 25 :

 525 Figure 5.25: Mean fluctuating velocity v ′ 2 field for the test case U40.

Figure 5 . 26 : 1 ]

 5261 Figure 5.26: Shear stress profiles for the test case U40 at the stations x D = 2, 4, 5 upstream of the cavity

Figure 5 .

 5 Figure 5.28 shows the shear stress profiles for test case U5.8 at the stations x D = 2, 4, 5 and repeat the profiles from a flat plate boundary layer. The mean Reynolds shear

Figure 5 . 27 :

 527 Figure 5.27: Shear stress profiles for the test case U40 at the stations x D = 6, 7, 8 in the cavity

Figure 5 . 28 :

 528 Figure 5.28: Shear stress profiles for test case U5.8 at the stations x D = 2, 4, 5 upstream of the cavity

Figure 5 . 29 := 8 .

 5298 Figure 5.29: Shear stress profiles for the test case U5.8 at the stations x D = 6, 7, 8 in the cavity

domain of size 0 Figure 5 . 30 :L l i n j 2 Figure 5 .Figure 5 . 31 := 4 .Figure 5 . 32 :Figure 5 . 33 :

 0530255314532533 Figure 5.30: Sound Pressure Level (dB) in the domain for the test case U5.8.

ATable 5 . 9 :Figure 5 . 34 :

 59534 Figure 5.34: Iso contours of vorticity ω z of a three-dimensional cavity

A

  numerical study on two-dimensional and three-dimensional shallow cavities with the aspect ratio L D = 4 was done with incoming turbulent boundary layers. Large eddy simulation with (3 rd order in space and time) has been performed to simulate the cavity flow. Influence of boundary layer on the mode selection was studied by varying thickness of boundary layer and the Mach number of the flow. The low velocity test case corresponds to experiments of Haigermoser. Three test cases with different Mach number (all less than 0.3) were simulated and analysed and compared with the experiments and existing literature. Among the three test cases conducted, one test case (U5.8) was performed with incoming turbulent boundary layer which was generated from the equilibrium turbulent boundary layer approach. The test cases with stream wise velocities u ∞ = 20 m/s and u ∞ = 40 m/s falls in the wake regime where as the other test case u ∞ = 5.8 m/s operates in shear mode. This clearly shows the influence of thickness of the boundary layer which introduces strong shear flow to oscillate at a dominant frequency, and the Mach number selected for the configuration.

  law and log law in the inner region has been proposed.A new blending function of the mixing length model has been proposed to improve the agreement between flat plate zero and adverse pressure gradient turbulent flows. A parameter of the blending function has been determined from experimental data and modified in comparison with the Direct Numerical Simulation of adverse pressure gradient boundary layer flow. There is a good agreement between zero pressure gradient cases and DNS experiments. In the case of an adverse pressure gradient flow, some discrepancies exist, some of which can be explained as problems with the boundary conditions implemented in the simulations (DNS), the other being the limit of the asymptotic approach, especially at low Reynolds number and high pressure gradient flow. Large eddy simulation with Two step Taylor Galerkin Colin scheme (3 rd order in space and time) has been performed to simulate the cavity flow. Classic Smagorinsky and filtered Smagorinsky are the turbulence models used in this work. Characteristic boundary conditions are used at the non-reflecting boundaries and at the walls. Two-dimensional cavity flows were simulated with incoming thick turbulent boundary layer, with an constant aspect ratio of 4. The thickness of the boundary layer introduces a strong shear flow in the cavity which can oscillate at a dominant frequency, with a shear layer mode at low velocity and wake mode at moderate velocities. The low velocity test case corresponds to the experiments, with a shear layer mode providing a validation for the numerical simulation. Both modes were analysed for the mean, instantaneous turbulent velocity and vorticity fields. The iso contours of turbulent fluctuations depict regions of high turbulent activities. The levels of Intensity fields obtained from the experiments agree with the values obtained from the numerical simulation, al though the fields are different.

  DOCTORAT DE L'UNIVERSITE DE TOULOUSEDélivré par : Institut National Polytechnique de Toulouse École doctorale : Mécanique, Energétique, Génie Civil, Procédés (MEGeP) Spécialité : Dynamique des fluides Soutenance : le 10 Novembre 2009

Table 2 .

 2

	Rockwell and Knisely [125] observed three-dimensional pattern in a water channel ex-
	periment for a wide rectangular cavity layer upstream. A hydrogen bubble technique was used to visualise the spanwise wavy L D = 1.08 and W = 3.76 with laminar boundary D
	structure emerging in the shear layer near the cavity trailing edge. Ahuja and Men-

doza

[START_REF] Ahuja | Effects of cavity dimensions, boundary layer and temperature on cavity noise with emphasis on benchmark data to validate computational aeroaoucstic codes[END_REF]

conducted experiments on the effect of cavity dimensions, boundary layer, and temperature on cavity noise for subsonic flows with turbulent boundary layer upstream of the cavity. They determined that the ratio L W the cavity length to width ratio, pro-1: Observations of Hiwada et al

Table 2 .

 2 Boundary layer thickness δ 99 17 mm 2: Details of experiments on cylindrical cavity by Marsden et al

	2.4. Direct Numerical Simulation

and details are given in the table 2.2. The tubulent boundary layer thickness at the incoming boundary condition was chosen smaller than that observed experimentally, and adjusted empir-

  Note that the mass heat capacities at constant pressure c p,k and volume c v,k are supposed constant between T i and T i+1 = T i + 100. They are defined as the slope of the sensible

					51	(2.9)
	enthalpy C p,k =	∂h s,k ∂T	and sensible energy C v,k =	∂e s,k ∂T	. The sensible energy

henceforth varies continuously with the temperature and is obtained by using a linear interpolation:

ref = 10 -12 W m -2 . The reference intensity level I ref is related to the reference pressure p ref by the relationship valid for propagating plane waves:

		2.8. Aeroacoustics
	applications. It causes human discomfort. The most notorious flow noise is that from
	aircraft jet engines, which continues to be an area of intense investigations in response
	to tightening regulation of airport noise. The flow in the cavities is unsteady and, at
	typical landing speeds, may features large-scale instabilities. In a civil aircraft, the high
	lift systems and the landing gear are the most acoustically active airframe components.
	The study of fluid flow over cavities is also relevant for a wide range of applications
	car sunroof, turbomachinery etc. Aeroacoustics, the study of air flow-induced noise
	is concerned with the sound generated by turbulent and/or unsteady vortical flows
	including the effects of any solid boundaries in the flow. The importance of aeroacoustics
	in vehicle and aerospace industry has increased during this decade. In these vehicule
	applications, the Mach number flows are typically small and the flows are often heavily
	separated due to the complex geometries present. The pressure perturbations p ′ ( p ′ =
	p -p 0 ) which propagate as waves and which can be detected by the human ear. For harmonic pressure fluctuations the audio range is:
	20 Hz ≤ f ≤ 20 kHz	(2.41)
	The Sound Pressure Level(SPL) measured in decibel (dB) is defined by:
	SP L = 20 log 10	p ′ rms p ref	(2.42)
	where p ref = 2×10 -5 P a for sound propagating in gases The sound intensity I = I • n is defined as the time averaged energy flux associated to the acoustic wave, propagating
	in direction n. The intensity level(IL) measured in dB is given by:
	IL = 10 log 10	I I ref	(2.43)
	where in air I

  the Lighthill stress tensor and a ∞ is the speed of sound. The equation (2.47) is exact and includes all physics as no assumption is made in deriving it from the governing equations. By assuming ρ ∼ ρ ∞ in term T ij , equation 2.47 becomes explicit. With this assumption, influence of acoustics on the fluid dynamics is not found in the Lighthill's equation. The Lighthill equation 2.47 is the most widely used acoustic

  Ce qui explique la difficulté a validé correctement des calculs de couches limites turbulentes(issus de la théorie, de simulations) et des valeurs expérimentales. Quoi qu'il en soit, l'accord entre le modèle asymptotique et les DNS ou expériences est excellent sur le profil de vitesse. L'accord est moins bon sur le profil de contrainte turbulente, spécialement pour des nombres de Reynolds R τ très faibles. Il est dû au fait que la zone logarith-

	3. Inflow conditions and asymptotic modelling
	donc permettre de proposer des profils de vitesse u + pour une large gamme de nombre
	de Reynolds Re τ ou Re θ basés sur la vitesse de frottement ou bien sur l'épaisseur de
	quantité de mouvement.
	L'approche asymptotique ( section 3.4 ) suit les premiers travaux de Mellor & Gib-
	son [97] et a été reprise et améliorée par de nombreux auteurs, dont Cousteix &
	Mauss [25]. Le calcul de la contrainte turbulente est basée sur une viscosité turbulente
	définie par une longeur de mélange et une fonction d'amortissement au voisinage de la
	3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . .	40
	3.2 Boundary Layer . . . . . . . . . . . . . . . . . . . . . . . . . .	40
	3.3 Analytical method . . . . . . . . . . . . . . . . . . . . . . . . .	44
	3.4 Successive Complementary Expansion Method . . . . . . . .	46
	3.5 Zero pressure gradient boundary layer . . . . . . . . . . . . .	55
	3.6 Adverse pressure gradient boundary layer . . . . . . . . . .	65
	3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .	73
	Résumé étendu en français
	Conditions d'entrée et approche asymptotique
	Ce chapitre a pour objet de définir des profils de couches limites turbulentes pouvant
	servir de conditions d'entrée dans les simulations numériques, en s'appuyant sur une
	approche asymptotique de résolution de la couche limite turbulente déficitaire en con-
	dition d'équilibre (l'épaississement est supposé nul lorsque le gradient longitudinal de
	pression est nul).
	Après quelques brefs rappels sur la couche limite turbulente, est décrite dans la section
	3.3, l'approche analytique avec un profil de vitesse moyenne longitudinale déficitaire
	auto-similaire (eq. 3.11) dans la couche externe contenant une loi de sillage de type
	Coles ( eq. 3.12 ). Cette dernière corrige la loi logarithmique u + = la zone intermédiare dite logarithmique. Une loi de sillage cubique a été implémentée 1 ln y + + B de κ
	par Rona et al [128] et validée pour une couche limite sans gradient de pression par
	rapport à des expériences ou des simulations numériques directes. Cette approche peut

paroi, dans la sous-couche visqueuse. Dans l'approche, on calcule le profil déficitaire auto-similaire en résolvant une équation de similitude non linéaire. Au voisinage de la paroi où la contrainte visqueuse domine, le profil des vitesses u + est intégré simplement et numériquement (eq.3.24) pour converger vers la loi logarithmique, dès lors qu'on dépasse la sous-couche visqueuse. En superposant les profils des vitesses internes et externes dans la zone logarithmique, on calcule le coefficient de frottement pariétal et on en déduit toutes les grandeurs classiques d'une couche limite. Habituellement le modèle de longueur de mélange employé dans cette approche est celui de Michel et al. [99]. Ici est proposé un nouveau modèle (eq. 3.21) fonction d'une paramètre n qui a pour but d'améliorer les comparaisons dans le cas de couche limites avec gradient de pression. L'intégration de l'équation de similitude, non linéaire, du profil de vitesse déficitaire n'est pas évidente car cette é equation dégenère au voisinage de la paroi et de la limite haute de la couche limite. Tous les détails numériques, avec approximation analytique au voisinage des dégénérescences sont donnés dans la section 3.4.7. Sont données aussi les relations analytiques des quantités sans dimensions utiles dans les validations. Dans la section 3.5 est abordée la validation de l'approche asymptotique avec la nouvelle longueur de mélange pour une couche limite sans gradient de pression. Nous comparons nos résultats sur les vitesses u + à des expériences et des simulations numériques directes (DNS de Skote et al) ainsi que la contrainte turbulente sans dimension τ + . L'effet du paramètre n est analysé et n est calibré à n = 4 pour obtenir une meilleure validation sur la viscosité turbulente ou la longueur de mélange normalisée déterminées dans les expériences. n = 2.7 correspond à la valeur donnant des résultats identiques au modèle de Michel. On relève la grande sensibilité des grandeurs de la couche limite aux valeurs des nombres de Reynolds Re τ ou Re θ ou bien à la valeur de u + e = u e /u τ . mique n'existe plus à faible Reynolds et que la méthode pour déterminer le coefficient de frottement turbulent peut provoquer une discontinuité forte sur la dérivée du profil de contrainte turbulente. C'est un artifact qui montre la limite de l'approche.

Table 3 .

 3 1: Experimental velocity profiles. Rona et al[START_REF] Rona | Generation of a turbulent boundary layer inflow for RANS simulations[END_REF].

	3.5. Zero pressure gradient boundary layer

  Turbulent boundary layer profiles fitted to eq. 3.16. Symbols as in table 3.1 Rona et al

	Re θ • 300 * 697 △ 1003 • 1430 ⊳ 2900
	× 3654 ⊲ 5200
	12633
	♦ 13000
	+ 22845
	⋆ 31000
	10
	Figure 3.5:

Table 3 . 3

 33

	3. Inflow conditions and asymptotic modelling
		Testcase	β c		Re θ	H		G
		ZPG	0	350 -525 1.60 -1.57	-
		APG1	0.24 390 -620 1.62 -1.57	7 -6
		APG2	0.65 430 -690 1.64 -1.63 8 -8.3
		Table 3.2: Description of Skote's testcase [148]
	Testcase 'Skote' data Re τ	u e u τ	F 1	H	G	Re θ Re τ	Re θ
	ZPG1	u200		222 19.54 3.0 1.59	7.2	1.9	422
	ZPG2	u350		272 20.45 3.3 1.54	7.2	2.3	588
	APG1	u350		251	20.6 3.8 1.58	7.5	2.4	606
	APG2	u335		251	21.7 4.4 1.625 8.35 2.7	681

: Analysis of the Skote's data Michel's model.

Table 3 .

 3 .3), and R θ have been calculated. The figures are rounded off to 1 or 2 digits. All the results are given in table 3.3. 4: Skote's testcase, asymptotic analysis of the numerical data produce coherent values of the mean parameters. Since all the numerical data were not treated, one can state that, in the paper[START_REF] Skote | Direct numerical simulation of self-similar turbulent boundary layers in adverse pressure gradients[END_REF], the range in R

	By comparing the both tables 3.2 and 3.3, one can observe that the post-treatment

θ in table 3.2 is a little bit under-estimated and that the shape factor in the both case ZPG2 and APG2 are over-estimated.

By maintaining constant pressure gradient parameter and Reynolds number Re τ , corresponding to a mean value of the Reynolds number Re θ given in Skote's paper, the four testcases have been carried out with the asymptotic analysis. The results are given in table 3.4.

For very low values of Reynolds number Re τ , the asymptotic method never converges.

It has been explained in a previous section that it is not possible to join the internal and external region in an intermediate log-law region. DNS are not restricted by the Reynolds number, naturally.

  in 2009 and have been done in collaboration with Dr Aldo Rona, from Leicester University, UK. The AVBP solver . . . . . . . . . . . . . . . . . . . . . . . . . 90

	Chapter 4
	Numerical simulation and LES
	models
	Contents
	4.1 de l'adaptation de maillage et d'améliorer la précision des solutions. C'est un code par-
	allèle qui résout les équations de Navier-Stokes compressibles en régime laminaire et
	turbulent, en 2D et 3D. Des cas stationnaires et instationnaires peuvent être traités.
	Pour le calcul des cas instationnaires turbulent plusieurs modèles de sous-maille pour la

4.2 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . . 91 4.3 Large Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . 104 4.4 Governing equations for LES . . . . . . . . . . . . . . . . . . 105 4.5 Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . 113 4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Résumé étendu en français Simulation numérique et modèles LES Le code de résolution AVBP Ce chapitre est consacré au code de calcul utilisé pour simuler l'écoulement de cavité en explicitant la méthode numérique, les conditions aux limites basées sur les caractéristiques et la modélisation LES utilisée. AVBP est un code dévelopé au CERFACS capable de simuler des écoulements sur des maillages de tout type. L'emploi de maillages hybrides est faite dans l'objectif de profiter de l'efficacité des maillages non-structurés, LES ont été implémentés. Bien qu'étant au départ dédié à l'aérodynamique externe, il a été étendu aux configurations internes et mêmes pour des écoulements réactifs. La loi d'Arrhenius permet d'étudier la combustion dans des configurations complexes.

  [START_REF] Ahuja | Effects of cavity dimensions, boundary layer and temperature on cavity noise with emphasis on benchmark data to validate computational aeroaoucstic codes[END_REF]

	Cette approximation est applicable à tout type de cellule et donc au maillage hybride. Le
	résidu 4.2 est tout d'abord calculé pour chaque élément en faisant une simple intégration
	sur les faces. Quelque soit le maillage on essaie d'obtenir des triangles (déjà existants
	ou crées par découpage). La valeur du flux est obtenue en moyennant sur quatre trian-
	gles (deux divisions suivant la diagonale). Cette technique 'linear preservation property
	'permet, dans l'algorithme, de préserver la précision sur un maillage irrégulier. Sous

forme discrétisée l'équation 4.2 sur un volume arbitraire s'écrit

  5) où D k Ω j est une matrice de distribution qui transfert le résidu des centres des cellules Ω j au noeud k (scatter operation), et V k est le volume de contrôle associé à chaque noeud. La conservation est assuré si k∈Ω j D k Ω j = I. Dans le présent contexte, l'équation 4.5 est résolue pour obtenir la solution d'état stationnaire en utilisant le pas temporel Euler ou Runge-Kutta.

	La famille des schémas concernée utilise la définition suivante pour la matrice de
	distribution:

  La discrétisation temporelle est explicite pour tous les schémas numériques dans AVBP. Numerical simulation and LES models où u est la vitesse de propagation d'une perturbation dans l'écoulement, a ∞ la vitesse du son, ∆x est la longueur de la maille et CFL est le nombre de Courant-Friedrichs-Lewy. La valeur CFL nécessaire pour la stabilité dépend du schéma choisi. Dans AVBP, il est fixé à 0.7.

	Le schéma de Lax-Wendroff
	Les principaux schémas convectifs sont le schéma de Lax-Wendroff (LW) de Lax &
	Wendroff [81, 82] en volume fini et le schéma à deux pas de Taylor-Galerkin (TTGC)
	de Colin & Rudgyard [19] en éléments finis. Ces deux schémas sont respectivement de
	second et de troisième ordre en temps. Le schéma diffusif est typiquement un schéma
	compact de second ordre. Les éléments utilisés par AVBP sont des triangles et quadran-

L'implémentation de ce type d'approche est aisée et le temps de calcul par itération est faible. Le schéma explicite a cependant un pas de temps ∆t limité par le critère de stabilité: ∆t < CF L min (∆x) max | u | +a ∞ (4.10) 4. gles en 2D et des tetrahèdres, prismes, pyramides et hexahèdres en 3D. L'intégration temporelle est explicite pour assurer la précision. Le schéma de Lax-Wendroff (précis au second ordre en temps et en espace) est basé sur l'expansion de Taylor expansion en temps pour la solution de U .

  .34)

	En substituant 4.33 et 4.34 dans 4.31:

  Numerical simulation and LES models qui vérifient si la viscosité artificielle est nécessaire ou pas. On a donc un paramètre de calibrage ζ Ω j pour chaque cellule Ω j égal à zéro ou un. ζ Ω j = 0 quand la solution est bien résolue et donc pas d'utilisation de viscosité artificielle et ζ Ω j = 1 quand il Pour chaque cellule Ω j , le capteur de Jameson ζ J Ω j est le maximum sur tous les noeuds de la cellule du capteur par noeud de Jameson ζ k :

	faut l'utiliser. Dans AVBP il y a deux capteurs "capteur de Jameson"ζ J Ω j [69] et celui
	"capteur de Colin"ζ C Ω j [18].
	Capteur de Jameson
	Les méthodes de discrétisation spatiales dans AVBP sont centrées. Il est connu que ces
	schémas présentent des petites oscillations autour de la solution. Il est d'usage d'ajouter

de la viscosité artificielle aux équations discrétisées pour lisser les très forts gardients.

Les modèles utilisés dans AVBP sont basés sur la combinaison des termes de capture de choc (second ordre) et un terme de dissipation (quatrième ordre). Il y a des capteurs 4.

Table 4 .

 4

	99)

Pour plus de détails et comparaisons entre les différentes conditions voir Nicoud &

Poinsot

[START_REF] Nicoud | Boundary conditions for compressible unsteady flows[END_REF]

et le manuel de AVBP

[10]

.

Conditions aux limites dans AVBP

Plusieurs types de conditions aux limites existent dans AVBP en raison de son utilisation pour plusieurs applications. Table 4.1 donne les conditions aux limites utilisées dans cette étude. 1: Conditions aux limites.

Aux parois on utilise des conditions d'adhérence. Pour la condition thermique, la paroi peut être adiabtique ou isotherme. en entrée et sortie les relations caractéristiques (entrantes et sortantes) (NSCBC). En sortie une relxation est appliquée à la pression.

  4.2. Numerical method the scaling by n d . It should be noticed that no weighting is required for the LW term because it is computed on the dual cell. Substituting equation 4.104 and 4.120 into equation 4.118 leads to:

  SensorsA sensor ζ Ω j is a scaled parameter which is defined for every cell Ω j of the domain that takes values from zero to one. ζ Ω j = 0 means that the solution is well resolved and that no artificial viscosity should be applied while ζ Ω j = 1 signifies that the solution has strong local variations and that artificial viscosity must be applied. This sensor is obtained by comparing different evaluations (on different stencils) of the gradient of a given scalar (pressure, total energy, mass fractions, etc. ). If these gradients are identical, then the solution is locally linear and the sensor is zero. On the contrary, if these two estimations are different, local non-linearities are present, and the sensor is activated. The key point is to find a suitable sensor-function that is non-zero only at places where stability problems occur. Two sensors in AVBP which are used in this work are "Jameson-sensor"ζ J Ω j from Jameson et al

	steep solution variations. It is a common practice to add artificial viscosity (AV) term
	to the discrete equations to avoid spurious modes and in order to smooth very strong
	gradients. we describe here the different artificial viscosity methods used in AVBP. The
	different artificial viscosity models which are used in this work are characterized by the
	linear preserving property which leartificial viscosityes unmodified a linear solution on
	any type of element. The models are based on a combination of a "shock capturing"term
	(called 2 nd order artificial viscosity) and a background dissipation term (called 4

th order artificial viscosity). In AVBP, adding artificial viscosity is done in two steps. Initial a sensor detects if artificial viscosity is necessary, as a function of the given flow characteristics. Then a certain amount of 2 nd and 4 th artificial viscosity is applied, depending 4.2. Numerical method on the sensor value and on user-defined parameters.

Table 4 .

 4 3: Inflow boundary : Correspondence between the ∂W notation and the strength array in the AVBP implementation (in 3D) corresponding to this boundary condition is thus :

	The set of equations
	123

Table 4 .

 4 Numerical simulation and LES models 4: Outflow boundary : Correspondence between the ∂W notation and the strength array in the AVBP implementation (in 3D)

		Outflow boundary
	Type	Way	AVBP	in derivations [103]
	acoustic wave	out	strength(4)	∂W 2
	entropy wave	in	strength(1)	∂W s
	transverse shear	in	strength(2)	∂W 3
	transverse shear	in	strength(3)	∂W 4
	acoustic wave	in	strength(5)	∂W 1
	species waves	in	strength(5 + k)	∂W 4+k

  Contents 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 5.2 Two-dimensional cavity . . . . . . . . . . . . . . . . . . . . . 136 5.3 Three-dimensional rectangular cavity . . . . . . . . . . . . . 161 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164 Des conditions aux limites de type caractéristiques sont introduites sur les frontières du domaine, hors paroi avec des paramètres de relaxation pour éviter les influences nonphysiques des quantités entrantes et sortantes. Les conditions sont décrites dans le chapitre précédent, mais les conditions et les paramètres, pour respectivement les deux premiers cas et le dernier cas tests, sont fournis dans les tableaux 5.6 et 5.7. Re τ et Re θ ). Il s'avère, comme on pouvait s'y attendre, qu'en dehors de la couche interne, la couche limite produite par l'écoulement n'est pas une couche limite turbulente en équilibre. On le remarque en particulier sur les valeurs de la vitesse adimensionnelle u + e largement sous évaluées par le modèle asymptotique. Les différences entre les solutions asymptotiques et les solutions numériques s'accroissent avec la vitesse infinie amont. L'existence de la cavité modifie aussi le gradient de pression longitudinal dans la couche limite amont, l'influence est négligable pour un mode de cisaillement et plus importante pour un mode de sillage. Ceci semble se confirmer en observant l'évolution de la contrainte de frottement pariétale pour le cas U5.8 (figure 5.8). La contrainte τ w la profondeur D de la cavité et la vitesse infinie amont. Les points de mesure se trouvent répartis longitudinalement juste au-dessus du niveau de la cavité.Les comportements sont similaires pour les deux cas U20 et U40 bien que les phénomènes soient toujours accentués lorsque la vitesse augmente. Ainsi la vitesse u est parfairement périodique après une certain temps T * On observe des pics importants réguliers sur la vitesse verticale au voisinage du coin aval de la cavité. Ils représentent les signatures des structures tourbillonnaires de fortes intensités qui impactent le coin aval tout en sortant de la cavité. Pour ces deux cas ce sont typiquement des modes de sillage 2D (wake mode) qui apparaissent. Cela semble en accord avec les travaux de Rowley[START_REF] Rowley | simulation, and control of cavity flow oscillations[END_REF] où l'influence de la valeur relative de l'épaisseur de quantité de mouvement θ/L est l'un des deux 5. Analysis of the cavity flows paramètres déterminant (le second est le nombre de Mach). Pour la cas à faible vitesse (figure 5.11 ), l'évolution temporelle de la vorticité ω z montre une périodicité dans la dynamique de l'écoulement arrivant beaucoup plus tôt. On remarque que la vorticité demeure négative, comme dans une couche limite, au point notée 3 sur la figure, ce qui indique l'absence d'un décollement. La vorticié positive sur la paroi verticale aval de la cavité indique l'existence de décollements, de cisaillements forts et la naissance d'un tourbillon contrarotatif. Enfin pour le point de mesure situé à la limite verticale de la cavité et proche de la partie aval, on observe que la vorticité est essentiellement négative tout en ayant des très faibles valeurs positives, ce qui indique l'existence d'un mode de cisaillement, contrairement au cas à plus grandes vitesses. Les courbes ne semblent pas montrer de changement de mode (de cisaillement vers le mode de sillage ou inversement) au cours du calcul, et pour les trois cas tests, comme cela est parfois observé dans certaines simulations. Sur la Figure 5.12 est présentée l'évolution de la vorticité instannée ω z au cours d'une période et pour le cas U40. Elle est typique du mode de sillage. Un tourbillon négatif se forme dès le début de la couche de cisaillement, sur le bord amont de la cavité, donnant naturellement naissance à un second tourbillon de vorticité opposée au fond de la cavité. Larsson et al [80] ont reporté une fréquence fondamentale de St L = 0.245. Colonius et al [21] ont trouvé similairement St L = 0.248. Shieh & Morris [144] ont donné St L = 0.216 pour une cavité identique mais pour un nombre de Mach M = 0.6 dans un écoulement turbulent.Les courbes iso-rms des fluctuations turbulentes u ′ 2 et u ′ 2 , pour le cas U5.8 sont montrées respectivement sur les figures 5.21 et 5.22. On peut les comparer aux mesures expérimentales, les variances u ′ u ′ et v ′ v ′ , fournies par Haigermoser et al[START_REF] Haigermoser | A timeresolved particle image velocimetry investigation of a cavity flow with a thick incoming turbulent boundary layer[END_REF].Les figures5.26, 5.27, 5.28 et 5.29 présentent des profils des fluctuations turbulentesu ′ u ′ , u ′ v ′ et v ′ v ′ ,dans la couche limite en amont de la cavité et dans la cavité, pour les

	5. Analysis of the cavity flows
	Ecoulement de cavité 2D Les figures 5.14(a) et 5.14(b) présentent les profils de vitesse longitudinale moyenne dynamique. tions acoustiques ne semblent pas être propagées par l'écoulement en aval de la cavité.
	dans trois sections de la cavité et pour les deux cas de vitesse aval U20 et respectivement Similairement, dans les expériences de Haigermoser [60], les lignes de niveau étaient
	Présentation des cas tests et du calcul U5.8. On observe pour tous les cas, les fortes variations de vitesses au sein de la concentriques par rapport à la cavité, sans indiquer de directivité, et le maximum était
	La géométrie de la cavité et le domaine de calcul sont décrits sur la figure 5.1 et cavité, indiquant une forte activité tourbillonnaire en moyenne, ainsi qu'une zone où cas U5.8 et U40. On remarque que dans la zone de couche limite, les différents profils aussi de 92dB.
	dans le tableau 5.1. Le rapport d'aspect de la cavité est L/D = 4. Le maillage est Evolution de la couche limite cette activité est plus réduite, lorsque sur les courbes, la pente est autour de zéro. On sont typiques, mais que l'intensité est beaucoup plus faible pour le cas U5.8 que pour le Pour le cas U40 (figure 5.31) Le maximum de bruit est situé sensiblement à la
	dessiné en partie sur la figure 5.2 et est composé d'environ 225 000 points. Le maillage correspondant à la cavité est de 221 × 121 cellules. Des maillages plus grossiers ont été testés mais il s'avère que la complexité de la dynamique de l'écoulement, avec l'apparition de très forts cisaillements ou gradients aux voisinages des parois de la cavité ou des parois horizontales amont et aval, a nécessité une densification des mailles et donc un accroissement de points dans ces zones. La première maille dans l'écoulement est fixée à une distance d'environ y + = 2 des parois. Il n'existe pas de zone tampon pour amortir les remarque que le cisaillement est naturellement beaucoup plus faible dans le cas de mode cas U40, de pratiquement deux ordres de grandeurs. Dans la cavité les évolutions sont même position, confirmant finalement ce qu'on savait déjà, mais l'amplitude a augmenté Naturellement, avant d'atteindre la cavité (x/D ≤ 5), l'écoulement se comporte comme de cisaillement que du mode de sillage. comme attendues, très différentes entre le cas avec un mode de sillage et celui avec un jusqu'à 134dB. On remarque que les forts tourbillons présents dans l'écoulement (cavité une couche limite turbulente. Nous avons donc chercher à qualifier cette couche limite. Une analyse spectrale a été menée dans le but de valider les résultats obtenues. Les mode de cisaillement. On retrouve des ordres de grandeurs corrects. Les fluctuations et en aval); associés au mode de sillage perturbent fortement la progapation des pertur-Pour les deux cas où la vitesse infinie amont est de 5.8 m/s et de 20m/s, sont signaux sont mesurés au voisinage du coin aval de la cavité, là où les phénomènes pour le mode de cisaillement sont qualitativement cohérentes avec les résultats de Bertier bations acoustiques, tout en les amplifiant. Une faible directivité vers l'amont semble représentés respectivement le profil de vitesse adimesnionnelle u + dans cette région de physiques sont les plus intenses. Les spectres de l'énergie associée à la perturbation et al [5] et sont typiques des écoulements décollés. Dans toutes ses courbes, des pics existée. Ahuja & Mendoza [2] ont rapporté dans leur expérience une très faible direc-l'écoulement sur les figures 5.4 et 5.5. Des comparaisons ont eu lieu avec des profils de de vitesse longitudinale, et pour les deux premiers cas de vitesse bien que relative-locaux peuvent être associés aux évènements rencontrés dans la cavités (décollement, tivité (à la perpendiculaire de la cavité) pour la même géométrie. Rowley et al [132], couche limite en équilibre dans des conditions équivalentes (même nombres de Reynolds ment chaotiques, permettent de retrouver approximativement la pente en -5/3 pour les forts cisaillements, centre de rotation). dans des simulations numériques directes a trouvé une directivité autour de 135 • , mais
	instabilités qui pourraient apparaître sur les frontières ouvertes du domaine de calcul. Les hautes fréquences. Deux pics, respectivement pour un nombre de Strouhal de 0.194 et pour une cavité de rapport d'aspect de 2, et pour un nombre de Mach de 0.6.
	Résumé étendu en français Analyse des écoulements cavités Dans les écoulements compressibles à faibles nombres de Mach (subsoniques), l'amplitude des perturbations acoustiques est d'un ou deux ordres de grandeurs plus faibles que les amplitudes des instabilités ou des fluctuations d'origine hydrodynamique. Ceci ex-différents problèmes numériques rencontrés associés à une sous définition du maillage, de 0.205, et pour les cas U20 et respectivement U40, sont observables. Pour un cas test Ecoulement de cavité 3D à des instabilités numériques, à un manque de précision dans le calcul des forts gradients et au choix des conditions aux limites avec des réflections d'ondes non physiques ne seront pas développés en détail dans le document, bien qu'ayant pris une part très importante du temps de travail. Les équations modèles utilisées (LES) et toutes les méthodes numériques employées par le code AVBP, qui est notre outil de base, ont déjà été décrites avec précision dans les chapitres précécents. Le tableau 5.2 rappelle les paramètres numériques importants des simulations associés aux 3 principaux cas tests nommés U20, U40 and U5.8. Le nombre dans le nom indique la vitesse infini amont du cas en m/s. Rappelons que les conditions à l'entrée du domaine ne comportent pas de champ fluctuant turbulent. Seulement le profil moyen d'une couche limite turbulente est pre-scrit, correspondant à une épaisseur de quantité de mouvement θ fixée par les cas tests, tout en sachant que le rapport entre cette épaisseur θ et la profondeur et la longueur de la cavité joue un rôle crucial dans la dynamique de l'écoulement et dans la génération du bruit. Le dernier cas test à basse vitesse correspond à des paramètres très proche des décroît normalement avec la distance longitudinale, puis semble 'diverger' au voisinage du début de la cavité. Des comparaisons avec la couche limite aval à la cavité ne sont pas possibles car celle-ci est fortement décollée et instationnaire. Ecoulement autour et dans la cavité L'écoulement dans la cavité est fortement instationnaire, mais devient périodique après un certain temps, comme l'atteste les figures 5.9, 5.10 où sont montrées les évolutions temporelles des deux composantes de la vitesse (u, v) en fonction d'un temps adimen-l'aval. Il finit par s'éclater sur le coin aval. On observe au cours du cycle des zones à très fort cisaillement dans la cavité, zones qui numériquement nécessite de la viscosité artificielle pour rester stable. Juste après le coin aval, le tourbillon positif initialement dans la cavité est éjectée de celle-ci pour aller se mélanger avec le tourbillon négatif déjà éclaté, ce qui entraîne, sur la paroi aval, une très forte instationnarité de la couche limite et de forts décollements, suite au passage des différentes structures tourbillonnaires. Ces évolutions complexes vont générer du bruit, lorsqu'un tourbillon impacte le coin aval, et des fortes augmentations de traînée à cause des décollements et des forts gradients de pression favorables puis adverses qui apparaîssent cycliquement. une sorte de langue de vorticité négative typique. La longeur d'onde des instabilités Le calcul en trois dimensions a pour ojectif de se rapprocher de cas d'étude réels. Ainsi dans certaines configurations 2D, on observe des modes de sillage dans la cavité qui finalement n'existent plus lorsqu'on effectue un calcul en 3D. Ces modes de sillage sont Conclusion Une étude numérique d'un écoulement turbulent de cavité 2D et 3D a été menée, pour un rapport d'aspect de 4, et pour 3 cas de vitesse amont de l'écoulement. Le modèle très similaire au cas U40, Les spectres des vitesses longitudinales instantannées, pour les mêmes cas et points sont tracés sur les figures 5.18 et 5.19. On observe plus disctinctement la présence des modes fondamentaux mais aussi l'existence d'harmoniques. Pour le cas à 40m/s il y a un signal à basse fréquence qui semble correspondre à une résonance acoustique de tube, associée à la longueur du domaine de calcul. On remarque que les zones de forte activité de la turbulence sont différentes entre la du coin aval de la cavité alors que dans la simulation, à part pour un pic sur le coin, d'ailleurs difficilement observables dans les expériences quand ils existent. D'un point de vue numérique, le domaine bidimensionnel a été étendu par translation de plan à la troisième dimension en diminuant la résolution 2D. Ceci est un problème important en terme de qualité des résultats obtenus. Les conditions aux limites dans la troisième direction sont celle d'une symétrie suivant la dénomination dans le code AVBP. Sur la figure 5.34 où est montrée la vorticité instantannée ω z , on observe l'existence d'un mode de cisaillement. Une vue plus détaillée montrerait les structures turbulentes tridimensionnelles dans la zone en aval de la cavité. Ce travail est malheureusement trop préliminaire et mériterait d'être approfondi. ations turbulentes au sein de la cavité sont propagées dans tout l'écoulement. A partir LES dans des cas à très faibles nombres de Mach a été utilisé. Différentes conditions d'entrée (profil de couche limite turbulente) ont aussi été introduites, dont l'un basé sur la solution asymptotique de la couche limite d'équilibre. L'écoulement turbulent, en valeur moyenne et en termes de fluctuations a été analysé. Le mode de cisaillement a été obtenu pour le cas correspondant à la plus faible vitesse, choisi pour correspondre au cas expérimental de Haigermoser. Pour les autres cas, seul le mode de sillage a été trouvé en 2D. Il disparaît en 3D conformément à la littérature sur le sujet. Les niveaux de pression acoustique ont été calculés pour deux cas présentant soit le mode de cisaillement, soit le mode de sillage. Globalement une faible directivité vers l'amont a été montrée. La vitesse est un facteur qui accroit les niveaux maximums, et les lignes de sionnel basé sur Le premier tourbillon capte de l'énergie au fluide, s'amplifie tout en étant convecté vers La figure 5.13 présente la même vorticité au cours d'un cycle mais pour le cas simulation et les expériences, par contre les valeurs des intensités sont très proches Aéroacoustique pour les cas 2D niveaux sont fortement perturbés lorsque l'écoulement turbulent est décollé ou chahuté à faible vitesse U5.8. L'existence du mode de sillage est très visible, puisqu'on observe (attention à la racine carrée). Dans l'expérience l'activité est concentrée au voisinage Les ondes de pression générées par la dynamique tourbillonnaire intense et les fluctu-par des structures tourbillonaires relativement intenses.
	plique qu'en général les outils (numériques, expérimentaux) pour déterminer les champs expériences de Haigermoser [60] qui a participé au programme AeroTraNet à Turin (voir (oscillations longitudinales) est aussi visualisé. La dynamique tourbillonnaire est cette l'activité turbulente est plus proche du fond de la cavité. Finalement un accord sur les des simulations LES, et en appliquant l'analogie de Lighthill-Curle, il a été possible de
	moyens et fluctuants et les perturbations acoustiques soient différents. Ainsi dans ce les paramètres dans le tableau 5.5). La forte épaisseur relative de la couche limite par fois ci concentrée au voisinage de la paroi verticale aval. Essentiellement deux tourbillons niveaux est obtenu malgré les énormes différences existantes entre les conditions de la déterminer le niveau de bruit, mesuré en SPL dedans mais surtout autour de la cavité.
	chapitre, l'écoulement est résolu à l'aide de simulations numériques de grandes échelles rapport à la hauteur de la cavité est illustrée, pour le premier cas test U20, par la figure de signes opposés restent prisonniers de la cavité, tout en déstabilisant la couche de simulation et l'expérience ainsi que les nombreuses incertitudes sur les conditions exactes Dans l'analogie, seules les intégrales de surface ont été prises en compte, les intégrales
	(LES) alors que le champ de pression acoustique instationnaire est déterminé par l'analogie 5.3 (δ/D = 2.2). Pour les deux premiers cas tests, le profil moyen turbulent est fourni cisaillement et en l'alimentant en vorticité positive. La vorticité positive est générée de l'expérience. de volumes étant souvent négligées pour les écoulements à faible vitesse.
	acoustique de Lighthill-Curle. par une loi en puissance alors que dans le dernier cas test, nous avons introduit un par des instabilités centrifuges induites par la rotation forcèe du fluide. La couche de Les figures 5.24 et 5.25 montrent les mêmes quantités pour le cas test U40. L'existence Les iso-contours de niveau de pression, en SPL, sont dessinés sur la figure 5.30 et
	Sont présentés successivement, l'étude numérique des écoulements de cavité bidimen-profil de couche limite en équilibre, calculé très précisement par l'approche asymptotique cisaillement vient ensuite s'éclater sur la coin amont, moins vigoureusement que pour les du mode de sillage introduit une localisation différente de l'activité turbulente. Le coin pour le cas U5.8. Globalement, loin de la cavité, les lignes de niveau sont concentriques
	sionnels puis tridimensionnels et l'étude aéroacoustique. Une conclusion clotûre ce décrite dans le troisième chapitre. cas à plus grande vitesse. Les zones de cisaillement dans la cavité semblent aussi moins aval est toujours un point avec un maximum local, mais on oberve aussi plusieurs zones par rapport à la cavité avec une légère directivité (orientabilité) vers l'amont, confirmant
	chapitre. Notons que l'étude 3D, incomplète, a juste pour objectif de montrer la lim-importantes. En aval de la cavité, la couche limite reste attachée tout en convectant les d'agitation : au sein de la cavité, près de la paroi du fond, dans les zones décollées et que la source sonore se situe bien au voisinage du coin aval de cette cavité. Le maximum
	itation de l'étude 2D. instabilités issues de la couche de cisaillement et énergétisées dans la cavité. sur la paroi en aval de la cavité où les tourbillons éjectés de la cavité créent leur propre est de l'ordre de 92dB, situé légèrement dans la cavité près du mur aval. Ces perturba-

Table 5 .

 5 2: Numerical Parameters of the two-dimensional test cases U20, U40 and U5.8

	5. Analysis of the cavity flows
		U20	U40	U5.8
	Turbulence model	Filtered Smagorinsky Smagorinsky Filtered Smagorinsky
	Integration scheme TTGC	TTGC	TTGC
	Time integration	Runge-Kutta	Runge-Kutta Runge-Kutta
	Time step (△t) CFL number	2.90 × 10 -8 s 0.7	2.88 × 10 -8 s 2.88 × 10 -8 s 0.7 0.7
	Fourier number	0.3	0.3	0.3
	Artificial viscosity SLK sensor	SLK sensor	SLK sensor
		Artificial viscosity coefficients
	4 th order	0.05	0.05	0.05
	2 nd order	0.2	0.2	0.2

  2 , O 2 , CO 2 , CH 4 , etc). In this work, neither combustion nor chemical reaction between species is studied. Therefore only two gases Nitrogen (N 2 ) and Oxygen (O 2 ) are handled and its combination results in air which is approximately 80% nitrogen and 20% oxygen by volume. The numerical methods which are used in the computations are the Lax Wendroff scheme from section 4.2.5 and the Two step Taylor Galerkin Colin (TTGC) scheme 4.2.6. Classic Smagorinsky model (with model constant C s = 0.18) and filtered Smagorinsky model (with model constant C S F = 0.37) are used to determine the turbulent viscosity ν t . More details of which can be obtained from chapter 4. The table 5.2 summaries the numerical parameters followed in the two-dimensional test cases: U20, U40 and U5.8. These three test cases are named with respect to the stream wise velocity values u ∞ = 20 m/s, u ∞ = 40 m/s, u ∞ = 5.8 m/s.The AVBP is a parallelised solver. All two-dimensional simulations are performed on the super computers which are given in the table 5.3.

Table 5 .

 5 4: Flow parameters of the test cases conducted.

	5.2. Two-dimensional cavity

Table 5 .

 5 5: Flow parameters of the test case carried by Haigermoser[START_REF] Haigermoser | Investigation of cavity flows using advanced optical methods[END_REF]

	was observed by Colonius et al [22]. For cavity flows, the ratio determining the mode (shear mode or wake mode). The flow parameters relevant to the L plays a major role in θ
	three cases in this work are given in the table 5.4.
	The inlet conditions of the first two test cases are imposed with mean turbulent
	boundary layers where as for the test case U5.8, the boundary layer profile is generated
	using equilibrium turbulent boundary layer approach (ETBL) (see section 3.4). This
	test case is similar to the Particle Image Velocimetry (PIV) experiment in a water
	cavity flow carried out by Haigermoser [60]. The characteristics are summarised in
	the table 5.5. The turbulent boundary layer in the experiment is thick and contains
	naturally turbulent quantities where as in the present numerical work, the turbulent
	boundary layer which is imposed at the inlet of the computational domain does not
	carry any turbulent quantities.
	A mean stream wise velocity profile for the thick turbulent boundary layer of the test
	case U20 is shown in the figure 5.3. For this test case U20, the boundary layer thickness
	δ is 22 mm which is greater than the depth of the cavity (D = 10 mm). It should
	be noted that the ratio layer in non-dimensional form and the ratios corresponding to the test cases given in L is a measure used to express the thickness of the boundary θ
	the table 5.4.

Table 5 .

 5 6: Boundary conditions and corresponding values for the test cases U40 and U20

	5.2. Two-dimensional cavity

Details in table 5.6 relate the location of the boundaries and its respective boundary conditions with the parameters for the test cases U20 and U40 where as the table 5.7

Table 5 . 7

 57

: Boundary conditions and corresponding values for the test case U5.8 carries the boundary condition details of the test case U5.8.

  ≤ T * ≤ 174

				0.03 0.03
	0.03 0.03 0.03 0.03	0.03 0.03	0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03 0.03	0.03 0.03	0.03 0.03 0.03 0.03 0.15 0.15	0.03 0.03 0.03 0.03 0.09 0.09	0.03 0.03 0.03 0.03	0.03 0.03 0.03 0.03	0.03 0.03 0.03 0.03	0.0000 0.0375 0.0750 0.1125 0.1500 q 0.03 u ′2 0.03 0.03 0.03

  Figure 5.24: Mean fluctuating velocity u ′ 2 field for the test case U40.

			38 0.38 0.38 0.38
	0.38 0.38	0.38 0.38	0.38 0.38
	0.38 0.38	0.38 0.38	0.38 0.38
			0.500 0.375 0.250 0.125 0.000

  Cette thèse s'intéresse à l'étude par simulation numérique de l'écoulement d'une couche limite compressible bidmensionnelle arrivant sur une cavité en se focalisant sur les aspects à la fois dynamique et acoustique. Elle s'intègre dans l'action de formation initiale européenne Marie Curie appelée AeroTraNet.Dans notre cas la couche limite turbulente amont a été initialement définie avec un profil de vitesse moyenne en loi puissance. Une extension de l'approche assymptotique usuelle basée sur la loi déficitaire dans la région externe de la couche et sur la loi linéaire et la loi log dans la région interne a été proposée dans le cas d'une couche limite en équilibre. Spécifiquement, une nouvelle fonction du modèle de longueur de mélange, qui améliore l'accord avec les écoulements de plaque plane avec ou sans gradient de pression adverse, a été proposée. L'unique paramètre de cette fonction a été déterminé à l'aide des données expérimentales et a été modifié à l'aide de la simulation numériquedirecte d'une couche limite avec gradient de pression. L'accord est très bon avec le cas sans gradient de pression. Avec un gradient de pression adverse quelques désaccords existent qui peuvent s'expliquer en partie par le problème des conditions aux limites dans les calculs DNS, en partie par les limites de l'approche asymptotique surtout à faibles nombres de Reynolds et grands gradients de pression longitudinaux de l'écoulement.Le conditions aux limites basées sur les caractéristiques sont utilisées comme conditions non réfléchissantes sur les bords du domaines et aux parois. Les écoulements bidimensionnels de cavité sont simulés avec la couche limite turbulente amont épaisse. le rapport d'aspect de la cavité est de 4. L'épaisseur de la couche limite introduit un fort écoulement cisaillé au-dessus et dans la cavité qui peut osciller à la fréquence dominante, avec un mode de cisaillemnt à faible vitesse et au mode de sillage à vitesse modérée. Le cas test à faible vitesse qui est celui correspondant à l'expérience de Haigermoser [61] avec mode de cisaillement a permis de valider nos résultats. Les deux modes ont été analysés en examinant l'écoulement moyen, l'écoulement fluctuant et le champ de vorticité. La zone de turbulence forte a été mise en évidence grâce aux iso-valeurs de la moyenne des fluctuations turbulentes. La comparaison avec le champ des intensités turbulentes obtenues par Haigermoser a permis de monter l'accord parfait entre les calculs et l'expérience. L'intensité turbulente est liée au type de mode qui dépend du rapport de l'épaisseur de quantité de mouvement de la couche limite à la profondeur de la cavité et du nombre de Mach de l'écoulement. Passant d'une vitesse de 5.8 à 40m/s fait augmenter l'intensité turbulente d'un ordre de grandeur. Un large ensemble de champs de pression instantannée a été déterminé par simulations de grandes échelles et ont servi à alimenter l'analogie de Curle qui permet de calculer la pression acoustique dans la zone visée. Le niveau des émissions sonores est lié à la vitesse amont et au type de mode dans la cavité. Dans tous les cas une faible directivité (orientabilité) dans la direction amont a été observée. Le mode sillage avec une dynamique tourbillonnaire complexe impacte localement et directement sur le niveau des émissions sonores et sur l'orientabilité. Conclusions dre en compte la physique des écoulements turbulents complexes. Pour une grande plage du nombre de Reynolds, les grandeurs moyennes et les fluctuations turbulentes de couches limites turbulentes épaisses, déterminées expérimentalement peuvent être ajoutées à celles existantes pour fournir des conditions d'entrée plus réalistes en 2D et 3D. L'étude des cavités avec différents rapports d'aspect peut être faite avec des couches limites épaisses à différents nombres de Reynolds pour analyser les parmètres qui influencent les mécanismes d'instabilité et les modes.Dans le calcul des niveaux de pression sonore à faibles nombres de Mach, l'intégrale de volume a été négligée dans le code Matlab R . Pour améliorer la précision de l'évaluation des niveaux de pression sonore, l'intégrale de volume devra être ajoutée à celle de surface pour prédire l'acoustique à nombre de Mach élevé et pour des écoulements fortement compressibles. L'influence de l'épaisseur de la couche limite sur les niveaux de pression sonore et sur l'orientabilité de la propagation serait intéressante à poursuivre.

	Perspectives
	Avec l'augmentation des resources de calcul (de IDRIS, Paris et de CALMIP, Toulouse),

L'approche numérique est basée sur la simulation de grandes échelles dont les équations sont résolues par le code AVBP du CERFACS. L'analyse des sources et des propagations acoustiques est basée sur la théorie de Lighthill-Curle en collaboration avec C. Haigermoser. Couramment, l'écoulement sur une cavité peut être décrit en regardant trois parties: la couche limite turbulente amont qui se développe et croît avant d'atteindre la cavité, l'écoulement turbulent dans et sur la cavité qui peut être soit un mode de sillage ou un mode de cisaillement, et l'écoulement en val la cavité qui est une couche limite instationnaire décollée ou non et qui convecte les structures éjectées par la cavité. La simulation de Grandes Echelles (LES) avec le schéma de Taylor Galerkin Colin à deux pas (3 eme ordre en espace et en temps) a été utilisée pour simuler l'écoulement de cavité. Le modèle de Smagorinsky et Smagorinsky filtré ont été adoptés dans ce travail. A 20m/s, le cas bidimensionnel a été étendu au cas 3D pour montrer que le mode de sillage est un pur artefact de la simulation numérique 2D. En 3D seul le mode de cisaillement est observé. Dans le cas tridimenionnel, initialement l'écoulement reste uniforme en envergure et oscille suivant le mode sillage. Avec le temps le mode sillage disparaît et il est suivi du mode de cisaillement. La couche limite turbulente amont conduit à un détachement tourbillonnaire irrégulier provenant du coin amont de la cavité et à un battement basse fréquence de la zone cisaillée. L'explication est que les structures de la turbulence de cette couche limite amont continuent à exister dans la zone cisaillée de la cavité et impactent le coin aval pour générer les mécanismes d'instabilité. L'étude de l'aéroacoustique qui a suivi a consisté en l'analyse des émissions sonores et des sources. La simulation de grandes échelles et la simulation directe de l'écoulement de cavité pourront être réalisées avec de grands domaines de calcul et des maillages fins. L'amélioration des modèles comme WALE, Smagorinsky filtré et dynamique permettra de mieux pren-

Courtesy:Paul Dopson, APG Photography.

For bilinear and trilinear elements (quads, hexahedra and pyramids for example) the gradient of the shape function over the element is no more constant. This difficulty is overcome by adding a correction to the residual computed as for linear element.

fournis, par lui-même et les thésards travaillant avec AVBP. This research project AeroTraNet has been supported by a Marie Curie Early Stage Training Fellowship of the European Community's Sixth Framework Programme under contract number MEST CT 2005020301. Thanks to European Commission. Thanks to IDRIS, Paris and CALMIP, Toulouse for the computing facilities. Thanks to Laia and Kaushik, my colleagues from AeroTraNet project at IMFT for the discussion, friendship, support, encouragement, motivation and help throughout my stay in Toulouse. Thanks to all the members from AeroTraNet project at Politecnico di Torino, Università di Roma Tre and University of Leicester during the project meetings. This AeroTraNet is an Early Stage research Training (EST) in the European Research Area (

the two-dimensional cases. In the region above the cavity, the grid is stretched using exponential function. Like the two-dimensional case, the grid is refined near the walls (to resolve the gradients).

Numerical schemes and LES Model

Large eddy simulation was performed on a three-dimensional test case 3DU20 for the velocity of u ∞ = 20m/s and Re D = 13.68 × 10 3 . A velocity profile generated by power law (see subsection 3.2.3) was imposed at the inlet of the domain. Parameters followed in this three-dimensional test case are given the table 5.8.

Boundary conditions

Results

As it is mentioned earlier that the purpose of conducting the three dimensional simulation is to demonstrate the wake mode in the three-dimensional cavity. It should be Conclusions integral of the Curle's analogy was neglected in the code.Apart from improving the accuracy in determining the sound pressure level, inclusion of the volume integral with the existing surface integral could help to capture the acoustics for high Mach number and compressible flows. Influence of boundary layer thickness on the sound pressure level and the directivity of the propagation of the sound could be an interesting study to continue.