Madam Podvin

Mr Madam Nore

Cordier

Firstly, I would like to express

Résumé

L'amélioration des performances des avions de transport est incontournable pour réduire les émissions polluantes à la source. L'un des moyens vise à réduire la traînée totale en agissant sur les aspects aérodynamiques. Dans ce contexte, l'étude des couches limites turbulentes est importante car la traînée de frottement peut représenter jusqu'à plus de 50% de la traînée totale en vol de croisière. La simulation numérique directe (DNS) de couches limites turbulentes nécessite une résolution spatiale élevée, en particulier près des parois pour capturer l'ensemble des échelles de longueur inhérentes à la turbulence. Il a été démontré que le nombre de points de grille (N ) devrait varier N ∼ Re 13/7 dans le cas d'une simulation des grandes échelles (LES) bien résolue et en N ∼ Re 37/14 dans le cas d'une DNS. Il est donc particulièrement difficile, voire impossible d'effectuer des simulations bien résolues à des nombres de Reynolds très élevés. Pour réduire le coût de calcul lié à la haute résolution spatiale nécessaire dans la région de proche paroi, des techniques qui éliminent la région de paroi ont été étudiées. Les solutions classiques sont généralement basées sur l'introduction d'un modèle de paroi. L'idée poursuivie ici est différente et il s'agit un modèle implicite à partir des POD.

Dans cette approche, l'écoulement est calculé dans un domaine numérique réduit, ce qui exclut la région coûteuse de la paroi. Une condition limite instationnaire est définie à chaque pas de temps sur les variables conservatrices, et appliquée aux frontières du domaine. Cette procédure a été mise en oeuvre dans le cadre incompressible, où le champ de vitesse doit être imposé à toute la frontière. Podvin et Fraigneau ont proposé une condition de paroi basée sur la décomposition orthogonale aux valeurs propres (POD) pour générer le champ synthétique. L'objectif de ce travail est d'étendre au cadre compressible l'approche développée en POD incompressible. Un aspect clé de la formulation compressible est que, contrairement au cas incompressible (elliptique), les ondes peuvent se propager à travers les frontières, en particulier pour la densité et la pression, de sorte que les conditions limites doivent être adaptées, en utilisant la méthode des caractéristiques. Certaines variables sur la frontière peuvent être calculées directement à partir de l'intérieur du domaine, tandis que d'autres devront être imposées. L'idée est d'utiliser la POD pour construire les champs imposés sur la frontière. La POD consiste à représenter le champ comme la superposition de modes spatiaux dont l'amplitude est fonction du temps. La forme des modes spatiaux est supposée connue a priori. Le défi consiste alors à déterminer l'amplitude temporelle correcte des modes. Des DNS de l'écoule-ment turbulent du canal sont d'abord effectuées sur deux configurations d'écoulements pour générer des bases de données et calculer les modes spatiaux POD : un écoulement subsonique du canal à un nombre de Mach de M a = 0.5 et un nombre de Reynolds Re τ = 180, et un écoulement supersonique à M a = 1.5 et Re τ = 220. Pour évaluer les conditions aux limites synthétiques aux bords du domaine réduit, on construit les coefficients temporels POD en projetant les variables de l'écoulement sur les modes spatiaux. Les bases POD peuvent être calculées relativement aux variables primitives ou conservatives. Nous avons constaté que le meilleur résultat est obtenu à partir des variables conservatives parce que ce choix impose la conservation du flux massique à travers le plan où sont appliquées les conditions synthétiques. Une procédure de redimensionnement a cependant été introduite pour s'assurer que les rms des variables compressibles reconstruites correspondent aux valeurs rms de référence du plan reconstruit. Des calculs de l'écoulement de canal réduit ont ensuite été effectués à l'aide des conditions limites synthétiques basées sur la POD et mises en oeuvre à une hauteur spécifique au-dessus des parois. Les conditions aux limites synthétiques ont été testées pour deux hauteurs différentes. Les résultats obtenus dans la configuration de canal réduit ont été comparés à ceux des DNS de canal complet pris comme référence. Une bonne concordance est observée pour le profil de vitesse moyenne de Van-Driest, la température et la contrainte de cisaillement turbulente dans le cas subsonique. La procédure s'est avérée robuste compte tenu de la hauteur à laquelle est appliquée la condition synthétique. Dans l'écoulement de canal supersonique, les intensités turbulentes, la contrainte de cisaillement turbulente et les profils moyens de température, de densité et de vitesse longitudinale sont reproduits avec une bonne qualité. Une réduction du coût de calcul est obtenue dans le canal réduit pour les simulations d'écoulements subsoniques et supersoniques ; le temps CPU est réduit de 60% pour h + = 18 et 70% pour h + = 54.

2.5

Profiles of the Reynolds stress components in the incompressible channel flow coming from [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]: rms values of the velocity components (< u i u i >),

on the left, and the turbulent shear stress (< u w >), on the right. . . .

2.6

Mean profiles in the supersonic channel flow: mean density (top left), temperature (top right) and mean streamwise velocity component (bottom).

2.7 Profiles of the Reynolds stress components in the supersonic channel flow (< u i u i > and -< u w >), normalized by the friction velocity u τ :

present DNS results, on the left, and, on the right, reference profiles from incompressible DNS [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] (red line), and from supersonic DNS (M a =

List of Tables

Wall turbulence

Many geophysical and industrial flows are characterized by turbulence, which is associated with high friction drag. Friction drag can represent up to 50 % of the full drag for a commercial airplane with consequences for fuel consumption, larger operational range, higher achievable speeds and greater endurance. This makes skin friction manipulation a crucial objective for drag reduction.

A considerable effort has been devoted to wall turbulence control ([39; 50; 57; 81; 84]).

However, despite the success achieved with a variety of approaches, our understanding of wall turbulence remains incomplete. Numerical simulations of wall-bounded flows represent an invaluable tool to obtain information about the flow statistics and provide insight into its physics.

Turbulent flows are characterized by a high Reynolds number Re = U H/ν which measures the ratio between the advective time scale of the flow H/U , where H and U are respectively outer length and velocity scales, and the viscous time scale of the flow H 2 /ν, where ν = µ ρ is the kinematic viscosity. Viscous effects are predominant in a very thin layer close to the wall, which constitutes the main source of vorticity for the turbulence.

Further away from the wall, the flow is dominated by the outer scales. This gives rise to a very wide range of scales, which makes simulation of real-life applications particularly difficult, despite the increase in computing power. As an example, the range of scales on a commercial plane flying at a cruise speed of 900km/hr goes from O(10)m to a few microns. This means that the ratio between the largest and the smallest scales should be on the order of 10 6 , yielding 10 18 in three dimensions.

Mean-flow kinetic energy transferred to the wall layer by the Reynolds stresses is converted into turbulence kinetic energy (through turbulence production) and into heat (through viscous dissipation). A relevant measure of the turbulence of the flow is given by the dissipation . The smallest scale that can be expected to be found in the flow depends on the dissipation and the dynamical viscosity ν:

η = (ν 3 / ) 1/4 (1.1)
It is called the Kolmogorov scale.

To determine the intensity of the turbulent fluctuations a relevant scale can be derived from the wall shear

ρu 2 τ = µ dU dz | wall (1.2)
Defining the wall units based on the friction velocity u τ and viscosity ν, one can make the wall-normal distance nondimensional as follows

z + = zν/u τ ,
and the wall layer can be divided into three regions [68] • the viscous sublayer: 0 < z + < 5 : U = z +

• the buffer layer: 5 < z + < 30: overlap region

• inertial sublayer 30 < z + < 0.2Re τ : U = 1 κ log(z + ) + G where κ is the Von Karman constant (κ ∼ 0.4 and G ∼ 5). The production of turbulent energy is maximal at z + = 15. One can see that strong interactions take place in the wall layer. Due to the small scales present in the flow, a very high resolution is necessary. The number of grid points necessary to resolve the flow down from the largest scale to the smallest scale, which is the Kolmogorov scale, can be shown to increase like Re 3/4 . For a three-dimensional direct numerical simulation of a turbulent flow, this means that the total number of grid points should increase like Re 9/4 . This number has been recently reevaluated by Choi and Moin [START_REF] Choi | Grid-point requirements for Large Eddy Simulations: Chapmans estimates revisited[END_REF] to Re 37/14 . An important feature of turbulent flows is the presence of coherent structures, as reviewed by Dennis [START_REF] Djc | Coherent structures in wall-bounded turbulence[END_REF]. Several types of coherent structures include:

These region can be identified in

• an alternation of high and low speed longitudinal streaks (shown in Figure 1.2),

• cross-stream vortices whose axis is roughly aligned with the mean shear,

• hairpin or horseshoe vortices (shown in Figure 1.2),

• ejections of low-speed fluid away from the wall,

• sweeping motions of high-speed fluid away from the wall. These structures are connected and exist over a wide range of scales.

One important question which has not been settled so far is the relationship between the coherent structures and the different sublayers. On the one hand, Jimenez and Pinelli [START_REF] Jiménez | The autonomous cycle of near-wall turbulence[END_REF] showed that the dynamics of the wall region were essentially independent of the outer layer. On the other hand, Mathis et al. [START_REF] Mathis | Large-scale amplitude modulation of the small-scale structures in turbulent boundary layers[END_REF] established that the influence of the outer region can be felt very close to the wall. These apparent contradictions highlight the complexity of the wall region. We have seen that direct numerical simulation or DNS requires a very high resolution in space in the near-wall region, which in turns constrains the time step. The computing cost of even a moderate Reynolds number flow is very high.

Simulating turbulence: DNS, LES and RANS

A reduction of complexity can be reached with Large Eddy Simulation or LES (see [START_REF] Tenaud | Introduction à la Simulation des Grandes Echelles pour les écoulements de fluide compressible[END_REF] for a review). The basic idea is to solve accurately for the motions of all turbulence scales above a certain cut-off, so as to capture the large-scale interaction that occur in many unsteady, non-equilibrium processes. It is cheaper than DNS as it is not necessary resolve the flow down to the Kolmogorov or dissipative scales. For wall-bounded flows, the question is whether the wall region is resolved or modelled for LES, i.e. whether the near-wall eddies are resolved or modelled on the computational grid. Chapman

[15] estimated that the total number of grid points for wall-resolved LES (WRLES) is proportional to Re 9/5 , whereas that for wall-modelled LES (WMLES) increases like Re 2/5 . These estimates have been recently revised by Choi and Moin [START_REF] Choi | Grid-point requirements for Large Eddy Simulations: Chapmans estimates revisited[END_REF] to Re 13/7 for the case of wall-resolved LES and Re in the case of wall-modelled LES.

A still simpler representation of the flow is constituted by the Reynolds Average Navier-Stokes or RANS equations, which are the governing equations based on the Reynolds decomposition. These equations govern the evolution of time-averaged quantities and require the introduction of a closure scheme or turbulence model to account for the nonlinear effect of the fluctuations [68].

The RANS turbulence models can be divided into several categories: the zero-equation models [16; 68; 83] (algebraic model), the one-equation models [41; 74], and the twoequation models [47; 52; 53; 82]. Their main disadvantage is that they can only provide time-averaged information on turbulence, which may not be sufficient enough for prediction and engineering design [7]. They are essentially based on the definition of a turbulent viscosity, which may limit their relevance for transitional flows, swirling flows and turbomachinery [68].

Wall models

We can see that the flow structure in the near-wall region is very complex (see the turbulent boundary layer in Figure 1.4). Figure 1.4: Iso-surfaces of the second invariant of the velocity gradient tensor in the near wall region. The iso-surfaces are coloured with the local streamwise velocity [START_REF] Naqavi | Direct Numerical Simulation of a wall jet: flow physics[END_REF].

Several strategies were proposed to avoid the huge cost for the simulation in near wall region. A recent review can be found in Bose and Park [START_REF] Bose | Wall-modeled Large-Eddy Simulation for complex turbulent flows[END_REF]. Piomelli and Balaras [START_REF] Piomelli | Wall-layer models for Large-Eddy Simulations[END_REF] and Piomelli [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] proposed that wall modelling approaches could be divided into three types (which are all loosely connected):

• Equilibrium-stress models

The goal of these models is to predict the wall friction. They are based on the assumption that a constant-stress layer exists near the wall, which implies that the velocity at the first point in the outer layer obeys a logarithmic profile. Early examples include Deardoff [START_REF] Deardorff | A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers[END_REF] and Schumann [START_REF] Schumann | Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli[END_REF]. It is successful for attached flows but is not appropriate to describe flows in which the boundary layer is not at equilibrium. More recently, a generalization of the method based on the integral method of Von Karman and Polhausen was proposed by [88] . It is called iWM-LES and allows the inclusion of more physics while preserving the economy of the equilibrium models. The cost of a simulation based on iWMLES is estimated by Piomelli [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF] to increase like Re 0.6 .

• Zonal approaches

In this approach, the thin boundary layer equations are solved in a region very close to the wall and weakly coupled to a large-eddy simulation in the rest of the domain. This is a two-layer approach where the outer layer imposes its scales on the inner layer. It was used successfully by Balaras et al. [6] for a rotating channel flow. The cost of a zonal approach is supposed to be about 20% higher than a standard equilibrium-stress model.

• Hybrid RANS/LES methods

In contrast with the previous approach, the inner layer has its own time and length scales determined by a URANS simulation. The first instance of such a method is Detached Eddy Simulation (DES) [75; 79]. The cost of this approach is roughly estimated to increase like Re 0.6 log(Re) [START_REF] Piomelli | Wall-layer models for large-eddy simulations[END_REF]. An issue is the correct treatment of the interface region, as the lack of resolved eddies in the interface region will typically result in an upward displacement of the logarithmic layer in the LES region, which is known as the logarithmic law mismatch (LLM) [START_REF] Baggett | On the feasibility of merging LES with RANS for the near-wall region of attached turbulent flows[END_REF]. Such methods are most successful for flows with a source of destabilization (such as strongly separated flows) which accelerates the generation and growth of these eddies. Most of the error is present in the first cell, where the stress-carrying eddies the size of which scales with the wall distance are underresolved. Different attempts have been made to correct the LLM. Park and Moin [START_REF] Park | An improved dynamic non-equilibrium wall-model for Large Eddy Simulation[END_REF] implemented a dynamic eddy viscosity model that corrects the effect of the resolved Reynolds stress on the wall friction (see Figure 1.5).

Compressibility effects

Several efforts have been made to extend these approaches to compressible flows. Compressibility effects introduce additional physics such as viscous heating near the wall and shock/boundary layers interactions. This requires the formulation of thermal boundary conditions and modelling for the heat flux. An important question then is whether the strong Reynolds analogy holds, i.e. the local skin friction and heat transfer are proportional. It is well known that this analogy fails for separated flows, such as those submitted to strong adverse pressure gradients, as the skin friction goes to zero while the heat transfer does not.

• Equilibrium-stress models Catchirayer et al. [14] extended the iWMLES [88] method developed for incompressible flows to compressible flows. Parametrization of the density and temperature variations was carried out, and the integral compressible boundary layer energy equation was included into the formulation. The method was tested for plane channel flow over a range of Mach numbers. In both subsonic and supersonic cases the mean profiles were in good agreement with direct numerical simulation. Bermejo-Moreno et al. [START_REF] Bermejo-Moreno | Confinement effects in shock wave/turbulent boundary layer interactions through wall-modelled Large-Eddy Simulations[END_REF] used an equilibrium wall model to describe an oblique shock/boundary layer interaction inside a rectangular cross-sectioned duct.

• Zonal Approaches Extension of the TBLE (Thin Boundary Layer Equation) approach to compressible flows was considered by Kawai and Larsson [START_REF] Kawai | Dynamic non-equilibrium wall-modeling for Large Eddy Simulations at high Reynolds numbers[END_REF]. It requires resolution of the full compressible Navier-Stokes equations in the near-wall region as the conservation of mass and energy is no longer enforced automatically.

Approximate boundary conditions 1.5.1 Slip boundary condition

Instead of solving the TBL equations, Chung and Pullin [START_REF] Chung | Large-Eddy Simulation and wall modelling of turbulent channel flow[END_REF] integrated them up to some distance from the wall and an equation for the wall stress was derived using wall-parallel filters and assuming a local inner scaling (i.e. based on wall units). From the wall stress, a slip velocity at some lifted virtual wall height can be determined by assuming a logarithmic velocity profile for the streamwise velocity and a constant-stress layer. The local, Karman-like constant is determined using a subgrid scale model based on near-wall vortices stretched in the streamwise direction. This amounts to setting a Dirichlet-like condition at the interface. In the original derivation of the model [START_REF] Chung | Large-Eddy Simulation and wall modelling of turbulent channel flow[END_REF], the velocity at the lifted wall was only defined in the streamwise direction, but cross-stream components were added in later editions of the model [START_REF] Cheng | Large-Eddy Simulation of separation and reattachment of a flat plate turbulent boundary layer[END_REF].

Control-based strategies

Control theory has been used to prescribe the wall-parallel shear stresses in order to produce a pre-defined target velocity profile [77]. They used an adjoint formulation for solving the optimal controlled problems. Although these approaches are cumbersome, they shed insight into the influence of the wall stress on the LES field and can allow precise characterization of the LLM issue associated with hybrid methods. The adjoint formulation also makes it possible to evaluate the sensitivity of the wall model with respect to a priori coefficients or the location of the interface.

Synthetic wall boundary conditions

In this part we present an alternative approach where a synthetic wall boundary condition at the edge of the wall layer is derived and plugged into the simulation. The idea is to bypass entirely the computationally expensive wall region and to solve the equations in a reduced domain where the boundary is no longer the wall, but a virtual plane located at a given height within the wall layer. The approach was developed by Podvin and Fraigneau The difference with the approximate boundary conditions described in the previous sections is that the near-wall region is treated in a completely implicit manner through the POD eigenfunctions. It therefore avoids extra modelling assumptions about the near-wall region, as the condition is exclusively based on the second-order flow statistics, which are assumed to be known a priori for instance from a precursor simulation.

The challenge is to provide adequate temporal amplitudes for each POD spatial structure so that the reconstructed field mimics the behavior of the turbulent flow within the layer.

Different variants of the method were implemented for an incompressible turbulent plane channel flow over a range of Reynolds numbers [66; 67]. It was found in particular that identifying the phase of the structures, i.e the locations of the boundary planes in the wall region were essential for the success of the procedure. Results shown in Figure 1.8 indicate that the turbulent statistics in the reduced domain are in good agreement with those of the reference channel.

We note that a simpler rescaling procedure was used by Mizuno and Jimenez [START_REF] Mizuno | Wall turbulence without walls[END_REF] and applied to the velocity field in Fourier space. The field on the boundary was extracted from a field on a higher plane and each of its components in Fourier space was rescaled with a factor that depends on the horizontal wavenumbers. As a result, the wall shear stress was incorrectly estimated and the virtual origin of the flow field was displaced, unlike the POD-based reconstruction procedure.

Outline of the thesis

The goal of this thesis is to extend the approach presented above, which was developed for incompressible turbulence, to compressible boundary layers. In the incompressible formulation, the equations are elliptic and the boundary condition is prescribed on the entire boundary. In the compressible formulation, the equations are hyperbolic and the flow on the boundaries is determined from the Riemann invariants. Depending on the characteristic velocities, the boundary conditions are either propagated from within the domain or prescribed so as to mimic conveyance of information from outside the domain.

The thesis is organized as follows: in chapter 2, we present the numerical method. The reconstruction tools and the reconstruction method are described in chapter 3. The next chapters present results for turbulent channel flow. In chapter 4, the flow is subsonic and the synthetic boundary condition is applied on one wall only. In chapter 5, the flow is subsonic and both walls are replaced by synthetic boundary conditions. Results for the supersonic case are reported in chapter 6. We then present a conclusion and perspectives in chapter 7. This study deals with high speed turbulent flows in the continuous regime. The conventional governing equations are the compressible version of the Navier-Stokes equations which formalize the conservation of mass, of the momentum and the total energy. In the following, ρ is density, u is velocity vector, P is the pressure, E is total energy per unit of mass. All variables are dimensionless by the reference quantities (reference density ρ 0 , reference velocity U 0 , reference length H, and reference temperature T 0 )

The dimensionless equations of, respectively, the mass conservation, the conservation of the momentum and the total energy conservation read:

∂ρ ∂t + ∂ ∂x j (ρu j ) = 0; ∂ ∂t (ρu i ) + ∂ ∂x j (ρu i u j ) + ∂P ∂x i = ∂τ ij ∂x j ; (2.1) ∂ρE ∂t + ∂ ∂x j (ρE + P )u j = ∂ ∂x j τ ij u i - ∂ ∂x j q j .
In these equations, the Einstein summation rule is considered regarding the repeated indices.

All along the manuscript, the fluid is considered as an ideal gas. The pressure is then related to the internal energy by the following relationship that is here written by using the conservative variables [ρ, ρu, ρE]:

P = (γ -1) ρE - 1 2 ρu • ρu ρ . (2.
2)

The dimensionless viscous stress tensor (τ ) is expressed as:

τ ij = µ(T ) Re ( ∂u i ∂x j + ∂u j ∂x i - 2 3 δ ij ∂u k ∂x k ), (2.3)
and the dimensionless heat flux (q) is:

q j = - µ(T ) (γ -1)ReP rM 0 2 ∂T ∂x j . (2.4)
Following the equation of state written for dimensionless quantities, the temperature is related to the pressure and the density with:

T = γM 0 2 P ρ . (2.5)
The molecular dynamic viscosity (µ) is calculated once the temperature is known according to Sutherland's law following:

µ(T ) = T 3 2 1 + s T + s , (2.6)
with s = 110.4 T 0 .

Following this version of governing equations, the flow configurations regarded are defined once four constant values are prescribed:

• The specific heat ratio (γ): γ = C p C v
, where C v and C p are heat capacities, respectively at constant volume and constant pressure.

• The Reynolds number expressed as Re = ρ 0 U 0 H µ(T 0 ) , which is based on the reference values of the density (ρ 0 ), the reference velocity (U 0 ), the reference length scale H, and the viscosity (µ(T 0 )), based on the reference temperature (T 0 ). The Reynolds number compares the inertia forces with the viscous forces.

• The Mach number:

M 0 = U 0 C
, where C is the speed of sound related to the reference values by C 2 = γ P 0 ρ 0 .

• The Prandtl number: P r = Cpµ kt that depends on the fluid and has a specific value P r = 0.72 for air.

The numerical approach for solving the governing equations

Simulations conducted in this study were performed by using the in-house parallel (MPI) DNS solver, named CHORUS (Compressible High-Order Unsteady Simulation), that has been developed at LIMSI for unsteady compressible flow simulations. The ability of the CHORUS software to compute high Reynolds compressible flows has been demonstrated on various test-cases in previous studies (see for instance [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF]).

The resolution of the governing equations (2.1) is based on a finite volume approach. An operator splitting procedure is employed that splits the resolution into the Euler part and the viscous problem:

∂Q c ∂t + ∇ • F Euler = ∇ • F visc . (2.7)
Where Q c is conservative variable, the left side of (2.7) is the same as the left side in (2.1). The right side is divergence of viscous flux, which is the same as the right side in (2.1).

Here, the Euler part:

Q c (n+ * ) = Q c n -δt ∇ • F Euler (Q c n ), (2.8)
is discretized by means of a one-step high-order coupled time and space scheme, named OS7 scheme [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF] based on a Lax-Wendroff approach, which ensures a 7th-order accuracy in both time and space.

Besides, the space discretization of the viscous problem is obtained by means of a classical 2nd-order centered scheme for the viscous flux approximation, coupled to a 2nd-order Runge-Kutta time integration giving at last a second order accurate scheme in both time and space:

Q c (n+ * * ) = Q c (n+ * ) + δt ∇ • F visc (Q c (n+ * ) , ∇Q c (n+ * ) ), (2.9)
Q c (n+1) = 1 2 Q c (n+ * ) + Q c (n+ * * ) + δt 2 ∇ • F visc (Q c (n+ * * ) , ∇Q c (n+ * * ) ). (2.10)
The viscous terms are discretized with centered second order formulae that are expressed in Appendix A.

The multidimensional resolution is ensured by a directional splitting following the Strang procedure. All information relative to the numerical procedure can be found in [START_REF] Daru | High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations[END_REF].

Numerical configuration : the compressible turbulent channel flow

We consider a compressible turbulent channel flow configuration where the air flows in between two horizontal solid walls maintained at the same constant temperature that limit the channel at the lower and upper bounds. The configuration is shown in Figure 2.1. The streamwise, spanwise, and wall-normal directions of the flow are respectively denoted by x, y, and z. The dimensions of the channel are (

L x × L y × L z ) = (2π × 4 3 π × 2)H
, where H is the half of the space between the two horizontal solid walls and is taken as the reference length scale (H = 1). The corresponding velocity components in the space directions are respectively denoted u, v, and w. The mesh in the streamwise and spanwise direction has a constant grid spacing. On the opposite, in the wall-normal direction, the grid is tightened near the solid wall by using a hyperbolic tangent function to ensure that the first point above the wall satisfies the constraint z + < 1 (expressed in the wall units).

Boundary and initial conditions:

In the streamwise (x) and the spanwise (y) directions, we suppose periodic developments and periodic boundary conditions are prescribed. In wall normal direction z, a no slip boundary condition is prescribed on the solid walls and the temperature is prescribed. The density is then calculated by solving the continuity equation to ensure the conservation of mass.

At the initial state, the streamwise velocity profile in the normal to the wall direction is defined as:

u t=0 = 3 U 0 z(2H -z) 4H 2 , (2.11)
and the other components v, and w at set to zero. The initial temperature is related to the velocity profile following:

T = T wall - 3γP rU 2 0 ((1 - z H ) 4 -1) 4 , (2.12)
where T wall is the dimensionless temperature at the solid wall, considered as the reference temperature. The initial pressure is set constant and the density is calculated using the equation of state (2.5):

ρ = γM 2 0 P T . (2.13)
The reference time is defined as H/U 0 , and in the following, results will be presented in time units.

Treatment of the periodic boundary condition in the streamwise direction.

In our model (Figure 2.1), to simplify the simulation, the streamwise direction is assumed to be periodic. In fact, one must notice that the pressure does not evolve periodically, however the pressure gradient which appears in the Navier-Stokes equations (2.1) can be considered as periodic. To enforce a periodic flow motion and compensate looses of the mass, the momentum and the total energy due to viscous effect in the boundary layers, a macroscopic pressure gradient is added to the channel flow that could be viewed as an external force applied to the flow. The magnitude of this macroscopic pressure gradient must be calculated at each time step. To do so, we average the equations of the momentum components over an horizontal plane, in the x and y directions. For this, we define an operator, noted < • > xy applied to a variable Φ:

< Φ > xy = 1 L x L y Lx 0 Ly 0 Φ dx dy.
(2.14)

We then integrate the resulting momentum equation in the wall normal direction (z).

Assuming an upper-lower symmetry on the mean flow, we finally obtain the equation that drives the mean flow rate (Q m ):

L x ∂Q m ∂t = -L x L y L z f 1 - 2L x L y Re µ ∂ < u > xy ∂z wall , (2.15)
where L x , L y and L z are, respectively, the length, the width and the height of the channel. The extension to compressible flows of the algorithm proposed by Deschamps

[25] is here employed to update the driving force at each time step:

f n+1 1 = f n 1 + δt L y L z [α(Q n+1 -Q 0 ) + β(Q n -Q 0 )], (2.16)
with Q 0 is the target mass flow rate, Q n is the mass flow rate at time step n δt, and Q n+1 the one at time step (n + 1) δt. The forecasted mass flow rate (Q n+1 ) is then evaluated following:

Q n+1 = Q n -δt L y L z f n 1 + 2δt Ly Re µ ∂ < u n > xy ∂z | wall .
(2.17)

From numerical experiments [START_REF] Lenormand | Contribution à la simulation des grandes Échelles d' Écoulements turbulents compressibles en milieux confinés[END_REF], the define values of the constants in incompressible to the CHORUS Software and has already been validated on channel flow test-cases.

Results of the Direct Numerical Simulation

Statistical treatments of simulation data

We would like to describe here the statistical treatment we applied to extract information from the simulation data.

In the following, < • > represents the mean value in time and over the plane z = constant, expressed as:

< Φ >= 1 t end -t start t end tstart < Φ > xy dt,
where t start and t end are the dimensionless times corresponding to the beginning and the end of the mean procedure, and < • > xy is the average operator given by (2.14).

As the near wall region is dominated by the viscous shear stress, an adequate dimensioning of the quantities in this region is achieved by using the friction velocity u τ :

u τ = τ w ρ w ,
where ρ w is the density at the wall. τ w is the shear stress tensor evaluated at the wall:

τ w = µ w d < u > dz wall , (2.18)
with µ w the viscosity at the wall. Quantities are mostly presented by using the wall units based on u τ as the reference velocity, and the viscous length scale δ ν = µ w ρ w u τ as the reference length scale.

For instance, lengths are expressed in wall units as:

z + = Re H ρ w u τ z µ w , where Re H = ρ 0 U 0 H µ(T 0 )
. Velocity is expressed as u + = < u > u τ . One of the control parameter used in the numerical simulations of bounded flows is the Reynolds number based on the friction velocity:

Re τ = Re H ρ w u τ H µ w = Re H H δ ν . (2.19)
The development of the turbulent boundary layer on a horizontal wall are often judged on the quality of the normal to the wall profile of the mean streamwise velocity that must reach a classic profile. This profile usually presents several regions:

• The viscous sublayer that is situated very close to the wall. The wall damps the normal to the wall fluctuations and the flow is mainly dominated by the wall shear stress. As so, the wall shear stress (2.18) can be considered as constant in this layer, and consequently the velocity recovers a linear variation towards the direction normal to the wall:

u + = z +
• The outer layer is mainly scaled by the macroscopic quantities, let say the external velocity U e and the boundary layer thickness δ CL . Consequently, the streamwise velocity profile follows a velocity-defect law where < u > -U e u τ is a function of z δ CL .

• In between, an overlap layer exists, named the inertia region where an equilibrium is reached between the production of the turbulent kinetic energy and its dissipation.

This layer is characterized by a logarithmic velocity profile:

u + = 1 κ log z + + G, (2.20)
where κ is the Von-Karman constant, κ = 0.41. G is a constant determined through experiments; the generally admitted value is G = 5. In compressible flows, a more complex law is preferred to the classic log-law, that takes into account the variation of the Mach number and the temperature:

1 A arcsin 2A 2 U U e -B (B 2 + 4A 2 ) 1 2 + 1 A arcsin B (B 2 + 4A 2 ) 1/2 = u τ U e 1 κ log(z + ) + G , (2.21) with the A 2 = γ -1 2 M 2 0 T 0 T c
, and

B 2 = 1 + γ-1 2 M 2 0 T 0 T c -1.
T c is the temperature at the center of channel.

• A buffer layer also exists in between the viscous sublayer and the inertia region where the mean streamwise velocity profile switches from the linear profile to the logarithmic profile.

The generally admitted limits of these various regions are:

• the viscous sublayer: z + ≤ 5;

• the buffer layer: 5 < z + ≤ 30;

• the logarithmic region: z + > 30 Finally, if we assume that the velocity fluctuations (u ) can be evaluated as:

u = u-< u >,
the components of the Reynolds shear stress( < u ⊗ u >) are calculated as:

< u i u j >=< u i u j > -< u i >< u j > . (2.22)

DNS results of the subsonic channel flow

We consider a periodic plane channel flow bounded by two isothermal solid walls separated by a gap of L z = 2 H. The Mach number is M a = 0.5. The Reynolds number based on the bulk velocity and the channel half-height is Re H = 3000, corresponding to a friction Reynolds number equal to Re τ = 180. Present simulations have been conducted by using the unlimited OS7 scheme, described previously. This test-case have largely been studied through both DNS and LES. Although incompressible, the DNS reference results from J. Kim et al. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] are generally taken as a reference solution. The present results will be compared with these reference results.

The mesh of this test-case is 97 × 97 × 129 in the x × y × z directions, respectively. As the initial state corresponds to the solution for a laminar flow, a transition towards turbulence occurs during the simulation after a long time integration. After this transition, statistics are calculated over a very long time.

The mean streamwise velocity profile (< u >), non-dimensionalized by the friction velocity, is plotted versus the wall normal direction expressed in wall units (z + ) in Figure 2.3.

The convergence of the statistics of mean values in the present DNS can be judged on the perfect fit obtained between the profiles from the upper and the lower half of the channel.

A logarithmic law is clearly obtained in the present results which fits the classical log-law distribution (Figure 2.3). The general shape of the profiles is in good agreement with the reference that has been shown in Figure 2.5. The maximum values in the rms profiles and their locations seem to be rather well predicted. As far as the channel flow is fully developed and reaches the equilibrium state, the total shear stress (i.e.

-< u w > + 1 Re H ∂ < u > ∂z
) recovers a straight line. Dimensionless by the wall units and the friction velocity, the slope of this line is -1, [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]. The present computation recovers the right linear distribution of the rms rms of u rms of v rms of w Max value in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] 2.7 1.1 0.83 Position of max value in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] 0.1 0.18 0.26 Value on center in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] 0 Although the general shape of the profiles is in good agreement with the reference, there however exists some discrepancies compared with the reference values. The maximum values in the rms velocity profiles and their locations are reported in Table 2.1 and compared with the reference values from [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]. We can see that the maximum values are rather well predicted although some discrepancies are evidenced mainly attributed to the convergence of the statistics that is certainly not fully reached for the rms values.

The locations of the maximum seem to be rather well predicted.

In the following, theses results are considered as statistically converged enough to apply the Proper Orthogonal Decomposition.

DNS results for the supersonic channel flow.

We here consider the same geometrical configuration with however a higher Mach number (M a = 1.5). The Reynolds number based on the bulk velocity and the channel halfheight is Re H = 3000, corresponding to a slightly higher friction Reynolds number equal to Re τ = 220. This test-case has largely been studied through both DNS and LES, and the reference solution generally considered is from G. N. Coleman et al. [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF]. The present results will be compared with these reference results.

The mesh of the supersonic test-case is 257 × 161 × 161 in the x × y × z directions, respectively. The mesh dimenion is supersonic flow is larger than that in subsonic flow to well capture the sturucture of turbulence. In the supersonic case, the transition towards turbulence occurs after a very long integration time, much longer than in the supersonic test case. Once the turbulent flow is fully developed, statistics are calculated over an integration time equal to t endt start = 75 time units.

Mean profiles of the density, the static temperature and the streamwise velocity component are plotted versus the normal to the wall direction in Figure 2.6. The streamwise velocity component is dimensionless by the friction velocity and z is expressed in wall units, and the profiles of u in both bottom half channel and top half channel are the same.

A logarithmic law is clearly obtained in the present results which fits the classical log-law distribution (Figure 2.6, bottom); its extent is however weaker than in subsonic case.

Wall normal profiles of the temperature and the density (Figure 2.6) recover classical distributions. The mean density at the wall reaches 1.355 while its value at the channel center is 0.98. These values are in complete agreements with the reference values given in Table 2.2. The mean temperature is unity at the wall since it is the value taken as reference for the dimensionless variables, and reaches 1.37 at the channel center, which is weakly underestimated compared with the reference value (Table 2.2).

Table 2.2: Time averaged results from [START_REF] Coleman | A numerical study of turbulent supersonic isothermal-wall channel flow[END_REF]. † corresponds to the incompressible DNS results from [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]. [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF] are also reported to compare with.

Present DNS profiles compare rather well with the DNS results of [START_REF] Lechner | Turbulent supersonic channel flow[END_REF].

rms rms of u rms of v rms of w Max value from [START_REF] Lechner | Turbulent supersonic channel flow[END_REF] 2.8 0.95 0.7 Location of max value from [START_REF] Lechner | Turbulent supersonic channel flow[END_REF] 0.12 0.2 0.34 Value on the center from [START_REF] Lechner | Turbulent supersonic channel flow[END_REF] 0.9 0.72 0.72 Max value in present DNS To better compare with, we reported, in Table 2.3, the maximal values of the rms velocity profiles and their locations in the normal to the wall direction z. The values at the center of the channel are also reported in this table. Although some slight discrepancies are noticeable that might be attributed to a weak statistical convergence of the rms values, we can think that the present DNS give results that are validated and can favorably be used for model decomposition.

Comparison between results of subsonic flow and supersonic flow

To compare results of the subsonic channel flow with those of the supersonic channel flow, mean profiles of the density, the viscosity and the streamwise velocity are presented in Figure 2.8. The present results can be compared with results of [START_REF] Foysi | Compressibility effects and turbulence scalings in supersonic channel flow[END_REF] in To present the method, let us consider a fluctuating spatio-temporal field Q(x, t)

< Q(x, t) >= 0,
where < • > is a temporal average. The quantity Q can be a scalar or a vector, and can aggregate different physical quantities. Proper Orthogonal Decomposition of the field Q consists in writing the field as a denumerable superposition of spatial modes, the amplitude of which varies in time

Q(x, t) = ∞ n=1 a n (t)φ n (x). (3.1)
The modes φ are obtained as the eigenfunctions related to the following maximization problem

Ψ = M ax < (Q(x, t), φ(x)) > (φ(x), φ(x)) , (3.2)
where (., .) defines an inner product over the domain. In all that follows we will use the standard inner product method.

Direct Method

The solution of the maximization problem leads to the following eigenvalue problem (Fredholm integral of the first kind [40]):

(< Q(x, t), Q(x , t >, φ(x )) = λφ(x). (3.3)
We can see that < Q(x, t)Q(x , t) > is the autocorrelation tensor of Q at zero time lag.

The operator is symmetric positive definite, so that the eigenvalues are all positive and the eigenfunctions are orthogonal. Moreover, due to compactness of the operator, the eigenvalues are denumerable and can be ordered

λ 1 ≥ λ 2 ≥ . . . ≥ λ n .
The eigenvalues λ n represent the time-averaged contribution of the mode n to the flow energy -as defined by the inner product (., .).

By convention the eigenfunctions are normalized and the amplitudes a n (t) which can be expressed as In the case of homogeneous flow in one direction x, the spatial autocorrelation tensor depends only on the separation

a n (t) = (Q(x, t), φ n (x)), (3.4)
< Q(x, t), Q(x , t) >= f (x -x ), (3.7)
and in that case it is straightforward to show that the spatial POD modes are in that case Fourier mode:

φ n (x) = e inkx , (3.8)
Where k is wavenumber. This means that the problem is decoupled in the homogeneous direction and that the autocorrelation tensor can be considered directly in Fourier space.

If we consider a vector Q of p dimensions which depends on q spatial dimensions, the dimension of the eigenproblem is pq. The autocorrelation tensor C is typically computed from N samples:

C =< Q(x, t)Q(x , t) >≈ 1 N N n=1 Q(x, t n )Q(x , t n ). (3.9)
If we have a fully inhomogeneous problem, for a Cartesian grid with resolution (N x , N y , N z ), the number of spatial dimensions is q = N x N y N z . If the problem is homogeneous in the x direction, one can decouple the problem in Fourier space and it is only necessary to solve N x /2 problems of size pq where q = N y N z . In the case of the turbulent channel flow, the flow is homogeneous in two directions, and one needs to solve N x N y /2 problems of pN z dimensions.

We note if we use the standard inner product, the eigenproblem will be discretized as

< Q l (x i , t), Q m (x j , t > φ m (x j )w j = λφ l (x i ), (3.10)
where w j > 0 is the weight associated with the spatial discretization. In the case of an irregular mesh, this makes the problem non-symmetric. However, multiplying on each side by w 1/2 i makes the problem symmetric again since we have

w 1/2 i < Q l (x i , t)Q m (x j , t > φ m (x j )w 1/2 j w 1/2 j λφ l (x i )w 1/2 i . (3.11) If we let W = diag([w 1/2 1 , . . . , w 1/2
n ]), we can then solve for the eigenfunctions ψ n = W φ of the operator W CW .

Method of snapshots

In the case where the flow is not homogeneous, application of the direct method is very costly as the dimension of the problem in space is very high. Since the spatial autocorrelation tensor is computed from N snapshots, one can see that the limiting factor of the problem is the minimum of N and pq. In many cases N << pq, so that the direct method is unnecessarily cumbersome. The idea is therefore to look for the N modes φ as a combination of the N snapshots Q(x, t n ) [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part i: Coherent structures[END_REF]:

φ n (x) = A np Q(x, t p ).
(3.12)

Injecting this in (3.3), replacing the expression for the autocorrelation tensor with its sampled mean, and using the assumption that the N snapshots are linearly independent, leads to the new eigenproblem

1 N (Q(x, t n ), Q(x, t p ))A pq = λA nq . (3.13)
One can show that

A np = a n (t p ). (3.14)
This is equivalent to solving the problem in the space of snapshots instead of physical space. In this thesis, since we consider a plane turbulent channel flow, both approaches can be and have been implemented. We use the following definitions:

• For the method of snapshots, the standard inner product is defined as (Q(x, y, z, t), φ(x, y, z)) =

x y z Q(x, y, z, t)φ(x, y, z)dxdydz.

(3.15)

• For the direct method, we will use the Fourier transform of the field in the horizontal directions (the channel is periodic with length L x and L z in respectively the streamwise and the spanwise direction) so that The inner product is therefore defined for each wavenumber:

Q(x,
(Q kxky (z, t), φ kxky (z)) = z Q kxky (z, t)φ * kxky (z)dz, (3.17)
where * refers to the complex conjugate.

To alleviate notations and avoid confusion with the horizontal wavenumbers, we will indicate the index of the POD mode as a superscript in Fourier space and write

Q kxky (z, t) = ∞ n=1
a n kxky (t)φ n kxky (z).

(3.18)

Symmetry

The comparison of the two methods can be done in the context of symmetry. The turbulent channel flow satisfies on average the following symmetries (see Figure 3.1):

• translation symmetry in the streamwise direction x

• translation symmetry in the spanwise direction y

• reflection symmetry with respect to a vertical plane y = cst

• reflection symmetry with respect to the horizontal mid-plane z = 1 No symmetry is enforced through application of the method of snapshots to an arbitrary set of realizations. In contrast, applying the direct method in Fourier space means enforcing the first two symmetries. To enforce the following two symmetries, one can add snapshots which are the image of the original snapshots through the symmetry, thus enlarging the database.

In practice, we enforced the fourth symmetry in both approaches by adding to each snapshot of the database its image through the mid-plane. The domain of definition of POD was chosen to be the full boundary layer that is one half of the channel.

Convergence

A key issue is the convergence of POD, i.e its ability to approximate a field with a finite number of modes. Mathematically, the set of eigenfunctions associated with the continuous spatial autocorrelation tensor is infinite. Numerically, the set of eigenfunctions is finite and is typically determined by a finite set of fields. Its size corresponds to the numerical spatial discretization for the direct method and the number of samples for the method of snapshots. However it is important to realize that there is no guarantee that a field that is not part of the samples used to compute POD may lie entirely in the finite-dimensional space spanned by the POD basis as shown in Figure 3.2. The difference may not be small and will depend on the convergence of POD, i.e its ability to capture a significant fraction of energy with a limited number of modes. Since POD is energy optimal, we expect the energy of a field to be contained in a limited number of structures. However the convergence properties of POD will vary depending on the flow. As an example, 99.9 % of the controlled wake flow behind a cylinder at a Reynolds number of 200 [START_REF] Bergmann | Optimal control of the cylinder wake in the laminar regime by trust-region methods and POD reduced-order models[END_REF] can be represented with 14 POD modes. The convergence of POD is much slower in turbulent channel flow [START_REF] Keefe | The dimension of attractors underlying periodic turbulent poiseuille flow[END_REF].

One final issue is the inclusion or exclusion of the mean flow. What happens if we include the mean flow in the decomposition? In most cases the mean flow is orthogonal to the fluctuations in the sense of the defined inner product, and therefore constitutes an extra mode, which is typically the most energetic. However the temporal average of the first coefficient is no longer zero.

Results

POD can be applied to any combination of variables. We chose to apply it independently to the density ρ, the momentum ρu, ρv, ρw, and the energy ρe. We present spatial POD modes for both methods. We used 100 samples, separated by ∆t = 6.25H/U 0 .

We indicate between [ ] the variables over which POD is performed jointly. We show We can see that the dominant eigenfunction (φ ρu )(φ ρw ) < 0, which corresponds to the fact that the Reynolds stress < ρu w > < 0.

Linear Stochastic Estimation

General definition

A complementary technique to POD is constituted by Linear Stochastic Estimation.

Linear Stochastic Estimation of coherent structures was introduced by Adrian and Moin [START_REF] Adrian | Stochastic estimation of organized turbulent structure: homogeneous shear flow[END_REF] and reviewed in [START_REF] Adrian | Stochastic estimation of conditional structure: a review[END_REF].

It consists in estimating the flow structure given a data event E. The best mean square estimate of Q given the data E is the conditional average of Q given E, (Q|E). In Adrian's approach, Q is constituted by the velocity field u. If the elements of Q and E are joint normally distributed, it is well known that (u|E) is a linear function of E [START_REF] Adrian | Stochastic estimation of conditional structure: a review[END_REF]. If this is not the case (Q|E) can only be approximated as a linear function of the event data E (the approximation can be extended to include quadratic terms). To determine this function, we use observations of Q given the data E. Let us call Q the linear estimate of (Q|E). We have

Q i = M j=1 L ij E j , (3.19)
where M is the number of data events for which observations exist. The estimation coefficients L ij are chosen to minimize

M in|(Q|E) - M j=1 L ij E j | 2
This leads to the following system of M × M linear equations:

M j=1 < E j E k > L ij =< E k Q i > for all j, k = 1, . . . M, (3.20)
where < . > represents an ensemble average. The key is to identify the structure or eddy based on the event E.

Application

We now explain how this method can be applied to our purpose, which is to estimate the boundary condition Q 0 on a plane h 0 given the field Q in the upper layer h 0 < z ≤ 1

shown in Figure 3.5. We will consider the method of snapshots (physical space). The derivation in Fourier space is strictly similar. The idea is to solve the problem in the POD space, that is to determine the projection of Q 0 on a set of POD modes φ n , (n = 1, . . . , N ) corresponding to the full reference domain 0 ≤ z ≤ 1. We therefore look for Q 0 such that

Q 0 (x, y, z 0 , t) = N n=1 a n (t)φ n (x, y, z 0 ). (3.21)
There will therefore be many variables Q 0 as there are independent compositions. For instance, if POD is applied to all variables, there will be one field Q. If it is applied separately to the density, momentum, and internal energy, there will be three fields Q:

Q ρ , Q ρu and Q ρe . In the domain, the field Q is characterized by

Q(x, y, z, t) = N n=1 a n (t)φ n (x, y, z), (3.22)
where

a n (t) = Lx x=0 Ly y=0 1 z=0 Q(x, y, z, t)φ n (x, y, z)dxdydz. (3.23)
The question we wish to solve is the following: at an instant t, can we determine a n (t) from the knowledge of the flow field in the domain h 0 ≤ z ≤ 1? The procedure is similar to what has been in described in [START_REF] Podvin | POD-based wall boundary conditions for the numerical simulation of turbulent channel flows[END_REF].

If we let b n (t) = Lx x=0 Ly y=0 1 z=h 0 Q(x, y, z, t) • φ n (x, y, z)dxdydz, (3.24)
the best estimate for a n , (n = 1, . . . , N ) given b n , (n = 1, . . . , N ) is given by the system

b n (t) = Lx x=0 Ly y=0 1 z=h 0 N p=1 a p (t)φ p (x, y, z) • φ n (x, y, z)dxdydz, (3.25)
which can be rewritten as:

b n (t) = R np a p (t), (3.26)
where

R np = Lx x=0 Ly y=0 1 z=h 0 φ p (x, y, z) • φ n (x, y, z)dxdydz. (3.27)
Similarly, in Fourier space, one will solve [START_REF] Podvin | POD-based wall boundary conditions for the numerical simulation of turbulent channel flows[END_REF] b n kxky (t) = R np kxky a p kxky (t),

where

R np kxky = 1 z=h 0 φ p kxky (z) • φ n * kxky (z)dz. ( 3

.29)

If h 0 = 0, R = I (full domain). The matrix R represents the interaction of the structures on the restricted domain h 0 < z < 1. It is a linear problem of size N . However it is ill-conditioned due to the projection issues mentioned earlier. Solving the full matrix problem will typically result in large amplitudes, especially for high-order modes. One way to avoid this problem is to use regularization such as Tikhonov-based regularization [START_REF] Cordier | Calibration of pod reduced-order models using tikhonov regularization[END_REF]. We chose instead to directly filter the matrix based on the physical meaning of the matrix coefficients. If the value of h 0 is small, R should remain sufficiently close to the identity. This is true for the direct POD method in Fourier space where R is defined for each horizontal wavenumber pair.

Since R np =

Lx x=0

Ly y=0

1 z=h 0 φ p (x, y, z) • φ n (x, y, z)dxdydz, we expect R = R + ε. R np =            R 1,1 0 ... 0 0 R 2,1 R 2,2 ... 0 0 ... ... ... ... ... R n-1,1 0 ... R n-1,n-1 0 R n,1 0 ... 0 R n,n            +            0 ε 1,2 ... ε 1,n-1 ε 1,n 0 0 ... ε 2,n-1 ε 2,n ... ... ... ... ... 0 ε n-1,2 ... 0 ε n-1,n 0 ε n,2 ... ε n,n-1 0           
R can be approximated with its diagonal part. However in the case of the method of snapshots, the first mode corresponding to the mean flow has so much energy that its interaction with all other modes cannot be neglected. We therefore keep the first column and the diagonal of the matrix R, and we express these terms R. For the remaining term ε, we assume it is very small, so R ∼ R.

We noticed that for some higher-order modes, the diagonal term is no longer large with respect to the other terms in Figure 3.6. Since this is likely to bias the estimation strongly, we excluded the amplitudes of modes n from the reconstruction for which R nn was smaller than a threshold ts = 0.3 [START_REF] Podvin | POD-based wall boundary conditions for the numerical simulation of turbulent channel flows[END_REF]. The threshold level was chosen empirically.

Variations in the level did not affect the reconstruction significantly.

Results

With the estimation method, we could obtain estimated temporal coefficients a est n (t), the different fields can be reconstructed from the estimated amplitudes

Q P (x, y, z, t) = N n=1 a est n (t)φ n (x, y, z). (3.30)
The reconstruction error can be measured by the integral We continued the simulation and chose a field that does not belong to the set of snapshots used to compute POD (the time separation with last snapshot was 6 time units). The errors between the field and its projection are represented in Table 3.1. From Table 3.1, we can see that the discrepancies of density and energy on boundary plane are very small. The discrepancies of streamwise velocity u on boundary plane are limited, but could not be neglected and the discrepancies of spanwise velocity v and the wall-normal velocity w are not small. This means that the cross-stream components are not very well recovered by their projection. This may make the use of rescaling necessary.

e Q (t) = |Q P (x, y, h 0 , t) -Q(x, y, h 0 , t)| 2 dxdy |Q(x, y, h 0 , t)| 2 dxdy . ( 3

Reconstruction method

In this section we present the reconstruction procedure used to derive the synthetic wall boundary condition. We first give an overview of synthetic boundary conditions, focusing first on inlet boundary conditions. We then detail the reconstruction algorithm.

Inlet Synthetic boundary conditions: rescaling approaches

Determining appropriate boundary conditions is a key issue for the simulation of flows, as on the one hand it must describe the physics and on the other hand it must be defined in a suitable way for numerical resolution. In this section, we ask the question: how is it possible to define a boundary condition that mimics the physics of the flow?

The simplest answer would be to superpose random fluctuations on a desired mean profile, where Q =< Q > +Q R . This method has been implemented with some success. Lee et al [START_REF] Lee | Simulation of spatially evolving turbulence and the applicability of Taylors hypothesis in compressible flow[END_REF] used this method for direct numerical simulation of compressible turbulence.

Le and Moin [START_REF] Le | Direct Numerical Simulation of turbulent flow over a backward-facing step[END_REF] generated anisotropic turbulence for the generation of inlet boundary conditions for a backward facing step.

The amplitudes of the fluctuations can be set to satisfy some statistics such as the Reynolds stresses but typically these statistics are second-order and computed at one point only. It is very difficult to impose phase relationships between the fluctuations and higher-order correlations cannot be satisfied and the flow lacks turbulent structure.

As a result, the success of these methods has been limited. In particular such characteristics of the turbulent boundary layer as the momentum thickness or wall friction are not always recovered. This led Aksevoll and Moin [START_REF] Akselvoll | Large-Eddy Simulation of turbulent confined coannular jets[END_REF] to introduce a auxiliary or precursor simulation, from which they could select a location where these characteristics were close to the desired values. A simple auxiliary simulation is that of a periodic channel flow in both horizontal directions. Such a simulation was used by Kaltenbach [START_REF] Kaltenbach | Study of flow in a planar asymmetric diffuser using Large-Eddy Simulation[END_REF] to generate inlet flow conditions for LES of a plane diffuser. This idea can be extended to a turbulent boundary layer. However the boundary condition at the top of the domain has to be modified to reproduce that of a free-stream flow, which is done with a symmetry condition at the top of the channel (which ensures that the vertical velocity is zero). However, unlike the boundary layer, a velocity field in a channel is characterized by a non-zero mean advection. To address this issue, Spalart and Leonard [START_REF] Balaras | Two-layer approximate boundary conditions for Large-Eddy Simulations[END_REF] added source terms to the Navier-Stokes equations which led to an equilibrium spatially evolving boundary layer with a correct momentum thickness and wall friction. A considerable simplification of the method was provided by Lund et al [START_REF] Lund | Generation of turbulent inflow data for spatially-developing boundary layer simulations[END_REF] who only modified the boundary conditions of the equations. The idea is to use as an inlet boundary condition the velocity field extracted at the outlet and properly rescaled using Spalart and Leonard's ideas.

The mean flow is rescaled according to the law of the wall in the inner region and the defect law in the outer region. The velocity fluctuations are rescaled with the local turbulent intensity. In order to avoid over-determination of the problem, the momentum thickness is directly imposed at the inlet from the computation of the wall friction using a correlation similar to the Ludwig-Tillmann correlation [START_REF] Lund | Generation of turbulent inflow data for spatially-developing boundary layer simulations[END_REF]. This recycling method has been popular over the years [START_REF] Le | Direct Numerical Simulation of turbulent flow over a backward-facing step[END_REF] and several variants have been proposed [START_REF] Xiao | An LES turbulent inflow generator using a recycling and rescaling method[END_REF].

Inlet synthetic boundary conditions: Structure-based decompositions

A way to limit the lack of correlation between the fluctuations is to use structure-based decompositions, where the velocity field is represented as a superposition of coherent structures.

The synthetic eddy method (SEM)

In a series of papers [8; 31; 87], a method called the synthetic eddy method was developed.

The method is able to reproduce first and second order one point statistics as well as characteristic length and time scales. The shape of coherent structures is also correctly recovered.

We now briefly describe the principle of the method. Let us consider a shape function f σ defined over a one-dimensional support [-δ/2, δ/2] which satisfy the normalization

condition 1 δ δ/2 -δ/2 f 2 σ (x) dx = 1. (3.32)
This shape function is supposed to represent an eddy or a turbulent spot defined by its position x i and its length scale σ i .

Examples used include a tent function f j (r) = 1r/L where r = y 2 + z 2 + (U 0 t) 2

(see Figure 3.9). Figure 3.9: Tent like function defined by f j (r) = 1r/L , where y = 0, z = 0, L = 2, U 0 t = 0.

The contribution of the turbulent spot i to the velocity field is given by

u(x) = i f σ (x -x i ), (3.33)
where i is a random value which is either +1 or -1 and x i is a random position on the interval. This can be extended to three dimensions.

The method has been applied to turbulent isotropic turbulence [START_REF] Jarrin | A synthetic-eddymethod for generating inflow conditions for Large-Eddy Simulations[END_REF]. Benhamadouche et al [START_REF] Benhamadouche | Synthetic turbulent inflow conditions based on a vortex method for Large-Eddy Simulation[END_REF] applied the method to derive an inlet boundary condition for DNS and LES of turbulent channel flow. They used a 2D vortex method for the cross-stream components of the velocity and a Langevin equation for the streamwise velocity component to mimic a Reynolds stress transport model. A good agreement between the reduced model and DNS was observed. The method can be useful to generate synthetic turbulence on the interface of a hybrid RANS/LES flow solver. A comparison between the SEM method and recycling approaches has been carried out by [START_REF] Keating | A priori and a posteriori tests of inflow conditions for large-eddy simulation[END_REF]. Recently, the SEM method has been recently extended in order to include thermal fluctuations [START_REF] Oh | Extended synthetic eddy method to generate inflow data for turbulent thermal boundary layer[END_REF].

POD-based reconstructions

We have seen that POD provides an energy-optimal basis to decompose a field. POD modes therefore constitute natural candidates for structure-based reconstructions. The requirement is that these modes are determined a priori, which requires an early computation. If the field Q can be represented as

Q(x, t) = n a n (t)φ n (x), (3.34)
with a known basis φ n , the question is to determine the amplitudes of the modes a n (t).

The Poitiers group [26; 61; 62] has been instrumental in developing methods to provide amplitudes, ranging from random perturbations to coherent structures. Perret et al [START_REF] Perret | Generation of turbulent inflow conditions for Large Eddy Simulations from stereoscopic PIV measurements[END_REF] used a spectral transfer function to determine the spectral energy content of the temporal amplitudes to generate turbulent inflow conditions for the large-eddy simulation of a mixing layer as shown in Figure 3.10. This method was further refined in [START_REF] Perret | Turbulent inflow conditions for Large-Eddy Simulation based on low-order empirical model[END_REF], where the flow was separated into coherent and incoherent parts. The most energetic structures were represented with a low-dimensional model derived from the experimental database, while Gaussian amplitudes were used to represent the small-scale incoherent motion. A variant of the method was proposed by Druault et al. [START_REF] Druault | Generation of three-dimensional turbulent inlet conditions for large-eddy simulation[END_REF] who used linear stochastic estimation to identify the large scales of the flow from experimental measurements. More recently, Abéguilé et al [START_REF] Abéguilé | Générateur de conditions aux limites amont pour les simulations de type LES des écoulements de paroi[END_REF] have tested neural networks to generate inlet conditions for LES of wall-bounded flows.

Wall Synthetic boundary conditions 3.3.3.1 Current approaches

Bergmann and Iollo proposed a domain decomposition method which couples a high and a low-fidelity model to reduce the computational cost of flow simulation [START_REF] Bergmann | A zonal galerkin-free POD model for incompressible flows[END_REF]. The idea is to identify parts of the domain which are not strongly sensitive to the choice of the design parameters. In these areas a low-fidelity simulation is sufficient and can be achieved from a POD-based projection -the spatial POD modes are pre-determined from an earlier computation with a high fidelity simulation on the full domain. They divide the field into a small portion where it is solved by high-fidelity model, and the field in other regions is described by a Galerkin-free Proper Orthogonal Decomposition. They proposed a prediction error indicator for splitting the computational domain between the high-fidelity and low-fidelity regions and coupling the two solutions. POD reconstruction is then used in the low-fidelity regions while the high-fidelity solver is applied only in the region where the flow field depends strongly on the choice of the design parameters.

In channel flow turbulence it is quite challenging to define a synthetic boundary condition near the wall since the wall region is crucial for the generation of turbulence. A large part of the turbulence production occurs in the layer z + < 50 [START_REF] Jiménez | Near-wall turbulence[END_REF]. Special care must therefore be given to the synthetic boundary condition so that the relevant pattern dynamics can be reproduced. In incompressible flow, a reconstruction procedure was developed to reconstruct the field on the boundary of a reduced domain which excludes a portion of the wall layer. Due the incompressible formulation, only the velocity field was considered and the velocity field was prescribed on the entire domain. In the next section, we extend this procedure by adapting it to the compressible formulation, that is: i) multiple physical quantities, and ii) in the context of wave propagation.

The reconstruction procedure

A pre-requisite for the procedure (step 0) is the POD basis for the reference flow. Which means the POD basis is computed from reference simulation with snapshot methods.

Then at each time step, the reconstruction procedure can be divided into several stages:

1. Estimation of the amplitudes using linear estimation 2. Reconstruction of the velocity field 3. Rescaling of the field

Determination of the characteristics and construction of the boundary condition

We will focus on the following aspects:

• the choice of the POD variables (step 0):

-primitive or conservative -selection of variables to aggregate

• the definition of the rescaling (step 3)

• the association of the variables with characteristics (step 4)

We now describe the last two steps.

Step 3: Rescaling

In the previous chapter we mentioned the convergence issue of the truncated POD basis.

Especially in the case of the method of snapshots, owing to the complexity of the flow, it is almost certain that a field that does not belong to the original set of snapshots will not coincide with its projection on the POD basis. The projected flow will therefore lack energetic content.

Another reason for an incorrect energetic content has to do with the approximation we use for R ∼ R, which results in an estimation error. It is not possible to control the magnitude of these errors. There is therefore no guarantee on the estimated amplitudes of the POD modes.

Since a correct energetic content in the fluctuations is essential for turbulence to be sustained, a rescaling factor was therefore introduced in order to ensure that the energetic content of the fluctuations is sufficient. The idea is to impose that the rms level of the reconstructed fluctuating field matches the reference level on the boundary plane z = h 0 .

For a given quantity Q, we define a factor f Q which compares the rms level of the quantity Q on the plane in the reference simulation (full channel) to that of its estimation

f Q = Q rms ref erence Q rms estimated . (3.35)
Rescaling can be done individually for each component of the vector Q. The estimated field could be rescaled as follows:

Q Est i (x, y, z 0 , t (m+1) ) = (Q P i (x, y, z 0 , t (m+1) ) -Q i mean ) × f Q + Q i mean , (3.36)
where i is the indice of quantity, Q P i is estimated field with POD reconstruction, Q i mean is spatial mean field of estimated field

Q P i .
In this thesis, the rescaling method is the compensation, they are the same things.

Step 4: Implementation of the reconstruction

We now detail how we determine the boundary condition from our estimation and the flow characteristics.

There, to solve the Navier-Stokes equations, we can treat the viscous term as external term S and the equation to solve is :

∂Q c ∂t + ∇ • F (Q c ) = S. (3.37)
The left side of (3.37) is the Euler equation and the right side is divergence of the viscous term which could be expressed as ∇ • F vis (Q c , ∇Q c ). This equation can be written as:

∂Q c ∂t + ∂F (Q c ) ∂Q c • ∇Q c = ∇ • F visc (Q c , ∇Q c ). (3.38)
where Q c represents the conservative variables [ρ ρu ρv ρw ρE ] T . ∂F (Qc) ∂Qc is the Jacobian matrice:

∂F (Q c ) ∂Q c = R s ΛR -1
s that has real eigenvalues λ, which could be expressed as (u

• n, u • n, u • n, u • n+ n C , u • n-n C
), and n is the direction normal to the cell face considered. Right eigenvectors are expressed as follows [74]:

R s =            u -C 0 0 u u + C u -C 0 0 u u + C w 0 1 w w H -uC v w u 2 +v 2 +w 2 2 H + uC 1 0 0 1 1            , (3.39)
and the left eigenvectors are the rows of

L s = R -1 s = C 2 2            -2v/C 0 2/C 0 0 -2w/C 0 0 2/C 0 2(1 -a 2 ) 2a 1 u 2a 1 v 2a 1 w -2a 1 a 2 -u C -(a 1 u -1 C ) -a 1 v -a 1 w a 1 a 2 + u C a 1 u + 1 C -a 1 v -a 1 w a 1            , (3.40)
where

a 1 = γ-1 C 2 and a 2 = 1 2 (u 2 + v 2 + w 2 )a 1 .
When (3.38) are projected onto the left eigenvectors R -1 s , we obtain:

R -1 s ∂(Q c ) ∂t + ΛR -1 s ∇ • Q c = R -1 s DQ c Dt = R -1 s ∇ • F visc (Q c , ∇Q c ). (3.41)
The right term could be considered as the source term R -1

s S = R -1 s ∇ • F visc (Q c , ∇Q c
), and it is known, so the main work is to reslove Euler equations: In Naiver Stokes equation, the Riemann invariants could be expressed as:

R -1 s ∂(Q c ) ∂t + ΛR -1 s ∇ • Q c = 0. ( 3
dω = R -1
s dQ c -δtR -1 s S. Then we applied this method in our channel for wall parallel plane. Figure 3.12: Boundary condition on the wall parallel plane.

There, in Figure 3.12, we could see that if characteristics wave is outwards, the eigenvalue λ is negative (λ < 0), then it corresponds to the numrical condition, which means to calculate the field at time t n at position z 0 , it could be calculated from field at time t n-1 in resolved region. And if characteristics wave is inwards, the eigenvalue λ is positive (λ ≥ 0), it corresponds to the physical condition, which means to calculate the field at time t n at position z 0 , we need to estimated the field in implicit region at time t n-1 , in this case, POD method is applied for estimation.

In more details, if characteristics wave is outwards, the eigenvalue λ is negative (λ < 0), and we need to solve the Navier-Stokes equation with Riemann invariants in order to avoid the instability. The approximation of flux in (3.38) is second order accurate, which means (3.41) is solved as:

Q c n+1 -Q c n δt = - ∂(F n (Q c ) -F n visc (Q c , ∇Q c )) ∂x = -( 3F n i-1 (Q c ) -4F n i (Q c ) + F n i+1 (Q c ) 2δx - 3F n visc(i-1) (Q c , ∇Q c ) -4F n visc(i) (Q c , ∇Q c ) + F n visc(i+1) (Q c , ∇Q c ) 2δx
).

(3.43)

This equation provides solutions Q c n+1 when λ < 0 .

If the characteristics wave is entering the domain, which means λ > 0, we have to prescribe the corresponding Riemann invariant whith predicted field in implicit region:

P n+1 = P * ⇒ ∂P ∂Q c ∂Q c ∂t δt = P * -P n . (3.44)
Where P * is predicted field in implicit region. When we consider the variables [ρ u v w T ] T . ∂P ∂Q can be calculated using :

∂P ∂Q c (Q n+1 c -Q n c ) = P * -P n . (3.45) ∂P ∂Q c =               1 0 0 0 0 - u ρ 1 ρ 0 0 0 - v ρ 0 1 ρ 0 0 - w ρ 0 0 1 ρ 0 u 2 + v 2 + w 2 -E C v ρ - u C v ρ - v C v ρ - w C v ρ 1 C v ρ               . (3.46)
When using the variables [ρ ρu ρv ρw T ] T on boundary condition:

∂P ∂Q c =             1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 u 2 + v 2 + w 2 -E C v ρ - u C v ρ - v C v ρ - w C v ρ 1 C v ρ             . ( 3.47)
The variables we consider are • For the three first eigenvalues u • n (n is inward wall normal direction):

[Q 1 , Q 2 , Q 3 , Q 4 , Q 5 ] = [ρ, (ρ)u, (ρ)v, (
when u(x, y, z 0 , t m )•n > 0.:

         (ρ)u(x, y, z 0 , t (m+1) ) = Q Est 2 (
x, y, z 0 , t (m+1) ); (ρ)v(x, y, z 0 , t (m+1) ) = Q Est 3 (x, y, z 0 , t (m+1) ); T (x, y, z 0 , t (m+1) ) = Q Est 5 (x, y, z 0 , t (m+1) ); -otherwise, the equation of the three first Riemann invariants are solved with resolved region.

• For the fourth eigenvalue u • n+ n C (C is the local speed of sound):

when (u(x, y, z 0 , t m )•n+ n C) > 0.: (ρ)w(x, y, z 0 , t (m+1) ) = Q Est 4 (x, y, z 0 , t (m+1) );

otherwise, the equation of the fourth Riemann invariant is solved with resolved region.

• For the fifth eigenvalue u • n-n C:

when (u(x, y, z 0 , t m )•n-n C) > 0.: ρ(x, y, z 0 , t (m+1) ) = Q Est 1 (x, y, z 0 , t (m+1) );

otherwise, the equation of the fifth Riemann invariant is solved with resolved region.

First test: Reduced simulation using reference flow fields as boundary conditions

As a first test, we collect fields from the reference simulation on a boundary plane within the wall layer and use these fields as synthetic boundary conditions in the reduced domain. This is done for a single wall.

We chose to associate the variables[T u v w ρ] with the characteristic eigenvalues

[u • n, u • n, u • n, u • n+ n C, u • n-n C],
where n is the direction normal to the wall: wise, spanwise, and wall normal directions of the flow will be respectively denoted by

u • n-n C is
x, y, and z; the corresponding components of the velocity are u, v and w. Periodic boundary conditions are applied in both the streamwise and the spanwise directions.

The solid wall at the upper part of the domain is maintained at a prescribed temperature (T w ), and the lower boundary plane is determined by synthetic boundary condition.

The dimensions of the reduced channel are (L x , L y , L z ) = (2π, 4π 3 , 1.9)H. The friction Reynolds numbers is Re τ = 180. The Mach number is M a = 0.5.

We choose a longitudinal section x = π to compare the DNS and reduced channel. In the simulation, the boundary region of 0 < z < 0.1 was taken off, and we used the original field from full channel as a boundary condition. The simulation was run for 500 time steps corresponding to 0.3U 0 /H in both full channel and reduced channel.

Computational basis

Here, we give details for the simulation in full channel for both subsonic flow and supersonic flow in Table 3 The profiles of full channel will be compared with results in reduced channel in later.

Chapter 4

Synthetic boundary condition on one wall

We firstly test the synthetic boundary condition on one wall of the channel, represented in The vertical extent of the reduced channel is defined by z ∈ [0. [START_REF] Abéguilé | Générateur de conditions aux limites amont pour les simulations de type LES des écoulements de paroi[END_REF][START_REF] Adrian | Stochastic estimation of conditional structure: a review[END_REF], where h = 0.1 corresponds to h + 0 = 18. The region where the turbulent production is maximal z + ∼ 15 is therefore excluded from the simulation. The mesh dimension is reduced by 17% compared to the one in full channel.

POD was applied independently to three variables: the density (ρ), the velocity field (u) and the total energy (E). As discussed in the previous chapter, POD-based reconstruction with a limited number of modes tends to underestimate the magnitude of the field in particular the small scales, and therefore to underestimate the wall-normal velocity field. In order to compensate for this effect, a rescaling technique is implemented. It consists in imposing a global rms value for the reconstructed fluctuations based on the reference level measured in the full channel. As we recall from Chapter 3 (see Equation 3.36) , for any physical quantity Q,

Q Est (x, y, z 0 , t (m+1) ) = (Q P (x, y, z 0 , t (m+1) ) -Q mean ) × f Q + Q mean . ( 4
Q DN S mean (z) = Σ x,y,t Q DN S (x, y, z, t) N t • N x • N y . (4.
2)

Q P mean (z, t) = Σ x,y Q P (x, y, z, t) N x • N y . (4.3)
In all that follows we compare the reference channel to

• the reduced channel with one synthetic boundary condition without rescaling,

• the reduced channel with one synthetic rescaled boundary condition where both definitions of the mean were tested. Mean profiles of velocities and density are shown in Figure 4.2. We can see that the mean profile of density matches very well the reference profile. Again the introduction of rescaling step through rescaling improves considerably the results in particular for the mean streamwise profile, for which the logarithmic law appears to be well recovered. We can see in Figure 4.3-4.4 that the statistics in the reduced channel are closer to those in the reference channel for the rescaled boundary condition, in particular for the wall-normal velocity which has a significant small-scale content.

Results are shown in

We note that with the choices made for the boundary condition, the wall-normal component corresponds to the eigenvalue u • n+ n C = w + C > 0, which is always positive, so that the wall-normal component is always determined from the POD reconstruction at all points on the boundary. In contrast the streamwise and spanwise components are associated with the eigenvalue λ = u•n = w, which means that the reconstruction is only used on cells with a positive (inwards) wall-normal velocity w > 0. For the cells with a velocity directed towards the wall, the boundary condition is evaluated from the fluxes calculated from the inner domain. This discrepancy between the different components of the velocity field may lead to errors that are difficult to control.

When defining the rescaling factor f Q , the mean value can be evaluated either locallyi.e from a horizontal spatial average at each instant -or globally; using both a time and a horizontal spatial average. Comparison between these two methods of rescaling step is presented in Figures 4.5 To evaluate the performance of the method, we computed the relative error between statistical quantities in the internal region (z + ∈ [START_REF] Mizuno | Wall turbulence without walls[END_REF]180], where z ∈ [0.3, 1]) between different conditions in reduced channel and DNS. The relative error is calculated as

e Q = 1 z 2 -z 1 z 2 =1 z 1 =0.3 |Q Reduced stat -Q DN S stat | |Q DN S stat | dz, (4.4)
where Q Reduced stat and Q DN S stat correspond to time and horizontal averages in the reduced channel and in the reference channel respectively. We consider the internal region from z = 0.3 to z = 1 to leave some distances for adjustment of boundary conditions. Relative errors for rms of u, v, w, and the mean of density and streamwise velocity u are reported in Table 4.1. For the method with rescaling step, we compare the two definitions of the mean for the rescaling procedure. The difference between the two rescaled methods is slight, with a small advantage for the definition based on the local value, particularly for the density.

Consequently, for the simulation with primitive variables, rescaling step with local mean value is the best choice.

Results at altitude h + 0 = 18 with conservative variables

The comparison was also carried out for a POD based on conservative variables [ρ], [ρu, ρv, ρw], [ρE]. As in the primitive variable case, POD was applied independently to the density [ρ], the momentum [ρu, ρv, ρw], and the total energy [ρE]. We tested:

• the reconstruction without rescaling step,

• the reconstruction with rescaling step using the mean value which is calculated from the reference simulation (4.2),

• the reconstruction with rescaling step using the local mean value from the reduced simulation (4.3).

Firstly, we compare the results without rescaling step and with rescaling step with global mean value. The errors defined as in (4.4) are reported in Table 4.2. We can see that using the local mean value in the rescaling procedure leads to smaller errors than when using the reference value.

Comparison between primitive and conservative variables in reduced channel

As seen in Figures 4. 13-4.15, there is no significant difference between the method with primitive variables and conservative variables, at least in the conditions of our numerical experiment. Small discrepancies are however recorded in the wall normal velocity component (w). This can be explained by the fact that, as ρw is prescribed in the conservative approach, the mass flux across the boundary plane is better controlled in compressible approach than when using the primitive variables.

Instantaneous flow fields

We compare instantaneous flow fields on a cross-stream plane at the location x = π at the end of the short simulation time that is 6 time units. These fields are obtained There is not evident difference between fields with conservative variables and primitive variables except the field of density.

The density field is different for the reference simulation and the reduced channel with conservative variables or primitive variables. This discrepancy can be explained by the fact that it is calculated by solving the equation on the first Riemann invariant with a second-order upwind approximated flux from the the reduced channel versus a seventhorder discretization in the full simulation. However the magnitude of the discrepancy is small.

The large scales of the streamwise velocity field are well recovered in the reduced channel.

Small discrepancies can be noticed in the cross-stream components such as the spanwise velocity and the wall-normal velocity, but the overall agreement is good. Ejections and sweeps corresponding to relatively large values of w are correctly captured in both two reduced channels. These results indicate a good agreement between the reference simulation and the reduced channel after a limited amount of time.

Results at height h

+ 0 = 54 (h 0 = 0.3)
The same procedure was implemented with a boundary plane located at a height of h + 0 = 54 corresponding to h = 0.3. The mesh size in the reduced channel therefore decreases by more than 29% compared to the full channel. POD was applied independently to three primitive variables:

• the density [ρ]
• the velocity field [u, v, w]

• the total energy [E] is observed for density. The logarithmic region for u appears to be recovered (we note that its extent is small at this relatively moderate Reynolds number).

As already observed for the smaller height h + 0 = 18, rescaling step through rescaling improves the agreement significantly for turbulent intensities and shear stress.

The relative errors defined as (4.4) are reported in Table 4.3. We calculated the relative errors over a region z ∈ [0.3, 1.0], which contains all the reduced region above h + 0 = 54 (h = 0.3) in bottom half channel. We can see that the error of rms values between reference and reduced channel is small. We also tested the two definitions of the mean value for the rescaling factor f Q (see Figures 4.23 Examination of errors in Table 4.3 shows that the relative errors of density are very small in the three methods. The error on the mean velocity is also very small. The largest error is observed for the wall normal velocity component w. Using rescaling step decreases the error on all variables except the density (which remains very small). The error is divided by 40% for < u >, by a factor of 2 for rms of w, by a factor of 3 for rms of v and almost by 5 for rms of u. In the simulation, we need some distance to adjust the turbulent intensities, and the reduced channel statistics are virtually identical with the reference channel at a height of z + = 60 for u, z + = 100 for v, and z + = 140 for w. again, using the local mean value seems to yield slightly better results (here for the mean density). This is confirmed by the quantitative measure given in Table 4

Quantity error

Comparison between reduced-channel simulations based on POD with primitive variables and with conservative variables

We also compare the results between POD reconstruction with primitive variables and conservative variables using a rescaling factor based on the local mean value. There is no real difference between the two approaches, primitive or conservative variables, except the profile of the mean velocity in the spanwise and normal the wall directions.

Comparison of Tables 4.3 and 4.4 shows that the errors are very close. However the mean profile of w is closer to the reference for the conservative variables. Because with conservative variables, we can better control the mass flux in reduced channel.

Summary

In this chapter we consider a reduced channel where one wall is replaced by a virtual plane in the wall layer and compare its statistics with those of a reference channel.

The streamwise and spanwise velocity components are rather well recovered. However some discrepancies can be observed for the wall-normal component. We examined the influence of the basis used to compute POD and found that it did not affect the results.

We found that rescaling the boundary condition so as to adjust its energy content led to a significant improvement of the statistics, in particular for the wall-normal component.

Two types of rescaling were considered, depending on the definition of the mean. We found that using a local mean value led to a slightly closer agreement with the statistics.

Using conservative or primitive variables did not change the results significantly. This

shows the robustness of the approach.

The spatial autocorrelation tensor is reconstructed for each Fourier mode from 200 samples extracted from 100 snapshots (as explained in the previous chapters, one snapshot corresponds to two samples, as the span of each sample is the half-channel). Each sample was therefore transformed in Fourier space and the POD problem in horizontal Fourier space was solved only in the inhomogeneous direction z. In tensor notation, the eigenproblem can be expressed as

z R ij kxky (z, z )φ * jn kxky (z )dz = λ n kxky φ in lk (z), (5.1)
where n refers to the POD mode index and the spatial autocorrelation tensor is defined by

R ij kxky (z, z ) =< Q i kxky (z, t)Q * j kxky (z , t) >, (5.2)
with < . > a spatial average over wall-normal direction, and Q i kxky represents the horizontal Fourier transform of the i-th component of the Q field

Q i kxky (z, t) = x y Q(x, y, z, t)e 2 √ -1π( kx Lx + ky Ly ) dxdy, (5.3)
which is evaluated numerically as

Q i kxky (z, t) = 1 N x N y Nx/2 m=-Nx/2 Ny/2 p=-Ny/2 Q(x m , y p , z, t)e 2 √ -1π( kx Lx + ky Ly ) δ x δ y . (5.4)
δ x = Lx
Nx-1 and δ y = Ly Ny-1 are spatial intervals in streamwise and spanwise directions respectively. The dominant eigenvalue for each Fourier wavenumbers pair is represented in Figure 5.1 left, where the POD reconstruction is applied on [ρ], [ρu, ρv, ρw] and [ρe] respectively.

Synthetic boundary conditions at h

+ 0 = 18 (h 0 = 0.1)
As indicated in Chapter 3, we solve only for modes n for which the diagonal part R nn is larger than a threshold ts -in practice we use ts = 0.3. This means that only modes

Unrescaled boundary conditions

Firstly, the initial field is reduced region which is taken from full channel, we apply synthetic boundary conditions for the field and compute turbulent statistics over a total period of 6 time units. In reduced channel, we did not apply rescaling step on variables, which means no compensation. There, the profile of full channel is calculated for the same duration as in reduced channel. We can see that the streamwise intensity of u is well recovered in Figure 5.3, as well as the cross-stream components v and w. The Reynolds stress is also well recovered.

Although the mean velocity profile in the reduced channel appears to be relatively close to the reference one, it turns out that the computation becomes unstable after 8 convective time units. We have not elucidated the reason for this. In all that follows we will present results for the method of snapshots. For each decomposition, the boundary condition was reconstructed and implemented on each wall for the reduced domain at a height h + 0 = 18 i.e. a distance of 0.1 from the walls in nondimensional length (see Figure 5.9). In all that follows unless specified otherwise results correspond to statistics obtained over a period of 36 convective time units, which is equal to 36 H U 0 . And the statistic duration in full channel is 120 time units for the following results in this chapter.

Results without rescaling

We first compare the four POD methods defined above in Figures 5.10 -5.12. [ρe] is more adequate. We conclude that results are not satisfactory when the reconstruction is performed without rescaling.

Results with rescaling

We now consider rescaled boundary conditions. Based on results from the previous section, we did not consider the decompositions based on primitive variables. For the decomposition of [ρ], [ρu, ρv, ρw], [ρe], which yields the best results, two types of rescaling were considered. In one case, rescaling was applied to each variable, while in the other case rescaling was only applied to the momentum -note that the rescaling factor was defined independently for each velocity component. Comparison of the simulation in the reduced domain with the reference case is shown in Figures 5. [START_REF] Bose | Wall-modeled Large-Eddy Simulation for complex turbulent flows[END_REF] Figure 5.13 presents the mean velocity profiles. There are relatively few differences between the decompositions. A major result is that for all cases the logarithmic law is recovered but tends to be slightly underpredicted (2%) in the center of the channel. A summary of the performance of the synthetic boundary conditions computed with different POD is given in Table 5.1. As in the previous chapter, the error on a statistical profile Q stat (z) where stat refers to a mean value or a rms, is computed as

e Q = 1 z 2 -z 1 z 2 z 1 |Q Reduced stat -Q DN S stat | |Q DN S stat | dz, (5.5)
where z 1 = 0.3 and z 2 = 2 -0.3, which corresponds to a distance to the wall of about 55 wall units. Results are given for rescaled and unrescaled conditions.

Even if h + 0 = 18 (h 0 = 0.1), we use a distance of z 1 = 0.3 in order to leave some height over which the flow can adjust to the boundary conditions. • For the mean quantities, the rescaled conditions give better results than the unrescaled ones. The error for the mean velocity is on the order of 5% for unrescaled conditions, and around 2% for scaled ones. While the amplitude of density and momentum oscillate near constant values, a slow decrease is observed for the energy, which tends to stabilize for the last 240000 iterations (last 25 fields in Figure 5.16, 144 time units). This suggests that the simulation has just reached statistical convergence. Our original snapshot basis for POD (corresponding to was obtained over a period of non-stationary convergence. To examine if this affects the synthetic boundary condition, we constructed another POD basis from 100 snapshots corresponding to the period labeled "new samples" in Figure 5.17 This choice is due to memory limitations for the storage of the POD basis of momentum.

Some choices are possible in Riemann invariants. Moreover, when prescribing the velocity field (for instance when u • n > 0) it is possible to set one of its components equal to zero.

Various combinations were tested, A summary is provided in Table 5.

Correction step in the estimation procedure of the POD amplitudes

In this paragraph, we explore a way to address the issue mentioned earlier, namely that the variables associated with the characteristic w are either propagated from within the simulation or prescribed by the POD reconstruction, which may result in discrepancies on the boundary plane and the generation of energetic small scales. To reduce these discrepancies, the idea is to use the boundary condition once it is determined on the limit plane to improve the determination of the POD amplitudes. This means solving an additional linear problem similar to that of the standard reconstruction method, but only limited to the plane of reconstruction.

Once the field Q(x, y, z 0 , t) is determined with the standard reconstruction procedure described in Chapter 3, we look for the best estimate A n (t) such that

Q(x, y, z 0 , t) = N n=1
A n (t)φ n (x, y, z 0 ). (5.6) This leads to the normal equations

A n (t) = Ψ -1 nm γ m , (5.7)
where

γ m =
x y φ m (x, y, z 0 )Q(x, y, z 0 , t)dxdy,

Ψ nm =
x y φ m (x, y, z 0 )φ n (x, y, z 0 )dxdy.

(5.8)

From figure 5.24, we can see that the diagonal terms in Ψ nm are much larger than others terms, we except the pattern to the well-conditioned, and it could be better for resolving the linear system. We then applied the adjustment condition to calculate new estimated temporal coefficient (condition C). The reason for these relatively disappointing results might be the absence of a correction on the ghost cell planes, which will affect the computation of the fluxes.

5.5 Spectra in the reduced channel at h + 0 = 18

Based on the statistic of first order and second order in section 5. Rescaling of the momentum components is applied to the synthetic boundary condition.

We continued the simulation to 72 time units and present the results in Figures 5. 28-5.31. We can see that the turbulent intensities tend to be overpredicted compared to the results obtained for 36 time units. We point out that no global rescaling of the full flow fluctuations was performed at any time during the computation. We observe that the streamwise turbulent intensity (< u 2 > 1/2 ) is overpredicted by 26.75 %, the spanwise turbulent intensity (< v 2 > 1/2 ) is overpredicted by 31.52 %, and the wall normal turbulent intensity (< w 2 > 1/2 ) is largely overpredicted by 70 %.

However the shear stress (Figure 5.31) is slightly underestimated, which means that the 5.5 Spectra in the reduced channel at h + 0 = 18 115

The integration time was 72 time units. The general distribution of the kinetic energy is rather well captured by the reduced channel simulation. We can however notice a shift of the kinetic energy towards greater wave numbers in the reduced channel, in both streamwise and spanwise directions. This suggests that POD reconstructed synthetic boundary conditions contain more kinetic energy in smaller length scales than DNS.

The shift is clearly noticeable in the spectra of the streamwise (Figure 5.32) as well as the wall-normal component (Figure 5.34) of the velocity. It explains the higher than expected turbulent intensities. We speculate that it could be created by the mismatch between the synthetic boundary condition and the flow in the inner domain. We note that the size of the inertial region is small, due to the relatively low value of the Reynolds number.

Results at h

+ 0 = 54 (h 0 = 0.3)
The synthetic condition considered above was applied at a larger height h + 0 = 54. The test was run over a total time of 36 convective time units for the statistically converged POD basis (obtained with the set of snapshots corresponding to the "new samples").

For POD reconstruction, all modes are used for density ρ and internal energy ρe, 120 modes are used for momentum. And the rescaling step is applied just on momentum.

Conclusion

In this chapter we have applied synthetic wall boundary conditions to both walls of a channel. The influence of various parameters was examined. We first examined a Fourier-based reconstruction, but we were not able to obtain a stable formulation. The reason for this is unclear as the approach was successful in the incompressible case.

We then tested different decompositions for the variables. We found that the best choice was a decomposition on [ρ], [ρu, ρv, ρw], [ρe]. The mean velocity profile was well recovered and the logarithmic region was correctly reproduced. However it seems that without further rescaling the turbulent intensities diverge slightly over very long time scales.

We also examined the influence of the snapshot basis, the Riemann invariants, and we tried to improve the estimation of the POD amplitudes. We found that examination of the spectra provides useful insight to evaluate the synthetic boundary condition. Two different heights were tested h + 0 = 18 and h + 0 = 54. Results were found to be largely insensitive to the height of the boundary plane. However we found that the procedure was sensitive to the choice of the variables for the decomposition. The reason for this is not entirely clear and will require further investigation. By selecting the variables in an appropriate manner, we were able to recover the mean velocity profile within 2%. The Reynolds shear stress was also correctly predicted. However, the turbulent intensities tended to be overpredicted. Spectral analysis showed that this overprediction was linked to the excess energy at higher streamwise wavenumbers than in the DNS.

The potential advantage of using synthetic wall boundary conditions is to significantly reduce the computational cost of a simulation, through • a reduction in the mesh size, due to the exclusion of the thin but highly resolved wall region

• a reduction in the computational time step, due to the increase in the size of the smallest cell However this computational gain is limited by the extra cost required to compute the synthetic boundary condition. Table 5.4 summarizes the characteristics of the reference and of the reduced simulations for the two boundary heights. We note that no effort was made to optimize the implementation of the boundary condition, which may explain why the total reduction of the computational cost is relatively small compared to the mesh and time step reduction. In order to reduce the memory requirements, we carried out an interpolation procedure in order to compute POD modes on a reduced channel. The condition was then reconstructed on the full mesh. Therefore, an interpolation procedure on a smaller grid has been employed to reduce the dimension of the problem [START_REF] Hazewinkel | Spline interpolation[END_REF]. The mesh size was reduced by 2 in each horizontal direction. All the points in the normal to the wall direction have however been retained leading to a new mesh with N x new ×N y new ×N z = 129×81×161 grid points.

Boundary plane height

The spatial modes are then calculated on the new grid. As far as the reconstructed synthetic conditions on the boundary planes require the knowledge of the POD spatial modes at each points of this plane, an interpolation procedure is needed to estimate the conservative variables at each points of the original mesh on boundary planes. For simplicity, a bi-linear interpolation is employed in each plane (x × y). If x and y are the coordinates of the point in the original mesh, and x 0 , y 0 , x 1 , and y 1 are the coordinates in the new mesh (see Figure 6.1), a function (f ) known at each point of the new mesh is interpolated in the original mesh following: in the supersonic case. We note that in contrast, there are differences for density in Figure 6.6, which is always computed from information within the channel. The density therefore appears to be sensitive to the reconstruction procedure.

f (x, y) = c 0 (1 -y d ) + c 1 y d , with c 0 = f (x 0 , y 0 ) × (1 -x d ) + f (x 1 , y 0 ) × x d , c 1 = f (x 0 , y 1 ) × (1 -x d ) + f (x 1 , y 1 ) × x d , and
x d = x -x 0 x 1 -x 0 , y d = y -y 0 y 1 -y 0 .

Statistics in reduced channel

Simulation with POD reconstruction of first 35 samples

In Figures 6.7-6.9, we present the results in reduced channel averaged over 26 time units, the DNS results are obtained from 75 time units, and the rescaling step is applied just on momentum. We compare the results for the simulation without rescaling step and with rescaling step and with the full channel data. From the results, we could see that the results with rescaling step are much better than without rescaling step. We calculated the relative error (Table 6.1) of quantities in the center region (z 1 < z < z 2 ) with z 1 = 0.3 and z 2 = 1.7: We could see that the rms values are rather well recovered in the center region of the channel. The rms value of w is however slightly overestimated. Besides the rms values, the shear stress < u w > is very well recovered in the center region.

e Q = 1 z 2 -z 1 z 2 z 1 |Q Reduced stat -Q DN S stat | |Q DN S stat | dz. ( 6
Relative errors are then calculated (6.1) for the solution over 40 time units; they are reported in Table 6.2. While relative errors on the rms velocity profiles have increased from the simulation at 26 time units, they stayed rather weak and lower than 0.1. We can notice that these relative errors are much smaller than in the subsonic flow (see the previous Chapter), indicating that the synthetic boundary conditions seem to work better in supersonic flows. The sonic line M a = 1 was found to be at z = 0.08, which is below the boundary plane.

This means that the local Mach number for the synthetic boundary condition is higher than 1, so that information is enduringly propagated downstream. However in subsonic flow, the propagated direction could be either downstream or upstream as Figure 6.14 shows. It is possible this phenomenon could be helpful for the stability of calculation for supersonic flow in reduced channel. .

Simulation with POD reconstruction using new 30 samples

In the previous section, results were obtained with POD calculated on only 35 time samples. To obtain more samples as in the subsonic flows, we resumed the DNS for 990,000 supplementary time steps (100 samples). We performed 3 POD's for [ρ], [ρu, ρv, ρw], [ρe] on 100 samples , separated by t = 2.3 H U 0 . On Figure 6.15, we plot amplitudes of the first POD mode of density, momentum and internal energy.

We easily notice that time variations of the first mode amplitudes of the density and the momentum are in opposition phase between the upper and lower half part of the channel. We also notice that these time evolutions do not show any drift over time (Figure 6.15). This is not the case for the internal energy where the amplitude of the first POD mode clearly shows a drift over the first 70 samples that distinctly decreases over the remaining 30 last samples. We then employed the last 30 samples, corresponding to set B of As for the computational time saving in Table 6.4, the mesh dimension in reduced channel is 50% of that in full channel, and the time step is 1.7 time of that in full channel, which means than the efficiency what we could obtain should be around 29%, however, because POD reconstruction also takes the time, so finally the real efficiency is around 42 %, and the gain is around 2.5. We can see that the spectrm of energy of u is rather well recovered in reduced channel. A shift towards larger wave numbers in streamwise direction (k x ) could however be observed. The energy spectrm of v in the reduced channel coincides rather well with the reference. Regarding the energy spectra of w, a shift is also observed towards small wavenumbers in the spanwise direction (k y ).

Spectra in the supersonic flow

In conclusion, the spectra of velocities are rather well recovered in reduced channel, and they are better than those in subsonic flows, the reason could be that in supersonic simulation, when the boundary plane is higher than sonic line, that information is enduringly propagated downstream. Which could be helpful for stability of simulation.

Chapter 7

Conclusions and perspectives

Conclusion

The objective of this thesis is to explore how synthetic wall boundary conditions could be developed for the compressible simulation of wall turbulence. The idea is to bypass entirely the simulation of the wall region which requires a very fine spatial resolution near the wall, consequently reduce the computational time step, and to replace it with a boundary condition located on a plane within the wall layer which mimics the features of full wall turbulence at that height. The computational domain is thus reduced, as we only simulate the flow at some distance from the wall.

The originality of the thesis is to extend the approach traditionally used for inlet boundary conditions to wall conditions, and to rely on Proper Orthogonal Decomposition modes to construct the boundary condition. The challenge of the compressible formulation is that the nature of the boundary condition depends on the wave propagation characteristics. This entails that the success of the procedure depends on the adequacy of the boundary condition prescribed with Proper Orthogonal decomposition (POD) with the solution propagated along the Riemann invariants.

Using POD, the fields on the boundary which need to be prescribed are written as a superposition of spatial modes, which are assumed to be known before-hand. The amplitudes of the modes need to be determined at each time step. This is done using linear stochastic estimation based on the information available in the reduced computational domain at the previous time step.

In the course of this thesis different aspects of the procedure were examined and tested.

The goal is (i) to examine whether the procedure is robust to the different choices that can be made, (ii) to determine the parameters that lead to the best results (iii) to provide physical or numerical arguments for these choices. We emphasize that our approach remains exploratory at this stage.

The procedure was tested in a turbulent channel flow with periodic boundary conditions in horizontal directions. The Reynolds number based on the bulk velocity and channel half-height was Re = 3000. The temperature remains constant at the walls. One subsonic case (Chapters 4 and 5) and one supersonic case (Chapter 6) are considered. They respectively correspond to the cases M a = 0.5 and M a = 1.5. The synthetic boundary conditions were first implemented on one wall (Chapter 4) then on both walls (Chapters 5 and 6). The statistics were computed over a total period of several decades of convective time units. Two heights within the wall layer were considered h + 0 ∼ 20, and h + 0 ∼ 50. Results were found relatively robust with respect to the boundary height.

In satisfactory cases, above some recovery height, a very good agreement was observed between the statistics of the reference channel and those of the reduced computational domain with the synthetic boundary condition, namely:

• The reference mean velocity was computed within 1-2% in the computational domain and the logarithmic law was recovered.

• The Reynolds stress was correctly estimated.

The turbulent intensities were overall well predicted, but the cross-stream turbulent intensities were overestimated. Spectral examination showed that this corresponded to a higher than expected energetic content in the small scales.

An important question is to define how POD should be applied to the different variables.

It was found that the conservative formulation was generally preferable to the primitive formulation. It was also found that it was best to use three independent decompositions for the density, the momentum and either the energy or the temperature. This may be because these quantities are dominated by different scales, so that splitting the decomposition improves the convergence of the representation. Among all the choices possible for the variable associated with energy, tests were conducted for the total energy, internal energy and the temperature. Best results were observed for the internal energy.

We checked that the results were largely insensitive to the choice of the set of snapshots used to perform POD, once statistical convergence has been established. Agreement with the reference statistics was observed for a few tens of POD modes. The maximal number of POD modes included in the reconstruction was 200.

In the reconstruction procedure, it is possible to include a rescaling step in order to ensure that the global energetic content of the fluctuations reaches a certain level. We checked that a local definition of the variance was sufficient to define appropriate rescaling factors.

The rescaling step was found to improve results considerably. This can be understood by the fact that our reconstruction is based on projection, so that even if the amplitudes were perfectly estimated, the reconstructed field will lack some energy. Moreover, it was found most important to rescale the momentum.

The implementation of the boundary condition depends on the definition of the Riemann invariants for the different variables. The adequacy of the boundary condition was found to depend strongly on the choice of the Riemann invariants. The reasons for this were not made entirely clear.

Although good results over the simulation time considered were obtained for both subsonic and supersonic case, the agreement between the reference channel and the reduced computational domain was better in the supersonic case. We speculate that this could be due to the fact that in the supersonic the flow velocity on the boundary plane is already supersonic, so that all information is propagated downstream. The computational time saving is achieved with a gain around 2.5 in reduced channel.

There, we can see POD is costly, but a significant cost reduction is achieved.

Perspectives

The work presented here remains largely exploratory. We believe that two main directions of investigation are necessary at this stage:

1. The validity of the approach should be established over longer simulation times that is several hundreds of time units. It was not possible to carry out this study in the present thesis as our focus was to determine all the different parameter choices that could be made.

2.

A significant effort should be made to understand the results associated with some parameter choices. Again, we did not have the time to analyze the results in detail. This should help establish the robustness of the procedure and define its best operating conditions. These two directions are complementary, as establishing the validity of the conditions over long periods may require adaptation of the procedure, such as the introduction of additional filtering techniques. This requires a thorough understanding of the effect of these techniques on the flow computation.

Once these two issues are addressed, optimizing the reconstruction procedure is essential in order to capitalize on the reduction of the spatial mesh and the increase in the computational time step.

It will then be interesting to consider higher Reynolds numbers, in which an inertial region can be more clearly defined. Results for incompressible flow suggests that increasing the Reynolds number tends to improve the procedure, as the separation between the wall region and the outer layer (or core region, in the case of a channel) increases. Increasing the Reynolds number will also allow us to test higher boundary planes, and therefore to obtain further reduction.

And then, we could obtain: The extension to compressible flows of the algorithm proposed by Deschamps is here employed to update the driving force at each time step: Direct Numerical Simulation of turbulent boundary layers requires a high spatial resolution specially near the wall to capture the different length scales inherent of turbulence. To reduce the computing cost related to the high spatial resolution, we investigate techniques that bypass the wall region by using an artificial boundary condition at the bound of a reduced computational domain. A reconstruction strategy based on POD is implemented in order to ge-nerate synthetic boundary conditions such that the correct spatial characteristics and dynamics of the turbulence are recovered in the core region of the channel. The boundary condition was tested at two different heights. Statistics in the reduced channel were compared with those of a reference case. A good agreement is observed for the mean velocity profile and the turbulent shear stress. The procedure was found to be robust with the boundary height and results in a drastic computational cost reduction.

f n+1 1 = f n 1 + ∆t L y L z [α(Q n+1 -Q 0 ) + β(Q n m -Q 0 m )], (

Universit é Paris-Saclay

Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

1 w 1 . 1

 111 Wall turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Simulating turbulence: DNS, LES and RANS . . . . . . . . . . . . . . . 1.3 Wall models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Compressibility effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Approximate boundary conditions . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Slip boundary condition . . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Control-based strategies . . . . . . . . . . . . . . . . . . . . . . . 1.5.3 Synthetic wall boundary conditions . . . . . . . . . . . . . . . . . 1.6 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 DNS of the compressible channel flow 2.1 The governing equations for compressible flows. . . . . . . . . . . . . . . 2.2 The numerical approach for solving the governing equations . . . . . . . 2.3 Numerical configuration : the compressible turbulent channel flow . . . . 2.3.1 Boundary and initial conditions: . . . . . . . . . . . . . . . . . . . 2.3.2 Treatment of the periodic boundary condition in the streamwise direction. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Results of the Direct Numerical Simulation . . . . . . . . . . . . . . . . . 2.4.1 Statistical treatments of simulation data . . . . . . . . . . . . . . vii CONTENTS 2.4.2 DNS results of the subsonic channel flow . . . . . . . . . . . . . . 2.4.3 DNS results for the supersonic channel flow. . . . . . . . . . . . . 2.4.4 Comparison between results of subsonic flow and supersonic flow .3 Reconstruction of synthetic boundary conditions 3.1 Proper Orthogonal Decomposition . . . . . . . . . . . . . . . . . . . . . 3.1.1 Direct Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Method of snapshots . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.5 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Linear Stochastic Estimation . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 General definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Application . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Reconstruction method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Inlet Synthetic boundary conditions: rescaling approaches . . . . 3.3.2 Inlet synthetic boundary conditions: Structure-based decompositions 3.3.2.1 The synthetic eddy method (SEM) . . . . . . . . . . . . 3.3.2.2 POD-based reconstructions . . . . . . . . . . . . . . . . 3.3.3 Wall Synthetic boundary conditions . . . . . . . . . . . . . . . . . 3.3.3.1 Current approaches . . . . . . . . . . . . . . . . . . . . . 3.3.3.2 The reconstruction procedure . . . . . . . . . . . . . . . 3.3.3.3 Step 3: Rescaling . . . . . . . . . . . . . . . . . . . . . . 3.3.3.4 Step 4: Implementation of the reconstruction . . . . . . 3.3.3.5 First test: Reduced simulation using reference flow fields as boundary conditions . . . . . . . . . . . . . . . . . . . Velocity in z direction H Half height, reference length U 0 Initial mean stremwise velocity, Reference velocity T 0 Temperature on the wall, Reference velocity ρ 0 Temperature on the wall, Reference density U 0 Reference velocity Velocity expressed with friction velocity ts Threshold α n (t), n = 1, 2 Number of mode α n (t), t = 1, 2 Number of sample xiii CONTENTS α Temporal coefficient φ n (x, y, z), x, y, z Position in spatial mode φ n (x, y, z), n = 1, 2 Number of mode in spatial mode Main velocity profile from Orlandi et al [58]. . . . . . . . . . . . . . . . . 1.2 Horseshoe vortex in top (Figure from Theodorsen [80]), and streaks in bottom (Figure from Wang et al [85]). . . . . . . . . . . . . . . . . . . . 1.3 Classification of different methods for fluid simulation [70]. . . . . . . . . 1.4 Iso-surfaces of the second invariant of the velocity gradient tensor in the near wall region. The iso-surfaces are coloured with the local streamwise velocity [55]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Mean velocity profile U + in viscous wall units [60] at Re θ = 1250, 1840, 2530, 2990, and 3310 (from bottom to top). Profiles are shifted by U + = 3 along the ordinate for increasing Re θ . Solid lines, present WMLES; dashed line, log-law (κ = 0.41, G = 5.2); +, matching locations; , DNS of Sayadi et al. [71] at Re θ = 1250; •, DNS of Wu and Moin [86] at Re θ = 1840; experiment of Osterlund [59] at Re θ = 2530. . . . . . . . . . . 1.6 Plane channel flow [12]. Left: mean velocity profiles in wall coordinates: a)Re τ = 200 : -, DNS [38]; case 1: +, two-layer model, ,shifted model; and b) Re τ = 2000:-, resolved LES [63]; case 2: +, two-layer model, , shifted model; case 3: . . ., two-layer model. Right: turbulence intensities in outer coordinates for Re τ = 200:-,DNS [38]; case 1: +, two-layer model; and , shift model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xv 1.7 Results of interface validation [72]. Above: Axial velocity profiles. Below: Axial velocity fluctuations. Circles: experiments. Solid lines: LES with inflow from experimental data. Dashed lines: Integrated RANS-LES, RANS with inflow from experimental data, LES inflow derived from simultaneously running RANS solver. . . . . . . . . . . . . . . . . . . . . 1.8 Comparison between a full incompressible channel flow simulation and a reduced domain with a POD-based synthetic boundary conditions using phase estimation. [67]. Left) Mean velocity profile; Right) Turbulent intensities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Sketch of the computational domain of the channel flow configuration. . . 2.2 Mean streamwise velocity profile in the wall normal direction: influence of the Reynolds number on various regions (from [69]). . . . . . . . . . . 2.3 Profiles of the mean streamwise component of the velocity (u + ) versus the normal to the wall direction (z + ) in the upper and the lower half part of the subsonic channel flow. . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Profiles of the Reynolds stress components in the subsonic channel flow, normalized by the friction velocity u τ : rms values of the velocity components (< u i u i >), on the left, and the turbulent shear stress (< u w >), on the right. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1. 5 )

 5 [START_REF] Lechner | Turbulent supersonic channel flow[END_REF] (green line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi LIST OF FIGURES 2.8 Profiles of the mean quantities for the two Mach numbers: M a = 0.5 in black, and M a = 1.5 in red. Top left: mean density profile versus z. Top right, mean dynamic viscosity profile versus z. Bottom, mean streamwise velocity component normalized by the friction velocity (u + ) versus z + . . . 2.9 Effect of the compressibility on channel flow results coming from Foysi et al. [27]: On the left, mean density profiles (symbols) and mean viscosity profiles (lines). On the right, profiles of the Van Driest-transformed mean streamwise velocity. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Model of channel for symmetry. . . . . . . . . . . . . . . . . . . . . . . . 3.2 Schematic illustrating the POD convergence issue: a field that is not part of the original set of snapshots may differ from its projection onto the POD basis. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Left: POD spectrum for density, momentum, and internal energy. Right: eigenfunction φ 2 ρu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 POD dominant mode (n = 1) in Fourier space for the wavenumber pair (k x = 0, k y = 4); Top left: Dominant mode (n = 1) for φ 1 kxky (ρu). Top right:Dominant mode φ 1 kxky (ρv) . Bottom left: Dominant mode of φ 1 kxky (ρw). Bottom right: POD spectrum for velocity ρu, ρv, ρw. . . . . . 3.5 Model for region h 0 < z ≤ 1. . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Distribution of R np (h 0 = 0.1). . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Reconstruction error between the real field and its projection for a field in the POD snapshot set -Primitive variables in POD reconstruction and boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Reconstruction error of different physical quantities for the same field as in Figure 3.7-Primitive variables in POD reconstruction and boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Tent like function defined by f j (r) = 1r/L , where y = 0, z = 0, L = 2, U 0 t = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Two point correlation R ii (y, y * ) of the original velocity fields (left) and synthetic velocity fields (right): positive increment between contours of 0.5, -: negative increment between contours of -0.

4. 3 Figure 4 . 2 (

 342 Figure 4.2 (reconstruction with POD based on primitive variables, and average taken over 6 time units). . . . . . . . . . . . . . . . . . . . . . . 4.6 Influence of the definition of the mean in the factor f Q . Top Left: rms of u. Top right: rms of v. Bottom left: rms of w. Bottom right: shear stress < u w > -Legends as in Figure 4.2. . . . . . . . . . . . . . . . . . . 4.7 Mean profiles in the bottom half channel for 6 time units for POD based on conservative variables. Top left: density. Top right: streamwise velocity u. Bottom left: spanwise velocity v. Bottom right: wall-normal velocity w. 4.8 Turbulent intensities in the bottom half channel. Left: rms value of u. Right: rms value of v -Legends as in Figure 4.7. . . . . . . . . . . . . . . 4.9 Turbulent intensities in the bottom half channel. Left: rms value of w. Right: shear stress < u w > -Legends as in Figure 4.7. . . . . . . . . . . 4.10 Mean profile computed over 6 time units with POD based on conservative variables. Left: density ρ. Right: streamwise velocity u. . . . . . . . . . . 4.11 Mean profile computed over 6 time units with POD based on conservative variables. Left: spanwise velocity v. Right: wall-normal velocity w. . . . 4.12 Turbulent intensities computed over 6 time units with POD based on conservative variables. Top left: rms value of u. Top right: rms value of v. Bottom left: rms value of w. Bottom right: Reynolds stress < u w >. 4.13 Mean profiles of the bottom channel for 6 time units with DNS, conservative variables (rescaling step with local mean value), and primitive variables (rescaling step with local mean value). Top left: mean profile of density. Top right: mean profile of u. Bottom left: mean profile of v. Bottom right: mean profile of w. . . . . . . . . . . . . . . . . . . . . . . 4.14 Comparison between the reconstructions. Turbulent intensities. Left: rms value of u. Right: rms value of v -Legends as in Figure 4.13. . . . . . . .

2. 1

 1 Maximum values and their locations in the profiles of the rms values of the velocity components from our DNS results on the subsonic channel flow compared with the reference values from [38]. . . . . . . . . . . . . . 2.2 Time averaged results from [20]. † corresponds to the incompressible DNS results from [38]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Maximal values and their locations in the profiles of the rms values of the velocity components from the present DNS results on the supersonic channel flow compared with values obtained from [44]. . . . . . . . . . . 3.1 Errors e Q for field and its projection measured on the plane h 0 = 0.1. . .

3. 2 5 . 2

 252 Data basis for calculation in channel flow. . . . . . . . . . . . . . . . . . 4.1 Relative errors for simulations with primitive variables in reduced region 0.3 < z ≤ 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Relative errors for reduced-domain simulations with POD based on conservative variables in region 0.3 < z ≤ 1. . . . . . . . . . . . . . . . . . . 4.3 Relative errors for simulations with primitive variables in reduced channel 0.3 < z ≤ 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Relative errors for simulations with conservative variables in reduced channel 0.3 < z ≤ 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Statistical errors for unrescaled boundary conditions. . . . . . . . . . . . xxviii Statistical errors for rescaled boundary conditions. The words between parentheses indicate the type of scaling that was performed, depending on whether all components (global) or only the momentum were rescaled. 5.3 Comparison between different choices of Riemann invariants for the boundary conditions (eigenvalue expressed in bottom boundary plane). . . . . . 5.4 Computational requirements of the reference and reduced-channel simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.1 Relative error for simulations in reduced channel 0.3 < z ≤ 1.7. . . . . . . 6.2 Relative error for simulations in reduced channel 0.3 < z ≤ 1.7 obtained from 70 modes (Comparison with subsonic flow in Chapter 5). . . . . . . 6.3 Relative error for simulations in reduced channel 0.3 < z ≤ 1.7 obtained for 60 POD modes of Part B in Figure 6.16 (Comparison with subsonic flow in Chapter 5). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4 Computational time saving for supersonic flow. . . . . . . . . . . . . . .

Figure 1

 1

Figure 1 . 1 :

 11 Figure 1.1: Main velocity profile from Orlandi et al [58].

Figure 1 . 2 :

 12 Figure 1.2: Horseshoe vortex in top (Figure from Theodorsen [80]), and streaks in bottom (Figure from Wang et al [85]).

Figure 1 . 3 :

 13 Figure 1.3: Classification of different methods for fluid simulation [70].

Figure 1 . 5 :

 15 Figure 1.5: Mean velocity profile U + in viscous wall units [60] at Re θ = 1250, 1840, 2530, 2990, and 3310 (from bottom to top). Profiles are shifted by U + = 3 along the ordinate for increasing Re θ . Solid lines, present WMLES; dashed line, log-law (κ = 0.41, G = 5.2); +, matching locations; , DNS of Sayadi et al. [71] at Re θ = 1250; •, DNS of Wu and Moin [86] at Re θ = 1840; experiment of Osterlund [59] at Re θ = 2530.

Bocquet

  et al. [12] developed an extension of the TBL (Thin Boundary Layer) model for compressible LES of plane channel flow. The mean wall fluxes, primitive variable profiles and turbulent fluctuations were found to be in good agreement with direct numerical simulation (see Figure 1.6) .

Figure 1 .

 1 Figure 1.6: Plane channel flow [12]. Left: mean velocity profiles in wall coordinates: a)Re τ = 200 : -, DNS [38]; case 1: +, two-layer model, ,shifted model; and b) Re τ = 2000:-, resolved LES [63]; case 2: +, two-layer model, , shifted model; case 3: . . ., twolayer model. Right: turbulence intensities in outer coordinates for Re τ = 200:-,DNS [38]; case 1: +, two-layer model; and , shift model.

Figure 1 . 7 :

 17 Figure 1.7: Results of interface validation [72]. Above: Axial velocity profiles. Below: Axial velocity fluctuations. Circles: experiments. Solid lines: LES with inflow from experimental data. Dashed lines: Integrated RANS-LES, RANS with inflow from experimental data, LES inflow derived from simultaneously running RANS solver.

[

  66] for incompressible flows. It relies on the construction of a boundary condition based on Proper Orthogonal Decomposition (POD), which allows the representation of the condition as a superposition of structures which are extracted from the spatial autocorrelation tensor.

Figure 1 . 8 :

 18 Figure 1.8: Comparison between a full incompressible channel flow simulation and a reduced domain with a POD-based synthetic boundary conditions using phase estimation. [67]. Left) Mean velocity profile; Right) Turbulent intensities.

1

 1 The governing equations for compressible flows.

Figure 2 . 1 :

 21 Figure 2.1: Sketch of the computational domain of the channel flow configuration.

channel flows are α = 2

 2 δt , and β = -0.2 δt . Same values of α and β have also been used in the compressible channel flow. Detailed calculations to obtain the macroscopic force are reported in Appendix B for 3D flow. The external force calculation has been introduced

  Logarithmic and outer layer regions exist at sufficiently high Reynolds number. Their extent largely depends on the Reynolds number, as seen in Figure2.2.

Figure 2 . 2 :

 22 Figure 2.2: Mean streamwise velocity profile in the wall normal direction: influence of the Reynolds number on various regions (from [69]).

Figure 2 . 3 :

 23 Figure 2.3: Profiles of the mean streamwise component of the velocity (u + ) versus the normal to the wall direction (z + ) in the upper and the lower half part of the subsonic channel flow.

Figure 2 . 4 :

 24 Figure 2.4: Profiles of the Reynolds stress components in the subsonic channel flow, normalized by the friction velocity u τ : rms values of the velocity components (< u i u i >), on the left, and the turbulent shear stress (< u w >), on the right.

1 :

 1 Maximum values and their locations in the profiles of the rms values of the velocity components from our DNS results on the subsonic channel flow compared with the reference values from[START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]. turbulent shear stress in the center of the channel (Figure2.4) since the viscous diffusion is negligible in this region.

Figure 2 . 5 :

 25 Figure 2.5: Profiles of the Reynolds stress components in the incompressible channel flow coming from[START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF]: rms values of the velocity components (< u i u i >), on the left, and the turbulent shear stress (< u w >), on the right.

Figure 2 . 6 :

 26 Figure 2.6: Mean profiles in the supersonic channel flow: mean density (top left), temperature (top right) and mean streamwise velocity component (bottom).

  2.7 0.98 0.75 Location of max value in present DNS 0.1 0.22 0.35 Value on the center in present DNS 0.9 0.73 0.72 Table 2.3: Maximal values and their locations in the profiles of the rms values of the velocity components from the present DNS results on the supersonic channel flow compared with values obtained from [44].

Figure 2 .

 2 Figure 2.7: Profiles of the Reynolds stress components in the supersonic channel flow (< u i u i > and -< u w >), normalized by the friction velocity u τ : present DNS results, on the left, and, on the right, reference profiles from incompressible DNS [38] (red line), and from supersonic DNS (M a = 1.5) [44] (green line).

Figure 2 . 9 .

 29 who explored the compressibility effects on turbulence in channel flow by studying four configurations: M a = 0.3 with Re τ = 181; M a = 1.5 with Re τ = 221; M a = 3.0 with Re τ = 556; and M a = 3.0 with Re τ = 1030. The two first configurations approximately correspond to the present simulations. The present DNS results favorably compare with results coming from the literature and the effects of compressibility are very well reproduced.

Figure 2 . 8 :

 28 Figure 2.8: Profiles of the mean quantities for the two Mach numbers: M a = 0.5 in black, and M a = 1.5 in red. Top left: mean density profile versus z. Top right, mean dynamic viscosity profile versus z. Bottom, mean streamwise velocity component normalized by the friction velocity (u + ) versus z + .

Figure 2 . 9 :

 29 Figure 2.9: Effect of the compressibility on channel flow results coming from Foysi et al. [27]: On the left, mean density profiles (symbols) and mean viscosity profiles (lines). On the right, profiles of the Van Driest-transformed mean streamwise velocity.

  satisfy < a n (t) >= 0, (3.5) and < a n (t)a m (t) >= δ mn λ n . (3.6)

Figure 3 . 1 :

 31 Figure 3.1: Model of channel for symmetry.

Figure 3 . 2 :

 32 Figure 3.2: Schematic illustrating the POD convergence issue: a field that is not part of the original set of snapshots may differ from its projection onto the POD basis.

  Figure 3.3 for the method of snapshots. The eigenfunction φ n=2 ρu is shown also in Figure 3.3.

Figure 3 . 3 :

 33 Figure 3.3: Left: POD spectrum for density, momentum, and internal energy. Right: eigenfunction φ 2 ρu .

Figure 3 .

 3 Figure 3.3 shows the dominant fluctuation (n = 2) in the right (the mean corresponds to n = 1).

Figure 3 .

 3 Figure 3.4 shows the dominant mode in Fourier space for the wavenumbers (l = 0, k = 4).

Figure 3 . 4 :

 34 Figure 3.4: POD dominant mode (n = 1) in Fourier space for the wavenumber pair (k x = 0, k y = 4); Top left: Dominant mode (n = 1) for φ 1 kxky (ρu). Top right:Dominant mode φ 1 kxky (ρv) . Bottom left: Dominant mode of φ 1 kxky (ρw). Bottom right: POD spectrum for velocity ρu, ρv, ρw.

Figure 3 . 5 :

 35 Figure 3.5: Model for region h 0 < z ≤ 1.

Figure 3 . 6 :

 36 Figure 3.6: Distribution of R np (h 0 = 0.1).

. 31 )Figure 3 .

 313 Figure 3.7 shows the reconstruction error between a field at the plane h 0 and its reconstruction. The error is round-off error, as expected.

Figure 3 . 7 :

 37 Figure 3.7: Reconstruction error between the real field and its projection for a field in the POD snapshot set -Primitive variables in POD reconstruction and boundary conditions.

Figure 3 . 8 :

 38 Figure 3.8: Reconstruction error of different physical quantities for the same field as in Figure 3.7-Primitive variables in POD reconstruction and boundary conditions.

Figure 3 .

 3 Figure 3.8 shows the reconstruction between the real field at h 0 = 0.1 and the field estimated from the data z > 0.1, there, all 200 modes are used in reconstruction of ρ and E, 128 modes are used in reconstruction of [u, v, w].

2 Table 3 . 1 :

 231 Errors e Q for field and its projection measured on the plane h 0 = 0.1. Comparing Figure 3.7 and Figure 3.8 allows us to evaluate the performance of the esti-mation method.

Figure 3 . 10 :

 310 Figure 3.10: Two point correlation R ii (y, y * ) of the original velocity fields (left) and synthetic velocity fields (right): positive increment between contours of 0.5, -: negative increment between contours of -0.1 (from [61]).

. 42 )

 42 There, we could use Euler equation to describle the wave propagation. The Riemann invariant is defined as dω = R -1 s dQ c . If we replace the variables in Euler equations with Riemann invariants, then we have: ∂ω ∂t + Λ • ∇ω = 0 =⇒ dω dt = 0 along lines λ i = dz dt . Where ω are constant along characteristic lines and propagate at velocity λ i as shown in Figure 3.11.

Figure 3 . 11 :

 311 Figure 3.11: Wave propagation.

  ρ)w, T ] Which are calculated from POD reconstruction. The boundary conditions for these fields must fulfill the compatibility relationships based on the Riemann invariants. Considering a boundary normal to the z direction, the Navier-Stokes equations are projected onto the eigenvectors of the Jacobian of the Euler flux in the z -direction. Depending on the sign of the eigenvalues for this Euler flux Jacobian, either the quantity is prescribed when the eigenvalue is positive, or equations of the Riemann invariant are solved by up-winding the Euler flux discretization in the normal to the boundary plane z-direction when the eigenvalue is negative. More precisely, for each point (x, y) of the boundary plane z = z 0 :

  always negative in subsonic flow (C is speed of sound, which is much larger than u • n), so the density ρ is always determined by Riemann invariant; u • n+ n C is always positive, so the wall normal velocity w is determined by POD reconstruction;the others three terms depend on u • n, if u • n ≥ 0, they are determined by POD reconstruction, otherwise, they are determined from propagation of the quantities within the domain along the characteristics.

Figure 3 . 13 :

 313 Figure 3.13: Boundary condition on the boundary planes.

Figure 3 .

 3 Figure 3.14 presents the flow configuration and the computational domain. The stream-

Figure 3 .

 3 Figure 3.15: Field of Density on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel.

Figure 3 .

 3 Figure 3.16: Field of u on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel.

Figure 3 . 17 :

 317 Figure 3.17: Field of v on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel.

Figure 3 . 18 :

 318 Figure 3.18: Field of w on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel.

Figure 3 . 19 :

 319 Figure 3.19: Field of total energy E on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel.

Figures 3 .

 3 Figures 3.15 to Figure 3.19 show the instantaneous field in the full channel (left side) and in the reduced channel (right side). No difference was observed for all quantities that validate the boundary conditions and the procedure. It now remains to test the reconstructed boundary condition in the reduced channel. This is the object of the next chapter.

Figure 4 . 1 .

 41 Figure 4.1.The new boundary of the domain is a plane cutting across the wall layer. As seen in Chapter 2, the numerical method requires that physical quantities are defined not only on the boundary but also on planes of ghost cells outside the boundary in order to compute the gradients with accuracy. Four extra planes are required, which means that five planes have to be stored for the boundary condition. The simulation in the reduced domain was carried out for a time of ∆t = 6H/U 0 and the turbulent statistics compared with those of the full channel. As seen in Chapter 3, the reconstruction method is based on POD modes. For the results presented here, the method of snapshots was applied to 70 fields, separated by a time ∆t = 6H/U 0 . In order to evaluate the robustness of the method, snapshots of primitive variables and conservative variables were considered, and two different heights were considered: h + 0 = 18 and h + 0 = 54. The results of full channel and reduced channel are for the same statistic duration in this chapter.

Figure 4 . 1 :

 41 Figure 4.1: Model of reduced channel.

4. 1

 1 Results at height h + 0 = 18 (h 0 = 0.1) with primitive variables We can choose to apply the POD to primitive variables [ρ, u, v, w, E] or conservative variables [ρ, ρu, ρv, ρw, ρE].

  Figures 4.2-4.4 for the bottom half-channel. The red lines correspond to DNS, the purple lines correspond to the reduced channel simulation with a synthetic boundary condition, while the black lines correspond to the rescaled boundary condition.

4. 1 Figure 4 . 2 :

 142 Figure 4.2: Mean profiles averaged over time and expressed with horizontal directions in the bottom half channel over 6 time units for POD in primitive variables. Top left: mean profile of density; top right: mean profile of u; bottom left: mean profile of v ; bottom right: mean profile of w.

Figure 4 . 3 : 67 Figure 4 . 4 :

 436744 Figure 4.3: Turbulent intensities in the bottom half channel. Left: rms value of u. Right: rms value of v -Legends as in Figure 4.2.

68 Figure 4 . 5 :

 6845 Figure 4.5: Influence of the definition of the mean in the factor f Q on the mean velocity profile. Top left: density. Top right: streamwise velocity. Bottom left: Spanwise velocity. Bottom right: wall-normal velocity -conditions as in Figure 4.2 (reconstruction with POD based on primitive variables, and average taken over 6 time units).

4. 1 69 Figure 4 . 6 :

 16946 Figure 4.6: Influence of the definition of the mean in the factor f Q . Top Left: rms of u. Top right: rms of v. Bottom left: rms of w. Bottom right: shear stress < u w > -Legends as in Figure 4.2.

4. 2 71 Figure 4 . 7 :

 27147 Figure 4.7: Mean profiles in the bottom half channel for 6 time units for POD based on conservative variables. Top left: density. Top right: streamwise velocity u. Bottom left: spanwise velocity v. Bottom right: wall-normal velocity w.

Figure 4 . 8 :

 48 Figure 4.8: Turbulent intensities in the bottom half channel. Left: rms value of u. Right: rms value of v -Legends as in Figure 4.7.

Figure 4 . 9 :

 49 Figure 4.9: Turbulent intensities in the bottom half channel. Left: rms value of w. Right: shear stress < u w > -Legends as in Figure 4.7.

  Figures in 4.8-4.9 show that the turbulence intensities match well those of the reference simulation especially in the channel center. The effect of rescaling step is beneficial, as the rms value of the wall-normal velocity w matches better that of the reference simulation near boundary planes.

Figures 4. 10

 10 Figures 4.10, 4.11,[START_REF] Akselvoll | Large-Eddy Simulation of turbulent confined coannular jets[END_REF]12 show that there is no significant difference between the two

  Figures 4.10, 4.11,[START_REF] Akselvoll | Large-Eddy Simulation of turbulent confined coannular jets[END_REF]12 show that there is no significant difference between the two rescaling methods. Results show that the results are relatively insensitive to the choice of mean value in the rescaling factor.

4. 2 73 Figure 4 . 10 :

 273410 Figure 4.10: Mean profile computed over 6 time units with POD based on conservative variables. Left: density ρ. Right: streamwise velocity u.

Figure 4 . 11 :Figure 4 .

 4114 Figure 4.11: Mean profile computed over 6 time units with POD based on conservative variables. Left: spanwise velocity v. Right: wall-normal velocity w.

Figure 4 . 13 :

 413 Figure 4.13: Mean profiles of the bottom channel for 6 time units with DNS, conservative variables (rescaling step with local mean value), and primitive variables (rescaling step with local mean value). Top left: mean profile of density. Top right: mean profile of u. Bottom left: mean profile of v. Bottom right: mean profile of w.

Figure 4 .

 4 Figure 4.14: Comparison between the reconstructions. Turbulent intensities. Left: rms value of u. Right: rms value of v -Legends as in Figure 4.13.

Figure 4 .

 4 Figure 4.15: Comparison between the reconstructions. Turbulent intensities. Left: rms value of w. Right: shear stress< u w > -Legends as in Figure 4.13.

  from simulation of conservative variables and primitive variables with rescaling step of local mean values. To enhance the comparison, we only present the region 0.1 < z < 2 in the full channel. The comparison is presented in Figures 4.16

to 4. 20 .

 20 The left side corresponds to the full channel, the middle side corresponds to the reduced channel with conservative variables, and the right side corresponds to the reduced channel with primitive variables.

Figure 4 .

 4 Figure 4.16: Snapshot of density on middle vertical plane at 6 time units. Left: in full channel. Middle: in reduced channel with conservative variables. Right: in reduced channel with primitive variables.

Figure 4 . 17 :

 417 Figure 4.17: Snapshot of u on middle vertical plane at 6 time units. Left: in full channel. Middle: in reduced channel with conservative variables. Right: in reduced channel with primitive variables.

Figure 4 . 18 :

 418 Figure 4.18: Snapshot of v on middle vertical plane at 6 time units. Left: in full channel. Middle: in reduced channel with conservative variables. Right: in reduced channel with primitive variables.

Figure 4 . 19 :

 419 Figure 4.19: Snapshot of w on middle vertical plane at 6 time units Left: in full channel. Right: in reduced channel.

4. 4 78 Figure 4 . 20 :

 478420 Figure 4.20: Snapshot of E on middle vertical plane at 6 time units. Left: in full channel. Middle: in reduced channel with conservative variables. Right: in reduced channel with primitive variables.

4. 4 79 4. 4 . 1

 47941 Results at height h + 0 = 54 (h 0 = 0.3) Results at height h + 0 = 54 with primitive variables

Figure 4 . 21 :Figure 4 . 22 :

 421422 Figure 4.21: Mean profiles with POD based on primitive variables over a time of 6 time units. Top left: density. Top right: streamwise velocity u. Bottom left: spanwise velocity v. Bottom right: wall-normal velocity w.

Figures 4 .

 4 Figures 4.21-4.22 show that the turbulent intensities, shear stress and mean profiles in the center region are well recovered in the reduced channel simulation. A good agreement

  -4.25).

4. 4 81 Figure 4 . 23 :

 481423 Figure 4.23: Influence of the rescaling factor definition on mean profiles. Top left: mean profile of density. Top right: mean profile of u. Bottom left: mean profile of v. Bottom right: mean profile of w.

Figure 4 . 24 : 82 Figure 4 .

 424824 Figure 4.24: Influence of the rescaling factor definition on turbulent intensities. Left: rms of streamwise velocity u. Right: rms of spanwise velocity v.

Figures 4. 23

 23 Figures 4.23-4.25 indicate that the definition of the mean in the rescaling factor does not

  Figures 4.23-4.25 indicate that the definition of the mean in the rescaling factor does not affect the results.

4. 4

 4 Results at height h + 0 = 54 (h 0 = 0.3) 83 4.4.2 Results at height h + 0 = 54 for POD based on conservative variables In this section, the POD reconstruction was applied to conservative variables [ρ] [ρu, ρv, ρw], [ρE]. All other conditions are the same.

Figure 4 . 26 :

 426 Figure 4.26: Mean profiles with POD based on conservative variables. Left: density. Right: streamwise velocity u.

Figure 4 . 27 :

 427 Figure 4.27: Mean profiles with POD based on conservative variables. Left: spanwise velocity v. Right: wall-normal velocity w.

Figures 4. 27

 27 Figures 4.27and 4.28 show that rescaling step improves the statistics of the reduced

  Figures 4.27and 4.28 show that rescaling step improves the statistics of the reduced channel except the profile of w, because it contains more small scales which is more difficult to predict.

4. 4 84 Figure 4 . 28 :

 484428 Figure 4.28: Turbulent intensities computed over a time of 6 time units with POD in conservative variables. Top left: rms value of u. Top right: rms value of v. Bottom left: rms value of w. Bottom right: shear stress < u w >.

Figure 4 . 29 :Figure 4 . 30 :

 429430 Figure 4.29: Influence of the definition of the mean in the factor f Q on the mean velocity profile. Left: density ρ. Right: streamwise velocity u.

Figure 4 . 31 :

 431 Figure 4.31: Influence of the definition of the mean in the factor f Q . Top left: rms of u. Top right: rms of v. Bottom left: rms of w. Bottom right: shear stress < u w >conditions as in Figure 4.26 (reconstruction with POD based on conservative variables, and average taken over 6 time units).

  -4.31. Once 4.4 Results at height h + 0 = 54 (h 0 = 0.3) 86

Figure 4 . 32 :

 432 Figure 4.32: Mean profile in the bottom half channel for 6 time units with primitive variables (DNS, rescaling step with local mean value, rescaling step with global mean value). Left: Profile of v. Right: Profile of w.

Figure 5 . 1 :

 51 Figure 5.1: Left: Dominant eigenvalue λ 1 kxky for [ρu, ρv, ρw] as a function of the wavenuber. Right: Eigenvalues λ n 00 of [ρu, ρv, ρw].

Figure 5 . 2 :

 52 Figure 5.2: Mean profiles in reduced channel. Left: Temperature. Right: Density.

Figure 5 . 3 :

 53 Figure 5.3: Mean profiles in reduced channel. Left: streamwise velocity u in bottom half channel. Right: streamwise velocity u in top half channel.

Figure 5 . 4 :

 54 Figure 5.4: Turbulent intensity value for the different POD methods. Left: rms streamwise velocity u. Right: rms spanwise velocity v.

Figure 5 . 5 :

 55 Figure 5.5: Turbulent intensity value for the different POD methods. Left: rms wallnormal velocity w. Right: shear stress < u w >.

Figure 5 . 6 :

 56 Figure 5.6: Instantaneous density at z = 0.2 (z + = 36) after 6 time units. Left: full channel. Right: reduced channel.

Figure 5 . 7 :

 57 Figure 5.7: Instantaneous streamwise velocity field u at z = 0.2 (z + = 36) after 6 time units. Left: full channel. Right: reduced channel.

Figure 5 .

 5 Figure 5.2 shows discrepancies for both the density and the temperature, which is relatively surprising. Figures 5.4-5.5 present turbulent intensity values and the shear stress.

5. 2 18 :

 218 Reduced simulation at h + 0 = Definition of POD variables 94 5.2 Reduced simulation at h + 0 = 18: Definition of POD variables 5.2.1 Proper Orthogonal Decomposition Synthetic boundary conditions were defined for each channel plane. The channel is therefore divided into two independent boundary layers -each corresponding to one half of the channel. POD is then performed over each half-channel. To compute the temporal autocorrelation tensor, 100 instantaneous flow fields corresponding to the full channel are considered. Each field was split into two parts, each of which corresponds to a half-channel, therefore yielding two snapshots. POD therefore yield 200 modes out of which only 128 are used in the reconstruction due to the memory limitation for the spatial modes of momentum. This number of modes capture about 80% of the fluctuating turbulent kinetic energy, which means that by construction, the synthetic boundary condition will lack energy. However it will be compensated by the rescaling procedure.

Figure 5 . 8 : 18 : 2 is the internal energy. 4 .

 581824 Figure 5.8: POD spectrum for the following variables: density, momentum, and ρe.

Figure 5 . 9 :

 59 Figure 5.9: Reduced simulation domain.

5. 2 Figure 5 . 10 :

 2510 Figure 5.10: Mean profiles -Top left: streamwise velocity u in bottom half channel. Top right: streamwise velocity u in top half channel. Bottom: density. Red line is for DNS. Violet line if for ([ρ], [u, v, w], [E]). Black line is for ([ρ], [ρu, ρv, ρw], [ρE]). Blue line is for ([ρ], [ρu, ρv, ρw], [ρe]). Green line is for ([ρ, ρu, ρv, ρw, ρE]).

Figure 5 . 11 :

 511 Figure 5.11: Turbulent intensity values for the different POD methods. Left: rms streamwise velocity u. Right: rms spanwise velocity v. -Legends as in Figure 5.10.

Figure 5 .

 5 Figure 5.12: Turbulent intensity value and shear stress for the different POD methods. Left: rms wall-normal velocity w. Right: shear stress < u w >. -Legends as in Figure 5.10.

  -5.15. Statistics were integrated over 36 time units except for the method based on primitive variables, which diverged after 19 time units, because we can not control mass flux well in reduced 5.2 Reduced simulation at h + 0 = 18: Definition of POD variables 98 channel with primitive variables.

Figure 5 . 13 :

 513 Figure 5.13: Mean profiles -Top left: Profile of streamwise velocity u in bottom half channel. Top right: profile of streamwise velocity u in top half channel. Bottom: profile of density. Red line is for DNS. Violet line is for ([ρ], [ρu, ρv, ρw], [ρE]). Black line is for ([ρ], [ρu, ρv, ρw], [ρe]with momentum rescaling). Blue line is for ([ρ], [ρu, ρv, ρw], [ρe]with global rescaling). Green line is for ([ρ, ρu, ρv, ρw, ρE]).

Figure 5 .

 5 Figure 5.14: Turbulent intensity values and shear stresses for the different POD methods. Left: rms streamwise velocity u. Right: rms spanwise velocity v-Legends as in Figure 5.13.

5. 2 Figure 5 .

 25 Figure 5.15: Turbulent intensity values and shear stresses for the different POD methods. Left: rms wall-normal velocity w. Right: shear stress < u w >-Legends as in Figure 5.13.

Figure 5 . 18 :

 518 Figure 5.13 shows the mean density. Large differences can be observed between the different decompositions. The best profiles correspond to the rescaled decomposition based on [ρ], [ρu, ρv, ρw], [ρe] -independently of the type of rescaling used. The worst case corresponds to a unique decomposition for all variables. It is possible to interpret this result as a lack of degrees of freedom in the reconstruction. Figures 5.14-5.15 present the turbulent intensities and Reynolds shear stress. Large differences are observed between the different decompositions. The wall normal intensities are always over-predicted while the Reynolds stress is under-predicted, which shows a lack of correlation between the streamwise and wall-normal components of the fluctuations. For all second-order quantities, best results are observed for the decomposition [ρ], [ρu, ρv, ρw], [ρe] when rescaling is only applied to the momentum.

•

  For turbulent intensities, rescaled conditions are worse than unrescaled ones, except for the triple decomposition based on [ρ], [ρu, ρv, ρw], [ρe] (momentum rescaling), which means others decompositions compensate more energy for fluctuation. The overprediction is worst for the wall-normal turbulent intensity (over 60% for all other methods and only about 30% for the triplet decomposition [ρ], [ρu, ρv, ρw],[ρe]).• For the triple decomposition based on [ρ], [ρu, ρv, ρw], [ρe], it was best to perform rescaling only on the momentum components. Global rescaling led to a significant increase in the prediction of the turbulent intensities, in particular for the wallnormal component. The unrescaled condition yielded the best prediction for the wall-normal component, because it contains more small scales which are more difficult to control and predict.5.3 Influence of the snapshot basis5.3.1 Evolution of the amplitude of the dominant modeCarrying out POD for the full velocity field yields a dominant mode whose time average corresponds approximately to the mean profile. Examination of the temporal amplitude of the first mode therefore reveals how truly stationary the flow is. Figure5.16 shows the first temporal coefficient of the dominant mode for each POD variable (density, momentum and ρe), where the x axis is time scale, and y axis is amplitude.

Figure 5 .

 5 Figure 5.16: Temporal coefficient of first mode A 1p = α 1 (t p ) (t p is indice of time). Top left: for density[ρ]. Top right: for momentum [ρu, ρv,ρw], bottom: for [ρe].

  ). The new basis has an overlap of 25% with the original basis.

Figure 5 . 17 :

 517 Figure 5.17: Definition of the two snapshot bases for POD.

Figure 5 .

 5 Figure 5.18 represents the temporal coefficient of the first mode of ρe for the new 100 samples. We observe that the amplitude of temporal coefficient for ρE oscillates near a constant. It therefore appears to be a good idea to evaluate the boundary condition constructed from the new basis.

Figure 5 . 18 :

 518 Figure 5.18: Amplitude of the first mode A 1p of [ρe].

Figure 5 . 21 :

 521 Figure 5.21: Turbulent intensity value and shear stresses for [ρ], [ρu, ρv, ρw], [ρe] decomposition. Left: rms wall-normal velocity w. Right: shear stress < u w >.

Figures 5 . 5 . 4 . 1

 5541 Figures 5.19-5.21 show turbulent statistics computed over an integration time of 36 convective units. Here, all modes are used for density ρ and internal energy ρe, 120 modes are used for momentum.The mean streamwise velocity profile is very well recovered, especially in the inertial region characterized by the log-law (Figure5.19) in both parts of the channel. The density and temperature profile also match very well the DNS results. The streamwise turbulent intensity (< u 2 > 1/2 ), the spanwise turbulent intensity (< v 2 > 1/2 ) are well recovered, but the wall normal turbulent intensity (< w 2 > 1/2 ) is overpredicted. Because the boundary conditions for field of u and w are different, the field of w is always determined with POD reconstruction, the field of u could be either determined by POD reconstruction or calculated with resolved region, this results in a lack of correlation between the longitudinal and the wall-normal fluctuations, since the shear stress < u w > in Figure5.21 is well recovered in a large part of the core region.

3 .

 3 Figures 5.22 and 5.23 represent the results for the rescaled conditions A and B.

Figure 5 . 22 :

 522 Figure 5.22: Mean profiles -Top left: streamwise velocity u in bottom half channel. Top right: streamwise velocity u in top half channel. Bottom left: density. Bottom right: temperature. Red line is for DNS. Violet line is for boundary condition A. Black line is for boundary condition B.

Figure 5 . 23 : 22 Figures 5 .

 523225 Figure 5.23: Turbulent intensity values and shear stresses for the different POD methods. Top left: rms streamwise velocity u. Top right: rms spanwise velocity v. Bottom left: rms wall-normal velocity w. Bottom right: shear stress < u w >-Legends as in Figure 5.22

Figure 5 . 24 :

 524 Figure 5.24: From left to right: Matrix Ψ for [ρ],[ρu, ρv, ρw], [ρe].

Figure 5 .

 5 Figure 5.25: Mean profiles -Top left: streamwise velocity u in bottom half channel. Top right: streamwise velocity u in top half channel. Bottom left: density. Bottom right:temperature. Red line is for DNS. Violet line is for simulation of [ρ], [ρu, ρv, ρw], [ρe] with rescaling in momentum. Black line represents adjustment.

Figure 5 . 26 :

 526 Figure 5.26: Turbulent intensity values and shear stresses for the different POD methods. Left: rms streamwise velocity u. Right: rms spanwise velocity v.-Legends as in Figure 5.25.

Figure 5 . 27 :

 527 Figure 5.27: Turbulent intensity values and shear stresses for the different POD methods. Left: rms wall-normal velocity w. Right: shear stress < u w >.-Legends as in Figure 5.25.

2 .

 2 The best result for subsonic flow is from POD reconstruction of [ρ], [ρu, ρv, ρw], [ρe]. The associated Riemann invariants are λ 1 = λ 2 = λ 3 = u • n for the temperature, ρu and ρv, λ 4 = u • n+ n C for ρw, and λ 5 = u • n-n C for the density ρ in reduced channel.

Figure 5 . 28 :

 528 Figure 5.28: Mean profiles -Left: streamwise velocity u in bottom half channel. Right: streamwise velocity u in top half channel.

Figure 5 . 29 :

 529 Figure 5.29: Mean profiles -Left: density ρ. Right: temperature T .

5. 5

 5 Spectra in the reduced channel at h + 0 = 18 114 underprediction of the shear stress corresponds to a local lack of correlation between the streamwise and the wall-normal component. The fact that the shear stress is relatively close to the reference profile is consistent with a fair agreement for the mean velocity profiles shown in Figure 5.28.

Figure 5 . 30 :

 530 Figure 5.30: Turbulent intensities. Left: turbulent streamwise intensity < u 2 > 1/2 . Right: turbulent spanwise intensity < v 2 > 1/2 .

Figure 5 . 31 :

 531 Figure 5.31: Left: turbulent spanwise intensity < w 2 > 1/2 . Right: shear stress < u w >.

Figure 5 . 32 :

 532 Figure 5.32: Power spectra of density for the streamwise velocity u on a plane at an altitude z = 0.2 versus the wave numbers in x and y directions. Left: in full channel. Right: reduced channel.

Figure 5 . 33 :

 533 Figure 5.33: Power spectra of density for the spanwise velocity v on a plane at an altitude z = 0.2 versus the wave numbers in x and y directions. Left: full channel. Right: reduced channel.

Figure 5 . 34 :

 534 Figure 5.34: Power spectra of density for the wall-normal velocity w on a plane at an altitude z = 0.2 versus the wave numbers in x and y directions. Left in full channel. Right in the reduced channel.

Figures 5 .

 5 Figures 5.35-5.37 compare the one-dimensional streamwise spectra of the velocity components u, v, w at different heights z = 0.2, 0.6, 1.0 corresponding in wall units: z + = 36, 108, 180. The spectrum decay in the reduced domain appears quite similar to that in the reference domain. The main difference seems to be a smaller energetic content for the wall-normal component w which is noticeable in particular for the large scales.

Figure 5 . 35 :

 535 Figure 5.35: Power spectra of density for the streamwise velocity u on a plane at three different altitudes versus the wave numbers in x direction. Left : full or reference channel. Right: reduced channel. (Red line is for spectra at altitude z = 0.2. Violet line is for spectra at altitude z = 0.6. Black real line is for spectra at altitude z = 1.0. Black imaginary line is for inertial region).

5. 6 Figure 5 . 36 :

 6536 Figure 5.36: Power spectra of density for the spanwise velocity v on a plane at three different altitudes versus the wave numbers in x direction. Left: full or reference channel. Right: reduced channel.-Conditions as in Figure 5.35.

Figure 5 . 37 :

 537 Figure 5.37: Power spectra of density for the wall normal velocity w on a plane at three different altitudes versus the wave numbers in x direction. Left: full or reference channel. Right: reduced channel.-Conditions as in Figure 5.35.

Figures 5 .

 5  show respectively the mean profiles of the conservative variables and the different components of the shear stresses. The mean velocity profiles are well re-

Figure 5 . 38 :

 538 Figure 5.38: Mean profiles -Top left: streamwise velocity u in bottom half channel. Top right: streamwise velocity u in top half channel. Bottom left: density. Bottom right: temperature.

Figure 5 . 39 :

 539 Figure 5.39: Turbulent intensity values for the synthetic boundary condition at h + 0 = 54. Left: rms streamwise velocity u. Right: rms spanwise velocity v.

Figure 5 .

 5 Figure 5.40: Turbulent intensity value and shear stresses for the synthetic boundary condition at h + 0 = 54. Left: rms wall-normal velocity w. Right: shear stress < u w >.

  -5.40 appear to be overpredicted.The overestimation is about 20% for the streamwise component, while it reaches respectively 40% and 100% in the center of the channel for the wall-normal and the spanwise component. However the shear stress (see Figure5.40) is well predicted. The very low values of the Reynolds stress on the boundary plane can be interpreted as a lack of correlation between the streamwise and the wall-normal component. This might be due to the fact that part of the boundary condition is prescribed, while the rest is computed from the simulation.

1

 1 Mesh interpolation for PODNow we consider the supersonic case M a = 1.5. The supersonic simulation was carried out with mesh dimension N x × N y × N z = 257 × 161 × 161, as mentioned in Chapter 2.

Figure 6 . 1 :6. 2

 612 Figure 6.1: Sketch of the interpolation procedure on the refined mesh.

Figure 6 . 2 :

 62 Figure 6.2: Instantaneous field of streamwise velocity u at altitude z = 0.2. Left: full channel at 4.3 time units. Right: reduced channel at 4.3 time units.

Figure 6 . 3 :

 63 Figure 6.3: Instantaneous field of spanwise velocity v at altitude z = 0.2. Left: full channel at 4.3 time units. Right: reduced channel at 4.3 time units.

Figure 6 . 4 :

 64 Figure 6.4: Instantaneous field of wall normal direction velocity w at altitude z = 0.2. Left: full channel in 4.3 time units. Right: reduced channel at 4.3 time units.

Figure 6 . 5 :

 65 Figure 6.5: Instantaneous field of total energy at altitude z = 0.2. Left: full channel at 4.3 time units;. Right: reduced channel at 4.3 time units.

Figure 6 . 6 :

 66 Figure 6.6: Instantaneous field of density at altitude z = 0.2. Left: full channel at 4.3 time units. Right: reduced channel at 4.3 time units.

Figure 6 . 7 :

 67 Figure 6.7: Mean profile of quantities for simulation in 26 time units. Top left: mean profile of u for bottom half channel. Top right: mean profile of u for top half channel. Bottom left: mean profile of density. Bottom right: mean profile of temperature. (Red line is for DNS. Violet line is for simulation without rescaling step. Black line is for simulation with rescaling step).

Figure 6 . 8 :

 68 Figure 6.8: Left: rms value of u. Right: rms value of v.-Conditions as in figure 6.7.

Figure 6 . 9 :

 69 Figure 6.9: Left: rms value of w. Right: shear stress < u w >.-Conditions as in Figure 6.7.

Figure 6 . 10 :

 610 Figure 6.10: Mean profile of quantities for simulation in 40 time units. Left: mean profile of u for bottom half channel. Right: mean profile of u for top half channel.

Figure 6 . 11 :

 611 Figure 6.11: Mean profile of quantities for simulation in 40 time units. Left: mean profile of density. Right: mean profile of temperature.

Figure 6 .

 6 Figure 6.12: rms values for simulation in 40 time units. Left: rms value of u. Right: rms value of v.

Figure 6 . 13 :

 613 Figure 6.13: rms value and shear stress for simulation in 40 time units. Left: rms value of w. Right: shear stress< u w >.

Quantity error time units rms u rms v rms 2 :

 2 Relative error for simulations in reduced channel 0.3 < z ≤ 1.7 obtained from 70 modes (Comparison with subsonic flow in Chapter 5).

Figure 6 .

 6 Figure 6.14: Wave propagetion in streamwise direction. Left: subsonic flow. Right: supersonic flow.

  Figure 6.15: Amplitudes of the first mode for 100 samples (corresponding to the mean field). Top left:[ρ]. Top right:[ρu, ρv, ρw]; bottom: [ρe].

Figure 6 .

 6 Figure 6.16: Procedure of simulation for DNS.

  -6.20.

Figure 6 . 17 :

 617 Figure 6.17: Mean profile of quantities for simulation in 40 time units. Left: mean profile of u for bottom half channel. Right: mean profile of u for top half channel.

Figure 6 . 18 :

 618 Figure 6.18: Mean profile of quantities for simulation in 40 time units. Left: mean profile of density. Right: mean profile of temperature.

Figure 6 . 19 :

 619 Figure 6.19: rms values for simulation in 40 time units. Left: rms value of u. Right: rms value of v.

Figure 6 . 20 :

 620 Figure 6.20: rms value and shear stress for simulation in 40 time units. Left: rms value of w. Right: shear stress< u w >.

  Figure 6.21: Spectra of energy of u at altitude z = 0.2. Left: full channel. Right: reduced channel.

3 )

 3 We integrate in plane yz, then we have:

  B.7)Normally with experimental experience, α = 2/∆t and β = -0.2/∆t. With this definition, we can use macroscopic pressure graident to maintain the flow dynamics. [88] X. Yang, J. Sadique, R. Mittal, and C. Meneveau. Integral wall model for Large Eddy Simulations of wall-bounded turbulent flows. Physics of Fluids, 27(2):025112, 2015. 7, 8 Titre : Condition aux limites synth étiques pour la turbulence compressible pr ès de la paroi avec reconstruction du POD Mots cl és : Simulation, turbulence, r éduction des co ûts de calcul R ésum é : La simulation num érique directe des couches limites turbulentes n écessite une r ésolution spatiale élev ée, en particulier proche des parois afin de capturer les diff érentes échelles de longueur inh érentes à la turbulence. Pour r éduire le co ût de calcul li é à la haute r ésolution spatiale, nous étudions les techniques qui éliminent la r égion de paroi en utilisant une condition aux limites artificielle aux bords d'un domaine de calcul r éduit. Une strat égie de reconstruction bas ée sur la POD est mise en oeuvre afin de g én érer des conditions aux limites synth étiques permettant de r écup érer les caract éristiques spatiales et dynamiques correctes de la turbulence dans la r égion centrale du canal. La condition aux limites a ét é test ée à deux hauteurs diff érentes. Les statistiques du canal r éduit ont ét é compar ées à celles du canal complet prise comme valeur de r éf érence. Une bonne concordance est observ ée pour le profil de vitesse moyenne et la contrainte de cisaillement turbulente. La proc édure s'est av ér ée robuste à la hauteur à laquelle la condition aux limites est impos ée et permet une r éduction drastique des co ûts de calcul. Title : Synthetic boundary condition for compressible near wall turbulence Keywords : Simulation, turbulence, cost reduction Abstract :

  1 (from[START_REF] Perret | Generation of turbulent inflow conditions for Large Eddy Simulations from stereoscopic PIV measurements[END_REF]). . . . . . 3.11 Wave propagation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Boundary condition on the wall parallel plane. . . . . . . . . . . . . . . . 3.13 Boundary condition on the boundary planes. . . . . . . . . . . . . . . . . 3.14 Model of simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.15 Field of Density on middle plane x=π at 0.3 time units. Left image is for

DNS. Right image is for reduced channel. . . . . . . . . . . . . . . . . . . 3.16 Field of u on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel. . . . . . . . . . . . . . . . . . . . . . 3.17 Field of v on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel. . . . . . . . . . . . . . . . . . . . . . 3.18 Field of w on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel. . . . . . . . . . . . . . . . . . . . . . 3.19 Field of total energy E on middle plane x=π at 0.3 time units. Left image is for DNS. Right image is for reduced channel. . . . . . . . . . . . . . . 4.1 Model of reduced channel. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Mean profiles averaged over time and expressed with horizontal directions in the bottom half channel over 6 time units for POD in primitive variables. Top left: mean profile of density; top right: mean profile of u; bottom left: mean profile of v ; bottom right: mean profile of w. . . . . .

Table 2 .

 2

	.8	0.65 0.65

  y, z, t) =

	Q kxky (z, t)e 2iπ kx Lx e 2iπ	ky Ly .	(3.16)
	kx ky

Table 3 . 2

 32 .2. Which is applied as parameters in the later chapters.

	Subsonic	Ma Re	Duration	Statistic	Number	Number	Duration
			for tran-	duration	of sam-	of POD	for Sam-
			sition		ples	modes	ples
		0.5 3000	700 time	120 time	100	200	600 time
			units	units			units
	Supersonic Ma Re	Duration	Statistic	Number	Number	Duration
			for tran-	duration	of sam-	of POD	for Sam-
			sition		ples	modes	ples
		1.5 3000	1500	75 time	30	60	75 time
			time	units			units
			units

: Data basis for calculation in channel flow.

  .1) 4.1 Results at height h + 0 = 18 (h 0 = 0.1) with primitive variables 65 where f Q is the rescaling factor determined in Chapter 3 (see Equation 3.35), and Q mean is the mean value which can be computed either as a time and spatial average (global mean value) or simply as a spatial average (local mean value):

Table 4 .

 4

	Quantity error Without rescaling step Rescaling step with	Rescaling step with
			global mean value	local mean value
	rms u	0.1284	0.0078	0.071
	rms v	0.096	0.0874	0.0599
	rms w	0.1608	0.1133	0.1131
	< ρ >	6E-4	2.7E-4	1.9E-4
	< u >	0.0144	0.009	0.0065

2: Relative errors for reduced-domain simulations with POD based on conservative variables in region 0.3 < z ≤ 1.

  .4.

	Quantity error Without rescaling step Rescaling step with	Rescaling step with
			global mean value	local mean value
	rms u	0.0948	0.0249	0.0248
	rms v	0.1240	0.0446	0.0443
	rms w	0.2066	0.1062	0.1063
	< ρ >	1.73E-4	3.63E-4	1.19E-4
	< u >	0.0072	0.0047	0.0046
	Table 4.4: Relative errors for simulations with conservative variables in reduced channel
	0.3 < z ≤ 1.

Table 5 .

 5 2: Statistical errors for rescaled boundary conditions. The words between parentheses indicate the type of scaling that was performed, depending on whether all components (global) or only the momentum were rescaled.

	unrescaled POD	Time units < u > < ρ >
	[ρ] [ρu ρv ρw] [ρE]	36	0.0402 0.0048
	[ρ ρu ρv ρw ρE]	36	0.0386 0.0039
	[ρ] [ρu ρv ρw] [ρe]	36	0.0452 7.5238E-4
	[ρ] [u v w] [E]	36	0.0528 0.007
	unrescaled POD	Time units rms u rms v	rms w
	[ρ] [ρu ρv ρw] [ρE]	36	0.1239 0.1348	0.3150
	[ρ ρu ρv ρw ρE]	36	0.1076 0.1324	0.3933
	[ρ] [ρu ρv ρw] [ρe]	36	0.1844 0.152	0.158
	[ρ] [u v w] [E]	36	0.1826 0.1508	0.2093
	Table 5.1: Statistical errors for unrescaled boundary conditions.

Table 5 .

 5 [START_REF] Adrian | Stochastic estimation of organized turbulent structure: homogeneous shear flow[END_REF] shows that it is important to define the variables on the boundary with the Riemann invariants otherwise the computation in the reduced domain simulation does not appear to be stable. Prescribing a velocity component to be zero in order to avoid the introduction of the small scales, in particular for w, did not seem to yield good results. The best choice appeared to associate ρ, ρw, ρv, ρu, T with the eigenvalues w -C, w + C, w, w, w.

	Eigenvalue w	w w	w + C	w -C	Stabilization in 36 time units
	1	ρu	ρw T	ρv	ρ	No
	1	ρu	ρv ρw	ρ	T	No
	3	T	ρ ρv	ρw	ρu	No
	4	ρu	ρw ρ	T	ρv	No
	5	ρu	ρw ρv	T	ρ	No
	6	ρ	ρw ρv	ρu	T	No
	7	ρu	ρw d(ρv) = 0	T	ρ	No
	8	T	ρw d(ρv) = 0	ρu	ρ	No
	9	ρ	ρw d(ρv) = 0	ρu	T	No
	10	ρu	ρv ρ	T	ρw	No
	11	ρu	ρv d(ρv) = 0	T	ρw	No
	12	ρu d(ρv) = 0 ρ	T	ρw	No
	13	ρ	ρu ρv	ρw	T	Yes	(Discrepancies
						for mean density and
						temperature)
	14	T	ρu d(ρv) = 0	ρw	ρ	Yes (Discrepancies for
						spanwise velocity)
	15	T	ρu ρv	ρw	ρ	Yes (Best result)

Table 5 .

 5 3: Comparison between different choices of Riemann invariants for the boundary conditions (eigenvalue expressed in bottom boundary plane).

Table 5 .

 5

		Mesh dimension	δt CPU time (36
			convective units)
	Full channel h + 0 = 0 Reduced channel h + 0 = 18	1213761 799765	6E-4 48 hours 1.8E-3 20 hours
	Ratio Reduced/Full	66%	300% 42%
	Reduced channel h + 0 = 54	498677	3.6E-3 14 hours
	Ratio Reduced/Full	41%	600% 29%

4: Computational requirements of the reference and reduced-channel simulations.

Table 6 .

 6 4: Computational time saving for supersonic flow.

	Ma=1.5 (h + 0 ∼ 20)	Mesh dimension Time step Total simulated time (40 H/U 0 )
	Reduced channel	3.35E6	4.3E-4 121 h
	Full channel	6.7E6	2.5E-4 288 h
	Ratio Reduced/Full	50%	170% 42% = (50%/1.7)

fields sampled every 10000 time steps indicated as "old samples" in Figure5.17)

Acknowledgements

Chapter 5

Synthetic boundary conditions on both walls

In this chapter we replace both walls of the channel with synthetic boundary conditions.

An important consequence is that since there are no more walls, there is no direct way to evaluate the friction and therefore the pressure gradient. The external force and the flow rate were therefore set at a constant value determined from the reference simulation.

The goal of this chapter is to determine the influence of various parameters of the reconstruction procedure:

• Fourier-based POD decomposition is tested (section 5.1),

• Different choices for the variables of decomposition are examined (section 5.2),

• The representativity of the snapshot set is examined (section 5.3),

• The influence of the Riemann invariants is tested (section 5.4),

• Different heights are tested (section 5.6).

Fourier-based reconstruction

In this section we construct synthetic boundary conditions using Fourier-based POD. Relative errors defined as in (6. Appendix A

Viscous flux

For the viscous stress, we have the following expression based on discretization:

We defined that Π is identity matrix. For the discretization in x direction:

For the discretization in y direction:

For the discretization in z direction:

For better understanding of discretization, we define P art 1 ,P art 2 ,P art 3 as following (P art 1 ,P art 2 ,P art 3 are discretizations of τ is x, y, z directions respectively):

And we have the following results as:

)δx i δy j (A.9)

)]δx i δy j (A.10)

As for the term U τ :

For better understanding of discrezation, we define P art 4 , P art 5 , P art 6 as following for better understanding:

δx i δy j (A.12)

And we have the following results as:

)δx i δy j (A.16)

)]δx i δy j (A.17)

For the part µ∇T M 0 2 (γ-1)P rRe : Then, we integrate between 0 and L x : For the pressure, we suppose that the pressure gradient ∂P ∂x could be expressed as the combination of macroscopic pressure gradient f 1 and periodic pressure gradient ∂p ∂x ., so we have ∂P ∂x = f 1 + ∂p ∂x .