Jean Zay

Rapid increase in the early 2000s Examples in the french industry 52 Z3

Keywords: URANS, SDM, OES, VLES, DES, LNS, PANS, XLES, PITM, TPITM, FSM, SBES, etc 2000 Rapid development RANS

� Large-eddy simulation (LES) � Unsteady statistical methods: URANS, OES, SAS � Zonal and continuous hybrid RANS/LES methods 7 8 � Moore's law : Moore (1975) : The number of transistors of a processor double every second year (10 years=factor of 32) source: http://www.tablette-tactile.net 39 � Supposed to be valid until about 2015 ("quantum wall") � However, the clock frequency increase slowed down after 2003 and in particular after 2006: "thermal wall" = the heat produced cannot be removed anymore source : http://www.gotw.ca/publications ⇒ Evolution toward multi-core processors for personal computers 40 3.5.2. Supercomputers � Analogy: building of a wall (adapted from V. Perrier) Task Working time � How to reduce the working time (computational time)? 41 � Task decomposition ❀ Principle � � � � Does not work: take care of the scheduling of the tasks! � Parallel machines Association of several processors or cores (from 2 to millions) Laptop � 6 cores � 122 GFlops (Giga = 10 9 ) � 16 GB Cluster Univ. Pau

Numerical simulation

Results: wall temperature Structured mesh (hexahedra)

Better stability and accuracy properties. Difficulties for meshing complex geometries Unstructured mesh (tetrahedra).

Meshing complex geometries is easy. Problems in near-wall regions ⇒ Mixed meshes: hexahedra in the boundary layer, tetrahedra elsewhere

Turbulent scales

Vizualisation of a mixing layer. From [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF]. � Often very specific: weakly flexible but very efficient (confined to simple geometries, to a specific physical phenomenon, a particular method, etc.)

τ ij = u * i u * j -u * i u * j �
� There are many of them: at least 10 at the institute P'

Collaborative codes (French centered)

� Some research centers (CEA, ONERA, CERFACS, IFPEN, etc.) or companies (Airbus, Dassault, EDF, etc.) develop and use collaborative codes (several groups or several research centers, many users)

� This development implies a significant financial effort � But is beneficial for: � They often propose both RANS and LES models

� Many have now hybrid models (DES, SAS)

� There is a rise in Lattice-Boltzmann methods (Powerflow, Xflow, ProLB)

Conclusion of this general introduction

� Many things can be done with CFD, but the limiting factor is the computing power, because of turbulence � One has to model turbulence � For industrial applications (design), the standard method is RANS modeling � Many parameters have an influence on the results: � The experience gained increases the knowledge of the range of application of computations (repeatability, influential parameters, recurrent modeling problems, experiment/computation discrepancies, etc.)

❀
� This know-how must be sustained:

❀ models evolve ❀ computing power increases very fast

Introduction

� Aim of turbulence modeling: to replace the Navier-Stokes equations by a model.

❀ The numerical resolution of the equations (the model) must be as cheap as possible: overnight computations are necessary for parametric studies and much shorter for optimization.

❀ The model must be predictive: only the parameters of the flow must be necessary ⇒ no a priori knowledge of the solution.

❀ The model must represent at best the physics of the flow.

75

❀ The model must give the useful quantities:

✄ At least global quantities: aerodynamic forces (drag, lift, etc.), heat transfer between fluids and solids, mixing of a pollutant, etc.;

✄ very often: variations of these quantities with parameters (velocities, temperature differences, shape parameters), in order to optimize a system;

✄ but also: separation location of boundary layers, pressure field at the wall, flow structure, shocks location, acoustic sources, etc.;

✄ and more and more: response to a control strategy (blowing, suction, MEMS, etc..).

⇒ Turbulence modeling, is the science that consists in building these models, for the purpose of their use by engineers.

This science is seeking an impossible goal: provide a simple, cheap, robust model, for all situations.

� Einstein's notation is used: implicit summation over indices repeated in a term.

In simpler words:

∂φ i ∂t + u k ∂φ i ∂x k = 0
is a compact expression standing for

∂φ 1 ∂t + u 1 ∂φ 1 ∂x 1 + u 2 ∂φ 1 ∂x 2 + u 3 ∂φ 1 ∂x 3 = 0 ∂φ 2 ∂t + u 1 ∂φ 2 ∂x 1 + u 2 ∂φ 2 ∂x 2 + u 3 ∂φ 2 ∂x 3 = 0 ∂φ 3 ∂t + u 1 ∂φ 3 ∂x 1 + u 2 ∂φ 3 ∂x 2 + u 3 ∂φ 3 ∂x 3 = 0
There is no summation over i since it is not repeated in a single term.

This two-term equation actually represents 3 equations containing 4 terms. It is compact!

Be careful:

∂ 2 φ ∂x 2 k (no repetition) is different from ∂ 2 φ ∂x k ∂x k 79
� Although the basic equations used in fluid mechanics are well known, it is worth recalling them here since the way they are derived is very similar to what will be made below in turbulence modeling. The derivation can be summarized in three steps:

❀ writing of the conservation equations; ❀ closure of the system using constitutive relations, which introduce state variables that describe the macroscopic state of the fluid at a given point and a given time;

❀ writing of equations of state, which describe the evolution of these state variables.

Conservation equations:

� Mass conservation (continuity):

∂ρ * ∂t + ∂ρ * u * i ∂x i = 0 � Momentum conservation: ∂ρ * u * i ∂t + ∂ρ * u * i u * j ∂x j = ∂σ * ij ∂x j + ρ * g i � Energy conservation: ∂ρ * e * ∂t + ∂ρ * u * i e * ∂x i = σ * ij s * ij - ∂ ∂x i γ i
� These equations describe the behavior of the fluid at macroscopic scale, whatever the type of fluid. However, they involve 14 unknowns:

❀ density ρ * ;
❀ the three components of velocity u * i ;

❀ the six independent components of the stress tensor σ * ij (symmetric tensor) ;

❀ internal energy e * ;

❀ and the three components of the heat flux γ i ;

for 5 equations only.

� In order to close the system, constitutive equations must be introduced, which describe the properties of the fluid (and thus are dependent on the type of fluid). 81

Constitutive relations:

A fluid is Newtonian if it can be described by the following linear laws giving:

� the stress tensor σ * ij as a function of the strain tensor s * ij (Newton's law):

σ * ij = � -p * + λ ∂u * l ∂x l � δ ij + 2µs * ij ( 1
)
where

s * ij = 1 2 � ∂u * i ∂x j + ∂u * j ∂x i �
, and δ ij is the identity tensor (Kronecker's symbol),

� and the heat flux γ i as a function of the temperature gradient (Fourier's law):

γ i = -k ∂T * ∂x i (2)
� We have provided the 9 equations which were missing.

� The properties of the fluid (dynamic viscosity µ and volume viscosity λ, conductivity k) are known.

� But two state variables (pressure p * and temperature T * ) have been introduced

which describe the state of the fluid at a given point and a given time (i.e., they describe at the macroscopic scale the properties of the molecular agitation).

⇒ We need additional relations, equations of state to obtain the evolution of these state variables.

83

Equations of state:

� For a perfect gas, we can write:

p * = ρ * rT * (3)
and

de * = C v dT * (4)
� The first relation is replaced by the incompressibility constraint

∂u * i ∂x i = 0 (5)
in the case of an incompressible flow (thus introducing a famous paradox: the equation of state does not involve pressure).

� The solutions of this system are very complex, and contain all the physical phenomena observed in a Newtonian fluid flow: vortices (in particular turbulence), shocks, acoustic waves, transformation of kinetic energy into internal heat, etc.

� This system can sometimes by simplified, but one must first answer these questions:

❀ Which assumptions make theses simplifications possible?

❀ Is it sufficient to address the problem?

� Example: it is often assumed that µ, λ, k and C v are independent of temperature.

Warning: this is not valid for large temperature differences.

� Example: incompressible flow (div u * i = 0). Incompressible fluids do not exist! However, assuming that the effects of compressibility are negligible is justified in numerous cases (low Mach numbers).

Warning: this excludes acoustic waves and shocks from the solution.

� Low-Mach-number approximation:

❀ For low Mach numbers, it can be assumed that density does not depend on pressure, but only on temperature

ρ * = f (P * , T * )
❀ For a perfect gas, we have

p * = ρ * rT * such that dρ * = ∂ρ * ∂P * � � � T * dP * + ∂ρ * ∂T * � � � P * dT * = - ρ * T * dT * and ρ * = ρ * 0 T * 0 T *
❀ The flow is thus considered incompressible, but density varies as a function of the inverse of the temperature: the fluid is dilatable.

❀ The equations of motion can be derived using asymptotic expansions at the limit of small Mach numbers.

� A very common approximation, when the temperature differences are weak, is the Boussinesq approximation: the density variations can be negelected in the Navier-Stokes equations (ρ * = ρ 0 ), except in the buoyancy term ρ * g i (buoyant force).

� Usual additional simplifications:

❀ µ, k and C v are assumed to be independent of temperature ❀ A Taylor series expansion at first order is made of the density around the reference temperature T * 0 :

ρ * = ρ * 0 -βρ * 0 (T * -T * 0 ) with β = - 1 ρ * 0 ∂ρ * ∂T (Isobar compressibility coefficient)
such that the buoyancy term reads:

ρ * g i = ρ * 0 g i -ρ * 0 g i β(T * -T * 0 )

87

� When the temperature variations are weak and the buoyancy effects are negligible (forced convection), the density ρ * can be considered as a constant (simply denoted by ρ).

The Navier-Stokes equation then reduce to

∂u * i ∂x i = 0 (6) ∂u * i ∂t + u * j ∂u * i ∂x j = - 1 ρ ∂p * ∂x i + g i + ν ∂ 2 u * i ∂x j ∂x j ( 7
)
∂T * ∂t + u * i ∂T * ∂x i = α ∂ 2 T * ∂x i ∂x i + 2µ ρC v s * ij s * ij ( 8
)
� Remarks:

❀ the volume viscosity λ has disappeared;

❀ the kinematic viscosity ν = µ/ρ and the thermal diffusivity α = k/ρC v naturally appear;

❀ the dynamics does not depend on the thermal field ⇒ the dynamics can be resolved, Eqs. ( 6) et [START_REF] Bentaleb | A hybrid Temporal LES/RANS formulation based on a two-equation subfilter model[END_REF], and then the velocity field used for solving the temperature equation ⇒ the temperature is then a passive scalar ;

❀ if the flow is isothermal (T * = constant), equation ( 8) disappears.

89

� If the gravitational field is oriented along -z:

g i =     0 0 -g     one can define: p * * = p * + ρgz ⇒ - 1 ρ ∂p * ∂x i + g i = - 1 ρ ∂p * * ∂x i � It is
seen that when the buoyancy forces can be neglected, the gravity has no influence on the velocity field: gravity only influences the pressure.

� Most of the time, in CFD codes, when the buoyancy forces can be neglected, p * * is used, and then g i disappears from the equations. Warning: this does not mean that the gravity is neglected, but rather that the hydrostatic pressure is included in p * * . ⇒ The computation will not provide the pressure, but p * * ! p * can be obtained from p * *ρgz.

The Reynolds decomposition

� In RANS modeling, the first modeling consist in assuming that the instantaneous quantities are random variables: ⇒ a random behavior of the flow is assumed, while the Navier-Stokes equations are obviously deterministic. This is the chaotic character of the system that makes this assumption possible.

⇒ The deterministic character of turbulence is not accounted for: in particular, the presence of coherent structures (large-scale eddies, with, very often, a pseudo-periodic character) can play a significant role in the flow (see the figure below).

⇒ The flow is thus presumed fully turbulent: the problem of transition to turbulence is outside of the scope of the present course. � The flow is decomposed into mean and fluctuating parts by using the Reynolds average: in RANS modeling, this average is denoted by .

Mean velocities:

U i = u * i Fluctuating velocities: u i = u * i -U i
Mean pressure: P = p * * (the hydrostatic pressure is thus included)

Fluctuating pressure: p = p * * -P

93

� The Reynolds average of the turbulent (random) variable f * is its statistical average, equivalent to the ensemble average (law of large numbers), defined by:

f * (x, t) = lim N →∞ � 1 N N � n=1 f * n (x, t) � (9)
i.e., evaluated at each point and each time by performing a large number of experiments.

Realizations Ensemble average

In many cases, the ensemble average can be expressed in a different form which is simpler to measure:

� Statistical stationarity = statistic quantities are independent of time:

f * (x, t) = f * (x).
In that case, the ensemble average is equivalent to a temporal average:

f * (x) = lim T →∞ 1 T � T 0 f * (x, t)dt (10)
� Statistical homogeneity in one (or several) direction = statistical quantities are independent of the location (x hereafter): f * (x, y, z, t) = f * (y, z, t). In that case, the ensemble average is equivalent to a spatial average:

f * (y, z, t) = lim L→∞ 1 L � L 0 f * (x, y, z, t)dx ( 11
)
� Statistical periodicity = statistical quantities are periodic with a period τ :

f * (x, t) = f * (x, φ)
where φ = 2πt/τ modulo 2π. In that case, the ensemble average is equivalent to a phase average:

f (x, t) = lim N →∞ 1 N + 1 N � n=0 f (x, t + nτ ) ( 12
)
95

The closure problem

� The introduction of the Reynolds decomposition leads to the averaged Navier-Stokes equations:

∂U i ∂x i = 0 ( 13
)
∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j ( 14
)
� There are 4 equations for 10 unknowns:

P , U , V , W , u 2 , v 2 , w 2 , uv, uw, vw.
⇒ the system is open.

� Writing the transport equations for the Reynolds stress tensor u i u j :

∂u i u j ∂t + U k ∂u i u j ∂x k � �� � C ij = ν ∂ 2 u i u j ∂x k ∂x k � �� � D ν ij - ∂u i u j u k ∂x k � �� � D T ij - 1 ρ u i ∂p ∂x j - 1 ρ u j ∂p ∂x i � �� � φ * ij -u i u k ∂U j ∂x k -u j u k ∂U i ∂x k � �� � P ij -2ν ∂u i ∂x k ∂u j ∂x k � �� � ε ij ( 15
)
leads to 6 new equations... ... but unfortunately to 34 new unknowns:

u i u j u k , u i ∂p ∂x j , ∂u i ∂x k ∂u j ∂x k ⇒ the problem is still open.
� One could write equations for these unknowns, but the process always produces more unknowns than equations.

� One must decide when to stop.

97

� First moment closure (or eddy-viscosity modeling):

❀ The 4 equations ( 13) and ( 14) are solved.

⇒ the first moments are obtained: U , V , W , P ❀ It is necessary to "invent" a relation between the variables that are not resolved as functions of the variables that are resolved: this relation providing the second moments (the u i u j 's) as functions of the first moments is called a first order model or, more widely, an eddy-viscosity model. ❀ This is the equivalent of a constitutive relation for a material: for instance, for a Newtonian, incompressible fluid, the stress tensor is related to the strain tensor by

σ * ij = -p * δ ij + 2µs * ij
This is a model for the fluid, the model of Newton. This model closes of the system (the Navier-Stokes equations).

❀ An eddy-viscosity model is the equivalent of the model of Newton for the "turbulent fluid": this is a relation between the Reynolds stress tensor (which exactly play the role of a stress in [START_REF] Craft | A Reynolds stress closure designed for complex geometries[END_REF] to the mean pressure and the mean strain tensor.

Moreover, the simplest models are based on the Boussinesq relation, similar to the Newton relation:

-ρu i u j = - 2 3 ρkδ ij + 2µ t S ij where S ij = 1 2 � ∂U i ∂x j + ∂U j ∂x i �
where k = 1 2 u i u i is the turbulent energy.

99

� Second moment closure, SMC (or Reynolds stress models, RSM):

❀ The 4 equations ( 13) and ( 14) are solved.

⇒ to obtain the first moments: U , V , W , P ❀ The 6 transport equations for the Reynolds stress tensor [START_REF] Craft | Development and application of a cubic eddy-viscosity model of turbulence[END_REF] are also solved ⇒ to obtain the second moments: u i u j ❀ It is necessary to "invent" relations giving the variables that are not solved as functions of the variables that are solved: such relations giving the unknowns

� u i u j u k , u i ∂p ∂x j , ∂u i ∂x k ∂u j ∂x k �
as functions of the first and second moments form a second moment closure.

❀ These are also equivalent to constitutive relations, but very complex ones.

Example: the Daly-Harlow model for the triple correlations is

u i u j u k = -C s k ε u k u l ∂u i u j ∂x l
❀ A Reynolds stress model solves the transport equations for the Reynolds stresses, while eddy-viscosity models evaluate them from the mean velocities and pressure by a simple algebraic relation: Reynolds stress models contain obviously "more physics".

❀ In particular, the production terms does not require modeling, since they only involve moments of order less or equal to two:

P ij = -u i u k ∂U j ∂x k -u j u k ∂U i ∂x k
This terms have a fundamental influence. � But it is not sufficient: everything is dependent on the type of flow! ❀ For attached flows (airfoil at low angle of attack), a simple model can be sufficient to obtain the drag (mixing length model).

❀ For more complex flows, a more sophisticated model can be necessary to obtain the drag! ❀ As a general rule, the understanding of the physical mechanisms that play a role in a flow and a good knowledge of the potential of the models to reproduce these mechanisms is necessary to find the relevant adequacy.

❀ A "universal" model, that is able to predict everything in any situation does not exist.

103

� Example: the Ahmed body (simplified car)

From [START_REF] Manceau | Computation of the flow around a simplified car using the rescaled v 2f model[END_REF] From [START_REF] Ahmed | Some salient features of the time-averaged ground vehicle wake[END_REF] ❀ The flow structure is very complex.

❀ The pressure at the surface of the body is very dependent on the separated regions.

⇒ A very refined model is necessary, even to simply obtain the drag (actually, nowadays, none of the available RANS models give satisfactory results!)

∂u * i ∂t + u * j ∂u * i ∂x j = 1 ρ ∂σ * ij ∂x j = - 1 ρ ∂p * ∂x i + ν ∂ 2 u * i ∂x j ∂x j ∂ρC v T * ∂t + u * i ∂ρC v T * ∂x i = ∂γ i ∂x i + σ * ij s * ij = α ∂ 2 ρC v T * ∂x i ∂x i + 2ρνs * ij s * ij
Modelling the stress tensor σ * ij introducing a molecular viscosity leads to the appearance of two effects:

❀ A diffusion in the momentum equation ν ∂ 2 u * i ∂x j ∂x j ❀ A dissipation of energy into heat 2ρνs * ij s * ij 107
� The idea of Boussinesq can be expressed as:

❀ The enhancement of mixing by turbulence can be modelled as a mean diffusion effect ❀ The resistance to the flow can be modelled as a mean dissipation effect ⇒ The flow behaves as if a strong additional viscosity was introduced � Mean momentum equation

∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j
❀ -ρu i u j plays exactly the same role as σ * ij ⇒ it is denoted as the Reynolds stress tensor ❀ It is thus modelled by introducing a turbulent viscosity:

ν is replaced by ν + ν t where ν t � ν ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j = ν ∂ 2 U i ∂x j ∂x j + ∂ ∂x j � ν t ∂U i ∂x j � = ∂ ∂x j � (ν + ν t ) ∂U i ∂x j � 109
❀ We know nowadays that the model writes

u i u j = -2ν t S ij + 2 3 kδ ij where S ij = 1 2 � ∂U i ∂x j + ∂U j ∂x i

�

This is the famous Boussinesq relation ❀ Note that the half-trace of the tensor is correct:

1 2 u i u i = -ν t S ii + 1 3 kδ ii = k since S ii = 0
where k = 1 2 u i u i is the turbulent energy 110 ❀ The introduction of the turbulent viscosity indeed leads to:

✄ A strong diffusion in the mean momentum equation: 13) and [START_REF] Craft | A Reynolds stress closure designed for complex geometries[END_REF].

∂ ∂x j � (ν + ν t ) ∂U i ∂x j � ✄ A strong

Digression:

The introduction of the Boussinesq relation in the mean velocity equations

∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j leads to ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂ ∂x i � P + 2 3 ρk � + ∂ ∂x j � (ν + ν t ) � ∂U i ∂x j + ∂U j ∂x i �� Defining: P * = P + 2 3 ρk ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P * ∂x i + ∂ ∂x j � (ν + ν t ) � ∂U i ∂x j + ∂U j ∂x i �� ( 16
)
⇒ It is seen that, similar to gravity in the Navier-Stokes equations, the turbulent energy in the averaged equations do not modify the velocities, but only the pressure ⇒ the knowledge of ν t is sufficient to evaluate the mean velocities.

Warning: in many CFD codes, it is indeed equation ( 16) that is solved ⇒ P is not provided by the computation, but rather P * .

113

� The mixing is mainly due to the large eddies (integral scale) ⇒ ν t is linked to their characteristics � Indeed, if we consider, in order to simplify, the mixing of a passive scalar (temperature, concentration, etc.):

❀ For instance, at some time, let us introduce a drop of dye (instantaneous concentration=φ * ) in a turbulent flow with zero mean velocity ❀ Turbulent mixing is due to the advection of fluid particles by the fluctuating velocity (the eddies)

❀ The dye is mixed by convection at the instantaneous level of description.

� What is observed on average? (i.e., repeating the same experiment a large number of times)

� The dye spot gradually spreads and the concentration decreases, due to mixing.

⇒ Everything, on average, is similar to a standard diffusion process.

115

⇒ Modelling the average effect of turbulent mixing by an additional diffusion (turbulent diffusion) is legitimate.

❀ This averaged effect of convection by turbulent agitation (turbulent transport) is then called turbulent diffusion.

❀ Again, the similarity with continuum mechanics is clear: what is called molecular diffusion is the effect, at macroscopic scale, of the mixing due to molecular agitation.

� The diffusive model must be calibrated in order to reproduce the average effect of turbulent mixing.

❀ Equation of pure convection:

dφ * dt = ∂φ * ∂t + u * k ∂φ * ∂x k = 0
❀ Let us denote τ the time necessary for the large-scale eddies, of size �, to transport the dye over the distance �.

❀ Order of magnitudes:

φ * τ � u φ * � ⇒ τ � � u
where u and � are the integral scales (velocity and length), i.e., the characteristic scales of the largest eddies.

117

❀ Decomposition: φ * = Φ + φ and u * k = u k (mean velocity= 0)

❀ ∂φ * ∂t + u * k ∂φ * ∂x k = ∂Φ ∂t + ∂φ ∂t + u k ∂Φ ∂x k + u k ∂φ ∂x k = 0
❀ On average:

✄ ∂Φ ∂t + ∂φ ∂t + u k ∂Φ ∂x k + u k ∂φ ∂x k = 0 ✄ ∂Φ ∂t + ∂φ ∂t + u k ∂Φ ∂x k + u k ∂φ ∂x k = 0 (linearity) ✄ ∂Φ ∂t + ∂φ ∂t + ∂u k Φ ∂x k + ∂u k φ ∂x k = 0 because ∂u k φ ∂x k = u k ∂φ ∂x k + φ ∂u k ∂x k and ∂u k ∂x k = 0 (incompressibility) ✄ ∂Φ ∂t + ∂φ ∂t + ∂u k Φ ∂x k + ∂u k φ ∂x k = 0 (commutation) ✄ ∂Φ ∂t + ∂u k φ ∂x k = 0 because Φ = Φ, φ = 0, u k Φ = u k Φ = 0
❀ It is seen that the averaged effect of turbulent convection is reflected by the presence of the second moment u k φ.

❀ Since the eddies are not solved, we try to represent the mean mixing effect by a diffusion:

❀ On average, we want the diffusion model to spread the dye spot over the same distance � during the same time τ as the turbulent mixing does.

119

❀ Equation of pure diffusion:

✄ ∂Φ ∂t = ∂ ∂x k � α t ∂Φ ∂x k � ✄ Order of magnitudes: Φ τ � α t Φ � 2 ⇒ τ � � 2 α t
⇒ It can be seen that, for the diffusion model to correctly represent the mean mixing effect, we must have α t � �u ❀ Remark: in this case of a passive scalar, the turbulent diffusivity only depends on the characteristics of the turbulent velocity field.

� Applying these arguments to the momentum (per unit volume) ρu * i rather than the scalar φ * , the same conclusions are drawn:

❀ The mixing of momentum is due to the mean effect of the convection of momentum by the fluctuating velocity:

u j ∂ρu i ∂x j = ∂ρu i u j ∂x j
❀ It is the term involved in the mean momentum equation (Reynolds-stress tensor).

❀ We model this mean effect of fluctuating convection by a diffusion:

∂ ∂x k � ν t ∂U i ∂x k � ❀
U ∂U ∂x + V ∂U ∂y = - 1 ρ ∂P ∂x + ν ∂ 2 U ∂y 2 -
) is then v = v * -V = v * > 0
❀ It is assumed that it remains unchanged for a while and then is mixed with the ambient fluid after some distance � m called the mixing length.

In other words, its momentum is preserved, and after a distance � m , the particle exchanges its momentum with the ambient fluid.

123

❀ Just before the moment when it is mixed, its instantaneous velocity is then unchanged: u * = U (y)

❀ Its fluctuating velocity (with respect to the local mean velocity) is then:

u = u * -U (y + � m ) = U (y) -U (y + � m ) � -� m ∂U ∂y < 0
❀ u is negative when the particle moves upwards (v > 0), because ∂U ∂y > 0. It can be easily seen that u > 0 when it moves downwards (v < 0). Therefore, uv < 0, and, on average, uv < 0 (more accurately, uv is always of a sign opposite to that of ∂U ∂y ).

❀ We have seen above that the turbulent viscosity is ν t = �u, where � is the size of the largest eddies. It can be observed experimentally that � m � � � κy, where κ is called Karman's constant. Therefore

ν t = � 2 m � � � ∂U ∂y � � � ⇒ uv = -� 2 m � � � ∂U ∂y � � � ∂U ∂y ( 17
)
❀ These arguments are similar to those of the kinetic theory of gazes. The mixing length is analog to the mean free path.

� To be valid whatever the orientation of the shear, this relation must be generalized as

ν t = � 2 m � 2S ij S ij
� This idea has proven very useful not only for boundary layers but also for free shear flows (jets, mixing layers, wakes), for which the mixing length can be related to the thickness of the shear layer. ⇒ This is also an eddy-viscosity model. But this is the first 1-equation model.

It is still necessary to prescribe the length scale l.

❀ Kolmogorov additionally proposes to evaluate the length scale by

� = √ k ω
where ω is a characteristic frequency for the large eddies.

He proposes to solve a transport equation for ω ⇒ this is still an eddy-viscosity model. But this is the first 2-equation model.

It is based on the scales k and ω.

1945 Chou writes the transport equations for the Reynolds stresses and shows that the pressure term can be decomposed into a rapid and a slow part

φ * ij = - 1 ρ u i ∂p ∂x j - 1 ρ u j ∂p ∂x i = B ij ���� φ 1 ij + (A ijml + A jiml ) ∂U l ∂x m � �� � φ 2 ij
⇒ this is the first step towards Reynolds stress models.

129 1951 Rotta proposes the first model of the Reynolds-stress transport equations.

❀ In particular, he proposes the following model for the slow term, using a simple "return to isotropy" model:

φ l ij = -C 1 εa ij where a ij = u i u j k - 2 3
δ ij is the anisotropy tensor.

❀ Based on the physical role of this term observed in the return-to-isotropy experiment of an initially anisotropic turbulence:

When the strain disappears, experimental observations show that the turbulence returns to an isotropic state (u i u j = 2 3 kδ ij , i.e., a ij = 0), under the effect of the only non-zero source terms (slow term and dissipation term).

❀ In order to reproduce this return-to-isotropy effect, the simplest way is to model the slow term such a way that it depends linearly on the anisotropy. The ε factor is necessary for dimensional reasons.

❀ This model for the slow term, known as the Rotta model, is still used today.

❀ In association to models for the rapid term, it is the basis of the Reynolds-stress models available in industrial codes.

Jones and Launder

❀ Using previous work of Davidov (1961), Harlow and Nakayama (1967) and [START_REF] Hanjalić | Two-dimensional asymetrical turbulent flow in ducts[END_REF]: two-equation model using ε as the second variable (k-ε model).

❀ The final form of the model with the set of coefficients will be given by [START_REF] Launder | The numerical computation of turbulent flows[END_REF]. (making use of the previous work of Donaldson (1968), [START_REF] Daly | Transport Equations in Turbulence[END_REF], [START_REF] Hanjalić | A Reynolds stress model of turbulence and its application to thin shear flows[END_REF] and Naot, Shavit and Wolfshtein (1973))

❀ This is a 7-equation model (u i u j + ε).

⇒ The foundations are laid: the rest of the story consists in a continuous improvement of these approaches.

General principles for modeling

� Some principles can be applied in modeling (=development of turbulence models).

The use of such principles has very significant advantages:

❀ They provide a mathematical framework for modeling.

❀ They give physical constraints that enable the modeler to make relevant choices.

❀ They thus provide a method.

� However:

❀ They are not sufficient for ensuring that the models will "work properly": that is the reason why modeling is still an active field of research.

❀ They are not necessarily compatible with each other! ❀ Some of them can be considered more important than others, but it is impossible to rank them: different modelers would give different rankings.

Here is a list of the most common principles, that will be detailed in the following of this chapter: ❀ this principle simply requires that a model must consist in a closed system of equations ⇒ there must be as many equations as unknowns.

Dimensional homogeneity

� This principle means that the model used for a given quantity must be of the same dimension as the quantity itself.

� Example:

ν t = � 2 m � 2S ij S ij
Both sides of this relation are indeed of the dimension

L 2 T -1 .
� This simple example is obvious, but when formalizing a bit more, this basic principle can be easily forgotten.

R = [u i u j ]: ρ R = F (µ, ρ, U, ∇U)
❀ The objectivity principle implies that U cannot be involved in this relation.

❀ The fact that the Reynolds number is high implies that µ is negligible and cannot be involved either. ❀ These scales play exactly the same role as the state variables (P , T , U , S) in thermodynamics: they characterize of the state of the turbulence at a given point and a given time.

139

� Moreover, dimensional analysis provides information about the shape of the relation.

❀ Coming back to the previous example: let us introduce the state variables k and ε in the constitutive relation:

ρ R = F (ρ, ∇U, k, ε) (18)
❀ One of the fundamental principles of physics consists in the independence of the relation of the system of units.

❀ Let us choose the following system of units:

L = � ; T = τ ; M = ρ � 3
where � and τ are respectively the integral scale and the turnover time of the largest eddies. They are linked to k and ε by � = k 3/2 /ε and τ = k/ε. ❀ In this system, the relation must remain the same:

ρ * R * = F (ρ * , ∇ * U * , k * , ε * )
where the * indicate that the variables are expressed in the new system of units.

❀ This yields

ρ * = ρ M L -3 = 1 ; R * = R L 2 T -2 = R k ; ∇ * U * = ∇U T -1 = k∇U ε ; k * = k L 2 T -2 = 1 ; ε * = ε L 2 T -3 = 1 which gives R k = F � 1, k∇U ε , 1, 1 � ❀ Let us define the function G(X) = F (1, X, 1, 1)
❀ It is seen that relation ( 18) is reduced to a relation between two non-dimensional quantities

R k = G � k∇U ε � (this is Vaschy's proof of Vaschy-Buckingham's Π theorem).
⇒ Dimensional analysis thus constraints the form of the constitutive relations.

Completeness

A model is complete if it does not require any a priori knowledge of the flow. In other words: the formulation of the model must be the same for every flow.

� Example: the mixing length model is not complete.

❀ Indeed, we have seen that in order to close the system we had to prescribe the mixing length as a function of the location in the flow.

❀ The mixing length is dependent on the local characteristics of the flow: it is necessary to have information about the flow.

� Similarly, the one-equation model of Prandtl giving the turbulent viscosity as ❀ Reynolds-stress models (second-moment closures)

ν t = C µ √ k � is
3. For Reynolds-stress models, k can be evaluated from the Reynolds stress. It is only necessary to solve a transport equation for a second scale (ε, ω, etc.).

Objectivity

� Reminder about continuum mechanics:

Objective quantity = it does not depend on the observer (i.e., if it is independent on the reference frame in which it is measured)

� More formally: change of reference frame = translation + rotation, both functions of time:

x i (t) = R ij (t)x � j (t) + T i (t)
x i (t): coordinate vector of a point M in the frame of the first observer

x � i (t): coordinate vector of the same point in the frame of the second observer � A quantity is objective if and only if it satisfies the following relations:

❀ For a scalar s : s = s � ❀ For a vector V :

V i = R ij V � j ❀ For a tensor M : M ij = R ik M � kl R -1 lj 145
� Examples:

❀ The mean velocity U is not objective: for a second observer in translation compared to the first observer (no rotation), the velocity is

U i = dx i dt = dx � i dt + dT dt = U � i + dT dt
❀ The mean-velocity-gradient tensor ∇U is not objective.

❀ The strain rate tensor S = 1 2 (∇U + t ∇U) is objective.

❀ The rotation rate tensor Ω = 1 2 (∇U -t ∇U) is not objective.

❀ The fluctuating velocities are objective:

u i = u * i -U i .
❀ As a consequence, the Reynolds-stress tensor R = u i u j is objective.

� Consequence for modeling: the model representing an objective quantity must be objective.

� Example:

We have used in section 4.2, about Wang's constitutive relation, the fact that the mean velocity vector U cannot enter the constitutive relation.

Indeed, since U is not objective, the Reynolds stress tensor would not be objective if U was involved.

� On the contrary, the symmetric part S of the velocity gradient can be involved, since it is objective. ❀ Similarly, by changing the origin of the spatial coordinates, it can be shown that the constitutive relation does not explicitly depend on x.

❀ The properties of a material do not depend either on the Galilean or non-Galilean character of the reference frame: they are not sensitive to acceleration or rotation.

149

� Example: the model of Newton for (Newtonian) fluids is the constitutive relation which gives the stress tensor as a function of the strain tensor:

σ * ij (x, t) = -p * (x, t) δ ij + 2µ s * ij (x, t)
This relation satisfies the material frame indifference: it does not involve explicitly the time t, the point x or the rotation of the reference frame. This is equivalent to say that ❀ the fluid does not deteriorate with time (its properties are constant with time);

❀ the fluid is homogeneous (its properties are constant in space);

❀ the properties of the fluid do not depend on the rotation:

✄ This is true because the molecular agitation has a very short characteristic length scale (mean free path � 70 nm) compared to the scale of the problem.

✄ Therefore, the Coriolis force cannot play a role.

✄ At macroscopic scale, the properties of the fluid do not depend on the rotation of the reference frame.

� On the contrary, as concerns the constitutive relation for the Reynolds stress, i.e., for the "turbulent material", this principle cannot be applied.

❀ Indeed, the transport equations for the Reynolds stresses are modified in a non-Galilean reference frame: the effect of the Coriolis forces is significant (for instance in meteorology, turbomachinery, etc.).

❀ This effect cannot be neglected because there is no scale separation: the large turbulent structure are at the same scale as the problem (scale of the mean flow variations).

Large eddies are at the same scale as the mean recirculation region. There is a continuum of scales from large to small scales ⇒ no scale separation.

151

⇒ The constitutive relation for the Reynolds stresses must therefore explicitly involve the rotation rate of the reference frame with respect to a Galilean reference frame.

System rotation vector ω i

System rotation tensor Ω * ij = ε mji ω m ε ijk is the Levy-Civita symbol:

ε ijk = � +1 if (i, j, k) = (1, 2, 3), (2, 3, 1) or (3, 1, 2) -1 if (i, j, k) = (3, 2, 1), (2, 1, 3) or (1, 3, 2) 0 otherwise
� Question: how can we model an objective quantity (the Reynolds stress tensor R ij = u i u j ) using the non-objective quantity ω i ? � A long controversy has opposed specialists about this issue. A consensus has only been reached recently (1999).

� For instance: there is no way to build an objective quantity using S and ω only.

� ∇U + t ∇U � S ij = 1 2 � ∂U i ∂x j + ∂U j ∂x i � Skew-symmetric part (non-objective): Ω = 1 2 � ∇U -t ∇U � Ω ij = 1 2 � ∂U i ∂x j - ∂U j ∂x i � �
The relation is objective if and only if it can be written in the form

R = f (S, W)
where

W ij = Ω ij + ε mji ω m is

Realizability

The moments of order 1 (e.g., U i ), of order 2 (e.g., u i u j ), etc., have some mathematical properties.

� We have considered that the instantaneous velocities are random variables.

� We are trying to model the moments of these random variables (particularly of order 2, the Reynolds stress tensor R = u i u j ).

� A model can represent these moments of a random variable (i.e., is realizable) only if satisfies some properties. The Reynolds stress tensor being composed of the second moment of the velocity u * i , it satisfies the realizability conditions. In particular, the Reynolds stress tensor is a positive-semidefinite tensor.

∀ n i , � n 1 n 2 n 3 �   u 2 uv uw uv v 2 vw uw vw w 2     n 1 n 2 n 3   = R ij n i n j ≥ 0
It is often what is meant by realizability in turbulence modeling.

However, numerous properties can be applied to higher moments.

Realizability of the Reynolds stress tensor

The positivity of the Reynolds stress tensor can be written in different ways:

� The 3 eigenvalues are positive.

� The principal minors are positive:

u 2 ≥ 0 uv 2 ≤ u 2 v 2 u 2 v 2 w 2 + 2 uv uw vw -v 2 uw 2 -w 2 uv 2 -u 2 vw 2 ≥ 0
(others can be obtained by permutation:

v 2 ≥ 0, w 2 ≥ 0, uw 2 ≤ u 2 w 2 , vw 2 ≤ v 2 w 2 )
� The 3 principal invariants are positive (coefficients of the characteristic equation):

I R 1 ≥ 0 ; I R 2 ≥ 0 ; I R 3 ≥ 0 157
� These properties can be reformulated in terms of the invariants of the anisotropy tensor.

� The main interest of this approach is that the invariants of the anisotropy tensor can be more easily related to the structure of turbulence.

� Two slightly different definitions are used:

❀ the American school defines the anisotropy tensor as:

b ij = u i u j 2k - 1 3 δ ij
❀ and the European school:

a ij = u i u j k - 2 3 δ ij
� The principal invariants of the tensor b are denoted as I, II and III (A 1 , A 2 and A 3 for the European school); {} denotes the trace � The curves H = -1 and H = 1 correspond to an axisymmetric turbulence (2 eigenvalues are equal). H = -1 when the third eigenvalue is less than the other two and H = 1 when the third eigenvalue is greater than the other two. � At the center of the channel, turbulence is close to an isotropic state.

I = {b} = 0 ; II = 1 2 ({b} 2 - � b 2 � ) = - � b 2 � 2 ; III = 1 6 � {b} 3 -3 {b} � b 2 � + 2 � b 3 �� = � b 3 � � The
� When the wall is approached, close to the peak of production, the u 2 component is much larger than the other two ⇒ the curve then approaches the tip of the triangle, which corresponds to the one-component limit.

161

� At the wall, the fluctuating velocity component normal to the wall (v here) tends to zero faster than the other two components (two-component limit)

⇒ when y + → 0, the curve goes towards the two-component limite line, i.e., the line G = 0.

� The parameter G = 1 9 + 3 III + II is an invariant which characterizes of a 2-component limit of turbulence. It is called Lumley's flatness parameter.

� This parameter can be used in modeling to sensitize the model to the presence of a wall or a free surface.

� The invariant H can be used to characterize an axisymmetric state of turbulence.

du i u j dt = P ij + D ν ij + D T ij + φ * ij -ε ij
� For a model to be realizable, u 2 ≥ 0 must for instance be ensured under all circumstances.

❀ This constraint implies that u 2 = 0 must be a minimum of u 2

❀ If u 2 =
❀ If u 2 = 0 then du 2 dt > 0 (sufficient condition).
❀ The drawback of this constraint is that it prevents the turbulence to approach a two-component state.

❀ This condition should be called strong realizability condition, since the condition is actually restrictive.

� The reasoning can be applied to any quantity that must remain positive, as for instance u 2 v 2uv 2 , or the invariants of the Reynolds stress, etc.

As a practical matter, writing transport equations for these quantities can be very complicated.

Conclusions about realizability

� Realizability gives a number of relations that guide the modeling.

� However, ensuring realizability under all circumstances is very difficult.

� As a practical matter, it can only be imposed in some simplified situation: for instance, by considering homogeneous turbulence, the transport terms vanish and one can focus on the terms φ * ij and ε ij .

165

Concerning the importance of imposing the realizability, there are two schools:

� The first school considers that realizability is crucial:

❀ Realizability ensures that no non-physical solutions, like negative energies, can be observed.

This risk exists in situations when one of the components is considerably damped: close to a wall, in the case of a strong stratification (stable) or of strong rotation.

❀ Even in situations when such solutions are not observed, realizability is useful for numerical stability:

Computations being often initialized by completely arbitrary values, it is possible to go through negative energies, which is a source of divergence of the computation.

Realizability is then a numerical robustness factor.

❀ Realizability is a tool that helps guiding the modelers.

❀ This is mathematically more elegant.

� The second school considers that realizability cannot be the main guide for modeling:

❀ Realizability leads to complex models, often nonlinear ⇒ the gain in terms of numerical robustness is compensated by the loss of linearity.

❀ Situations where non-realizable solutions are observed are extreme situations, being almost pathological cases.

❀ The appearance of non-realizable solutions can be avoided by a posteriori corrections of the models. For instance:

Universality

Ultimate aim of modeling: write a model that can be applied to all flows (or at least to a wide class of flows: for instance, single-phase flows of incompressible fluids).

� Utopian: Reynolds averaging = huge loss of information ❀ trying to build a model able to correctly reproduce all the phenomena, in all situations, is hopeless.

❀ Consensus today: modeling will never reach this goal.

❀ What can be expected from a model?

✄ Being applicable to all the situations.

✄ Giving results that are satisfactory for the widest possible range of flows.

✄ Having limitations that are known.

� Many research programs, in particular European programs, aim at applying numerous models to numerous test cases, in order to evaluate their potential, their shortcomings and possibly evaluate uncertainties.

� A recent and fashionable research topic is UQ (Uncertainty Quantification): the uncertainty due to the model is difficult to evaluate.

� However, a weak form of universality must be satisfied:

❀ The coefficients of a model must be fixed after a calibration procedure in simple flows and must thereafter be considered as "universal", i.e., applicable to any flow.

❀ Modifying the coefficients of a model to fit experimental results is not admissible.

This kind of tuning has been very detrimental to the reputation of CFD. Historically, first models (mixing length models) were developed for boundary layer flows, in particular for aeronautics.

� New models are expected to provide boundary layer results at least as good as those provided by existing models.

U + = y + U + = 1 κ ln y + + C DNS of Moser et al. at Re τ = 590.
� There are not many theories in turbulence: the log layer theory (Prandtl, Izakson, Millikan) is one the most solid of them.

� It will be seen that the Karman constant κ (slope in the lin-log plot) is one of the constraints used to calibrate the models.

Numerical robustness

� As introduced at the beginning of this course, the turbulence models are not supposed to stay confined in the research field, but to be used in industry.

� The numerical difficulties faced when trying to integrate the system of equations (difficulties to converge, and even divergence) are often underestimated by modelers.

� A compromise must be found between two very important constraints:

❀ a good representation of the physics;

❀ numerical robustness and ease of use.

� Example: in order to model the slow part of the redistribution term, Speziale, Sarkar and Gatski (1991) propose the following nonlinear form:

φ l ij = - � C 1 + C * 1 P ε � εb ij + C 2 ε � b ik b kj - 1 3 b kl b kl δ ij � ❀
The nonlinear term is necessary to correctly reproduce the return to isotropy experiment of an initially anisotropic, homogeneous turbulence (for both strong and weak anisotropies).

❀ This term has however a limited influence in most of the common flows.

❀ Since it introduces a nonlinearity which make the numerical resolution more difficult, this term is often dropped.

� Warning: It is often possible to improve the numerical scheme ⇒ simplifying the models in not always the better solution.

∂U i ∂x i = 0 ; ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j
� A first moment closure = constitutive relation (model) which expresses the second moments (u i u j ) as a function of the first moments.

� We are seeking a relation in the form:

R = f (ν, ρ, P, U, β n )
where f is any functional, which can involve variables at different points and different times.

� The scalar quantities β n are state variables which globally characterize the state of the turbulence (k, ε, ...).

General assumptions

Different remarks and assumptions make simplifications possible:

� The mean pressure does not appear in the exact Reynolds-stress equations ⇒ there is no reason to take it into account in the model.

� One can assume that f only involve the present instant t (assumption of instantaneity):

⇒ arbitrarily, the possibility of a memory effect of turbulence (which exists), is excluded.

� At high Reynolds number, the effects of viscosity are negligible: the global diffusive properties of turbulence are linked to the larger eddies of the flow (energetic structures) and, when the Reynolds number ν t /ν is sufficiently high, are independent of the molecular viscosity.

⇒ ν can then be removed from the relation for the moment. One must reconsider this issue for a low-Reynolds number turbulence.

175

� As usual in continuum mechanics, assumption of locality:

❀ The relation only involves the neighborhood of point x.

❀ For instance, in the constitutive relation, the dependance on a the scalar field

β i (x � ), ∀x � ∈ Ω
is replaced by a dependance on the derivatives of the variable evaluated at point x

∂ k β n ∂x n k (x), k ≥ 0
� Most of the time, derivatives of order ≥ 2 are neglected (standard assumption of a materially simple medium)

⇒ one considers that the variables and their first derivatives are sufficient to describe the properties of turbulence.

� In particular: the mean velocity gradient ∂U i ∂x j must be retained in the models.

� Most of the time, gradients of scalar quantities β n (k, ε, ...) are also neglected.

� Objectivity constraint (the modeled R must be objective) ⇒ U cannot enter the relation.

� As explained before: objectivity + sensitivity to system rotation ❀ the constitutive relation cannot be a funtion of ∂U i /∂x j (not objective)

❀ rather a function of:

✄ the strain rate tensor:

S ij = 1 2 � ∂U i ∂x j + ∂U j ∂x i �
✄ the absolute rotation rate tensor W:

W ij = 1 2 � ∂U i ∂x j - ∂U j ∂x i � + ε mji ω m
This constraints and assumptions lead to the following formulation:

R(x, t) = f � ρ, S(x, t), W(x, t), β n (x, t) � 177

Dimensional analysis

� We know that the diffusive properties of turbulence must be described by two characteristic scales

ν t = �u = � 2 τ -1
� � and τ are the integral scale of turbulence and the turnover time of the large eddies.

� Consequently, in order to correctly represent the diffusive effect of turbulence, � and τ must be introduced as state variables.

R = f (ρ, S, W, �, τ )
� Dimensional analysis imposes that this relation must be valid whatever the system of units.

� Let us arbitrarily choose the following system of units:

L = � ; T = τ ; M = ρ � 3
� In this new system of units, the relation is

R * = f (1, S * , W * , 1 , 1)
where

R * = R � 2 τ -2 , S * = τ S and W * = τ W. � Defining g(x, y) = f (1, x, y, 1, 1) this relation reduces to R * = g (S * , W * ) ( 19
)
R * = g 0 I + g 1 S * + g 2 (S * W * -W * S * ) + g 3 � S * 2 - 1 3 � S * 2 � I � + g 4 � W * 2 - 1 3 � W * 2 � I � + g 5 � W * S * 2 -S * 2 W * � + g 6 � W * 2 S * + S * W * 2 - 2 3 � S * W * 2 � I � + g 7 � W * S * W * 2 -W * 2 S * W * � + g 8 � S * W * S * 2 -S * 2 W * S * � + g 9 � W * 2 S * 2 + S * 2 W * 2 - 2 3 � S * 2 W * 2 � I � + g 10 � W * S * 2 W * 2 -W * 2 S * 2 W * �
� The coefficients g i are not constant, but rather functions that can only depend on the six independent invariants formed from the tensors S * and W * :

� S * 2 � , � W * 2 � , � S * 3 � , � W * 2 S * � , � S * 2 W * 2 � , � S * W * S * 2 W * 2 �
, where {.} denotes the trace.

� Despite the numerous (very restrictive) hypotheses (instantaneity, locality, etc..) the relation is still very complex, involving fifth order nonlinear terms.

� Fortunately, the Cayley-Hamilton theorem limits the degree of each tensors in this relation to 2.

� Such a complex relation is almost inapplicable:

❀ The model would be very difficult to numerically integrate.

❀ The coefficients g i are not provided by the theory. Even though they are often considered constant, the calibration process would be very difficult: how to find values of the 10 coefficients from simple flows?

⇒ the nonlinear terms of highest degree must be "sacrificed".

181

� However, there is a very recent renewal of nonlinear models � Therefore, in some commercial codes, cubic models are proposed.

❀
� The constitutive relation is then truncated by taking g 7 to g 10 equal to zero.

183

Quadratic models

� In a statistically 2D flow (averaged quantities independent of z, for instance), the last six term (g 5 -g 10 ) can be shown dependent on the first five terms (g 0 -g 4 ).

� The relation thus simplifies as

R * = g 0 I + g 1 S * + g 2 (S * W * -W * S * ) + g 3 � S * 2 - 1 3 � S * 2 � I �
� This relation is the exact solution for Eq. ( 19) in 2D flows.

� It is also used in 3D flows, as simplified models.

� We have

u i u j = g 0 � 2 τ 2 δ ij +g 1 � 2 τ S ij +g 2 � 2 (S ik W kj -W ik S kj )+g 3 � 2 � S ik S kj - 1 3 S kl S kl δ ij �
� This type of models is also included in some commercial codes.

Linear models

� Finally, one can decide, for simplicity and robustness reasons, to only keep the linear terms, which leads to

u i u j = g 0 � 2 τ 2 δ ij + g 1 � 2 τ S ij � The trace is 2k = 3g 0 � 2 τ 2 since k = 1 2 u i u i , δ ii = 3 and S ii = 0 such that g 0 = 2 3 and � 2 τ 2 = k are used � Remember that ν t ∝ �u = � 2 τ It is usual to note the proportionality coefficient C µ with C µ = - g 1 2 such that we have u i u j = 2 3 kδ ij -2ν t S ij
� This is the Boussinesq relation.

185

� It is seen that this relation (which was initially introduced by a simple phenomenological reasoning) is the consequence of a long series of simplifying assumptions:

❀ instantaneity: no influence of the history of the strain and the turbulence;

❀ locality: the turbulence is only influenced by its vicinity;

❀ materially simple medium;

❀ linearity of the constitutive relation.

� It is thus clear that linear eddy-viscosity models cannot perform well in all the situations.

Comments about the linear models � The constitutive relation can also be written as

a = -C µ k ε S
⇒ proportionality between the anisotropy tensor and the strain tensor ⇒ the eigenaxes of the anisotropy tensor are aligned with those of the strain tensor ⇒ their eigenvalues are proportional. All this properties are wrong in general. 189

[u i u j ] =                               � u i u j Bouss � =                              
❀ Boussinesq constitutive relation: well adapted to reproducing the shear stresses uv and vw.

These stresses are produced by the wake and the mixing layer, respectively:

P 12 = -v 2 ∂U ∂y ; P 23 = -v 2 ∂W ∂y ( 20
)
❀ On the contrary: not able to distinguish among the three diagonal components.

Indeed, when the plane strain is zero, it predicts:

u 2 = v 2 = w 2 = 2 3 k (21)
❀ The uw component is produced by the wake/mixing layer interaction

P 13 = -uv ∂W ∂y -vw ∂U ∂y ( 22
)
✄ It is simply zero with the Boussinesq relation.

✄ A nonlinear constitutive relation is then necessary.

❀ Second moment closures, which have exact production terms, are of course able to reproduce the correct Reynolds stress profiles.

❀ The following figure shows the comparison between results given by several eddy-viscosity models: linear models and a cubic model.

⇒ only the nonlinear eddy-viscosity model is able to reproduce the turbulence anisotropy and the component due to the wake/mixing layer interaction.

191

Comparison between linear models and the cubic model of Merci and Dick.

Test case: wake/mixing layer interaction (Experiment from Béharelle, Delville, Bonnet, 1997).

Reproduced from [START_REF] Manceau | Elliptic Blending Model: A New Near-Wall Reynolds-Stress Turbulence Closure[END_REF].

Contra-rotating jets: PROFILES

Eddy-viscosity models

Location: 0004 (x = 150 mm)

• Experiments LEA/k-ε Ghent-U/Low-Re-k-ε-Yang-Shih EC-Nantes/k-ω-BSL EC-Nantes/k-ω-SRSST Ghent-U/cubic-k-ε -0,22 -0,11 0

One-equation models

� The constitutive relations that have been derived, either linear or nonlinear, involve two independant turbulent scales, characterizing the energetic eddies.

� The particular choice of the scales does not matter a priori:

❀ The integral length scale � ❀ The eddy-turnover time τ ❀ The turbulent energy k = � 2 τ 2 ❀ The dissipation rate ε = � 2 τ 3 ❀ The eddy-viscosity ν t = � 2 τ ❀ The specific dissipation ω = ε k = 1 τ ❀ Any other combination ϕ = � a τ b
∂U i ∂x i = 0 ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P * ∂x i + ∂ ∂x j � � ν + ν t � � ∂U i ∂x j + ∂U j ∂x i

��

where

P * = P + 2 3 ρ k .
❀ Solving this system of equations gives U i and P * .

⇒ k is not needed to obtain the mean velocity.

⇒ A complete model can be obtained by writing a transport equation for ν t .

❀ P * is NOT the pressure, but since k = 0 at solid boundaries, P * = P at the surface of the obstacles ⇒ the forces (drag, lift) can be computed.

� The only model using a transport equation for ν t used in industrial applications is the Spalart-Allmaras model (1992).

❀ ν t is given by ν t = f v1 � ν, and � ν is obtained through a transport equation:

∂� ν ∂t + U j ∂� ν ∂x j = c b1 � S� ν -c w1 f w � � ν d � 2 + 1 σ ∂ ∂x k � � ν + � ν � ∂� ν ∂x k � + c b2 σ ∂� ν ∂x k ∂� ν ∂x k
❀ The details about the coefficients c ? and the functions f ? can be found in the original article or, for instance, in the book of Wilcox (1998).

❀ The functions f ? are in particular introduced such that the model is valid in the near-wall region.

195

❀ The model is used a lot in aeronautics, for different reasons:

✄ It is easy to integrate numerically.

✄ For attached flows, it produces results as good as those given by zero equations models (Cebeci-Smith, Baldwin-Lomax).

✄ For detached flows, it gives a much better description of the flow field than zero equation models.

❀ However, the model is too simple (only one equation) to be valid in a wide range of flows.

❀ For instance, it give wrong results in jet flows (see Wilcox, 1998), overestimating the expansion rate by more than 40%.

k-ε models

� Two reasons for using two scales and not only one:

❀ One wants to get rid of the limitations of the Spalart-Allmaras model.

❀ One uses a nonlinear constitutive relation:

✄ For instance, the quadratic constitutive relation is:

u i u j = 2 3 � 2 τ 2 δ ij + g 1 � 2 τ S ij + g 2 � 2 (S ik W kj -W ik S kj ) + g 3 � 2 � S ik S kj - 1 3 S kl S kl δ ij � ✄
The two scales � and τ cannot be combined in only one scale.

� As said earlier, various scales can be used instead of � and τ .

197

� The most standard choice nowadays is the determine � and τ from the turbulent energy k and its dissipation rate ε (k-ε models):

� = k 3/2 ε ; τ = k ε
� Another rather common choice: turbulent energy k and specific dissipation rate

ω = ε/k (k-ω models): � = k 1/2 ω ; τ = 1 ω 5.3.1. Transport equations � The exact transport equation for k is ∂k ∂t + U k ∂k ∂x k = -u i u k ∂U i ∂x k -ε + ∂ ∂x k � ν ∂k ∂x k � + ∂ ∂x k � 1 ρ pu k - 1 2 u i u i u k � �
The right hand side involves:

❀ sources terms (production and dissipation, the latter being actually always a sink)

❀ transport terms: viscous diffusion, pressure "diffusion" and turbulent transport.

199

Source terms and transport terms � The transport equations are all in the symbolic form

∂f ∂t + ∂U k f ∂x k = ∂ ∂x k F k + s
where f can be a scalar, a vector or a second order tensor.

� The terms to be modeled are all on the right hand side.

200

� Integrating this equation over a control volume yields

∂ ∂t � Ω f dV = - � ∂Ω U k f n k dS + � ∂Ω F k n k dS + � Ω s dV
⇒ The total variation of f over Ω is due to two types of terms:

❀ Transport terms through the surface ∂Ω by the fluxes U k f and F k :

These fluxes are of convective or diffusive type:

✄ they do not create nor destroy the quantity f inside the volume Ω, contrary to source terms ✄ but rather transport f through the surface Modeling of the source terms of k

Ω U k f F k
∂k ∂t + U k ∂k ∂x k = -u i u k ∂U i ∂x k -ε + ∂ ∂x k � ν ∂k ∂x k � + ∂ ∂x k � 1 ρ pu k - 1 2 u i u i u k � � Production term P = -u i u k ∂U i ∂x k :
❀ The model for P is simply obtained by introducing the constitutive relation (model for u i u j ).

❀ For instance, linear constitutive relation (Boussinesq relation)

u i u j = -2ν t S ij + 2 3 k δ ij ⇒ P = 2ν t S ij S ij
� Two terms require modeling:

❀ ∂ ∂x k � 1 ρ pu k - 1 2 u i u i u k � (so-called turbulent diffusion) ❀ ε (dissipation) 203
Modeling of the transport terms � The two turbulent terms pu k and -u i u i u k are generally modeled as a whole.

� It is assumed that the effect of these terms is diffusive ⇒ similar to the viscous term: 1

ρ pu k -u i u i u k = ν t σ k ∂k ∂x k where σ k is a Prandtl number.

This model is called the Simple Gradient Diffusion Hypothesis (SGDH).

� This type of model would be perfectly justified:

❀ if the turbulent agitation had a behavior of "pure agitation" as the molecular agitation ❀ if this agitation were at a scale very small compared to the mean flow scale (in other words, if turbulence had characteristics similar to those of the molecular agitation).

� Actually, this hypothesis can be criticized:

❀ There is no scale separation in turbulence ❀ There are coherent structures (not a pure agitation) ❀ Pressure "diffusion" does not always behave the same way as the turbulent diffusion:

✄ In some cases pressure diffusion transport energy in the opposite direction ⇒ It can be anti-diffusive ✄ Some authors have proposed models specific to this term.

� Using the SGDH model, the modeled k equation is

∂k ∂t + U k ∂k ∂x k = P -ε + ∂ ∂x k �� ν + ν t σ k � ∂k ∂x k � 205
The ε-equation

� ε is defined as ε = ν ∂u i ∂x k ∂u i ∂x k
� The exact transport equation for ε is very complex.

� Individual modeling of the different terms has never been successful.

� A phenomenological approach is retained:

The ε-equation must involve the same terms as the k-equation: convection, source term (production), sink term (dissipation), viscous diffusion, turbulent diffusion:

∂ε ∂t + U k ∂ε ∂x k = P ε -� ε + ∂ ∂x k �� ν + ν t σ ε � ∂ε ∂x k � 206
� The source/sink terms P ε and � ε must be modeled.

❀ It is proposed to related them to the corresponding terms in the k-equation, i.e.,

P and ε.

❀ A simple dimensional analysis provides a relation between the production and dissipation terms of ε to those of k:

P ε = C ε1 ε k P ; � ε = C ε2 ε k ε 207 �
The modeled equation for ε is finally written as

∂ε ∂t + U k ∂ε ∂x k = C ε1 ε k P -C ε2 ε 2 k + ∂ ∂x k �� ν + ν t σ ε � ∂ε ∂x k � �
u i u j = -2 ν t S ij + 2 3 kδ ij ν t = C µ k 2 ε ; P = 2ν t S ij S ij ∂k ∂t + U k ∂k ∂x k = P -ε + ∂ ∂x k �� ν + ν t σ k � ∂k ∂x k � ∂ε ∂t + U k ∂ε ∂x k = C ε1 ε k P -C ε2 ε 2 k + ∂ ∂x k �� ν + ν t σ ε � ∂ε ∂x k �
and calibrated the coefficients:

C µ = 0.09 ; σ k = 1 ; C ε1 = 1.44 ; C ε2 = 1.92 ; σ ε = 1.3
� These equations + these coefficients = the standard k-ε model. � Grid turbulence:

U 0 U 0 x y z x * y * z * x * y * z * x * = x -U 0 t
❀ Turbulence produced by a uniform flow through a grid and convected downstream ❀ No turbulence production downstream of the grid (S ij = 0) ⇒ turbulent energy decays along the x-axis (∂k/∂x < 0).

❀ Moving reference frame x * = x -U 0 t ⇒ ∂k/∂t < 0.
❀ Homogeneous turbulence is obtained locally (in a small region), i.e., all the spatial derivative ∂/∂x * i of turbulent quantities (statistical moments of order ≥ 2) are zero

❀ In this situation ∂k ∂t = -ε ; ∂ε ∂t = -C ε2 ε 2 k ❀ Solution: k = k 0 � 1 + t t 0 � -γ ε = ε 0 � 1 + t t 0 � -γ-1
where

t 0 = γ k 0 ε 0 and γ = 1 (C ε2 -1)
❀ Experimental data indicate that the turbulent energy decreases approximately like a power 1.1 of time ⇒ C ε2 = 1.92 is chosen.

� Homogeneous sheared turbulence

❀ A crucial property expected from a model is to reproduce shear flows, i.e., the response of turbulence to a mean shear rate.

❀ The simplest situation is sheared homogeneous turbulence.

U0 U (y) x y z x * y * z * x * y * z * x * = x -U0t
❀ Similar to standard grid turbulence, but the grid is not regular in the y-direction.

❀ Variation along y of pressure drop ⇒ streamwise velocity U is a function of y.

❀ In order to reproduce an homogeneous turbulence, a constant velocity gradient ∂U/∂y must be ensured, such that the production term P = -u i u j ∂U i ∂x j is independent of y.

❀ In that case, production is not zero, and the equations reduce to

∂k ∂t = P -ε ; ∂ε ∂t = C ε1 ε k P -C ε2 ε 2 k ❀ Define S = � 2S ij S ij ; η = Sk ε ; t * = tS
❀ Leads to the equation for η :

∂η ∂t * = 1 ε ∂k ∂t - k ε 2 ∂ε ∂t = -C µ (C ε1 -1) η 2 + (C ε2 -1)
❀ In the case P > ε, k and ε go to ∞ when t → ∞, but experiments show that P ε = C µ η 2 goes to a constant value (weak equilibrium).

❀ According to the equation for η, this value can only be

P ε = C µ η 2 = C ε2 -1 C ε1 -1 ( 23
)

213

❀ This result shows the strong link between this coefficient and the response of turbulence to a shear.

❀ Unfortunately, there is no unique calibration for this coefficient:

✄ Different experiments can be carried out for different intensities for the shear strain S.

✄ Different S correspond to different limiting values for P/ε.

✄ Each value of S then leads to a different calibrated coefficient C ε1 .

✄ Launder and Spalding actually chose the value C ε1 = 1.44, in order to reproduce the growth of a mixing layer.

✄ A slightly different value would be obtained using other type of flows for calibration (jets, wakes, ...).

✄ The case of C ε1 illustrates that calibration is, in general, a matter of compromise.

� Logarithmic region in channel flows and boundary layers ❀ Eliminating the velocity gradient between this two relations, and using the fact that in this region P = ε, leads to

U (y) U (x,
C µ = uv 2 k 2
❀ Experiments: in the log region uv/k � -0.3 ⇒ the standard value of C µ is 0.09 ❀ Moreover, using the relations valid in the log layer

P + = ε + = 1 κy + and k + = 1 C 1/2 µ uv + = 1 C 1/2 µ
the transport equation for ε gives � Choice of ε as a second scale: less obvious.

σ ε = κ 2 C 1/2 µ (C ε2 -C
❀ Actually, numerous models use a different scale:

k-ω, k-ω 2 , k-kl, k-ν t , k-ε/ √ k, etc.
❀ Any combination ϕ of � and τ can be chosen, as soon as ϕ is independent of k (it cannot be a power of k)

� If the same procedure is followed as for the k-ε model, a transport equation for ϕ is written

∂ϕ ∂t + U k ∂ϕ ∂x k = C ϕ1 ϕ k P -C ϕ2 ϕ k ε + ∂ ∂x k �� ν + ν t σ ϕ � ∂ϕ ∂x k � � Writing ϕ = k a ε b
, and applying the same calibration procedure as for the coefficients of the k-ε, the following relations are obtained:

C ϕ1 = a + b + b C µ γ η 2 ∞ ; C ϕ2 = a + b + b γ ; σ ϕ = b C 1/2 µ γ η 2 ∞ κ 2 C µ η 2 ∞ -1
� The parameters are given by the experiments:

C µ = 0.09 ; γ = 1.1 ; η ∞ = 4.8 ; κ = 0.41 219 � For instance, choosing the variable ϕ = ω 2 = ε 2 k 2 leads to C ϕ1 = 0.88 ; C ϕ2 = 1.82 ; σ ϕ = 2.38
� There is clearly a problem for some variables: for instance, considering

ϕ = � = k 3/2 ε , one gets C ϕ1 = 0.06 ; C ϕ2 = -0.41 ; σ ϕ = -1.19
The coefficient of the diffusion term is negative: the transport equation for � able to reproduce the three cases used for the calibration is counter-diffusive! � Among all the variables that lead to a positive σ ϕ coefficient, Launder and Spalding chose ε because the coefficient σ ε gives, without modifications, realistic expansion rates for free shear flows. This is not the case, for instance, of the variable ω 2 .

u i u j = -2 ν t S ij + 2 3 kδ ij ν t = k ω ; P = 2ν t S ij S ij ∂k ∂t + U k ∂k ∂x k = P -β � kω + ∂ ∂x k �� ν + ν t σ k � ∂k ∂x k � ∂ω ∂t + U k ∂ω ∂x k = α ω k P -βω 2 + ∂ ∂x k �� ν + ν t σ ω � ∂ω ∂x k � β � = 0.09 ; σ k = 2 ; α = 5/9 ; β = 3/40 ; σ ω = 2 221
� This model has some advantages compared to the standard k-ε model, but also major drawbacks.

� Advantages:

❀ The model is integrable down to the wall. This is a major advantage.

❀ It gives better results than the standard k-ε model in adverse pressure gradients situations ⇒ it predicts better the separation location.

❀ The ω equation is actually independent on k, since its production term reduces to

α ω k P = 2αS ij S ij
This feature is very interesting, since it uncouples the ω equation from the k equation, which is a factor of numerical robustness.

❀ It is possible to show that only the powers of ω satisfy this feature. The choice of ω is then optimal from a numerical point of view.

� Drawbacks:

❀ The ω variable tends to infinity at the wall

ω(y) ∼ y→0 6ν βy 2
which is the analytical solution of the ω equation in the near-wall region.

❀ As a practical matter, [START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] proposed the Dirichlet boundary condition ω = 10 6ν

βy 2 1
where y 1 is the distance between the center of the cell adjacent to the wall and the wall. ⇒ This model is seldom applied: however, it is used as a near-wall model associated to the k-ε model far from the wall (k-ω SST model, described later).

Advantages and limitations of eddy-viscosity models

Three parts:

� General indications concerning all the eddy-viscosity models.

� Focus on the standard k-ε model.

� Finally, on nonlinear eddy-viscosity models.

Eddy-viscosity models in general

� Advantages: can be described in one sentence:

They contain less differential equations than Reynolds stress models.

❀ This means that they are less expensive in terms of CPU time ❀ However, the difference is not always large (an implicit CFD code spends most of the time solving the velocity-pressure system).

❀ They are not necessarily more numerically stable than Reynolds stress models: stability depends on the linearity or nonlinearity of the model.

225

� Limitations:

❀ Eddy-viscosity models assume a relation between the Reynolds stress and the mean field (through S and W): this is not true in general ⇒ the prediction of complex flows, in particular 3D flows, can be problematic.

❀ This relation is assumed instantaneous ⇒ the response to a change in the strain is instantaneous: this is wrong.

Actually, turbulence is sensitive to the history of the strain (memory effect).

❀ The transport equation of the second scale (often ε) is very empirical: many people brought it into question, but very few have proposed other approaches.

❀ The models predict in general an expansion rate of axisymmetric jets larger than the expansion rate of plane jets, while experiments evidence exactly the opposite:

✄ this shortcoming is known as the round jet/plane jet anomaly.

✄ Corrections have been proposed (e.g., [START_REF] Pope | An explanation of the turbulent round-jet/plane-jet anomaly[END_REF], but are not widely used.

❀ Models based on the second scale ε have difficulties in reproducing the effects due to adverse pressure gradients.

In particular, this implies that separation points are not always well predicted. � Advantages:

❀ The model is easy to implement and numerically very stable.

❀ It correctly reproduces the shear stress in free shear flows (for instance uv produced by the shear ∂U/∂y).

❀ Its behavior is somehow predictable: the model makes always the same type of errors ⇒ without being able to give error bars, such systematic errors can be taken into account for the interpretation of the results.

❀ It is, despite its simplicity, applicable to much wider range of situation than the Spalart-Allmaras model.

� Limitations:

❀ Linearity of the model ⇒ virtually no chance to correctly reproduce complex flows, in particular 3D flows (see, for instance, the wake/mixing layer interaction).

❀ Proportionality between anisotropy and strain ⇒ the production is always predicted positive (P = 2ν t S ij S ij ).

❀ Not sensitive at all to the effects due to the streamline curvatures and a solid body rotation: : ❀ Finally, the standard k-ε model is not valid in near-wall regions. ❀ The stagnation point anomaly can be avoided.

✄
❀ Some nonlinear models are also valid in the near-wall region.

❀ Numerically less robust than linear models.

❀ The calibration process is more difficult: there are more coefficients ⇒ are difficult to isolate considering simple flows.

❀ They have a behavior much less predictable than the standard k-ε model: they can give very good results for some flows, and, for other flows, have a completely wrong behavior.

235

Correction for the stagnation point anomaly

Correcting this anomaly is crucial, such that numerous corrections were proposed.

� Durbin's realizability constraint [START_REF] Durbin | On the k-3 stagnation point anomaly[END_REF] remarks that the stagnation point anomaly can be related to a realizability problem.

❀ Eigenaxes of the strain rate S ij used as the coordinate axes.

The tensor S ij is then diagonal:

    λ 1 0 0 0 λ 2 0 0 0 λ 3     (λ i 's = eignevalues)
In this coordinate frame: the three diagonal components of the Reynolds stress (normal stresses), given by the Boussinesq relation are

∀α, u 2 α = -2C µ kT λ α + 2 3 k
where T = k/ε has been defined.

It can be shown that if T ≤ 1

C µ √ 6 S , i.e.
, if this bound is imposed on T :

T = min � k ε ; 1 C µ √ 6 S �
and the eddy viscosity is evaluated by

ν t = C µ kT instead of ν t = C µ k k ε , then
the realizability constraint u 2 α ≥ 0 is satisfied.

237

❀ This upper limit imposed to the turbulent time scale corrects the stagnation point anomaly: indeed, production becomes

P = min � 2C µ k 2 ε S 2 ; 2 √ 6 C µ k S �
The second term is linear in S: the strong overestimation of production in the regions of strong plane strain is thus avoided (see the figure of § 5.5.2).

Imposing a realizability constraint can then have surprising effects! ❀ The Cµ coefficient is a bit more complex, but as a summary, in a plane strain situation, it reduces to

Cµ = 1 A0 + A S k ε S
❀ When the strain is strong, the second term of the denominator is dominant, and production becomes

P = 2 √ 2 3 k S
which is in the same form as in the case of Durbin's correction.

❀ This model has become very popular since it was implemented in the commercial code Fluent, under the name

Realizable k-ε model.

❀ It is now available in other codes as well.

239

� The RNG k-ε model ❀ By a recursive procedure, it is possible to tend to the Reynolds decomposition

⇒ k-ε equations.
❀ Exactly the equations of the standard k-ε model, but:

✄ the coefficients are obtained directly from the theory.

✄ Unfortunately, it appeared that with that set of coefficients the model does not

give very good results in the log layer.

❀ Modification by [START_REF] Yakhot | Development of turbulence models for shear flows by a double expansion technique[END_REF]:

✄ Different coefficients ✄ Additional term (-R) in the ε equation: R = C µ η 3 (1 -η/η 0 ) 1 + βη 3 ε 2 k with η = √ 2S k/ε.
❀ This is the so-called RNG model, for Re-Normalization Group.

❀ Adding that term is strictly equivalent to introducing a variable C ε1 coefficient in the form: ❀ However, as shown by [START_REF] Laurence | Turbulence Modelling[END_REF], the improvement of the predictions in stagnation regions is due to the variable C ε1 coefficient, which strongly increases the dissipation, instead of diminishing the production ⇒ the model improves the predictions because of a compensation of errors.

C * ε1 = C ε1 - η(1 -η/η 0 ) 1 + βν 3

Reynolds-stress models (second moment closures)

Numerous limitations of the eddy-viscosity models ⇒ going to Reynolds-stress models have decisive advantages:

� Equations for the Reynolds stresses are solved:

❀ the models reproduce the memory effect of turbulence ❀ for instance, in the case of a sudden change in the mean strain, turbulence respond with a time delay.

� By increasing the modeling order, much is gained in the representation of the physics:

❀ instead of assuming that the Reynolds stresses have a given behavior (determined by the constitutive relation), their transport equations are solved ❀ these equations contain the main physical mechanisms that drive turbulence: production, redistribution, turbulent transport, viscous diffusion, dissipation.

� In particular, the production terms, that are sufficient to explain many phenomena, do not require modeling! ⇒ Reynolds-stress models represent a major breakthrough compared to eddy-viscosity models. 243

Reynolds-stress transport equation

∂u i u j ∂t + U k ∂u i u j ∂x k � �� � C ij = -u i u k ∂U j ∂x k -u j u k ∂U i ∂x k � �� � P ij + 1 ρ p � ∂u i ∂x j + ∂u j ∂x i � � �� � φ ij - 1 ρ ∂ ∂x k (u i p δ jk + u j p δ ik ) � �� � D p ij -2ν ∂u i ∂x k ∂u j ∂x k � �� � ε ij +ν ∂ 2 u i u j ∂x k ∂x k � �� � D ν ij - ∂u i u j u k ∂x k � �� � D T ij ( 24
)
C ij Convection (by the mean flow)

P ij Production � φ ij Pressure-strain or redistribution term � D p ij Pressure diffusion term � ε ij Dissipation D ν ij Molecular or viscous diffusion � D T ij

Turbulent transport or turbulent diffusion

Only the checked terms require modeling.

� The pressure term

φ * ij = - 1 ρ u i ∂p ∂x j - 1 ρ u j ∂p ∂x i has

been decomposed into two parts:

❀ Pressure-strain correlation (or redistribution):

φ ij = 1 ρ p � ∂u i ∂x j + ∂u j ∂x i � that is conveniently traceless: φ ii = 0
Since the equation for the turbulent energy k is obtained as half the trace of the equation for u i u j , this term disappears from the k-equation: it does not create and it does not destroy turbulent energy ⇒ it can only redistribute turbulent energy among the components of the Reynolds-stress tensor u i u j .

❀ Pressure diffusion: � The name "diffusion" can hide the fact that these terms have no reason to be purely diffusive, even if, in many cases, it is what is observed. 247

D p ij = - 1 ρ ∂ ∂x k (u i p δ jk + u j p δ ik ) that

The redistribution term

� Term that plays the most important role after production.

� Since the production does not require any modeling, it is essentially the redistribution term that characterizes the models.

Chou's theory (1945)

� Taking the divergence of the equation for the fluctuating velocity u i ⇒ Poisson equation:

∇ 2 p = ∂ 2 p ∂x k ∂x k = -2ρ ∂u k ∂x l ∂U l ∂x k -ρ ∂ 2 (u k u l -u k u l ) ∂x k ∂x l
� Theory of Chou: infinite domain (no wall).

� It will be seen later how to extend the theory to region close to the wall.

� Solution of the Poisson equation:

p(x) = - � IR 3 � -2ρ ∂u k ∂x l (x � ) ∂U l ∂x k (x � ) -ρ ∂ 2 (u k u l -u k u l ) ∂x k ∂x l (x � ) � 1 4π�x � -x� dV (x � ) � Involves the Green function 1 4π�x � -x� of the IR 3 space.
� This relation reminds us that the Navier-Stokes equations are integro-differential equations.

� Pressure fluctuations at a given point depends on the mean and fluctuating velocities in the whole domain ⇒ turbulence is non-local.

249

The redistribution term can then be written as

φ ij (x) = 1 ρ p(x) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � = � IR 3 ∂ 2 u k u l ∂x k ∂x l (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� � �� � φ 1 ij + � IR 3 2 ∂U l ∂x k (x � ) ∂u k ∂x l (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� � �� � φ 2 ij
� φ 2 ij depends on the mean velocity gradients.

❀ It responds instantaneously to their variations.

❀ It is called the rapid term.

� φ 1 ij is independent on the mean velocity.

❀ This term does not respond instantaneously to a sudden variation of the mean strain ❀ This is an turbulence-turbulence interaction term.

❀ Consequently it responds only when the turbulence has been modified by another mechanism.

❀ For this reason, it is called the slow term.

251

� Quasi-homogeneity hypothesis:

❀ In the rapid term φ 2 ij : (mean velocity gradient) × (two-point correlations).

❀ Two-point correlations go to zero when separation �x � -x� increases

L �x � -x� f (x, x � ) Two-point correlation function
❀ The correlation length scale L characterizes the distance over which the two-point correlations are "non-zero".

❀ If the velocity gradients evolve slowly compared to this correlation scale, it can be assumed that in the integral, the velocity gradients are constant.

⇒ This is called the quasi-homogeneity hypothesis.

❀ [START_REF] Bradshaw | On Local Approximations of the Pressure-Strain Term in Turbulence Models[END_REF] have shown, from DNS data in a channel flow, that this hypothesis is perfectly justified, except in the near-wall region.

❀ We will see later how to account for the near-wall effects.

❀ The quasi-homogeneity then leads to:

φ 2 ij (x) = ∂U l ∂x k (x) � IR 3 2 ∂u k ∂x l (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� = ∂U l ∂x k (x) � � IR 3 2 ∂u k ∂x l (x � ) ∂u i ∂x j (x) dV (x � ) 4π�x � -x� � �� � A ijkl (x) + � IR 3 2 ∂u k ∂x l (x � ) ∂u j ∂x i (x) dV (x � ) 4π�x � -x� � �� � A jikl (x)
� 253 � Locality hypothesis:

❀ We have now:

✄ A slow term, independant of the mean velocity: φ

1 ij (x) = B ij (x)
✄ A rapid term, function of the mean velocity:

φ 2 ij (x) = [A ijkl (x) + A jikl (x)] ∂U l ∂x k (x)
❀ The locality hypothesis consists in assuming that A ijml and B ij can be modeled by an algebraic relation which only depends on quantities evaluated at point x.

✄ According to their exact equations, these terms are actually non-local.

✄ This hypothesis has to be brought into question close to the wall, in particular. � Experiments show that the role of the slow term is essentially to make the turbulence tend to an isotropic state by redistributing the energy between the components (experiment: homogeneous, initially anisotropic, decaying turbulence)

� Indeed, this term being traceless, it does not play any role in the budget of turbulent energy k.

255

� We thus seek a model in the form

Φ 1 = f (R, β i )
� Similar to the case of eddy-viscosity models, we introduce 2 turbulent scales, � and τ

Φ 1 = f (R, �, τ )
� Dimensional analysis:

❀ This relation must be valid whatever the system of units.

❀ Let us chose

� = k 3/2 ε ; τ = k ε ❀ We get Φ 1 ε = f � R k , 1, 1 � ❀ This relation simplifies as Φ 1 ε = g � R k � by introducing g(X) = f (X, 1, 1)

257

� The theory of invariants (see for instance [START_REF] Lumley | Computational Modeling of Turbulent Flows[END_REF] shows that the most general relation writes

Φ 1 ε = α I + β b + γ b 2 where b = R 2k - 1 3
δ ij is the anistropy tensor (isotropic turbulence corresponds to b = 0). � α, β and γ are functions of are the principal invariants of b

I = {b} = 0 ; II = - � b 2 � 2 ; III = � b 3 � 3 �
The use of the Cayley-Hamilton theorem appears clearly here: all the powers of b higher than 2 can be reduced.

� trace of Φ 1 is zero ⇒ � Φ 1 � = 3α -2γII = 0
� The most general relation then is

Φ 1 = β εb + γ ε � b 2 + 2 3 II I � ( 25
)
� It is relatively simple.

� But it can still be simplified:

❀ Assuming γ = 0 and β = cste = -C 1 ⇒ Rotta's model.
This is the only linear model.

❀ Assuming γ = 0 and a variable β (function of the invariants): so-called quasi-linear model (linear tensorial form, but nonlinear coefficients)

❀ The nonlinear term can also be used and the coefficient considered constant:

these models are simply called non-linear models with constant coefficients.

259

� Example of quasi-linear model:

❀ Rotta's model Φ 1 = -C 1 εb is not realizable.
❀ In order to ensure realizability, it is sufficient to make the coefficient a function of the invariants.

❀ Fu, [START_REF] Fu | Modelling strongly swirling recirculating jet flow with Reynolds-stress transport closures[END_REF] propose:

C 1 = 1 + 6.2 √ -18 G II,
where G is Lumley's flatness parameter G = 1 9 + 3 III + II.

❀ Even though the form of (25) must be satisfied, there is wide flexibility in the choice of the form of the coefficients.

� Example of nonlinear model with constant coefficients:

❀ Return-to-isotropy experiments in homogeneous turbulence show that φ 1 ij does not vary linearly with a ij .

⇒ Some authors prefer calling φ 1 ij the resistance-to-strong-anisotropy term rather than the return-to-isotropy term.

❀ That is the reason why, for instance, Speziale, Sarkar and Gatski (1991) use the nonlinear model

β = -3.4 ; γ = 4.2
❀ The use of constant coefficients implies that the model is not realizable.

� Example of nonlinear model:

❀ Fully nonlinear models are seldom used. They are still confined to research. ❀ The most typical example of the nonlinear model of Craft and Launder (1996).

❀ The method followed to derive the model is:

✄ The starting point is the general relation [START_REF] Fadai-Ghotbi | Modélisation de la turbulence en situation instationnaire par approches URANS et hybrides RANS-LES. Prise en compte des effets de paroi par pondération elliptique[END_REF], associated with the general relation for the rapid term, which will be shown below.

✄ Realizability constraints are widely used: in particular, the model is derived in such a way that it behaves satisfactorily in two-component turbulence (near walls or near a free surface).

✄ That is the reason why the authors call it the TCL model, for Two-component Limit.

✄ The model has been derived for two decades ⇒ it has been tested against a wide range of flows, but has become so sophisticated and complex that it is cumbersome to handle.

❀ Here is the form of its slow term (expressed in terms of the "European' invariants):

Φ 1 = -(C 1 + √ A) ε a -C 1 C � 1 ε � a 2 - 1 3 A 2 I
� with:

C 1 = 3.1f A A 2d ; C � 1 = 1.1 f A =          A � 0.05 0.7 � 1/2 if A < 0.05 A 0.07 1/2 if 0.05 < A < 0.7 A if A > 0.7 A 2 = � a 2 � ; A 3 = � a 3 � ; A = 1 - 9 8 (A 2 -A 3 ) ; A 2d = min(A 2 ; 0.5) 263
6.2.3. Modeling of the rapid term � Theory of invariants (see [START_REF] Speziale | Modeling the pressure-strain correlation of turbulence: an invariant dynamical system approach[END_REF] shows that the general form is

φ 2 ij (x) = [A ijml (x) + A jiml (x)] ∂U l ∂x m (x) (26)
Φ 2 = β 1 kS + β 2 k � bS + Sb - 2 3 {bS} I � + β 3 k � b 2 S + Sb 2 - 2 3 � b 2 S � I � +β 4 k (bW + Wb) + β 5 k � b 2 W + Wb 2 � (27)
� Model of Speziale, Sarkar and Gatski (SSG)

❀ Speziale, Sarkar and Gatski (1991) studied solutions for short times (rapid distortion theory) and equilibrium solutions for long times (t → ∞) in homogeneous turbulent flows subjected to plane strains ❀ They showed these solutions can be reproduced using the quasi-linear formulation

Φ 2 = -g * 1 P b + � g 3 -g * 3 II 1/2 � kS +g 4 k � bS + Sb - 2 3 {bS} I � + g 5 k (bW + Wb)
❀ They associate this model with the slow term described in 6.2.2., which is

Φ 1 = -g 1 εb + g 2 ε � b 2 - 1 3 � b 2 � I � ❀ Coefficients: g 1 = 3.4 ; g 2 = 4.2 ; g * 1 = 1.8 ; g 3 = 0.8 ; g * 3 = 1.3 ; g 4 = 1.25 ; g 5 = 0.4

265

❀ Some comments about the SSG model:

✄ It has been applied to numerous case, with some success.

✄ In particular, it is well adapted to rotating flows.

✄ It is not valid in the near-wall region ⇒ it has to be used with wall functions.

✄ It is one of the most widely used Reynolds-stress model (with the Rotta+IP model).

✄ It is available in some industrial CFD codes.

❀ Remark: the nonlinearity of the slow term is often removed by taking g 2 = 0.

� Rotta+IP or LRR model ❀ Launder, Reece and Rodi (1975): associated the following model for the rapid term

φ 2 ij = -C 2 � P ij - 2 3 P δ ij � to Rotta's slow term φ 1 ij = -C 1 εa ij .
❀ P ij is the production tensor and P = 1 2 P ii is the production rate.

❀ The model is called the Rotta+IP model or LRR model.

❀ The model for the rapid term φ 2 ij (IP model, for Isotropization of Production):

✄ has a simple phenomenological interpretation.

✄ it can also be seen as a simplification of the SSG model: indeed, taking g 1 = 3.6 ; g 2 = 0 ;

g * 1 = 0 ; g 3 = 0.8 ; g * 3 = 0 ; g 4 = 1.2 ; g 5 = 1.2 gives the Rotta+IP model.

✄ Crow's constraint (rapid distortion theory) ⇒ C 2 = 3 5 .

✄ C 1 = 1.8 (calibration against the return-to-isotropy experiment).

267

❀ Phenomenological interpretation: tends to counteract the anisotropy of production.

❀ Consider the example of an isotropic turbulence subjected to a sudden strain (entrance of a distorted channel):

From Chassaing (2000)

❀ When the strain is applied at the entrance of the distorted channel, the φ 1 ij term is zero, since the turbulence is isotropic.

❀ On the contrary, since the production is suddenly non zero, the φ 2 ij term is activated.

❀ Exact equations:

∂u 2 ∂t = 0 + 2 ρ p ∂u ∂x -2ν ∂u ∂x k ∂u ∂x k ∂v 2 ∂t = - 4 3 D k + 2 ρ p ∂v ∂y -2ν ∂v ∂x k ∂v ∂x k ∂w 2 ∂t = 4 3 D k + 2 ρ p ∂w ∂z -2ν ∂w ∂x k ∂w ∂x k
(here we have used the fact that the turbulence is isotropic at the entrance of

the distorted channel ⇒ u 2 = v 2 = w 2 = 2 3 k).
❀ The Rotta+IP model then gives:

∂u 2 ∂t = 0 + 0 -ε 11 ∂v 2 ∂t = - 4 3 D k + 4 3 C 2 D k -ε 22 ∂w 2 ∂t = 4 3 D k - 4 3 C 2 D k -ε 33

269

❀ The IP model is active as soon as the turbulence is subjected to the strain.

❀ It instantaneously responds ⇒ it is indeed a "rapid" term.

❀ It contradicts the effect of production, which tends to make the turbulence anisotropic.

❀ The slow term is zero at the entrance of the distorted channel, but becomes significant later on, when the turbulence become anisotropic.

❀ At the end of the distorted channel, only the slow term can be responsible for the return to isotropy, since

S ij = 0 ⇒ φ 2 ij = 0.
❀ This is against this type of experiments that the C 1 coefficient is calibrated. ❀ These models are very nonlinear.

❀
❀ They have remained confined in the research field and are not implemented in industrial codes.

❀ Some of them, in particular the TCL model, have been successfully tested against a wide range of experimental data.

❀ A very good review of the existing Reynolds-stress models is given by Hanjalić and Launder (2011).

Problems in the near-wall region

� The SSG and Rotta+IP models are not valid in the near-wall region.

� The problem can be traced to the hypotheses used in the derivation of the model (locality, quasi-homogeneity) which must be brought in question in the near-wall region.

� The SSG model gives good results in the log layer, contrary to the Rotta+IP model. ❀ [START_REF] Gibson | Ground effects on pressure fluctuations in the atmospheric boundary layer[END_REF], (based on an idea due to [START_REF] Shir | A preliminary numerical study of atmospheric turbulent flows in the idealized planetary boundary layer[END_REF], proposed the inclusion of the following terms, active in the log layer:

φ 1p ij = C p 1 ε k � u k u m n k n m δ ij - 3 2 u k u i n k n j - 3 2 u k u j n k n i � f � L T d � φ 2p ij = C p 2 � φ 2 km n k n m δ ij - 3 2 φ 2 ki n k n j - 3 2 φ 2 kj n k n i � f � L T d � ❀ L T = k 3/2 /
ε is the turbulent length scale, d the distance between the point and the wall and n the unit vector normal to the wall.

❀ f is chosen such a way that f = 1 in the log layer and f = 0 in a free flow.

275

❀ These terms indeed correct the model in the log layer.

❀ They are generally proposed in CFD code in association with the Rotta+IP model.

❀ The original arguments associate theses terms to the wall echo (reflection of the fluctuating pressure on the wall) ⇒ these terms are called wall echo terms.

❀ Actually, these arguments are wrong (cf. [START_REF] Manceau | Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier-Stokes modelling[END_REF]): these terms actually represent the blocking effect.

Modeling of the turbulent transport terms

� The two transport terms due to turbulent fluctuations are:

❀ the triple correlation term (third order moments), also called turbulent diffusion

D T ij = - ∂u i u j u k ∂x k
❀ and the pressure diffusion term

D p ij = - 1 ρ ∂ ∂x k (u i p δ jk + u j p δ ik )
� Similar to the terms in the k-equation: it is standard to consider that this two terms have similar effects, of diffusive type, and to model them together, despite the remarks made above.

� The term to be modeled is then:

u i u j u k + 1 ρ (u i p δ jk + u j p δ ik )
� Historically, it has been considered that pressure diffusion is negligible, and the modeling has been based on measurements of u i u j u k and its exact equation.

u i u j u k = -C s k 2 ε ∂u i u j ∂x k
� This model can be criticized because it does not satisfy the interchangeability of the indices:

❀ i, j and k in the left hand side have symmetrical roles, contrary to the right hand side.

❀ However, the triple correlations only appear in the equations through their gradients. We have: � This model is too isotropic: when turbulence is anisotropic, the turbulent transport is enhanced in the direction in which the fluctuation are stronger.

∂u i u j u k ∂x k = - ∂ ∂x k � C s k 2 ε ∂u i u j ∂x k � ❀ In
u 2 v 2
� The model is used because it is numerically very stable. � It is often considered that the gain is not worth the numerical endeavor.

u i u j u k = -C s k ε u k u l ∂u i u j ∂x l � Model
u i u j u k = -C s k ε � u i u l ∂u j u k ∂x l + u j u l ∂u i u k ∂x l + u k u l ∂u i u j ∂x l � � Comes
� Other models, less complex [START_REF] Mellor | A survey of the mean turbulent field closure models[END_REF] or more complex [START_REF] Lumley | Computational Modeling of Turbulent Flows[END_REF][START_REF] Magnaudet | The modelling of inhomogeneous turbulence in the absence of mean velocity gradient[END_REF] exist, but are very seldom or never used. 281

Dissipation tensor modeling

� The dissipation of turbulent energy is due to the small scales.

� At high Reynolds number, the turbulent cascade theory shows that the small scales of the flow are rather isotropic.

� This remark suggests writing the dissipation tensor in an isotropic form, known as the Kolmogorov model:

ε ij = 2 3 εδ ij
where the 2/3 coefficient is introduced because ε = 1 2 ε ii and δ ii = 3 � It is only necessary to solve a transport equation for ε.

� The same equation as the standard k-ε model is used, except as regards the turbulent transport terms, in which a generalized gradient diffusion hypothesis is used:

∂ε ∂t + U k ∂ε ∂x k = C ε1 ε k P -C ε2 ε 2 k + ν ∂ 2 ε ∂x k ∂x k + ∂ ∂x k � C ε k ε u k u l ∂ε ∂x l � �
The isotropy of the small scales is not a valid hypothesis anymore when the turbulent Reynolds number Re t = ν t /ν is weak, in particular in the near wall region.

� Remark: similarly to the eddy-viscosity models, a transport equation for another quantity than ε can be used (in particular ω is sometimes used). � For complex flows, they have proved superior to eddy-viscosity models.

� However, commercial code developers still have to work to make these models as numerically stable as the eddy-viscosity models.

Explicit algebraic stress models

� As seen earlier, the Reynolds-stress models much better represent the physics than eddy-viscosity models.

� On the contrary, eddy-viscosity models have advantages as concerns simplicity and numerical robustness.

� From an industrial point of view, the best would be to combine the advantages of the two modeling levels, or at least to find a compromise.

� The method which is going to be presented is such a compromise.

� Using some hypotheses, it make the formal derivation of eddy-viscosity models from Reynolds-stress models possible.

� These models are called explicit algebraic (Reynolds) stress models.

du i u j dt = P ij + D ν ij + D T ij + D p ij � �� � D ij +φ ij -ε ij the transport equations of the anisotropy tensor b ij = u i u j 2k - 1 3 δ ij can be written: db ij dt = 1 2k du i u j dt - u i u j 2k 2 dk dt = 1 2k � P ij + D ij + φ ij -ε ij - u i u j k P - u i u j k D + u i u j k ε �
� First it is assumed that the anisotropy tensor is in equilibrium, i.e., that it is conserved on a streamline:

db ij dt = ∂b ij ∂t + U k ∂b ij ∂x k = 0
� This is wrong in general, but a very good approximation in some cases (boundary layer) and even exact in a channel flow.

� Then it is assumed that total diffusion (molecular, turbulent, pressure diffusion) have exactly the same anisotropy as the Reynolds stress:

D ij 2D - 1 3 δ ij = u i u j k - 1 3 δ ij
� This is a reasonable approximation, since diffusion is essentially due to turbulence fluctuations.

� Using these two approximations, the following, purely algebraic system (no partial derivatives) is obtained:

� P ij - u i u j k P � + φ ij - � ε ij - u i u j k ε � = 0 (28)
� One can use standard models for φ ij and ε ij , and obtain an algebraic system which directly gives the Reynolds stress.

� It is not necessary to solve differential equations for the Reynolds stress, only

for k and ε.

� This is what is called an algebraic model.

287

� Unfortunately, this tensorial equation (system of 6 independent equations) in nonlinear.

� Even using a linear Reynolds-stress model, since the Reynolds stresses appear in P , the product u i u j P is nonlinear.

� It is difficult to numerically solve this nonlinear system, which is very stiff.

� One can solve this equation implicitly: at each time step, the values of the u i u j at the previous time step can be used to simplify.

� Unfortunately, this idea also leads to very serious numerical problems, except if very small time steps are used.

Explicit algebraic methodology

� In order to avoid numerical problems, an analytical solution is sought.

� Example of the SSG model:

❀ Introducing the SSG model, and writing the algebraic relation [START_REF] Gibson | Ground effects on pressure fluctuations in the atmospheric boundary layer[END_REF] in a tensorial form gives

- 1 a 4 b -a 3 � bS + Sb - 2 3 {bS} I � + a 2 (bW -Wb) = a 1 S (29)
with

a 1 = 2 3 - 1 2 C 3 a 2 = 1 - C 5 2 a 3 = 1 - C 4 2 a 4 = gτ τ = k ε g = �� 1 + C * 1 2 � P ε - � 13 3 - C 1 2 � + 10 3 � -1
❀ The following identity has been used:

P = - 4 3 kS -2k (bS + Sb) + 2k (bW -Wb)

289

❀ We are going to express the anisotropy tensor in a tensorial basis which only depend on the mean flow.

❀ In our example, we choose a 3-term basis:

T 1 = S T 2 = SW -WS T 3 = S 2 -1 3 {S 2 } I
❀ In this basis, the anisotropy tensor can be written

b = α 1 T 1 + α 2 T 2 + α 3 T 3
where the α i 's are the projection coefficients.

❀ Remarks:

✄ The theory of invariants teaches us that ten tensors are necessary to exactly represent b.

✄ See in this course, the chapter about eddy-viscosity models: this are the ten tensors appearing in the most general form.

✄ Using a 10-term basis is intractable. We will restrict the basis to 3 tensors.

✄ Why 3? Because it is a basis in 2D flows. It is an acceptable approximation in 3D flows.

❀ A Galerkin projection of the algebraic system of equations ( 29) onto the tensorial basis is then performed, i.e.,

✄ The decomposition b = 3 � n=1 α n T n is introduced in the relation: 3 � n=1 α n � - 1 a 4 T n -a 3 � T n S + ST n - 2 3 {T n S} I � + a 2 (T n W -WT n ) � = a 1 S
✄ This relation is projected onto each tensor T m of the basis (the relation is multiplied by T m and the trace is taken).

✄ The three following relations (m = 1, 2, 3) are thus obtained:

3 � n=1 α n � - 1 a 4 {T n T m } -a 3 � {T n ST m } + {ST n T m } � +a 2 � {T n WT m } -{WT n T m } �� = a 1 {ST m }

291

� These 3 relations can be summarized in a matrix form

XA = Y where X =     α 1 α 2 α 3    
� Inverting this matrix analytically, an explicit expression of the α i 's is obtained, in the form

α i = f i � k, ε, P, � S 2 � , � W 2 ��
� This provides the anisotropy tensor in exactly the same general form as for a quadratic eddy-viscosity model

b = α 1 S + α 2 (SW -WS) + α 3 � S 2 - 1 3 � S 2 � I �
� This form is obtained because a 3-term basis was used. Other forms (linear, cubic, etc.) can be obtained with other bases.

� The coefficients are functions of the invariants k, ε, P , � S 2 � and � W 2 � .

Remarks

� This method does not require calibration: the expression for the α i 's is directly inherited from the underlying Reynolds-stress model.

� The variations of the coefficients as functions of the invariants preserves the reproduction of the same physical mechanisms as the underlying

Reynolds-stress model.

� But there are additional hypotheses: equilibrium of the anisotropy tensor and alignment of the diffusion tensor on the Reynolds-stress tensor, and the projection on a reduced tensorial basis is an approximation.

� In particular, the mechanisms linked to the production term are still present.

293

� It is seen that k and ε appear in the α i 's: their transport equations must be solved.

� If a 1-term basis had been used, S, we would have recovered a linear k-ε model:

b = - ν t k S.
The main difference is that ν t would have been a function of the invariants.

� The production P appears in the α i 's: since P = -2α 1 k � S 2 � , the relation for α 1 is actually nonlinear (at least cubic).

Conclusion � A two-equation model has been derived (k and ε).

� Its properties are close to those of the original Reynolds-stress model.

� This is a good compromise between representation of the physics and numerical robustness.

� This type of approaches is a very active field of research.

� If the underlying Reynolds-stress model is not valid in the near-wall region, the explicit algebraic stress model is not valid either (no miracle). � The near-wall region is a region where this is not true anymore � A large majority of the model presented so far is not valid in the near wall region.

� In this region, some physical phenomena take place such that the hypotheses upon which the models are based are not valid.

⇒ A particular treatment in the near-wall region is necessary. ❀ The wall-echo effect tends to enhance this phenomenon.

❀ Historically, it was erroneously considered that wall echo is at the origin of the decrease of redistribution in the near-wall region (this is actually due to wall blocking).

� Blocking effect:

❀ Any fluctuation in the direction of the wall generates a zone of high pressure that tends to slow down the flow in the wall-normal direction and to deviate the velocity towards the directions parallel to the wall.

❀ This effet leads to a selective damping of the wall-normal fluctuation.

❀ This blockage becomes sufficiently strong in the vicinity of the wall to invert the transfer of energy:

energy is transferred from the v 2 component to the others (essentially w 2 ) below

y + = 12.
High pressure � These two non-viscous effects are non-local (the wall is felt at distance by turbulence).

� The wall-echo effect is negligible compared to the blocking effect.

� All these effects in the vicinity of the wall bring into question:

❀ The quasi-homogeneity: very strong variations of the mean velocity gradient.

❀ The locality: the non-viscous effects are non-local.

❀ The high-Reynolds-number hypothesis: the lack of separation between the large and the small scales implies that viscosity is not negligible. � This hypotheses are very strict, they actually never apply. � However, the conclusions are fundamental: the velocity profile, the Reynolds stresses and the dissipation, made non-dimensional by the friction velocity and the viscosity, have a behavior which does not depend on the Reynolds number ⇒ they are somehow "universal".

U + = y + U + = 1 κ ln y + + C � In the viscous sub-layer: U u τ = yu τ ν (U + = y + ) � In the log layer: U u τ = 1 κ ln yu τ ν + C (U + = 1 κ ln y + + C) ε = u 3 τ κy ; uv = -u 2 τ ; u i u j =cste � For details, see Appendix A.
� In real flows (Re � = ∞, curve walls, pressure gradient, tridimensionality, separation, etc.) it has to be considered that ❀ either the "universal" behavior is perturbed ❀ or it is completely modified (at separation points, for instance) 303

Taylor series expansion at the wall

� The Taylor series expansion of the mean and fluctuating velocities at the wall

(y → 0) is                  U = A 1 y + A 2 y 2 + O(y 3 ) u = a 1 y + a 2 y 2 + O(y 3 ) v = b 2 y 2 + O(y 3 ) w = c 1 y + c 2 y 2 + O(y 3 ) p = p 0 + p 1 y + p 2 y 2 + O(y 3 )
� The simplifications have been applied:

       U (0) = 0 (no-slip condition) u i (0) (no-slip condition) ∂v ∂y (0) = 0 (incompressibility)
� Warning: these relations are only valid very close to the wall (typically, until y + � 5, i.e., in the viscous sublayer).

� The behavior of the Reynolds stresses can be deduced

                     u 2 = a 2 1 y 2 + 2a 1 a 2 y 3 + O(y 4 ) v 2 = b 2 2 y 4 + 2b 2 b 3 y 5 + O(y 6 ) w 2 = c 2 1 y 2 + 2c 1 c 2 y 3 + O(y 4 ) uv = a 1 b 2 y 3 + (a 2 b 2 + a 1 b 3 ) y 4 + O(y 5 ) uw = a 1 c 1 y 2 + (a 2 c 1 + a 1 c 2 ) y 3 + O(y 4 ) vw = c 1 b 2 y 3 + (c 2 b 2 + c 1 b 3 ) y 4 + O(y 5 ) (30)
� For k :

k = 1 2 u i u i = 1 2 (a 2 1 + c 2 1 ) y 2 + (a 1 a 2 + c 1 c 2 ) y 3 + O(y 4 ) ( 31
)
305 � Very important conclusions:

❀ Turbulence decreases as y 2 .

❀ The components involving v once decrease as y 3 (uv and vw).

❀ The component v 2 decreases as y 4 .

� Wall turbulence then tends to a very anisotropic state, called Two-component limit, i.e., v 2 << u 2 and v 2 << w 2 . � Even though the asymptotic behavior is only valid in the viscous sublayer, the non-locality of turbulence implies that this limiting state influences a much larger part of the near-wall region, up to the log layer.

Consequences for modeling

� Turbulence models are based on hypotheses which are not applicable in the near-wall region.

� There are several possibilities to circumvent this problem:

❀ using wall function which avoid the resolution of the near-wall region;

❀ introducing functions depending on the distance to the wall or the turbulent Reynolds number Re t in order to "force" the model to behave correctly ⇒ low-Reynolds number models;

❀ prescribing a turbulent scale in the near-wall region ⇒ two-layer models;

❀ reconsidering the hypotheses used in the derivation of the models, in particular the locality ⇒ elliptic relaxation; � Nowadays, the computing power is sufficient, but in a complex geometry, it is extremely difficult, if not impossible at all, to satisfy the condition y + � 1 everywhere (at best, generating the mesh is very long).

❀
� Wall function strategy:

❀ Any high-Reynolds number model is used (standard k-ε, RNG k-ε, Rotta+IP, SSG, etc.).
❀ The discretization point closest to the wall is placed in the log layer.

❀ The exact boundary conditions (U = 0, k = 0, u i u j = 0, ...) are not used, but rather wall functions, based on the "universal" laws obtained in the log layer.

❀ This artificial boundary conditions correspond to a slip wall with friction.

8.6.2. One-scale approach � The following relations can be used:

❀ Mean velocity (assumed parallel to the wall):

U u τ = 1 κ ln yu τ ν + C ( 32
)
❀ Turbulent energy and dissipation (example of a k-ε model):

ε = u 3 τ κy ; k = C -1/2 µ u 2 τ
� Warning: these relations are given for a wall located in y = 0. If the wall has a different orientation, a rotation of the axes is necessary.

311

� In these relations, the unknown u τ must be determined during the computation.

� It is obtained using an iterative method:

❀ In equation [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF], the value of U from the previous iteration at the discretization point closest to the wall is used.

❀ Equation ( 32) is solved by a iterative method (of the type of the Newton algorithm) to obtain the value of u τ .

� This u τ is used to evaluate the friction at the wall (in finite volumes for instance, this is what is necessary to solve the momentum equation):

µ ∂U ∂y = τ w = ρu 2 τ � This u τ is used in ε = u 3 τ κy and k = C -1/2 µ u 2 τ .
8.6.3. Two-scale approach � A major weakness of the previous method appears at separation or reattachement points (where U = 0):

❀ the friction velocity is zero;

⇒ ε and k are zero: this is not correct (turbulence does not disappear at this points).

� Two-scale approach: � They substantially reduce of the number of grid points.

❀
� They make the mesh generation process easier.

� They are quite robust: even though their range of applicability is in theory limited, they can be used all the time (this does not guarantee the quality of the results).

� It does not mean that the number of cells is reduced by a factor of 30.

❀ Indeed, meshes are progressive in the near-wall region: usually, a geometric expansion of the layers of cells is used:

Δ i = rΔ i-1
, where r is the expansion ratio.

❀ We want to mesh the near-wall region up to some distance d with a progressive mesh, and then the external region with a constant mesh.

Δ1 Δn

Near-wall mesh

Δext

External mesh 315 ❀ We want to choose the two thicknesses Δ 1 and Δ ext . The expansion ratio r and the number of cell layers n are imposed by the two relations

Δ ext = Δ 1 r n and d = Δ 1 r n -1 r -1
where d is the distance covered by the progessive mesh.

❀ Discretization point closest to the wall located at y + c = 1 means � The location of the discretization point closest to the wall must be very well controlled: it must be located in the log layer, i.e., between y + = 30 and about one tenth of the boundary layer thickness ⇒ the adequate size of the grid cells can be evaluated a priori, but must always be checked a posteriori (by evaluating u τ from the results).

Δ + 1 = Δ 1 u τ ν � 2 which implies Δ 1 � 2 ν u τ • y c Δ 1 ❀ Using the rule of thumb u τ � 5%U ∞ , we have Δ 1 � 2 ν 0.05U ∞ � d 0.025Re ❀ We have r = 1 + Δ ext d - 1
� There are elaborate wall functions which ❀ account for the influence of pressure gradients;

❀ ensure a more or less correct behavior if the discretization point closest to the wall is below ✄ If the mesh is locally coarse, AWF tend to standard wall functions.

y + =
✄ If the mesh is locally fine, AWF tend to exact boundary conditions.

✄ If the mesh is in between, AWF empirically blend these two limiting cases.

321

Low-Reynolds number models using damping functions

� In order to avoid the use of wall functions, the models must be integrated down to the wall.

� For this purpose, their incorrect behavior in the near-wall region must be corrected.

� The simplest way is to "force" the model to reproduce the data (experimental or from DNS) by introducing damping functions.

� For instance, in order to correct the standard k-ε model, a f µ function is generally introduced in the turbulent viscosity: 1981), for instance, propose:

ν t = C µ f µ k 2 ε � Lam & Bremhorst (
f µ = � 1 -exp � -0, 0165 √ ky ν �� 2 � 1 + 20, 5 Re t �
which corrects the prediction of uv in a channel flow. � Some models only depends on one of them:

❀ dependence on Re y = √ ky/ν ⇒ active only in the near-wall region.

❀ dependence on Re t ⇒ active in all the low-Reynolds number regions, including the near-wall region.

� There is a large number of models of this type.

� Damping functions are also necessary in the ε equation (in front of of C ε1 and

C ε2 ).
� The models using damping functions are simply (and ambiguously) called low-Reynolds number models.

� This approach avoids the use of wall functions.

� It requires a discretization point closest to the wall located at y + � 1.

� The damping functions are nonlinear (exponential) ⇒ sometimes, numerical difficulties.

� The damping functions are very empirical, chosen in order to match a particular type of flow (boundary layer).

⇒ lack of universality (but less than wall functions).

� This approach is also applicable to Reynolds-stress models.

f µ = 1 -exp � - Re y 50.5
� � A test is necessary to determine in which layer the point is: for instance the value f µ = 0.95 can be chosen as the limit of the application of the 1-equation model ⇒ the limit is in the log layer.

� These models require a discretization point closest to the wall located at y + � 1. � It appears that, in general, they are numerically more stable than low-Reynolds number models (with damping functions).

� Remark: the 1-equation model used in the near-wall region is not complete.

❀ It is however applicable because it is only used close to the wall where the value of � can be prescribed.

❀ It is thus assumed that in complex configurations the variation of the length scale � with the distance to the wall is the same as in a simple boundary layer.

� They are more empirical than low-Reynolds number models (prescription of �).

� This approach is also applicable to Reynolds-stress models. � The idea is then to write a set of equation that tends to the k-ω model in the near-wall region and to the k-ε model far from the wall.

� The transition between the two models is made by quite complex functions which depend on the distance to the wall ⇒ the principle is simple, but the details are complex.

� The model has been applied successfully to many configurations. This is a popular model in the aeronautics industry. � However, this is still a linear eddy-viscosity model ⇒ it suffers from the limitations of these models, but it inherits the good properties of the k-ω model in the near-wall region.

� Using the change of variables ω = �/(C µ k), the k-ε model can be written under the form of a k-ω model:

u i u j = -2 ν t S ij + 2 3 kδ ij ∂k ∂t + U k ∂k ∂x k = P -β * kω + ∂ ∂x k � (ν + σ k ν t ) ∂k ∂x k � ∂ω ∂t + U k ∂ω ∂x k = α P ν t -βω 2 + ∂ ∂x k � (ν + σ ω ν t ) ∂ω ∂x k � + 2(1 -F 1 )σ ω2 1 ω ∂k ∂x i ∂ω ∂x i
� These are not exactly the same equations as for the k-ω model: the coefficients in red are different; there is an additional cross-term 1

ω ∂k ∂x i ∂ω ∂x i � The model is driven by F 1 = tanh � � min � max � √ k β * ωy , 500ν y 2 ω � , 4ρσ ω2 k CD kω y 2 �� 4 � ❀ When F 1 → 1 (near-
C = C 1 F 1 + C 2 (1 -F 1 )
❀ k-ω coefficients:

α 1 =
ν t = a 1 k max(a 1 ω, SF 2 ) ; S = � 2S ij S ij
� To avoid the stagnation point anomaly, a bound is introduced in the production

P = min(ν t S 2 , 10β * kω)

Elliptic relaxation

� From the observation that most of the problems of Reynolds-stress models in the near-wall region come from the hypotheses of Chou's theory (locality, quasi-homogeneity), Durbin, in 1991, proposed a different theory.

� We have seen that the redistribution term can be written in the integral form:

φ ij (x) = � IR 3 ∇ 2 p(x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� ( 33
)

331

� The hypothesis of Durbin is that the two-point correlations decrease exponentially with the separation (r) between the two points

∇ 2 p(x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � = ∇ 2 p(x � ) � ∂u i ∂x j (x � ) + ∂u j ∂x i (x � ) � exp � - r L � r f (x, x � ) Two-point correlation function
� Using this hypothesis, it can be shown that [START_REF] Hanjalić | Natural convection in partitioned two-dimensional enclosures at higher Rayleigh numbers[END_REF] is the solution of the differential equation

φ ij -L 2 ∇ 2 φ ij = - L 2 ρ ∇ 2 p(x) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � (34)
� Far from the walls, in a quasi-homogeneous situation, L 2 ∇ 2 φ ij goes to zero, and the equation reduces to

φ ij = - L 2 ρ ∇ 2 p(x) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � (35)
� The analysis shows that pressure diffusion must be included in Eq. ( 34) ⇒ the discussion is actually about

φ ij + D p ij .
� Since far from the walls, Chou's hypotheses are valid, it can be considered that the standard models (Rotta+IP, SSG, etc.) can be applied.

� The right hand side of ( 35) can then be modeled by a standard model φ h ij .

333

� It can be deduced that the differential equation can be written as

φ ij -L 2 ∇ 2 φ ij = φ h ij �
This equation is known as the elliptic relaxation equation.

� It is seen that

❀ φ ij is not given by an algebraic relation, but rather by a differential equation (one per component).

❀ This approach can be used with any quasi-homogeneous model φ h ij (Rotta+IP, SSG, etc.).

❀ The locality and quasi-homogeneity hypotheses are not used, except for the right hand side φ h ij , which is active essentially far from the wall. ⇒ By reproducing accurately the redistribution term in the near-wall region, the model reproduces the blocking effect (nonlocal effect) ⇒ the two-component limit of turbulence is obtained.

� Drawbacks of the method

α -L 2 ∇ 2 α = 1 φ * ij = (1 -α 3 )φ w ij + α 3 φ h ij
where φ w ij is an asymptotically correct near-wall model.

❀
ν t = C µ k 2 ε by ν t = C µ v 2 k ε ( 36
)
❀ The v 2 scale (the Reynolds-stress component normal to the wall) is actually the correct scale. For instance, writing the Daly and Harlow model in a channel flow

shows that the only non-zero component of turbulent diffusion is � This was a daring hypothesis, but it appeared that the model gives good results in a many situations. � κ c is called the cutoff wavenumber.

∂u i u j v ∂y = - ∂ ∂y    C s v 2 k ε � �� � ν t ∂u i u j ∂y   
D t k = P -ε + ∇ • ((ν + ν t )∇k) D t ε = C � ε 1 P -C ε 2 ε T + ∇ • � (ν + ν t σ ε )∇ε � D t v 2 = k f - v 2 k ε + ∇ • � (ν + ν t )∇v 2 � f -L 2 ∇ 2 f = (C 1 -1) (2/3 -v 2 /k) T + C 2 P k ν t = C µ v 2 T ; P = 2ν t S ij S ij L = C L max � k 3/2 ε ; C η � ν 3 ε � 1/4 � ; T = max � k ε ; 6 � ν ε � 1/2 �
� The separation explicitly defines the computed scales and the scales called "subgrid scales" (although it corresponds better to the formalism, the term "subfilter scales" is seldom used).

� In LES, modeling concerns only a part of the turbulent motion, contrary to RANS computations.

Advantages and shortcomings

Compared to RANS computations, LES offers real advantages but also real shortcomings.

Advantages:

� The part to be modeled is smaller and more universal ⇒ simpler modeling and less crucial role of the model.

� Accounting for the unsteadiness.

353

Shortcomings:

� CPU cost much larger that RANS ❀ The reduction of the numbers of degrees of freedom is less significant, and a reduction of the number of dimensions of the problem is not possible (the large-scale motion remains 3D).

❀ The situation is even more difficult for wall bounded flows (close to the wall, the energetic turbulent structures are at very small scale) � Influence of the numerical scheme on the results

Filtering

� The objective is to compute only the so-called "large scale" component of each variable of the flow.

� This component is obtained by applying a low-pass spatial filter denoted by . , which eliminates from the instantaneous field any contribution of a scale smaller than the characteristic length scale of the filter.

Homogeneous case

� For homogeneous turbulence (statistical properties are independent of x), the filter can be also independent of x.

� In this case, the filtering operation corresponds in physical space to a convolution product of the form

f (x) = G � f = � G(x -x � )f (x � )dx � (37)
where the expression of the convolution kernel G depends on the type of filter.

The dependence on xx � indicates that the filter is assumed homogeneous.

� This convolution product can be expressed in spectral space in the form of a simple product

� f (κ) = � G(κ) � f (κ) (38)
where � G(κ) is the transfer function associated to the kernel G(x).

Fundamental properties of the filter

In the case of an homogeneous filter, three properties are to be satisfied

� Conservation of the constants � G(x -x � )dx � = 1 ⇒ c = c (39)
� Linearity

f + g = f + g (40)
� Commutativity with spatial differentiation

∂f ∂x i = ∂f ∂x i ( 41
)

357

� These three properties are also satisfied by the Reynolds-average operator.

� But, contrary to the Reynolds average, a filter is not necessarily idempotent.

Idempotence is the following property:

f = f ( 42
)
which consequence is

f � = 0 (43)
� As a practical matter, only the cutoff filter (see below) is idempotent.

Examples

Numerous possibilities exist as concerns the form of the filter kernel. The most common filters are (expressed here in 1D):

� the spectral cutoff filter � G(κ) = � 1 if |κ| ≤ κ c 0 otherwise (44)
whose kernel in physical space is

G(x) = sin(π x Δ ) π x Δ where Δ = π κ c ( 45
) Δ Spatial κ c Spectral 359 � the tophat filter G(x) =        1 Δ if |x| ≤ Δ 2 0 otherwise (46)
with the corresponding transfer functions

� G(κ) = sin( π 2 κ κ c ) π 2 κ κ c (47) Δ Spatial κ c Spectral � the Gaussian filter G(x) = � 6 πΔ 2 exp � -6 x 2 Δ 2 � (48)
with the corresponding transfer function

� G(κ) = exp � - π 2 24 κ 2 κ 2 c � (49)
Δ Spatial κ c Spectral 361

Inhomogeneous case

As a practical matter, the use of homogeneous filters are limited, for two reasons:

� Wall boundaries are not compatible with an homogeneous filter. Indeed, the width of the filter (necessarily non-local in physical space) must be modified � Inside the flow, the scales of the flow can be very different from one region to the other (for instance: reduction of the turbulent structures in the vicinity of the wall)

⇒ An optimal strategy is to use a filter width which adapts locally to the size of the energetic structures ⇒ The interest is similar to using an inhomogeneous mesh

An inhomogeneous filtering process can be written as

f (x) = � G(x, x � )f (x � )dx � (50)
� Contrary to the homogeneous case, it is not possible to express the filtering process in spectral space by a simple product of the type [START_REF] Izakson | On the formula for the velocity distribution near walls[END_REF].

� Inhomogeneous filters do not commute with spatial derivatives.

� The commutation error scales with the variations of the filter width Δ.

� It can only be kept in mind that moderating the variations of Δ limits of the error.

⇒ As a practical matter, in the very large majority of LES, the commutation error is ignored. � These variables satisfy the Navier-Stokes equations

∂u i ∂x i = 0 ( 51
)
∂u i ∂t + ∂ ∂x j (u i u j ) = - 1 ρ ∂p ∂x i + ∂ ∂x j � ν � ∂u i ∂x j + ∂u j ∂x i �� (52)
where the density ρ and the kinematic viscosity ν are assumed constant.

� To obtain the equations for the filtered variables u i and p: the filter is applied to the equations:

∂u i ∂x i = 0 ( 53
)
∂u i ∂t + ∂ ∂x j (u i u j ) = - 1 ρ ∂p ∂x i + ∂ ∂x j � ν � ∂u i ∂x j + ∂u j ∂x i �� (54)
⇒ the commutativity property ( 41) is used ⇒ these equations are exact only in the case of an homogeneous filter (commutation error is neglected). 365

The subgrid stress tensor

In order to obtain a transport equation for the filtered velocity, the equation is written as

∂u i ∂t + u j ∂u i ∂x j = - 1 ρ ∂p ∂x i + ∂ ∂x j � ν ∂u i ∂x j � - ∂τ ij ∂x j ( 55
)
where the subgrid stress tensor τ ij is introduced:

τ ij = u i u j -u i u j (56)
� Similarly to the RANS approach, the nonlinear term introduces new unknowns since u i u j cannot be expressed from u i and p.

� The τ ij tensor represents the global influence of the non resolved scales (called subgrid scales) on the dynamics on the large scales.

� It is necessary to close the system ⇒ to model the subgrid-scale tensor.

� A constitutive relation must be proposed, linking the unknown tensor τ ij to the resolved variables, u i and p.

� Since u i = u i + u � i , the subgrid scale tensor can decomposed into

τ ij = u i u j -u i u j = u i u j -u i u j � �� � L ij + u i u � j + u j u � i � �� � C ij + u � i u � j ���� R ij (57)
❀ L ij = Leonard's tensor (depends on explicitly resolved variables only): large scale-large scale interactions.

❀ C ij = cross-stress tensor: large scale-small scale interactions.

❀ R ij = subgrid Reynolds stress tensor: small scale-small scale interactions.

� Case of an idempotent filter (f = f , spectral cutoff filter only): L ij = 0 and

C ij = 0 ⇒ the subgrid stress tensor τ ij is exactly R ij .
� Non-idempotent filters: R ij remains the dominant term.

� Distinguishing among these three tensors highlights the different mechanisms.

� The three terms are sometimes considered separately. However, most of the time, a global modelling of τ ij is proposed. ❀ Functional models, which aim at correctly representing the divergence of the subgrid-scale tensor

∂τ ij ∂x j (58)
i.e., at reproducing the physical effect of the small scales on the large scales, without paying two much attention to the correct reproduction of the subgrid-scale tensor itself.

❀ Structural modelling, which aim at representing the subgrid-scale tensor τ ij itself.

� The most widely used models are those of the first type.

� "Mixed" models combine the two approaches.

Energy transfer between resolved and subgrid scales

� The equation of the kinetic energy associated to the resolved scales

K = 1 2 u i u i (59)
writes

∂K ∂t + u j ∂K ∂x j = - 1 ρ ∂(u i p) ∂x i + ν ∂ 2 K ∂x j x j -ν ∂u i ∂x j ∂u i ∂x j - ∂(τ ij u i ) ∂x j -(-τ ij ∂u i ∂x j ) � �� � ε ( 60
)
� The subgrid dissipation term ε play the most important role in the energy transfer from the resolved velocity field u i toward the subgrid velocity field u � i .

369

� It can be noted that the symmetry of the tensor τ ij also leads to

ε = -τ ij S ij (61)
where

S ij = 1 2 � ∂u i ∂x j + ∂u j ∂x i � (62)
is the large scale strain tensor.

� It is noted that, even though in a developed turbulence there is essentially a cascade of energy from the resolved scales to the subgrid scales (ε > 0), the transfer of energy can be locally negative (backscatter, ε < 0).

Schematic role of a subgrid-scale model

� A first try could consist in ignoring the role of the subgrid scales by using

τ ij = 0.
� Example: isotropic homogeneous turbulence ⇒ kinetic energy accumulates in the vicinity of the cutoff κ ≈ κ c (absence of dissipation).

ln k kc ln E(k) t 0 1 4 5 t t 3 t t t 2
� The final equilibrium state of the flow is described by a velocity field of which each component corresponds to a 3D white noise, i.e., E(κ) ∝ κ 2 .

371

� A LES of this type cannot lead to an acceptable result if it is free from numerical errors.

❀ Such "no model" computations are nowadays very frequent, under the generic name ILES (Implicit LES).

❀ A particular attention must (should) be paid to the dissipative properties of the numerical scheme.

� Physical interpretation: absence of small scales in the computation ⇒ strong underestimation of the dissipation, for which the small scales are responsible.

⇒ the most important for a subgrid-scale model is to represent the missing dissipation.

� Since the small scales are not resolved, the resolved dissipation can be considered negligible.

� So we want the average energy transfer from large to small scales be equal to the physical dissipation:

ε = ε
where ε is the statistical average of ε

ε =< ε > , ⇒
The subgrid-scale model must compensate for the absence of molecular dissipation by adjusting the modeled dissipation level.

The Smagorinsky model

� An additional viscosity is introduced in the system of equation, ν t (x, t), called the subgrid eddy-viscosity.

� The subgrid-scale tensor τ ij is evaluated using a Boussinesq-type relation, similar to the constitutive relation for Newtonian fluids

τ ij = -2ν t S ij + 1 3 τ kk δ ij (63)
❀ The subgrid viscosity plays the same role as the molecular viscosity.

❀ The subgrid energy 1 2 τ kk plays the same role as a pressure.

375

� Using this relation, the equations to be solved become

∂u i ∂t + u j ∂u j ∂x j = - 1 ρ ∂P ∂x i + ∂ ∂x j � (ν + ν t ) � ∂u i ∂x j + ∂u j ∂x i �� (64)
where the modified pressure P contains the trace of the subgrid-scale tensor

P = p + 1 3 ρ τ kk (65)
� Similarly to eddy-viscosity RANS modeling (Spalart-Allmaras model for instance), the subgrid scale energy does not require modeling, since it does not affect the velocity. It must be kept in mind that the pressure given by the computation does not correspond to the filtered pressure p.

� This model replaces the mixing due to the convection at subgrid scale by diffusion (as molecular diffusion is introduced in the Navier-Stokes equation to mimic the effect of molecular agitation).

As in RANS, it can then be shown that

ν t ∝ �u = � 2 τ ( 66
)
where � is the size of the largest subgrid eddies, proportional to the filter width (� = C s Δ), u is their characteristic velocity and τ their characteristic time-scale.

� It is assumed that the times scale of the largest subgrid eddies is the same as the time-scale of the smallest resolved eddies, τ = 1

S

, where we have defined

S = � 2S ij S ij .
� The proportionnality constant denoted by C s is called the Smagorinsky constant.

Finally, we have

ν t = C 2 s Δ 2 S ( 67
)

377

� Assuming a Kolmogorov spectrum

E(κ) = C K ε 2/3 κ -5/3 avec C K = 1.5
we get:

< S 2 >= 2 � κ c 0 κ 2 C K ε 2/3 κ -5/3 dκ = 3 2 C K ε 2/3 � π Δ � 4/3
The dissipation is thus related to the resolved velocity field by:

ε = C 2 S Δ 2 < S 2 > 3/2 (68)
� This relation yields [START_REF] Smagorinsky | General circulation experiments with the primitive equations. I. The basic experiment[END_REF] and applying the statistical averaging operator yields

C S = 1 π � 2 3C K � 3/4 � 0.18 � Substituting the Smagorinsky model τ ij = -2ν t S ij + 1 3 τ kk δ ij in the definition of ε ε = -τ ij S ij
ε =< ε >=< 2ν t S ij S ij >= C 2 s Δ 2 < S 3 > (70)
� Comparing with Eq. ( 68), the averaged modeled dissipation ε is indeed equal to the physical dissipation ε if we assume

< S 2 > 3 ≈< S 3 > 2 379
� As a practical matter, the value C s = 0.18 is adequate to compute an isotropic homogeneous turbulence (we have assumed a Kolmogorov spectrum), but must be empirically corrected for other flows.

� For instance: in a channel flow, C s = 0.065 is generally recommended ⇒ lack of universality of the Smagorinsky model.

The dynamic Smagorinsky model

� The dynamic procedure is a method aiming at evaluating the Smagorinsky constant from the knowledge of the large-scale motion (resolved) ⇒ without calibration.

� Let us introduce a second filter, called test filter, of the same type as the first one, but whose width is larger, denoted by � . :

� Δ > Δ (71)
A double-filtered variable is thus denoted by � f . The subgrid-scale tensor that appears in the double-filtered velocity equation is

T ij = � u i u j -� u i � u j (72)
� The principle of the dynamic procedure consists in assuming that the two tensors T ij and τ ij can be represented by the same model, using the same constant.

381

� For the Smagorinsky model:

τ ij - 1 3 τ kk δ ij = -2(C s Δ) 2 � 2S kl S kl S ij = -2Cβ ij T ij - 1 3
T kk δ ij = -2(C s � Δ) 2 � 2 � S kl � S kl � S ij = -2Cα ij � The tensor L ij is introduced L ij = T ij -� τ ij = � u i u j -� u i � u j - � � u i u j -� u i u j � = � u i u j -� u i � u j ( 73
)
which is only composed of explicit terms. Relation ( 73) is called Germano's identity.

� The identity (73) yields

L ij - 1 3 L kk δ ij = -2Cα ij + 2 � Cβ ij (74)
In this equation, the coefficient C is now a function of space and time.

� Problems in the estimation of C.

❀ First, ( 74) is of integral nature because of the application of the test filter upon Cβ ij . The circumvent this difficulty, it is common to approximate �

Cβ ij by C � β ij
⇒ C satisfies the following equation

E ij = L ij - 1 3
L kk δ ij + 2Cα ij -2C � β ij = 0 (75)
❀ A second difficulty lies in the necessity for C to satisfy the six independent relations [START_REF] Speziale | Turbulence modeling for time-dependent RANS and VLES: a review[END_REF]. Therefore, ( 75) is contracted into

E ij E ij = 0.
� The best way to obtain the solution is to seek the minimum as

∂E ij E ij ∂C = 0 (76)
� By defining

m ij = α ij -βij (77)
C is given by

C = - 1 2 m ij � L ij -1 3 L kk δ ij � m kl m kl ( 78
)

383

� Most the time, the field C(x, t) is too irregular, leading to numerical instabilities ⇒ it is regularized by averaging locally in space and/or time or by averaging in homogeneous directions.

� The dynamics procedure can give negative values of C. This is interpreted as backscatter (transfer of energy from the small scale to the large scale).

� The choice of Δ has a strong influence on C. A reasonable choice, recommended by many authors, is Δ = 2Δ.

� The dynamic procedure can be considered an automatic calibration of the model.

Scale-similarity models

� It is assumed that the properties of the subgrid scales are the same as the properties of the smallest resolved scales.

� The subgrid scale tensor associated to the smallest resolved scales can be computed by applying a second filter � . to the resolved velocity.

� It can be shown, using some simplification hypothesis, that

τ ij � � u i u j -� u i � u i (79)
� This approach is now to be classified in structural modelling, contrary the the dynamic Smagorinsky model, which corresponds to functional modelling.

385

� In the case of a double-filtering with the same filter � . = . (non-idempotent filter), we have exactly

u i -u i = u i + u � i -u i = u i + u � i -u i = u � i ( 80
)
� The assumption of scale similarity consists in approximating

R ij = u � i u � j � u � i u � j = (u i -u i )(u j -u j ) (81)
C ij = u i u � j + u j u � i � u i u � j + u j u � i = u i (u j -u j ) + u j (u i -u i ) (82)
� Such that the subgrid scale tensor

τ ij = R ij + C ij + L ij = u i u j -u i u j ( 83
)
can be evaluated from the resolved velocity field [START_REF] Bardina | Improved turbulence models based on large-eddy simulation of homogeneous, incompressible, turbulent flows[END_REF].

� From experiments or DNS data, it is possible to show that this approach provides a much higher level of local and instantaneous correlation with the exact τ ij than eddy-viscosity models (satisfactory anisotropy, local interaction in physical and spectral space).

� On the contrary, this type of models does not ensure the correct subgrid-scale dissipation.

τ ij - 1 3 τ kk δ ij = 1 2 � -2ν t S ij + L ij - 1 3 L kk � (84)
with

ν t = (C s Δ) 2 � 2S kl S kl (85)
and

L ij = u i u j -u i u i (86)

Deconvolution methods

� Other structural models consist in trying to reconstitutes the exact (unfiltered) velocity by applying the inverse of the filter (deconvolution). ❀ The resolved field is not the exact filtered field, such that the deconvolved field is not the exact field.

u i = G � u i ⇒ u i = G -1 � u i ⇒ τ ij = u i u j -u i u j = (G -1 � u i )(G -1 � u j ) -(G -1 � u i ) (G -1 � u j )
❀ The inverse of the filter is not known and would be very expensive to compute explicitly: an approximate deconvolution method is used, generally based on Taylor series expansions of the filter.

❀ The numerical discretization implies the presence of a projective filter: the scales smaller than the grid cells are not only damped by the filter, they are zero ⇒ the deconvolution cannot reconstitute these missing scales. ⇒ The subgrid scale tensor is not correctly reconstructed. In order to ensure a correct level of dissipation, the subgrid scale tensor must be corrected (regularization term).

Numerical viscosity and subgrid eddy-viscosity

� It is well known that some differencing schemes, in particular upwind schemes, introduce a diffusive error, equivalent to a numerical viscosity.

� In that case, it is common that the numerical viscosity is of the same order as, or even larger than, the subgrid eddy-viscosity.

⇒ The issue of numerical dissipation is crucial in LES.

� As mentioned above, many authors perform LES computations with such schemes, without any model, the subgrid dissipation being thus ensured by the numerical error (ILES = implicit LES).

� Some authors write numerical scheme in such a way that their effect is similar to a subgrid-scale model 391

Choice of the filter width

� The filter width Δ which enters the models must be adequately chosen.

� The choice must satisfy two constraints:

❀ The filter width must be sufficiently small to explicitly compute the energetic scales of the flow.

❀ It must be sufficiently large to reduce the cost of the LES compared to a DNS.

� As a practical matter, the choices of Δ and of the computational mesh are made in parallel:

❀ Indeed, the mesh must be sufficiently fine to solve the eddies of size Δ.

❀ In contrast, a mesh finer than Δ is useless and a loss of computational time.

� In the simple case of an homogeneous Cartesian mesh Δx = Δy = Δz, the following constraint should be satisfied

Δ ≥ Δx (87) � It is common to impose Δ = Δx (88)
which is optimal in terms of CPU cost.

� Relating the model to the grid step has a major consequence: numerical convergence cannot be obtained by reducing the grid step (Δx → 0) since the equations change with the grid step.

⇒ The notion of grid convergence no longer exists in LES! However, the sensitivity of the mean solution (statistics) to the mesh can be studied.

393

� For an homogeneous, anisotropic, Cartesian mesh Δx � = Δy � = Δz, it is usual to evaluate the filter width by

Δ = (ΔxΔyΔz) 1/3 (89)
but there are other possibilities, as, for instance,

Δ = � Δx 2 + Δy 2 + Δz 2 3 ( 90
)
or

Δ = max (Δx, Δy, Δz) (91)
� For complex, unstructured meshes, the filter width Δ is usually taken as the cubic root of the volume of the cell.

Wall turbulence

� In the near-wall region, there is no scale separation (the inertial range of the spectrum dissapears): it is not possible to apply a filter that separates energetic scales and dissipative scale ⇒ all the scales must be resolved � A sufficiently fine mesh must be used to actually perform a DNS in the near-wall region, but a LES elsewhere (Quasi-DNS).

� The standard subgrid-scale models, derived in homogeneous situations, are not valid in the near-wall region.

⇒ The model must be modified in the near-wall region. For instance, Van Driest's damping of the Smagorinsky constant:

C � s = C s � 1 -exp � - y + A �� ( 92
)

395

� In order to avoid the resolution of the near-wall region, wall functions can be applied.

❀ Writing valid wall functions for LES is much more difficult than for RANS, since they have to be valid for an unsteady resolved field.

❀ Example: Werner-Wengle wall functions � RANS modeling is developed to represent Reynolds-averaged quantities: this average is often equivalent to a temporal average, since many flows are statistically steady (stationary).

u + =    y + if y + < 11.
� However, many flows are non-stationary: for instance, in a pulsed flow, or a flow subjected to moving boundaries (e.g., in a car engine).

� Open questions:

❀ Can a RANS model be used in unsteady state?

❀ Or do we necessarily have to switch to the LES?

❀ Can intermediate approaches between RANS and LES be developed? ❀ In that case, the use of standard RANS models is justified.

❀ But, very often, it is not the case. From Carpy (2006).

⇒ Reynolds-stress modeling is the most appropriate level to account for this effect.

405

� Multi-scale modeling ❀ One of the most significant consequences of the turbulence equilibrium hypothesis is the evaluation of the characteristic length and time scales of the large structures:

� = k 3/2 ε ; τ = k ε
❀ This is based on the assumption that the rate of energy transfer from the large scales to the small scales is equal to ε.

⇒ It can be proposed to distinguish:

✄ the rate of dissipation ε that appears in the k equation;

✄ and the rate of transfer at large scale that appears in the evaluation of the length and time scales;

⇒ to resolve 2 transport equations (see, for instance, Schiestel, 1998).

❀ This issues are far from being solved and are still a topic of research.

URANS in statistically stationary flows

� There are numerous statistically steady cases in which, when a RANS computation is performed, it is not possible to obtain a steady solution.

In that case, it is often possible to obtain an unsteady solution if

❀ terms ∂/∂t are included in the equations;

❀ an unsteady numerical method is used.

� Example: wake of a triangular cylinder. ❀ However, experiments show that there is no clear periodicity:

✄ contrary to the laminar case, in the turbulent case, the vortex shedding is not exactly periodic ✄ an ensemble average (i.e., performing N experiments) indeed gives a time-independent solution (the "periodic" events are smoothed out because they never occur exactly after the same period).

409

❀ This problem is circumvented by considering that the vortex shedding is truly periodic and the periodic and turbulent part are separated by using a phase averaging.

< f (x, t) >= lim

N →∞ 1 N + 1 N � n=0 f (x, t + nT ) =< f (x, φ) >
where T is the period, and the phase is defined as φ = 2π t/T .

� The phase-average decomposition assumes that there is a periodic component, which leads to the following remarks:

❀ Before the computation, the period is not known, except in the case of a forced periodicity ⇒ the decomposition is implicit: the period is given afterward, by the solution!

❀ Is the solution always periodic? The answer is no! the solution can be:

✄ Steady (independent of time)! ✄ Pseudo-periodic (several periods exist).

413

� Case of a pseudo-periodic solution:

Motivation

� LES is a powerful tool, but is still to expensive for many industrial applications.

� In particular, in the near-wall region, LES requires a mesh of DNS type ⇒ quasi-DNS=very expensive (see the table at the beginning of the present course).

� Aim of the hybrid methods:

❀ defining a methodology that performs a LES in some regions and a RANS (or URANS) computation other regions.

❀ In particular, in the near-wall region, a RANS computation is required.

423

� Two approaches are possible: ❀ A first try consists in imposing the mean velocity from the RANS region to which a white noise is superimposed. This is absolutely not sufficient.

Zonal methods

LES

❀ Something similar to turbulence must be imposed (spatio-temporal coherence properties) ⇒ synthetic turbulence.

❀ In order to correctly couple with the RANS computation, the synthetic turbulence must satisfy, as much as possible, the mean properties given by the

VLES

VLES=Very Large Eddy Simulation [START_REF] Speziale | Turbulence modeling for time-dependent RANS and VLES: a review[END_REF] � The standard RANS equations are resolved with a Reynolds-stress model, which provided the Reynolds-stress tensor R ij , and, in the filtered velocity equations, the subgrid-stress tensor is modeled by

τ ij = � 1 -exp � -β Δ η �� n R ij
where η is the Kolmogorov scale and Δ a scale characteristic of the cell size.

� Thus, when the mesh is coarse (compared to Kolmogorov scale), the equations go to the RANS model equations.

� On the contrary, when the mesh is fine, the equations go to the Navier-Stokes equations (DNS).

� This model has not been applied directly.

LNS

LNS=Limited-Scale Simulation [START_REF] Batten | LNS-An approach towards embedded LES[END_REF] � This approach is directly inspired by the VLES of Speziale.

� It is based on an eddy-viscosity model.

� The eddy-viscosity given by the RANS model is used in the filtered velocity equations in the form:

ν t = α ν RANS t where α is α = min    C 2 s Δ 2 � 2 � S ij � S ij ν RANS t ; 1    (α = 1 ⇒ RANS ; α < 1 ⇒ LES).
τ ij = f s τ RANS ij + (1 -f s )τ LES ij
which reduces, in case both models are eddy-viscosity models, to

ν t = f s ν RANS t + (1 -f s )ν LES t �
Thus, this approach is not a new model, but a simple blend of two existing models.

� Note that the blending function f s is kept secret by ANSYS !

DES

DES=Detached-Eddy Simulation [START_REF] Spalart | Strategies for turbulence modelling and simulations[END_REF] � This approach is originally based on the Spalart-Allmaras RANS model.

∂� ν ∂t + U j ∂� ν ∂x j = c b1 � S� ν -c w1 f w � � ν d � 2 + 1 σ ∂ ∂x k � (ν + � ν) ∂� ν ∂x k � + c b2 σ ∂� ν ∂x k ∂� ν ∂x k
in which the length scale d is replaced by � This problem has motivated improved versions of DES (DDES: delayed detached-eddy simulation).

d DES = min

SAS

SAS=Scale-Adaptive Simulation [START_REF] Menter | SAS turbulence modelling of technical flows[END_REF] � This approach is based on the same principle as DES, but the dissipation is not modified as a function of the cell size, but rather as a function of the ratio:

L κS/U ��
where κS/U �� is the Von Karman scale: S = � 2S ij S ij and U �� = �∇ 2 U �.. � The idea is that the computation must adapt to the solution itself: typically, in the shear-dominated regions, where instabilities generate coherent structures, the model must turn to LES mode.

⇒ In contrast with DES, the model does not explicitely depend on the grid step: the model adapts itself to the solution observed during the computation.

� The underlying RANS model is a k-√ kL model.

PANS

PANS=Partially Averaged Navier-Stokes (Girimaji, 2003) � Here, the idea is different: the main parameter is the ratio

f k = k SGS k
Then, transport equations are derived for k SGS and ε.

� The equation which are obtained are formally identical to the standard k and ε equation, except for the coefficient C ε2 in the ε equation which is modified as

C * ε2 = C ε1 + f k (C ε2 -C ε1 )
� Thus, in order to perform a RANS computation in some region, the f k coefficient must tend to 1, and in order to perform a LES, f k must be < 1.

� Following the decrease of f k , ε increases and thus k decreases. ❀ All theses approaches are relatively empirical.

❀ They are essentially based on k-ε type models.

❀ Is there a mean to define things more rigorously and to use a Reynolds-stress model?

❀ A possibility: the PITM approach. From Fadai-Ghotbi (2007).

∂E ∂t = λ ij A ij + T -2νκ 2 E ⇒              du i dt = - 1 ρ ∂p ∂x i + ν ∂u i ∂x j ∂x j - ∂τ sgs ij ∂x j dτ sgs ij dt = P sgs ij + φ sgs ij -ε ij + D sgs ij dε dt = C ε1 ε k sgs P sgs - � C ε1 + C ε2 -C ε1 1 + βL 2/3 κ 2/3 c � ε 2 k sgs + D ε ⇒ LES
f (x, t) = � G(x -x � )f (x, t � )dx � (93)
is replaced by a temporal filter

f (x, t) = � G(t -t � )f (x, t � )dt � (94)
� The same methodology as for PITM is followed 443 � Eddy-viscosity models:

dk sfs dt = P + D T + D ν - k sfs T dε dt = standard � Second moment closures: dτ ijsfs dt = P ij + φ * ij + D T ij + D ν ij - k sfs T ε ij ε dε dt = standard � T = r 1 + � C ε2 C ε1 -1 � � 1 -r C ε1 C ε2 � T int where T int = k ε r = 1 β � U s √ k � 2/3 � ω c k ε � -2/3 ω c = min � π dt ; U s π Δ �
In the case of a flow along a plane plate, in the limit of high Reynolds numbers: it can be shown that the different quantities follow an "universal" behavior (theoretical result due to von [START_REF] Von Kármán | Mechanische Ähnlichkeit und Turbulenz[END_REF][START_REF] Izakson | On the formula for the velocity distribution near walls[END_REF][START_REF] Millikan | A critical discussion of turbulent flows in channels and circular tubes[END_REF] The equations below are given for the specific case of a fully developed channel flow, since the fact that the derivatives in the x-direction are zero makes the equations simpler. However, the same analysis is valid for a boundary layer. � This equation shows that the solutions must be sought as

Equilibrium

U = f (y, ν, u τ , h) uv = g(y, ν, u τ , h)
� We do not know any general analytical solution, but some relations can be found, valid in some regions, at the limit of infinite Reynolds numbers, using matched asymptotic expansions

Equations written with external variables

� The equations can be simplified by choosing the correct length and time units and by making Re τ go to ∞.

� A first possible choice is to choose the so-called external variables h and u τ to build the length and time units:

L = h ; T = h/u τ � One gets: U u τ = f � y h , ν hu τ , 1, 1 � ⇒ U u τ = F � η, 1 Re τ � where η = y h and uv u 2 τ = g � y h , ν hu τ , 1, 1 � ⇒ uv u 2 τ = G � η, 1 Re τ � 457
� Very simple formulations are obtained at infinite Reynolds number:

U u τ = F ∞ (η)
where

F ∞ (η) = F (η, 0) and uv u 2 τ = G ∞ (η) where G ∞ (η) = G(η, 0)
� It is seen that U and uv only depend on u τ and y/h, and not on the viscosity ν.

� Making (97) non-dimensional by u τ and h, gives:

- uv u 2 τ + 1 Re τ dU/u τ dη = 1 -η (98)
� When the friction Reynolds number Re τ goes to infinity, the viscous term is indeed negligible.

- uv profiles made non-dimensional using h and u τ . DNS data by [START_REF] Lee | Direct numerical simulation of turbulent channel flow up to Re τ ≈ 5200[END_REF].

uv u 2 τ = 1 -η ⇒ G ∞ (η) = η -1 (
� The region where (99) is not valid gets smaller and smaller when Re Total stress made non-dimensional using ν/u τ and u τ . DNS data by [START_REF] Lee | Direct numerical simulation of turbulent channel flow up to Re τ ≈ 5200[END_REF].

463

� Equations ( 98) and (100) are equivalent as long as the Reynolds number is not infinite.

� When Re τ → ∞, the first one is valid for η kept constant (such that y + → ∞)

and the second one for y + kept constant (such that η → 0). � In the limit Re τ → ∞, (98) becomes (99).

❀ The region of validity of (99), in which the viscous stress is negligible, is called central layer.

� In the limit Re τ → ∞, (100) becomes (101).

❀ This equation is only valid in the region where η → 0, which is called the wall layer.

❀ In this region, the flow is independent of the half-width h of the channel.

❀ In this region, the total shear stress, which is the sum of the turbulent and viscous stresses, is then a constant.

Matching

Wall layer (or internal region)

� (101) then shows that the solutions for U/u τ and uv/u 2 τ must be sought as functions of y + , (but not of h):

U + = f (y + ) ; uv + = g(y + )
� These relations are called the laws of the wall.

� Moreover, in the viscous sublayer, the first term of ( 101) is negligible, since the viscous effects are dominant, which gives

U + = y + 465

Central layer (or external region)

� (99) shows that the flow is independent of viscosity. It is convenient here to write a relation involving the velocity derivative:

dU + dη = h(η); uv + = m(η)
or, introducing H, the anti-derivative of h:

dU + dη = dH dη (102)
� Integrating (102) from the center of the channel, the velocity defect law is obtained:

U -U c u τ = H(η)
where U c is the velocity at the center of the channel. � Experiments show that κ � 0.41.

Logarithmic law

467

� The logarithmic law is deduced:

y + df dy + = 1 κ ⇒ dU + dy + = 1 κy + ⇒ U + = 1 κ ln y + + C (104)
� The region of validity of this law, which corresponds to y + → ∞ and η → 0 is called the logarithmic layer.

� It is thus the region of overlap of the wall layer and the central layer.

∂ρ * ∂t + ∂ρ * u * i ∂x i = 0 (106) � Momentum conservation ∂ρ * u * i ∂t + ∂ρ * u * i u * j ∂x j = - ∂p * ∂x i + ∂ ∂x j � µ � ∂u * i ∂x j + ∂u * j ∂x i �� + ρ * g i (107) � Energy conservation ρ * C v ∂T * ∂t + ρ * C v ∂u * i T * ∂x i = ∂ ∂x i � λ ∂T * ∂x i � + 2µs * ij s * ij ( 108
)
The terms in red (ρ * , µ, C v , λ and µ ) are functions of the temperature. ✄ For instance, Sutherland's law is often used to describe the evolution of the viscosity with temperature:

µ(T ) µ 0 = � T T 0 � 3/2 T 0 + S T + S ( 109
)
✄ The heat capacity C p and the Prandtl number Pr are often considered constant, such that

λ(T ) = µ(T )C p Pr = µ 0 C p Pr � T T 0 � 3/2 T 0 + S T + S ( 110
)
� Low-Mach-number approximation:

❀ For low Mach numbers, it can be assumed that density does not depend on pressure, but only on temperature

ρ * = f (P * , T * ) (111)
❀ For a perfect gas, we have

p * = ρ * rT * (112) such that dρ * = ∂ρ * ∂P * � � � T * dP * + ∂ρ * ∂T * � � � P * dT * = - ρ * T * dT * (113)
and

ρ * = ρ * 0 T * 0 T * (114)
❀ The flow is thus considered incompressible, but density varies as a function of the inverse of the temperature: the fluid is dilatable.

❀ The equations of motion can be derived using asymptotic expansions at the limit of small Mach numbers. 477

Boussinesq approximation

� Assuming that the temperature differences are not too large, the Boussinesq approximation can be applied:

❀ Density variations can be neglected (ρ * = ρ 0 )

❀ Except in the buoyancy term ρ * g i .

� Additional simplifications: in this case, it is standard to also assume that ❀ µ, k and C v are independent of the temperature;

❀ Density variations in the buoyancy term are linear in the temperature:

ρ * = ρ * 0 -βρ * 0 (T * -T * 0 ) with β = - 1 ρ * 0 ∂ρ * ∂T (115) such that ρ * g i = ρ * 0 g i -ρ * 0 g i β(T * -T 0 ) (116)
� The following equations are obtained

∂u * i ∂x i = 0 (117) ∂u * i ∂t + ∂u * i u * j ∂x j = - 1 ρ 0 ∂p * * ∂x i + ν ∂ 2 u * i ∂x j ∂x j -βg i (T * -T 0 ) (118) ∂T * ∂t + ∂u * i T * ∂x i = α ∂ 2 T * ∂x i ∂x i ( 119
)
where α = λ ρC p is the molecular diffusivity.

� Remarks:

❀ As usual, the hydrostatic pressure is contained in p * * :

p * * = p * -ρ 0 g j x i δ ij ( 120
)
❀ For low Mach number flows, the heat source due to viscous dissipation

2µs * ij s * ij is neglected.
* i ∂t + ∂u * i u * j ∂x j = - 1 ρ 0 ∂p * ∂x i + ν ∂ 2 u * i ∂x j ∂x j -βg i (T * -T 0 ) (121) ∂u * i u * j ∂x j � U 2 ref L ref ; βg i (T * -T 0 ) � βgΔT ⇒ Ri = βg ΔT L ref U 2 ref
� It can also be seen as a time-scale ratio ❀ The time scale of convective phenomena:

∂u * i ∂t � ∂u * i u * j ∂x j ⇒ τ conv = L ref U ref ❀ Time scale of buoyancy phenomena ⇒ Ri = τ conv τ buo (122)

481

� When the Richardson number is small compared to unity (Ri << 1), it is considered that the flow is in the forced convection regime.

� In this case, heat transfer has no influence on dynamics, the model for the Reynolds stress does not require modifications.

Reynolds decomposition

� Introducing the Reynolds decomposition, with T * = T + θ, leads to

∂U i ∂x i = 0 (123) ∂U i ∂t + U j ∂U i ∂x j = - 1 ρ ∂P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂u i u j ∂x j -βg i (T -T 0 ) (124) ∂T ∂t + U i ∂T ∂x i = α ∂ 2 T ∂x i ∂x i - ∂u i θ ∂x i (125)
� Compared to the isothermal case, three additional unknowns are involved, the turbulent heat fluxes u i θ.

483

� The turbulent heat fluxes play the same role for the mean temperature as the Reynolds stresses for the mean velocity: they increase the transport of internal energy (rather than momentum for the Reynolds stresses). As for the u i u j , first moment closures or second moment closures can be applied.

� Heat transfer modeling has received much less attention than Reynolds stress modeling. The main reason is that it is difficult to measure them. They have been investigated mainly since the 80's, using DNS. We will not go as much into details as for the Reynolds stresses. For a complete review, see for instance the article of Hanjalić [START_REF] Hanjalić | One-point closure models for buoyancy-driven turbulent flows[END_REF].

Forced convection

� When the Richardson number is small compared to unity (Ri << 1), it is considered that the flow is in the forced convection regime.

� In this case, heat transfer has no influence on dynamics, the model for the Reynolds stress does not require modifications.

� Similar to the case of the Reynolds stresses, the turbulent heat transfer can be modeled using algebraic relations (first moment closure).

Simple Gradient Diffusion Hypothesis (SGDH)

� The turbulent heat fluxes are simply written as

u i θ = -α t ∂T ∂x i (126)
� The turbulent diffusivity (or eddy diffusivity) is linked to the turbulent viscosity (or eddy viscosity)

α t = ν t /σ t . ( 127
)
� This this analogous to the Boussinesq relation for the Reynolds stresses (except that this is now a vector, rather than a tensor). It corresponds to the Fourier law for molecular diffusion.

� σ t is the turbulent Prandtl number. It can be taken as constant (in general, σ t = 0,9 or 1) or evaluated using an experimental correlation, as, for instance, the Kays and Crawford correlation

σ t = 1 0,5882 + 0,228(ν t /ν) -0,0441(ν t /ν) 2 [1 -exp(-5,165/(ν t /ν))] (128)
� Such a relation assumes that there is a direct link between turbulent diffusivity and turbulent viscosity, i.e., a similarity between the thermal and dynamic length and time scales.

� Assuming that the turbulent Prandtl number is of the order of 1 means that turbulence transports heat as rapidly as momentum. This assumption is called the Reynolds analogy. � Here, k θ is, similarly to the dynamic case, half the temperature variance

k θ = θ 2 /2 � ε θ is its dissipation rate.
� τ θ can be directly used to evaluate the turbulent diffusivity as

α t = kτ θ (129)
� However many authors suggest the use of a mixed time scale τ m , function of both τ θ and τ .

� For instance, Nagano, Tagawa, Tsuji (1991) suggest the relation

τ m = k ε � (2R) 2 + 3,4(2R) 1/2 � ν t ν � -3/4 � (130)
� It is seen that to evaluate the thermal-to-mechanical time-scale ratio R, both k θ and ε θ are required.

� Transport equations must be solved for these two variables, which leads to a 4

equation model: k-ε-k θ -ε θ , or, equivalently, k-ε-θ 2 -ε θ 2 .

489

� The exact transport equation for θ 2 reads

∂θ 2 ∂t + U k ∂θ 2 ∂x k = -2θu k ∂T ∂x k � �� � P θ -2α ∂θ ∂x k ∂θ ∂x k � �� � ε θ + ∂ ∂x k � α ∂θ 2 ∂x k � � �� � D α θ + ∂ ∂x k � -u k θ 2 � � �� � D t θ (131)
� As usual, a production term P θ , a dissipation term ε θ , a molecular diffusion term D α θ and a turbulent transport term D t θ . are recognized.

� For instance, D t θ can be modeled using a turbulent diffusivity

D t θ = ∂ ∂x k � α t ∂θ 2 ∂x k � (132)
� For ε θ , a transport equation is solved analogous to that used for ε, but the analysis of its exact equation suggest to include both dynamical and thermal production and dissipation terms

Dε θ Dt = C P 1 ε θ θ 2 P θ +C P 2 ε θ k P -C D1 ε 2 θ θ 2 -C D2 ε θ k ε+ ∂ ∂x j � � α + α t σ θ � ∂ε θ ∂x j � (133)

Generalized Gradient Diffusion Hypothesis (GGDH)

� The SGDH model is an isotropic eddy-diffusivity model.

� Analogous to the Daly-Harlow model for the triple correlations u i u j u k , an anisotropic formulation can be used,

u i θ = -C θ k ε u i u j ∂T ∂x j ( 134
)
� In this model, the u i u j 's are required, such that it is typically associated with a Reynolds stress model.

� It is also possible, as for the SGDH model, to use, rather than the dynamic time scale, k/ε, the thermal time scale τ θ or a mixed time scale, by solving additional transport equations for θ 2 and ε θ . � These models are detailed below, for the case of the mixed and natural convection regimes.

Mixed and natural convection

� When the Richardson number is not small compared to unity, the influence of buoyancy on the dynamics cannot be neglected � When the Richardson number is large compared to unity, the flow is considered in the natural convection regime (the flow is produced by the thermal field) � In the case intermediate between forced and natural convection, the regime is called mixed convection 493

Influence on the dynamics

� In order to understand the influence of the thermal field on the dynamics, it is necessary to write the exact transport equations for the Reynolds stresses � This case is the most simple case of a flow completely driven by the temperature field.

Du i u j Dt = P ij + D ν ij + D T ij + φ ij + D p ij -ε ij + G ij ( 135
G ij = -β � u i θg j + u j θg i � (
� The mean velocity and mean temperature equations reduce to � The vertical direction is a privileged direction.

0 = ν ∂ 2 U ∂z 2 - ∂uw ∂z + βg(T -T 0 ) (138) 0 = ν P r ∂ 2 T ∂z 2 - ∂wθ ∂z ( 139
� But it will be seen that the sign of the turbulent heat fluxes u i θ strongly depends on the thermal stratification, i.e., on the sign of ∂T /∂x, the vertical mean temperature gradient (zero here)

� Turbulence generates thermal agitation, which interacts with turbulence, by modifying the distribution of turbulent energy among the components. Such a complex effect typically requires a second moment closure for a correct reproduction.

Influence on redistribution

� However, buoyancy also has an indirect influence on some terms, which must be taken into account for their modeling ❀ Redistribution φ ij is modified: indeed, Chou's analysis is modified by the fact that the transport equation for the instantaneous velocity involves the buoyancy

φ ij = φ 1 ij + φ 2 ij + φ 3 ij ( 140
)
through a term denoted by φ 3 ij , which reads

φ 3 ij = � IR 3 βg k ∂θ ∂x k (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� (141)

501

✄ The term is absolutely similar to the rapid term.

✄ Rapid term φ 2 ij (reminder) φ 2 ij is a function of the mean velocity gradient and the dynamical turbulence

φ 2 ij (x) = � IR 3 2 ∂U l ∂x k (x � ) ∂u k ∂x l (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� (142)
Production:

P ij = -u i u k ∂U j ∂x k ⇒ IP model: φ 2 ij = -C 2 � P ij - 2 3 P δ ij � ✄ Buoyant term φ 3 ij : similar φ 3
ij is a function of βg and the coupling dynamical/thermal turbulence:

φ 3 ij (x) = � IR 3 βg k ∂θ ∂x k (x � ) � ∂u i ∂x j (x) + ∂u j ∂x i (x) � dV (x � ) 4π�x � -x� (143)
Production:

G ij = -βg i u j θ -βg j u i θ ⇒ IP model: φ 3 ij = -C 3 � G ij - 2 3 Gδ ij � 5.1.3.
Influence of buoyancy on dissipation � The dissipation term ε ij is also modified.

❀ This can be seen by writing the exact equation for dissipation that also involves buoyancy terms.

❀ Most of the time, the model for ε ij is not modified, but only the transport equation for ε, in which a buoyancy production term appears

G ε = -2βg i ν ∂u i ∂x k ∂θ ∂x k (144)
� A buoyant term must be added to the equation solved for ε:

dε dt = C ε1 ε k P k + G ε -C ε2 ε 2 k + D ε (145)

503

� There is no consensus in the literature: does buoyancy significantly influence the energy cascade ⇒ the dissipation? How?

❀ Negligible influence:

G ε = 0 ❀ Same influence as P k : G ε = C ε3 ε k G k with C ε3 = C ε1 or C ε3 < C ε1
❀ Unstable stratification: influence; Stable stratification: negligible influence:

G ε = C ε3 ε k max(G k ; 0)
❀ Influence function of the flux Richardson number:

C ε1 ε k (P k + G k )(1 + C ε3 Ri f ) with Ri f = -G k P k +G k

Transport equation for the turbulent heat fluxes

� Similar to the Reynolds stresses, the exact transport equation for the turbulent heat fluxes can be written � There are however important differences:

∂u i θ ∂t + U k ∂u i θ ∂x k = ∂ ∂x k � -u i u k θ � � �� � D t θi + ∂ ∂x k � 1 ρ p θδ ik � � �� � D p θi + ∂ ∂x k � νθ ∂u i ∂x k + αu i ∂θ ∂x k � � �� � D ν θi + 1 ρ p ∂θ ∂x i � �� � φ θi -βg i θ 2 � �� � G θi - � u i u k ∂T ∂x k + u k θ ∂U i ∂x k � � �� � P θi -(ν + α) ∂u i ∂x k ∂θ ∂x k � �� � ε θi ( 146
❀ Production P θi does not arise from the turbulence/mean velocity gradient interaction but rather from both the turbulence/mean temperature gradient and turbulent heat fluxes/mean velocity gradient interactions.

❀ The buoyancy production term G θi involves the temperature variance θ 2 .

❀ Molecular diffusion must be modeled, except for the case ν = α, i.e., when ❀ wθ, which is the only component that has an influence on the mean temperature (see above), is production by the turbulence/temperature gradient interaction.

Pr = 1, since in this case ∂ ∂x k � ν ∂u i θ ∂x k � . 507 5
❀ This remark justifies that gradient models (SGDH and GGDH) can work, despite the fact that they are oversimplified: the flux of energy from the hot wall to the cold wall follows the direction opposed to the temperature gradient and its intensity is driven by the turbulence intensity.

❀ There is no production of vθ: this component is only fed by the scrambling term.

513

� A positive temperature gradient produces a negative wθ négatif � that in turn implies a negative production of w 2 .

� A stable stratification thus damps the vertical fluctuations.

� Can lead to 2-component turbulence (in the atmosphere or the ocean, for instance).

⇒ This analysis of a simple case shows that thermal stratification influences the turbulence anisotropy. � It can be seen (Figure ) that the temperature gradient is zero at the center.

� So, in the central region, ∂T /∂z = 0, and turbulence is produced by the buoyancy term βgθ 2 only.

� If βgθ 2 is neglected (SGDH or GGDH) ⇒ laminar (linear) temperature profile. � Only the simplest model will be presented below, but more sophisticated models can be used.

In particular, the method based on the general formulation of tensorial relations can be applied.

� An analysis similar to that of Chou for the redistribution term shows that the term φ θi can also be separated in three terms:

❀ The slow term: similar to Rotta's model, Monin (1965) proposed

φ 1 θi = -C θ1 1 τ u i θ ( 147
)
where τ is the dynamic time-scale k/ε.

❀ But a mixed time-scale can also be used (cf. previous discussion about the SGDH model).

517

� Rapid term: it can be analytically shown (using a quasi-homogeneous assumption) that this term can be deduced from the model for φ 3 ij . If the model

φ 3 ij = -C 3 � G ij - 2
D ν θi = ∂ ∂x k � ν + α 2 ∂θu i ∂x k � (151)
� It has at least the merit of giving the exact term when ν = α. � Warning: for other fluids (liquid metals, for instance), this assumption is not valid.

Eddy-viscosity models

Reminder:

Boussinesq relation:

u i u j = -
u i u j = 2 3 kδ ij + k ε 1 -C 2 C 1 � P ij - 2 3 P k δ ij � � �� � Boussinesq part + k ε 1 -C 3 C 1 � G ij - 2 3 G k δ ij � � �� � Buoyancy extension (152)
⇒ Buoyancy extended model: (Davidson, 1990 [17]) (153)

u i u j = 2 3 kδ ij -2ν t S ij � �� � u i u j Bouss + C * θ τ � G ij - 2 3 G k δ ij � � �� � u i u j Buo
u i θ D θi = 1 2 � 1 k D k + 1 θ 2 D θ 2 �
� the algebraic equation is obtained

� P θi - u i θ 2k P k - u i θ 2θ 2 P θ 2 � + φ * θi - � ε θi - u i θ 2k ε - u i θ 2θ 2 ε θ 2 � = 0 (157)
� Introducing for instance the models described in section 5.4. yields

P U θi + P T θi + G i -u i θ � P k 2k + P θ 2 2θ 2 � -C θ1 ε k u i θ -C θ2 P U θi -C θ3 G i + u i θ � ε 2k + ε θ 2 2θ 2 � = 0 (158)

529

� It is deduced that

u i θ = P T θi + (1 -C θ2 )P U θi + (1 -C θ3 )G i C θ1 ε k + 1 2k (P + G -ε) + 1 2θ 2 � P θ 2 -ε θ 2 � (159)
� If equilibrium is assumed for the dynamic turbulence (P + Gε = 0) and the thermal turbulence (P θ 2ε θ 2 = 0), the following model is obtained

u i θ = -C θ k ε � u i u k ∂T ∂x k + (1 -C θ2 )u k θ ∂U i ∂x k + (1 -C θ3 )βg i θ 2 � ( 160
)
� This is the AFM (algbebraic flux model) (Hanjalić, Kenjeres, Durst, 1996 [33])

⇒ Second moment closure is the natural level to account for these phenomena: Reynolds stress model (RSM) + Differential flux model (DSM).

� Simplified models are also used: � We will not address this issue in this course, however, it is to be noted that, to obtain good predictions in the near-wall region, and in particular the heat transfer between the fluid and the wall, it is necessary (and sometimes sufficient) to reproduce correctly the near-wall dynamics, in particular the two-component limit of turbulence (the wall-normal fluctuations ensure the transport of energy).

functions must also be written (for the mean temperature, and, possibly, for θ 2 , ε θ , u i θ).

� Moreover, the universality of the dynamic wall functions is called into question when the thermal field has an influence. Despite many research efforts, none of the available wall functions are satisfactory in heat transfer.

� The extension to heat transfer of near-wall models is an active research topic.

One can go as far as for the dynamics as concerns the modeling sophistication (TCL model, elliptic relaxation, etc.).

535

� Contrary to a widespread belief, solving an engineering problem using a CFD code does not consists in safely clicking some buttons in the interface.

� Many parameters chosen by the user (model, mesh, boundary conditions, discretrization schemes, convergence criteria, etc.) have an influence on the results ⇒ results can be very user-dependent if the code is not correctly used.

� An educated use of CFD is necessary.

� The turbulence model is a major component, that have been extensively discussed in this course.

� This appendix provides an overview of the other important parameters and the way to chose them. ❀ Turbulent boundary layers develop on both sides.

❀ For the sake of simplicity, it is assumed that the boundary layers are turbulent from the leading edge. � For the code to give the correct results, here is a series of points that deserve attention:

❀ The equations.

❀ The near-wall mesh.

❀ The computational domain.

❀ The boundary conditions.

❀ The numerical schemes. � In the present 2D geometry, the mean variables, such as mean pressure, mean velocity, turbulent energy, dissipation, etc., are independent of z and t: φ(x, y, z, t) = φ(x, y)

� In such a 2D, stationary case:

❀ The mesh can be 2D. In industrial codes, it is usually a 3D mesh with only 1 cell in span with periodic or free-slip boundary conditions (symmetries).

❀ A steady-state algorithm is possible.

� For a RANS computation, we have to chose between "low-Reynolds number"

and "high-Reynolds number" strategies.

� As explained above, this choice has nothing to do with the value of the Reynolds number Re = U ∞ c/ν.

❀ We are talking about the resolution or not of the near-wall region (y + < 30).

❀ This is the region where the turbulent Reynolds number ν t /ν is low, i.e., not � 1. ❀ Let us just evaluate the number of cells necessary to cover the boundary layer (the computational domain is much larger).

❀ The boundary layer thickness δ is about 10 cm at the end of the wing (trailing edge).

❀ The number of layers of cells of constant size necessary to cover the boundary

layer is n = δ Δ 1 � 50000
� Therefore, a geometric expansion of the grid step is used

Δ i = rΔ i-1 (161)
As already seen in the RANS chapter, section 8.6.4:

� One wants to match a progressive near-wall mesh (Δ 1 → Δ n ) and the external mesh (Δ ext ).

Δ1 Δn

Near-wall mesh ❀ This is to cover the boundary layer. For the rest of the domain, the two strategies use the same meshes.

⇒ Using the low-Reynolds number strategy is possible nowadays.

❀ The problem is not the computing power, but rather the difficulty to mesh a complex domain with y + c � 1 ensured everywhere.

❀ That is why most of the CFD codes propose so-called adaptive wall functions (AWF) or enhanced wall treatment (see the RANS chapter, section 8.6.5).

551

� Remark:

❀ the grid step in the directions parallel to the wall is not dependent on the Reynolds number.

❀ It is to be adapted to the geometry (local curvature of the wall, geometrical details, etc.) and the expected flow structure (sudden variations of the velocity field). there is no mass flux through these boundaries.

Control by suction

� This is artificial: the flow domain is, in the real world, unbounded.

� Mass conservation implies

U inlet b S = U b s
where S is the inlet section and s the section at some location. U b is the bulk velocity (velocity averaged over the section).

⇒ there is an acceleration due to the restriction of the section (blockage effect).

� This artificial acceleration can affect the results.

553

� In order to minimize this effect, the acceleration must be as small as possible.

� The acceleration is linked to the blockage ratio

B = S -s S = 1 - s S
In the present 2D case, the blockage ratio reduces to

B = e H
where H is the height of the domain.

Location of the outlet boundary

� Outlet boundary condition also impose an artificial constraint of the flow computed.

� Oftentimes, the outlet conditions are just homogeneous Neuman conditions ∂φ ∂x n = 0 � Here, there is no help from the theory. The only available method is to compute the flow for several locations of the outlet and compare the results. � Sometimes, they can be rigorously obtained from an preliminary computation (e.g., fully developed pipe flow), but usually they must be infered from available information.

� Experimental data sometimes provide k, but not ε, that cannot be measured.

� Many codes propose entering the turbulence intensity I and the integral length scale � instead of k and ε, using the relations

k = 3 2 U 2 ∞ I 2 and ε = C 3/4 µ k 3/2 � ( 162
)
� This is very useful since � is sometimes available or can be inferred from the geometry.

� For instance, in fully developed flows in ducts (i.e., sufficiently far from the inlet of the duct), the turbulence intensity and integral length scale can be roughly estimated as

I � 7% ; � = 5%D H
where D H is the hydraulic diameter

D H = 4 A ζ
i.e., the ratio of the section A of the duct to the perimeter ζ of the duct (so-called wetted perimeter).

� Many codes automatically prescribe inlet conditions for turbulent variables if the user does not prescribe them.

❀ This is very dangerous! The code basically uses the values for a developed duct flow: they are only valid for duct flows, and only if the correct length scale is chosen.

❀ An educated user (i.e., a user who has read the present chapter) should never let the code choose the inlet boundary.

561

� In the present case (airfoil), the incoming flow is laminar � Boundary condition k = 0 and ε = 0 cannot be used, since they lead to the solution k = 0 and ε = 0 everywhere in the domain.

� The standard method consists in imposing a low turbulence intensity (e.g., I � 0.1%) and a low turbulent Reynolds number (e.g., R t � 10).

� The length scale is thus obtained using

R t = ν t ν = C µ ν k 2 ε
and relations (162), yielding

� = νR t � 3 2 C 1/4 µ U ∞ I 562

Initialization

� In the case of a steady-state computation, the initial conditions are not supposed to have an influence on the results.

� However:

❀ Some initial conditions can lead to a divergence of the computation ❀ The number of iterations necessary to reach convergence can significantly depend on the initial conditions ❀ Sometimes, there are several solutions to the system of equations. The solution obtained at the end of the computation can depend on the initial conditions [START_REF] Rumsey | Arbitrary steady-state solutions with the k-ε model[END_REF].

� The initialization with U i = 0, k = 0, ε = 0 is not always the best choice. Some constant, but non-zero values are sometimes preferable (same evaluations as for the boundary conditions).

� Some codes have automatic initialization procedures.

563

Numerics

Grid refinement studies

� There is no objective way to directly build a sufficiently fine mesh. The error due to spatial discretization cannot be estimated.

� The only safe procedure is to perform a grid refinement study: several computations using different meshes.

� The solutions using different meshes must be compared: as many quantities as possible must be plotted (global quantities, velocity profiles, turbulent energy profiles, streamlines, etc.).

� The quantities of interest must be plotted as a function of the number of cells n: for instance, if the main quantity of interest is C x (drag coefficient), the curve C x = f (n) must be plotted (with at least 3 points).

Convection scheme

� The most critical term in the equations is the nonlinear term (convection).

� Codes usually propose to choose the discretization scheme for this term among:

❀ The basic upwind differencing (UD) scheme: first-order accurate.

❀ A variety of second-order or third-order accurate schemes.

� The UD scheme is less accurate, but much more stable (thanks to numerical diffusion), in particular at the beginning of the computation (artificial initial conditions):

❀ It is usual to start the computation with UD until convergence and restart with a more accurate scheme.

❀ Sometimes, it is impossible to obtain a steady-state solution with second/third order accurate schemes (oscillations).

569

� Contrary to a common belief, it is fully acceptable to work with the UD scheme:

❀ The low order of accuracy just means that, for a given mesh, the solution is less accurate.

❀ A finer mesh must be used with UD to reduce the numerical error to the same level as with higher-order schemes.

❀ It means that the CPU time will be larger to obtain the same level of numerical error.

570 571

OutlineI. 9 1. Challenges of numerical simulation 1 . 1 .

 911 Introduction to CFD � Challenges � Fields of application � Different components of a computation � Computational cost � The different modeling approaches � Global picture of CFD codes 5 II. Industrial standard: RANS � Closure problem, similarities with continuum mechanics � Physical and mathematical principles � Eddy-viscosity models � Reynolds-stress models � Problems due to walls Introduction to computational fluid dynamics (CFD) for turbulent flows What is CFD? Configuration to be studied: Here, the flow above a wall with pins (cooling of turbine blades) solve them do iel=1,ncell pij = xprod(iii,jjj) phiij1 = -cvara ep(iel)* (cssgs1*xaniso(iii,jjj)+cssgs2*(aikakj-d1s3*deltij*aii)) phiij2 = -cssgr1*trprod*xaniso(iii,jjj) + trrij*xstrai(iii,jjj)*(cssgr2-cssgr3*sqrt(aii)) + cssgr4*trrij*(aiksjk-d2s3*deltij*aklskl) + cssgr5*trrij* aikrjk epsij = -d2s3*cvara ep(iel)*deltij w1(iel) = cromo(iel)*cell f vol(iel)*(pij+phiij1+phiij2+epsij) w2(iel) = cell f vol(iel)/trrij*crom(iel)*( cssgs1*cvara ep(iel) + cssgr1*max(trprod,0.d0) ) end do

15 1. 4 .❀⇒❀

 154 The cost issue� Strongly dependent on the particular case and the stage in the design cycle:❀ Some experiments are very expensive (cost of a large wind tunnel > 20 kE/day + scale models) or impossible (e.g., nuclear accident) ❀ The cost of a month of a CFD engineer is much lower. But software licenses are expensive (CFD codes, meshing software, post-processing software), computers, storage, ... In early stages of the design cycle, if several concepts are considered, CFD is attractive Scale models: very expensive. Takes a long time to manufacture. On the contrary, for aerodynamic optimization at the end of the design cycle: very rapid tests of small modifications are possible in a wind tunnel � CFD is attractive at the beginning of the design cycle ⇒ Complementarity 1

17

 17

1. 5 . 2 .�Difficult step 21 2. 2 .

 52212 Examples of industrial applications (RANS) What is to be meshed is not the elements, but the volume between the elements, contrary to structural analysis. � An external "skin" is needed � Only useful elements must be kept � Details must be removed (reduction of the size of the mesh) This step can consist in weeks of efforts = Volume discretization (mesh) � No analytical solution of the Navier-Stokes equations � Example: finite volume discretization: the fluid volume is divided into small control volumes (cells) in which the variables (3 velocity components, pressure) are considered constants N cells=N values for each variable � Approximation of the equations ⇒ 4 × N system of equations for 4 × N unknowns ⇒ can be solved by a computer � Fine mesh (large N ) ⇒ accurate solution (convergence toward the exact solution)

23 25 ����

 25 Refining the mesh in appropriate regions is crucial: � boundary layers; � around geometrical details; � in the wakes. Ahmed body Control by suction Local refinement: � keeps a reasonable number of cells; � ensures an accurate solution Meshing is then a crucial step Theory: grid convergence must be ensured (sufficiently fine mesh to obtain a good approximation of the exact solution) � In practice: not possible in industrial configurations Ex: external aerodynamics computation at PSA � Despite the large number of cells, solutions are still dependent on the mesh ⇒ a low quality mesh (not refined in crucial regions) is a significant source of errors = second difficult step 2.3. Physical modeling Real physics is very complex: what must be accounted for? � Compressibility (influence on the dynamics if the Mach number is M = U/c > 0, 3). Acoustics. Heat transfer (convection, radiation, conduction) Different phases: liquids, sprays, ice, particles, etc. Chemical reactions: combustion � Fluid-structure interaction: vibration, deformation of the geometry � Environment/variable geometry: wind (cars, buildings), transients, mobile elements (turbines, ...) � Turbulence A good knowledge of the physics is necessary to make the appropriate choices Many phenomena are difficult to represent ⇒ many physical approximations = Third difficult step, the most critical! 27 Similarly: the boundary condition issue: � the domain is artificially bounded � Hypotheses are necessary about the quantities at the boundaries � These hypotheses are also physical approximation

Flow 33 3. 3 .∝ Re 3 / 4 t 4 t 35 3. 3 . 2 . 1 )Smallest vortices � 15 µm 3 . 4 .� 1 ⇒ 4 �✄ 18 ✄ 6 ✄ 37 3. 5 . Available computer power 3 . 5 . 1 .❀❀ 2 ❀Flops 53 �� 1 � 55 ��field 63 ��

 33334435321341418637535125315563 around a square cylinder at Re = 22000 (DNS, Trias et al., 2015) Visualization of the magnitude of the pressure gradient Movie at: https://www.youtube.com/watch?v=c8zKWaxohng 31 Full domain Close-up (scale 2) Close-up (scale 4) Close-up (scale 8) DNS of homogeneous turbulence at Re λ = 732. From Yokokawa et al. (2002). Influence of the mesh � Mesh not sufficiently fine ⇒ the smallest eddies are not resolved Resolved eddies Unresolved eddies 3.3.1. Free turbulence (no walls) � Ratio of the largest vortices to the smallest vortices: l η Re t = ul ν where u and l are the characteristic scales of the largest eddies � Example: wake of a car (100 km h -1 ) Re t � 150000 ⇒ l η � 7500 ⇒ η � 100 µm � η determines the grid step for a direct numerical simulation of turbulence (DNS) � The number of grid cells is proportional to Re 9/Wall turbulence (boundary layer) � Ratio of the largest vortices to the smallest vortices ∝ Re τ Re τ = δu τ ν where δ is the boundary layer thickness and u τ the friction velocity � Number of grid cells ∝ Re 11/4 τ (less than 3, due to the expansion in the direction normal to the wall) ❀ Example: boundary layer of a car (100 km h -Evaluation of the computational cost � Number of floating point operations per iteration: Many numerical methods have a computational cost (number of operations) proportional to N ln N (where N is the number of cells) Time step linked to the grid step for stability and/or accuracy reasons: Courant-Friedrich-Levy (CFL) number: CFL = U Δt Δx � time-step Δt proportional to Re -3/For a good convergence of the statistics, the "fluid particles" must travel through the domain at least 10 times. � Example: car (100 km h -1 ) Number of grid cells � 10 Number of time steps � 10 Memory requirement � 10 20 bytes ✄ Computational cost � 10 27 floating point operations Performance of processors � Measured in Flops: FLoating-Point Operations Per Second = clock frequency × number of operations per clock cycle This unit is a measure of the available brute force � This theoretical performance is... theoretical ❀ Peak performance = number of operation the processor(s) (CPU) are able to reach ❀ Actual performance = performance really achieved during a computation � Can be much lower, in particular for parallel computer � Optimization of the computing time ⇒ find a high-performance scheduling ⇒ HPC=High-Performance Computing Designed for computing the same operation onto a vector (vectorial processor) rather than a real (scalar processor) ❀ Very efficient treatment of loops: DO I=1,10000 X(I)=Y(I)*Z(I) END DO Scalar processor: sequential treatment of the iterations Vectorial processors: N data treated simultaneously ❀ Historically: the Cray family Cray-2: first supercomputer above the GFlops in 1985 Many others, among which the Playstation Their cost is such that they are replaced today by parallel machines Evolution of the most powerfull computer from 1940 (adapted from T.B. Gatski, private communication) The memory necessary to compute the flow around a car at 100 km h -1 represents 11000 times the largest supercomputer today ⇒ The computation is far from possible � If the computation were possible, it would take 18 years! � Extrapolation: if the increase of the computing power continues at the same pace as today, DNS will be possible in car/aerospace engineering in 2080 (evaluation from Spalart, Boeing) Using available computer power in the car industry today, it would be possible to perform DNS of cars at only 1 km h -Using 100% of the largest super-computer: DNS of a car at 3 km h -1 Memory limitation. � Is such computing power desirable? Frontier (rank 1) Power: 21 MW (Electric consumption of a town of � 27000 inhabitants) The ecological footprint in terms of CO 2 emitted by one CPU hour in France has been estimated at 4.68 g (Berthoud et al., 2020) � For example, Perez Arroyo's (heroic) calculation (2021) of a fan-compressor-combustion chamber assembly of an aircraft engine in LES (with wall functions) required 31.6 million CPU hours on the Joliot-Curie national supercomputer, which corresponds to about 150 tons of CO 2 , i.e.: ❀ 1.2 million km for a car ❀ 3000 km of flight of an A320 ❀ The average consumption of a French person during fifteen years � The 18-year calculation cited above would result in 6.7 million tons of CO 2 (emissions of about 700,000 French people during one year) An operator . is applied to the equations of motion (Navier-Stokes equations) ❀ RANS: operator = statistical average ❀ LES: operator = low-pass convolution filter � Physical variables (velocity u * , pressure p * ) are decomposed into Navier-Stokes equations: ∂u * i ∂t = F(u * i , p * ) � Decomposition + Application of the operator ⇒ ∂u * i ∂t = F(u * i , p * ) -∂τ ij ∂x j Due to the nonlinearity of the equations:

5 .�

 5 Represents the influence of the non-resolved field on the resolved field (stress) ❀ RANS: τ ij = Reynolds-stress tensor ❀ LES: τ ij = subfilter-stress (ou subgrid-stress) tensor ⇒ τ ij requires modeling Global picture of the CFD codes � Codes can be classified into three main categories: Developed in research institutions by teams or even individuals

❀ 67 5. 3 .��❀❀ 69 5. 5 .

 673695 developing the internal skills and knowledge ❀ developing specific applications ❀ avoiding the "black-box" effect of commercial codes � Examples: ❀ TrioCFD (CEA) ❀ elsA (ONERA-Airbus-SAFRAN) ❀ CEDRE (ONERA) ❀ Code Saturne (EDF) ❀ AVBP (CERFACS, IFPEN) Commercial codes � Often very generalist ("able to address all the problems") � Licenses, in general expensive, must be bought (very expensive for parallel computing) They are black-boxes (no access to the source code) � They can be modified using user subroutines � Examples: The sources (the code, not only the executable files) of some codes are publicly available, freely distributed under an open-source license ❀ Code Saturne by EDF ❀ OpenFoam Incompac3d by Pprime/Imperial College � This strategy aims at: Enabling a development at a lower cost, by a community of users: Linux type strategy ❀ Enabling a community of users, exchanges, large-scale validation ❀ Sell associate services (advice, studies, etc.) : Google type strategy ❀ Demonstrate the skills and knowledge of the company (EDF case) Distinctive features of the codes � Beside the numerical methods, the type of meshes, etc., codes have their own features.

The mesh❀ 71 �

 71 Boundary conditions/size of the domain ❀ Discretization schemes � But the most influencial choice is the turbulence model: it is very important to understand the underlying hypotheses and the limitations The only rigorous way to ensure the quality of the results would be to test, for every case, the influence of all these parameters ⇒ perform a large number of computations � In practice, this is impossible � Know-how must be built-up � Oftentimes engineers follow a fixed methodology: recommended model, mesh type, domain size, numerical schemes, etc.

91

  Mixing layer downstream a thick plate Experiments by Perret and Delville (Institute Pprime, university of Poitiers/CNRS/ENSMA) RANS solution (k-ε low-Reynolds number model).

101�❀❀

  The question everybody wants to ask: what modeling level do we need? This question does unfortunately not have a simple answer, for different reasons: It depends on the effort one is ready to put in the study: Reynolds stress models are not much more expansive in terms of CPU cost, but slightly more difficult to correctly handle. It depends on the objectives: prediction of global quantities (drag, lift, global heat transfer) or of more precise details (separation location, flow structure, turbulent scales, anisotropies, etc.).

3. 3 . A quick historical review of the different types of models 3 . 3 . 1 .❀❀�❀�

 3331 Before the advent of computers: the pioneers 1877 Boussinesq invents the concept of eddy viscosity. Simple phenomenological reasoning: the "turbulent material" behaves like a Newtonian fluid, i.e., the mean stress is proportional to the mean strain ❀ This is the first turbulence model. The notion of eddy-viscosity is introduced. 105 What is the reasoning of Boussinesq? By observing river flows, he understood that: Turbulence strongly enhances mixing Reynolds (1883) experiment: flow in a pipe, transport of dye Laminar Turbulent ❀ Turbulence is at the origin of a resistance to the flow stronger than in a laminar flow Modern interpretation: Navier-Stokes equations

❀❀

  dissipation of the mean mechanical energy: 2ρ(ν + ν t )S ij S ij (energy is not directly dissipated into heat, but first transformed into turbulent energy, as will be seen later) 111 Comparing this relation with the constitutive relation for Newtonian fluids, it can be seen how this model is similar to the Newton model. ❀ The eddy-viscosity and the turbulent energy play roles equivalent to those of the molecular viscosity and the pressure, respectively. ❀ However, this model does not provide a relation to evaluate ν t and k, which strongly vary in the flow and from one flow to another. This type of models is called an algebraic model, or a zero-equation model. n-equation model = model involving n additional differential equations (compared to the averaged Navier-Stokes equations, equations (

1925�

  The turbulent viscosity must also satisfy: ν t � �u 121 Prandtl introduces the concept of mixing length to evaluates ν t .� We have seen that ν t � �u, but how to evaluate � and u? � In a boundary layer above a flat plate (idealized model for the boundary layer above an airfoil), one want to reproduce the streamwise velocity U . U Since variations along x are very slow compared to those along y, the only significant component of the Reynolds stress in the U -equation is uv. Keeping the dominant terms only, the U -equation reduces to

�❀❀

  Boussinesq's model yields uv = -ν t ∂U ∂y . How can one evaluate ν t ? � Displaced particle argument: Let us imagine a fluid particle located at y above the wall, that follows the mean flow (u * = U (y)) U We move it vertically upwards: in other words, we give it a velocity v * > 0. Its vertical fluctuating velocity (difference between the instantaneous velocity and the mean velocity

131 3 . 3 . 2 .

 332 After the advent of computers: the heroic age of CFD 1967 Cebeci and Smith ❀ Algebraic model (of mixing length type), improved later on by Baldwin and Lomax (1978) ❀ The variation of the mixing length with the distance to the wall is prescribed. ❀ These models have been and are still used intensively in aeronautics. They are not valid for separated flows.

❀

  This is the famous standard k-ε model. ⇒ It is the birth of modern RANS modeling. 1975 Launder, Reece and Rodi (1975) ❀ most well-known Reynolds-stress model (known as the LRR model, or the Rotta+IP model).

  133

137�

  Example: in an article published in 1997 in Journal of Fluid Mechanics, L. Wang proposes the following constitutive relation for the Reynolds stress

⇒⇒⇒

  These remarks lead to ρ R = F (ρ, ∇U) ❀ Now, there is no possible combination of ρ and ∇U = The constitutive relation introduced by Wang is not valid. It is necessary to introduce additional scales in this type of relation: the turbulent energy k, its dissipation rate ε, the turbulent time scale k/ε, etc. are used.

  141

147�❀✄✄

  Digression: the issue of the material frame indifference (MFI) � Noll's material frame indifference principle is one of the bases of the classical continuum mechanics. � Noll's principle: the properties of a material do not depend on the reference frame. ⇒ the form of the constitutive relation does not depend on the reference frame. No relation with objectivity: Objectivity = property of a physical quantity (temperature, force, etc.) MFI = property of a material Objectivity = purely descriptive, what observers see (or measure) MFI = how a material responds to some applied external constraint � Consequences: The constitutive relation for a material does not explicitly depend on time. The stress varies in time, but only because the strain varies in time: the relation between the stress and the strain is the same at all instant. This can be demonstrated from the material frame indifference by a change of the origin of time.

❀

  a constitutive relation of the form R = f (S, ω) is not correct (the modelled R is not objective) ❀ a constitutive relation of the form R = f (S) (linear eddy-viscosity model) is correct, but insensitive to system rotation 153 � A relation of the form R = f (S, Ω, ω) can be objective. ∇U = S + Ω Symmetric part (objective): S = 1 2

  called the absolute rotation rate tensor � Conclusion: modelers work with S and W.

155 4 . 5 . 1 .

 451 Use in modeling

159�

  The straight line G = 0 corresponds to the limiting case when one of the eigenvalues of the Reynolds stress is zero (two-component limit), going to a one-component turbulence (a second eigenvalue going to zero) at the top of the triangle.

4 .

 4 Consequences for the turbulence modelsThe transport equations for the Reynolds stresses (15) are

  167

169 4 . 7 .

 47 Consistency with the boundary layer theory� A model is expected to reproduce at least a standard boundary layer.

  et al.)

  171

173 5 . 5 . 1 .

 551 Eddy-viscosity models (first moment closures) Constitutive relation � Reynolds-averaged Navier-Stokes equations

�

  Theory of invariants (cf. Deville and Gatski, 2012): a symmetric tensorial function of two tensors (one symmetric, one skew-symmetric), is a polynomial of the form

❀

  The variable coefficients can be obtained by machine learning ❀ A neural network can give the values of the coefficients as a function of the invariants. In this case, the relation is not known explicitly ❀ Another approach consists in finding an explicit relation (function) using a genetic algorithm (Gene Expression Programming, GEP) In all cases, a large amount of data is required for the learning phase ❀ It is a very active research topic. It is not a mature technique ❀ It does not solve the numerical issues 5.1.4. Reduced models Cubic models � Craft, Launder and Suga (1996): necessary to retain the cubic terms to reproduce the curvature effects (sensitivity of the Reynolds stresses to the curvature of the streamlines).

12 ❀ 3 k

 123 figure below. ❀ The Boussinesq relation leads to the profiles shown in the right part of the figure.

Figure 4 :

 4 Figure 4: Eddy-viscosity models. Profiles at x/r0 = 1.25 (continued).

193��❀

  Reminder: a model is complete if no a priori knowledge about the flow is necessary ⇒ these scales must be obtained as solution of transport equations, not prescribed by the user. Particular case: linear constitutive relation (Boussinesq relation) R = -2ν t S + 2 3 k I The Reynolds-averaged Navier-Stokes equations become

s 201 ❀✄

 201 Source (or sink) terms s: generation (or destruction) terms for the quantity f . Examples: ✄ a chemical reaction can produce carbon oxide in the volume ⇒ source term in the transport equation for carbon oxide. ✄ a volume force (due to gravity, Coriolis, magnetic field, etc.) creates momentum (acceleration) ⇒ source terms in the transport equation for momentum. Turbulent production P is the generation term (or sometimes the destruction term) of turbulent energy k ⇒ source (or sink) term in the k-equation. ✄ Dissipation transforms turbulent kinetic energy into heat ⇒ sink term in the k-equation and source term in transport equation for heat.

209 5 .

 5 3.3. Calibration procedureThe coefficients are determined considering the following flows:

❀❀

  One of the most important constraint to impose to a model: correct reproduction of the logatithmic region of a channel or boundary layer, where U + = 1 κ ln y + + C et al.) U + = y + U + = 1 κ ln y + + C 215 Channel flow: production exactly reduces to P = -uv ∂U ∂y ❀ Boundary layer: this relation is not exact, but an excellent approximation ❀ Now, the k-ε model gives P = C µ k 2 ε � ∂U ∂y � 2

  ε1 ) where κ = 0.41 is the Karman constant. ❀ This relation leads to the standard choice σ ε = 1.3. 217 5.3.4. Why k and ε? � Choice of k: rather natural ❀ appears directly in the constitutive relation (however, it is obviously possible to evaluate k from other scales) ❀ k is mainly chosen because its exact transport equation is simple and the mechanisms are identified easily.

5. 4 .

 4 k-ω model � Wilcox (1988): back to the old Kolmogorov idea of using the characteristic frequency of the large eddies (also called specific dissipation) ω = ε C µ k as the second scale � Transport equation for ω in the same form as the ε equation � Linear constitutive relation

❀❀

  Remark: this condition is not consistent (the discrete solution does not converge towards the exact solution when the grid is refined)! ❀ It also implies that the correct near-wall behavior k ∼ y→0 y 2 cannot be reproduced 223 In external aerodynamics (e.g., a wing), the model predicts a boundary layer development sensitive to the turbulence level present outside of the boundary layer, which is usually fixed arbitrarily by the user ⇒ this model is not applicable to external aerodynamics.

227 5 . 5 . 2 .

 552 The standard k-ε model � This model is based on a linear constitutive relation, which has some advantages, and many limitations.

❀ 2 ⇒ 2 ✄✄

 22 It gives exactly the same solution in a channel flow with and without rotation. ⇒ It is completely irrelevant to turbomachinery problems, for instance. ✄ Rotation/curvature correction are often introduced. 229 In a simple shear flow: predicts equal diagonal components of the Reynolds stress u 2 = v 2 = w it is not able to reproduce the secondary flow due to anisotropy. For instance, in a square sectioned duct, the analysis of the transport equation of the mean vorticity Ω x = ∂W/∂y -∂V /∂z shows that the secondary eddies are produce by the term ∂ 2 (v 2w 2 ) ∂y∂z . The linear model gives Ω x = 0. After Speziale (1991). stagnation point): ✄ Exact production is in the form P = (w 2v 2 )D ✄ the Boussinesq relation leads to P = 4ν t D For strong strains (4ν t D is higher than w 2v 2 ): the Boussinesq relation leads to a strong overestimation of the production ⇒ This phenomenon is known as the stagnation point anomaly. 231 This anomaly is critical: -It leads to a strong overestimation of the turbulent energy in boundary layers downstream of stagnation points ⇒ It significantly delays separation. Turbulent energy contours. (a) Without the stagnation point anomaly. (b) With the stagnation point anomaly. From Durbin and Pettersson Reif (2001).

233 5 .❀

 5 5.3. Nonlinear modelsNonlinearities in the constitutive relation ⇒ present pros and cons compared to the standard k-ε model.� Advantages:❀ The constitutive relation is more subtle and can potentially reproduce complex flows. This relation can account for the effects the linear relation is not sensitive to: negative production, rotation effects, generation of secondary vortices, etc.

❀

  [START_REF] Yakhot | Renormalization group analysis of turbulence. 1. Basic theory[END_REF]: applying the theory of the renormalization group (from theoretical physics) to turbulence ⇒ one can theoretically derive the equations of the k-ε model.❀ Method based on a spectral partitioning of turbulence + simplifying assumptions concerning the interactions between the scale.

241❀❀

  The other coefficient are C µ = 0.085 ; C ε1 = 1.42 ; C ε2 = 1.68 ; σ k = 0.72 ; σ ε = 0.72 η 0 = 4.38 ; β = 0.012 The model is often proposed in commercial code and is rather popular in the industry, in particular because it corrects the stagnation point anomaly.

�

  write as the divergence of a flux ⇒ it is not a source term, rather a term that transports Reynolds stresses from one region to another region, as viscous and turbulent diffusion do: it is called pressure diffusion. � This decomposition make easier the interpretation of the role of the different terms. v e c t i o n M o l e c u l a r T u r b u l e n t Warning: this sketch must not hide the fact that: ❀ production terms are not always positive; ❀ the different components of production and dissipation are not equal; ❀ the different transport terms through the control surface can transport energy in both directions, depending on the situation. � The terms D T ij and D p ij are most of the time called "turbulent diffusion" and "pressure diffusion", respectively, because they are transport terms leading to effects analog to that of the molecular diffusion D ν ij . But they actually result from completely different mechanisms: ❀ D T ij results from the mean effect of the convection of energy by the fluctuating velocity. ❀ D p ij results from the mean effect of the fluctuating pressure forces.

6. 2 . 2 . 1 ij�

 221 Modeling of the slow term � The slow term Φ 1 = � φ , according to its equation, only depends on the state of turbulence ⇒ the mean field and density have no reason to enter the model.

�

  The fourth order tensor A ijml (x) is a function of velocity fluctuations ⇒ the Reynolds stress R and turbulent scales only (e.g., k, ε) enter the model. � Φ 2 is linear in the velocity gradient ∇U = S + W ⇒ the model must be a function of R, k, ε, S and W, linear in S and W.

  of the Rotta+IP and SSG models. DNS data from Moser et al. (1999).� So-called "wall echo" terms ❀ The figure above shows that the Rotta+IP model overestimates the redistribution in the log layer.

277 6 . 3 . 1 .

 631 Shir's model � Simplest model[START_REF] Shir | A preliminary numerical study of atmospheric turbulent flows in the idealized planetary boundary layer[END_REF]: considers the triple correlation as a simple diffusion of the Reynolds stresses by a turbulent viscosity (similar to the k or the ε equation), i.e., in using a simple gradient diffusion hypothesis (SGDH):

  this relation, k play indeed a role different of that of i and j.

279 6 . 3 . 2 .

 632 The Daly and Harlow model � Daly and Harlow (1970) introduced a generalized gradient diffusion hypothesis (GGDH):

���

  from the exact equations of the triple correlations with some drastic simplifications. Similar to the model of Daly and Harlow, but satisfies the interchangeability of the indices. The model works slightly better than the Daly and Harlow model. But it involves 27 terms for each Reynolds stress component instead of 9.

283 6 . 5 .��

 65 Conclusions about the Reynolds-stress models � It is very difficult to give as precise indications as for eddy-viscosity models about the relative performance of the Reynolds-stress models. � Indeed, the Reynolds-stress models have been much less used, in particular in industrial application: their implementation in commercial codes is recent and they are not used much. As a summary, only the two models have spread outside of the research domain: Rotta+IP model and SSG model. Both models are based on the Daly and Harlow model for the turbulent transport term and the Kolmogorov model for dissipation ⇒ they differ by their redistribution term.

285 7 . 1 .

 71 Basis hypotheses of algebraic modeling� From the transport equations of the Reynolds stresses

295 8 .

 8 The near-wall region � We have considered so far that molecular viscosity effects are negligible compared to turbulent effects since Re t = ν t ν � 1

8. 1 .u 2 ✷ v 2 � w 2 8 . 1 . 2 .

 1222812 Physical phenomena specific to the near-wall region 8.1.1. Viscous effects � strong mean velocity gradients (non-slip condition) ⇒ a turbulence production peak. � a damping of all the components of the fluctuating velocity. � a narrowing of the turbulence spectrum (vanishing of the inertial zone). * uv DNS data in channel flow at Re τ = 590 (Moser et al., 1999) 297 Non-viscous effects � Wall echo: ❀ It can be mathematically shown that the effect of the reflexion of the pressure fluctuation on the wall can be represented considering an image flow. ❀ The image flow generates pressure fluctuations which are added to the ones generated by the real flow.

p❀❀

  This effect leads to an increase of the pressure fluctuations and consequently to an increase of redistribution. In a boundary layer type of flow (flow parallel to the wall), the redistribution takes energy to the u 2 component and gives it to v 2 and w 2 .

  et al.)

�

  This laws have been obtained based on very strict hypotheses: ❀ flow above a plane plate (channel flow or boundary layer); ❀ very high Reynolds number (limit of infinite Re); ❀ no pressure gradient.

��

  The universal profile of the mean velocity in the near-wall region shows that, in order to solve the mean velocity down to the wall, the discretization point closest to the wall must be placed at y + � 1, since ∂U + ∂y + = 1. The condition is very strict: for instance, for a flow at a Reynolds number of 10 6 , with a characteristic length scale (geometry) of the order of one meter, the location y + = 1 corresponds approximately to y = 20 µm.⇒ In the past (70's-80's), such fine meshes were not affordable.

0 ❀❀ 8 . 6 . 4 .�

 0864 Instead of imposing the value of k, a homogeneous Neuman boundary condition is used: ∂k ∂y = Another velocity scale, u k , is evaluated from the turbulent energy by u 2 k = C 1/2 µ k where k is taken at the discretization point closest to the wall and at the previous iteration. ❀ u τ is evaluated from U The friction is imposed as µ ∂U ∂y = τ w = ρu τ u k 313 Advantages of the wall functions: They give good results in situations close to ideal conditions (very high Reynolds number, flow parallel to the wall, weak curvature of the wall, etc.).

❀❀319 8 . 6 . 5 .

 865 Example: near-wall meshes for Re = 100 000, with Δ ext = 10%d For a large-scale industrial case at large Reynolds number, a reduction by a factor of 2 or 3 of the total number of cells can be expected, which is already significant! Limitations of the wall functions � The relations does not hold in theory: at stagnation, separation, reattachment points, in recirculation regions, when the boundary layer is subjected to a pressure gradient, when the wall is curved, when the flow is 3D, etc.

323 �

 323 It can be seen that this function depends on both the distance to the wall and the turbulent Reynolds number.

325 8 . 8 .��

 88 Two-layer models � So-called two-layer models divide the flow in 2 regions: the near-wall region and the rest of the domain. � In the rest of the domain, the basic equations of the model are solved (for instance, the standard k-ε model). In the near-wall region, a one-equation model is used: the k equation is solved, but ε is obtained from an algebraic relation. � Example: the Norris and Reynolds (1975) model. A damping function is also used for the turbulent viscosity:

327 8 . 9 .

 89 The k-ω SST Model� To derive the SST model[START_REF] Menter | Two-equation eddy-viscosity turbulence models for engineering applications[END_REF] remarks that:❀ the standard k-ε model does not work in the near-wall region;❀ the k-ω model works well in the near-wall region, but is sensitive to the external turbulence level; ⇒ a "mixed" model would have all the advantages.

329 �

 329 wall region): the cross-term disappears and k-ω model is recovered. ❀ When F 1 → 0 (far from the wall): the cross-term is active and the equations are exactly the standard k-ε equations with the change of variable ω = ε/k. Coefficients are a blend of k-ω and k-ε coefficients

5 / 9 ;�

 59 β 1 = 3/40 ; σ k1 = 0.85 ; σ ω1 = 0.5 ❀ k-ε coefficients: α 2 = 0.44 ; β 2 = 0.0828 ; σ k2 = 1 ; σ ω2 = 0.856 To avoid an overestimation of the shear stress in adverse pressure gradient boundary layers (very important situation for aeronautics), a bound is introduced in the eddy-viscosity (hence the name SST for Shear Stress Transport)

  Effect of the elliptic relaxation on the models: Rotta+IP ; SSG. Comparison with DNS data in a channel flow at Re τ = 590 (Moser et al., 1999).

A

  priori test of uv given by: the standard eddy-viscosity; the v 2 -f eddy-viscosity. DNS data from[START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re τ = 590[END_REF] : channel flow at Re τ = 590 � The equations for v 2 and its redistribution term f 22 (simply called f ) are identical to those obtained in a channel flow (i.e., only the wall-normal component is considered in the model).� These equations, which are valid in a channel flow, are then considered valid in any configuration.

339 9 .�����

 9 Conclusion about RANS models� RANS modelling is a complex topic, which requires a lot of physical understanding and a bit of maths � There is no model that gives good results in all the situations � Nor is one model better than the others � This is the reason why there are many different models in the codes � In each application field (aerodynamics of planes, of cars, turbomachinery, heat transfer, ...), industries choose one model and often, they keep it for a very long time Research on RANS models has been active for 60 years! � There has been a big decline in research (not as concerns the use in industry) over the last 20 years: a lot of hope has been placed in unsteady approaches (LES, hybrid RANS/LES) There is a renewal based on new approaches. The data-driven approaches (AI, machine learning, data assimilation) make it possible:❀ to refine the models by using variable coefficients ❀ to consider a mixture of experts: AI make a choice or blends the models in the different regions of the flow ❀ The GEKO model of Ansys Fluent is a kind of mixture of experts but hand-made (regions chosen by the user) 341 Because of the cost and ecological impact of unsteady simulations, RANS will stay for a long time:❀ On its own as long as the results are satisfactory for a given application ❀ In combination with LES otherwise (hybrid RANS/LES) Objective: explicitly compute the energetic eddies (large scales) while modeling the dissipative eddies (small scales) on the left) and dissipation (axis of the right) spectra in log-log scale (left) and log-lin (right). Quantitites are made non-dimensional using k and ε.351 In spectral space: the objective of LES is to compute the dynamics of the structures corresponding to wavenumbers κ < κ c while modeling the effects of the structures such that κ > κ c .

�

  Notation: in LES, instantaneous variables are denoted by lowercase letters f , filtered variables f and the so-called subgrid variables f � = ff .

  355

363 4 .

 4 The filtered equations� In an incompressible flow, the unknown are the instantaneous velocities u i and pressure p (note that the notation is different from the one used in RANS modeling).

�

  Two model families are often distinguished:

ε

  Kinetic energy spectrum of a turbulence in equilibrium (P = ε = ε).

�

  Aim at combining the dissipative properties of eddy-viscosity models and the correct representation of the subgrid scale tensor of the scale-similarity models. � Example: the model of Bardina et al. (83) and the Smagorinsky model can be combined with equal weights

389 �

 389 Several issues are associated with this method.

81

 81

  flow. Werner-Wengle law compared to the mean velocity profile where y + = yu τ ν and u + = u u τ , with u τ the wall friction ρu 2 τ

399 1 .❀t❀❀❀

 1 RANS modeling in non-stationary flows� Let us consider a case in which the unsteadiness is due to a variation in time of the boundary conditions.❀ For instance, take-off of an airplane, flow inside a car engine, around an moving obstacle, in a pulsed duct (blood flow, for instance), pulsed jet (control). In this case, the Reynolds average (ensemble mean) is a function of space and time: f (x, t)Example of signals measured in transient or cyclic turbulent flows � In that case, a modeling issue can however be raised: are the models described in this course still valid? In particular, many modeling choices are based on the idea of a turbulence in equilibrium: constitutive relations in the eddy-viscosity models (no memory effect), use of ε in the length and time scales, etc. If the boundary conditions vary slowly: the turbulence can be considered quasi-steady. In that case, turbulence is close to equilibrium. This requires that the time scale of the external variations T ext is large compared to the time scale of turbulence τ .

401�

  Example: synthetic jet Synthetic jet issuing in a fluid at rest. Experiment of Yao et al. (2004). ❀ Solutions of a Reynolds-stress model (Rotta+IP) and the standard k-ε model. Rotta+IP model k-ε model Q isocontours (top) and vorticity (bottom) From Carpy (2006).

403❀⇒

  In that case, the main consequence of the unsteadiness is a misalignment of the anisotropy and strain tensors. θ=angle between the eigenvectors of b and S ❀ In the eigenaxes (e1, e2, e3) of b: Production: P = 2k b : S = 2kβ(λ1 -λ2) cos 2θ Rotta+IP model k-ε model Production and velocity vectors.

❀❀�❀

  This solution is puzzling: the Reynolds averaged is supposed, in this statistically steady flow, to be independent of time ❀ what is the meaning of what is obtained? 407 Even more puzzling: the Strouhal number is correct (St = f D/U ) and the time-averaged solution is better than a steady solution:From Durbin Many URANS computations are performed, in particular in the industry, because:❀ in many cases, it is difficult to obtain a steady solution ❀ the quality of the results is improved compared to steady RANS computation ❀ correct frequencies are obtained � Can we give a clear meaning to these solutions? A possibility, often proposed, is a change of point of view on the flow: There is a periodic component in the flow, by nature different from turbulence.

From

  [START_REF] Iaccarino | Unsteady 3D RANS simulations using the v 2f model[END_REF] ❀ The solution is absolutely not periodic ⇒ in that case, it is difficult to give a clear formal definition of the solution: the model induces an implicit filter.

�

  RANS computation: turbulent energy (kε) or Reynolds stresses (Reynolds-stress model); length and time scales of the large structures. ❀ This is a very active research topic. 427 From Druault et al. (2004) [19] � Full simulation of a mixing layer � Zonal simulations: (a) Inlet conditions extracted from the full simulation (b) White noise matching the energy level If κ c → 0, the model must go to a RANS model. � If κ c → ∞, the model must go to zero ⇒ Navier-Stokes equations (DNS). � In between: LES.

  429

�

  In a similar way, SBES simply blends the Reynolds stress given by a RANS model and the subgrid stress given by a LES model

���

  It is seen that the role of this modification is to increase the sink term in the equation of � ν when the mesh is sufficiently fine ⇒ the turbulent viscosity diminishes ⇒ tends to a subgrid-scale eddy-viscosity.� This model is very popular.433 It can be adapted to other models than the Spalart-Allmaras model, by increasing the dissipation rate in the k equation:ε DES = max � L C DES Δ , 1 � εwhere L is the integral scale of turbulence (for instance, SST-DES model) � A well identified problem is the gray zone: between RANS and LES, the subgrid model is not really valid. This problem can be at the origin of wrong solutions: for instance, if the grid is suddenly refined in a boundary layer, the computed flow can artificially separate (grid-induced separation).

437

  PANS simulation of the flow around a rudimentary landing gear (Q-isosurfaces colored by the velocity magnitude) From Krajnovic et al.

� 3

 3 zones defined in the spectrum:[0; κ c ] [κ c ; κ d ] [κ d ; ∞]� Integration in each zone of the transport equation for the energy spectrum:

  with transport equations for the subgrid stress tensor. � When κ c → 0: tends to a RANS model. � Example: adaptation of the elliptic-blending Reynolds-stress model to the PITM model ❀ Channel flow at Re τ = 395 ❀ Coarse mesh: 32x54x32 cells.

  uv. Resolved component (LES) and modeled component (SGS).

441� 3 . 9 .��

 39 Refinement of the mesh ⇒ tends a LES: Flow above a periodic hill (Q-isosurfaces colored by the velocity magnitude) HTLES / SRH turbulence model HTLES=Hybrid Temporal LES (Manceau, Friess, Gatski, 2011; Manceau 2018) SRH turbulence model=Scale-Resolving Hybrid turbulence model (name given by StarCCM+ [55]) The hybrid methods presented so far are either empirical or based on homogeneity assumptions (PITM) For inhomogeneous flows, replacing spatial filtering by temporal filtering makes it possible to derive a consistent hybrid RANS/LES formalism � The usual spatial filter of LES

�� 455 �❀

 455 In a channel flow, the equilibrium of a slice of fluid of the same width 2h as the channel requires that the friction at the wall τ w = µ dU dy Integrating between 0 and y the momentum equation in the x-direction gives 0 Using (95) and introducing the friction velocity defined by τ w = ρu 2 τ , This relation is exact: le total stress (turbulent + viscous) varies linearly from ρu 2 τ at the wall to 0 at the center of the channel.

99 )�

 99 This function very well describes the behavior of uv when Re τ → ∞,

�� 2 τ

 2 τ increases, and tends to vanish when Re τ → ∞ � But considering that it disappears does not allow us to satisfy the boundary condition uv = 0 at the wall. 459 By making Re τ → ∞ for a constant η, we have introduced a discontinuity at the wall. � Close to the wall, the viscous term is always necessary to describe the flow (no-slip condition). � It is then necessary to use another length scale such that the viscous term does not disappear at high Reynolds numbers. This relation is indeed satisfied by the flow for Re τ → ∞ in the near-wall region: + dU/uτ /dy + = 1

� 1 � 30 � 2 τ��

 1302 Channel flow at Re τ = 5200 (Lee & Moser, 2015 [46]) � Wall layer: η < 0.Central layer: y + > The dashed line is the log law � In the overlap region, (99) and (101) both give uv = u The turbulent shear stress is then a constant. � This result suggests that in this layer, the Reynolds stresses are constant, which implies that the turbulent transport is zero: the turbulence is said to be in local equilibrium. The transport equation of k then reduces to P = ε. Since P = -uv dU/dy, the equations (104) and (105) show that in this layer Back to the Navier-Stokes equations, we have � Mass conservation (continuity)

�

  Variation of the physical properties :Approximated dependence laws can be used for the physical properties of the fluid.

479 2 . The different flow regimes 2 �

 22 Relative weight of the buoyant term compared to the convective term = Richardson number ∂u

  485

487 4 . 2 .

 42 k-ε-θ 2 -ε θ models � In order to avoid using this restrictive hypothesis, the thermal-to-mechanical time-scale ratio R = τ θ /τ can be defined, where τ = k/ε and τ θ = k θ /ε θ .

491 4 . 4 .

 44 Second moment closure� A second moment closure can also be used, solving transport equations for u i θ.

) 5 . 1 . 1 .

 511 Modification of production� In comparison with the isothermal case, a buoyancy term appears

136 )� 2 G�

 1362 Transport equation for kDk Dt = P + Gε + D ν + D T + D p ii = -βg i u i θ Energy cascade:Flow between two infinite plates, separated by a distance h, generated by the temperature difference between the two walls: ∂/∂t = 0; ∂/∂x = 0; ∂/∂y = 0 V = 0; W = 0.

) 497 �❀��

 497 It can be seen that❀ As for the case of an isothermal channel flow (generated by a pressure gradient), only the Reynolds shear stress has an influence on the mean velocity (uv for a standard channel flow).❀ The buoyancy term plays in this flow the role of the driving force that was played by the pressure gradient in the standard channel flow case. Only the component wθ has an influence on the mean temperature profile. The production terms of the Reynolds stresses reduce to It can be seen that the temperature fluctuations interact with the velocity fluctuations to create (or suppress) turbulence (potential energy � turbulent energy).

) 505 �

 505 All the terms are already familiar, since they are similar to those appearing in the transport equations for the Reynolds stresses: ❀ turbulent transport D t θi ❀ pressure diffusion D p θi ❀ molecular diffusion D ν θi ❀ pressure-temperature gradient correlation, called scrambling term φ θi ❀ production P θi and G θi ❀ dissipation ε θi

. 3 .�

 3 Illustration of the role of the production terms 5.3.1. Differentially heated vertical channel It can be seen that ❀ The vertical direction is a privileged direction for the production of turbulent heat fluxes.

509 5 . 3 . 2 . 2 TFigure :

 5322: Figure: Flow between two infinite plates, separated by a distance h, generated by a pressure gradient ∂P/∂x, with a temperature difference between the two walls: ∂/∂t = 0; ∂/∂x = 0; ∂/∂y = 0; V = 0; W = 0.

5. 3 . 4 .

 34 Rayleigh-Bénard flow � Particular case of the flow without a pressure gradient (∂P/∂x = 0) with an unstable stratification .

5. 4 .

 4 Modeling of the unknown terms 5.4.1. Scrambling term � Similar to the redistribution term for the dynamics, this term is the most important one in the modeling process.

3 ) 5 . 4 . 2 .�

 3542 Gδ ij � is used, the model ofOwen (1973) is obtained for φ 2 θi φ 2 θi = C θ2 u k θ ∂U i ∂x k (148)� Buoyancy term: as for the dynamic buoyancy term, an isotropization of buoyancy production is used (Owen, 1973):φ 3 θi = C θ3 βg i θ 2 (149Turbulent transport term � The adopted model the model of Daly-Harlow, already used for the dynamics � The model represents both the triple correlations and the pressure diffusion by D t θi + D p θi = The following model is often used

��

  For the most common fluids, the Prandtl number is neither very small nor very large compared to unity � Most of the time, similar to the case of the Reynolds stress dissipation tensor, ε ij , isotropy is assumed for the dissipation vector of the turbulent heat fluxes ε θi This assumption implies ε θi = 0 (isotropic vector=0)

2 .�

 2 Buoyancy-extended Eddy-Viscosity ModelsReminder:Du i u j Dt = P ij + G ij + φ ijε ij + D ν ij + D T ij + D p ij Weak equilibrium + Equilibrium (P k + G k = ε):

� - 1 525

 1 Pr t extracted from a DNS of stably-stratified channel flow From Garcia-Villalba and del Alamo, 2011[START_REF] García-Villalba | Turbulence modification by stable stratification in channel flow[END_REF] � Modifying Pr t is a patch � At least, it should not be a constant but a function of the Richardson number � Variations must be modest � Using weak equilibrium assumptions

❀❀❀ 6 .

 6 RSM+Algebraic Flux Model (AFM) RSM+gradient models (SGDH, GGDH) Very often: Eddy-viscosity models (EVM), with very limited success 533 Final remarks: the near-wall region � The same difficulties as for the modeling of the dynamics are faced for the modeling of the heat transfer.

539�

  We consider the geometrically simple case of an airfoil in a uniform flow at U ∞ = 280 m s -1 .

U∞ = 280 m s - 1 e❀

 1 The geometry is assumed invariant by translation in the spanwise (z) direction.

541 1 .

 1 Choice of the equations and of the near-wall mesh� As explained above, it is not possible to perform a DNS of this high-Reynolds number flow. The flow is computed in RANS.

�

  Discretization point closest to the wall Discretization point closest to the wall Low-Re model High-Re model 543 The two methods have completely different requirements in terms of meshing of the near-wall region: ❀ Low-Reynolds number: fine mesh (y + c � 1) where y c is the distance to the wall of the discretization point closest to the wall (center of the control volume adjacent to the wall in the finite volume method).❀ High-Reynolds number: coarse mesh (y + c > 30).

• y c Δ 1 ���

 1 Since y + c = y c u τ ν , we need an estimate of u τ to evaluate the physical distance y c , in order to decide what is possible or not, and to build the corresponding mesh. An elaborate estimation is possible in the present particular case of a boundary layer, but not in general, complex configurations. Using the rule of thumbu τ = 5% U ∞ is sufficient in general. � Here, u τ = 0.05 U ∞ = 14 m s -1 ⇒ y c = ν u τ � 1.1 µm� For a low-Reynolds number strategy, the size Δ 1 of the cell adjacent to the wall must be as small as Δ 1 = 2y c � 2 µm 545 Using such small cells everywhere in the domain would lead to a huge number of cells.

Δext External mesh 547 ���r = 1 . 2 → n = 51 ��r = 1 . 1 → n = 54 r = 1 . 2 → n = 32 550�❀

 547125111541232 If one fixes the two thicknesses Δ 1 and Δ ext , the expansion ratio r and the number of cell layers n are imposed byΔ ext = Δ 1 r n and d = Δ 1 r n -1 r -1where d is the distance covered by the progressive mesh. For large Reynolds numbers,Δ 1 � Δ ext , such that r � 1 + Δ ext dThe number of layers necessary to cover the entire boundary layer δ (i.e., we use d = δ here) can then be easily evaluated:r = 1.01 → n = 625 r = 1.05 → n = 160 r = 1.1 → n = 89A reasonable number of layers is obtained using such meshes.549 Now, for a high-Reynolds number strategy, the size Δ 1 of the cell adjacent to the wall must be larger than Δ 1 = 30 × 2y c � 60 µm ❀ Layers of constant size necessary to cover the boundary layer: The number of layers is lower than for the low-Reynolds number strategy, but not much lower. The factor is only about 2.

552 2 . Domain 2 . 1 .

 221 Height of the domain� A free slip condition is generally applied at the top and bottom boundaries:

�

  Boundary conditions for turbulent variable (say k and ε) can have a strong influence.

Complementarity with wind tunnel experiments Experiments CFD

	Physical real-ism?	� Measurement of the real flow, but in a simplified environment (complex in-teractions with other objects, variable boundary conditions, perturbations, vi-brations, other unknown/uncontrolled phenomena	� Model: (equations) � Potentially: everything is possible (two-phase flows, heat transfer, acous-tics, ...) � But most of the time: simplified sit-uation
	Geometrical realism?	� Depends on the level of details in the scale model	� Potentially exact (CAD), but actu-ally always simplified
	Technical limi-tations?	� Limited access to physical quantities (some points/planes, not pressure and velocity at the same point, etc.) � Measurement uncertainties (a few percents)	� All the quantities are available, ev-erywhere � Error under control, depending on the available computing power

1.2. Numerical simulation: Why? 1.2.1. Because experiments are impossible Collision Milky Way/Andromeda 4 billion years from now Climate modelling 13 1.2.2. Because measurements are long/complex/costly/partial Aerodynamics of cars... ...and planes 1.3.

first moment, algebraic model (zero equation).

	1942 1945 Independently, Kolmogorov and Prandtl proposed to relate the turbulent velocity scale to the turbulent kinetic energy:
			u ∝	√	k
	❀ Introducing a coefficient C µ , we have
	Mixing layer	Plane jet	Round jet ν t = C µ √ k �	Plane wake
	relation known as the Prandtl-Kolmogorov relation.
					127
	❀ Two equations of state are thus necessary to evaluate these variables.
	� This approximation is very useful as a "rule of thumb". ❀ Prandtl proposes to solve a transport equation for k, based on the exact
	� Under more elaborate forms, it is still used in aeronautics for attached flows equation, in which the unknown terms are modeled.
	(see below Cebeci and Smith, Baldwin and Lomax).

From Viollet et al. (1998)

125

� This is also a

  not complete: it requires to prescribe the length scale �. � Historically, the first complete model is the Kolmogorov model, which, as seen

	earlier, evaluates the turbulent viscosity by
		ν t = C µ	k ω
	The model is complete because the two scales, k and ω, are given by equations
	of state, which are the transport equations:
	dk dt	= ...	dω dt	= ...
	⇒ This model does not require any a priori knowledge about the flow.
					143
	� There are only three possibilities for building a complete model:
	❀ Eddy-viscosity models
	1. Using the Boussinesq constitutive relation u i u j = -2ν t S ij + writing a transport equation for ν t .	2 3	kδ ij and
	2. Using any constitutive relation (including the Boussinesq relation) and
	evaluating ν t from two turbulence scales (k-ε, k-ω, etc.) for which transport
	equations are solved.

  analysis shows (cf. for instance[START_REF] Lumley | Computational Modeling of Turbulent Flows[END_REF] or Schiestel 1998) that the realizability of the Reynolds stress implies that the invariants II and III must remain inside the anisotropy triangle, also called Lumley's triangle.

	0.3						1 C l i m i t
					G=0
					2C turbulence
	0.2
	-II
						H=1
	0.1					axisymmetric turbulence
			H=-1
	axisymmetric turbulence		Isotropic turbulence
	0		-0.02	0	0.02	0.04	0.06
					III
	where G =	1 9	+ 3 III + II and H =	III 2 (-II/3) 3/2

✄

  The main problem is the estimation of the turbulent viscosity by the Prandtl-Kolmogorov relation ν t = C µ k 2 /ε:

	-1
	-2
	-3
	-4
	-5
	-6
	-7
	0	100	200	300	400	500	600

0

A priori test: comparison of uv and -2ν t S 12 . DNS data of

[START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re τ = 590[END_REF]

.

✄ A specific treatment in the near-wall region is necessary.

  [START_REF] Shih | A new k-ε eddy-viscosity model for high Reynolds number turbulent flows. Model development and validation[END_REF]: similar reasoning on the realizability of the Reynolds stress ❀ Variable Cµ coefficient, rather than an upper bound for the time scale (looking into details, this is more or less the same).

� Shih et al.'s correction ❀

  The model is completely linear.❀ It is less sophisticated than the SSG model, which does not mean that it gives less accurate results.

❀ This is the most widely used Reynolds-stress model. It is available in some industrial codes. 271 � Other models ❀ More complex models exist. ❀ In particular, they are built such a way to satisfy the realizability constraints (the Rotta+IP and SSG models are not realizable). ❀ It would be too long to present them here. Some examples are the models derived by ✄ Shih and Lumley (1985); ✄ Launder and Tselepidakis (1991);

✄

Craft and Launder (1996)

, TCL model, for Two-component limit (see its slow term in section 6.2.2.).

  is much more popular because it is anisotropic.� The model does not satisfy the interchangeability of the indices either.

	� The coefficient proposed by Daly and Harlow is C s = 0.22, but values between
	0.20 and 0.25 are used.

6.3.3. Hanjalić and Launder's model � Hanjalić and Launder (1972):

  It can be seen that for large Reynolds numbers, a ratio r of 1.1 means that the thickness Δ ext of the external layers is about 10% of d. The numbers of layers grows like ln Re.✄ The size of the cells in the directions parallel to the wall do not depend on the Reynolds number, but only of the geometry.

					0.025Re	-→ Re→∞	1 +	Δ ext d
	❀ Thus, for large Reynolds numbers,
				n =	ln Re ln r	+	ln (0.025(r -1)) ln r
	✄ n =	ln Re ln r	+	ln (0.025(r -1)/30) ln r		=	ln Re ln r	+	ln (0.025(r -1)) ln r	-	ln 30 ln r
	❀ Thus, for large Reynolds numbers, there is just a constant additional number of layers ln 30 ln r � 3.4 r -1 � 3.4 Δ ext /d to switch to the low-Reynolds number strategy.
	❀ This number only depends on the ratio	Δ ext

⇒ In total, the number of cells grows like ln Re, to be compared to Re 3 for DNS!

317

❀ With wall function, the thickness Δ 1 is larger by a factor of 30, such that d . For instance, for Δ ext = 10% d, there are about 34 additional layers, whatever the Reynolds number.

  30 (adaptive wall functions)� Remark: if this point is below the log layer (buffer layer, viscous sublayer), high-Reynolds number models are not valid ⇒ adaptive wall functions are necessarily used with low-Reynolds number models.

	⇒ Most of the CFD codes now propose such hybrid low-Re/high-Re strategies (adaptive wall
	functions (AWF), enhanced wall treatment, etc.): boundary conditions are used that blend wall
	functions and exact boundary conditions:

  The v 2 -f model is a drastic simplification of the Durbin's Reynolds-stress model using elliptic relaxation. It consists of replacing the eddy-viscosity

This model is available in StarCCM+, in Code Saturne and will be available soon in OpenFOAM. ❀ Eddy-viscosity model: v 2 -f model, proposed by Durbin in 1991, available in most of the commercial codes. ❀

�

  When Re τ is sufficiently large, there exists a region in which y + → ∞ and η → 0 are valid at the same time, i.e., In Eq. (103), the left hand side only depends on y + , while the right hand side only depends on η. Since they are not dependent on the same variables, and they are equal whatever the values of the variables, the two sides are necessarily constant.

	dU dy	=	u 2 τ ν	df dy + =	u τ h	dH dη
	or, multiplying by y/u τ	y + df dy + = η	dH dη	(103)
	�

� The constant is noted 1/κ, κ being called the von Kármán constant.

  Temperature profile in Rayleigh-Bénard convection. After

								515
	Comp. u 2 v 2	Prod. -2uw ∂U ∂z 0	Buo. prod. 0 0	Comp. uθ	U grad. prod. -wθ ∂U ∂z	T grad. prod. -uw ∂T ∂z	Buo. prod. 0
					vθ	0	0	0
	w 2	0		2βgwθ
	k	-uw	∂U ∂z	βgwθ	wθ	0	-w 2 ∂T ∂z	βgθ 2

  2ν t S ij + 2 3 kδ ij SGDH u i θ = -ν t Pr t ∂T ∂x i 5.6.1. Inconsistency of the model � Production term: G k = -βg i u i θ = βg i ν t Pr t ∂T ∂x i � Unstratified (or weakly stratified) flows:

  5.6.3. Variable turbulent Prandtl number?

	From Tennekes & Lumley [77]
	ν t Pr * t	∂T ∂y	= -	ν t Pr t	∂T ∂y	+ vθ Buo ⇒ Pr * t =	� 1 Pr t	+	1 ν t	vθ Buo ∂T /∂y

SGDH: u i θ = -ν t Pr t ∂T ∂x i with variable Pr t ? � Diffusion is due to mixing by large scales � These are the same scales for mechanical and thermal turbulence ⇒ Pr t = ν t κ t must be close to unity � Modifying Pr t for buoyant flows is a common practice (atmosphere/ocean) � Why? � Should have a buoyancy extension: vθ = - ν t Pr t ∂T ∂y + vθ Buo � To compensate: -

❀ It appears in the exact equation of the Reynolds stress tensor (production term).❀ It is crucial for describing the dynamics of turbulence (actually the most important parameter).

[START_REF] Bentaleb | A hybrid Temporal LES/RANS formulation based on a two-equation subfilter model[END_REF].

Statistical turbulence modeling (RANS)

Triple deomposition

URANS deomposition

411

� Introducing this decomposition, the equations are

� In the case of a statistically periodic flow (cyclo-stationary), as the case of the synthetic jet seen above, the Reynolds average and the phase average are identical ⇒ both lead to the same set of equations.

� This formal identity is the basis of the use in the URANS context of standard RANS models. � However, the turbulent stresses < u �� i u �� j > are subjected to an unsteady strain field ⇒ the equilibrium modeling can be brought into question.

SDM or OES modeling

SDM = Semi-Deterministic Modeling OES = Organized-Eddy Simulation � These two acronyms correspond to the same concept � The idea is to conceptually separate the contribution of the "organized structures" or "quasi-deterministic structures", from random turbulence.

From Ha Minh (1999) 415

� Formally speaking, a decomposition (using some filter) cannot be defined explicitly since, as indicated by the figure, the spectral signature of the structures we want to resolve is determined a posteriori by the solution.

� In the SDM (Ha Minh), as in the OES (Braza), the models are of k-ε type, in which the C µ coefficient is reduced:

to account for the misalignment of the anisotropy and strain tensors: � The difficulty lies in the conditions at the interface between the regions.

❀ The RANS model needs boundary conditions for the mean variables: the LES provides filtered variables depending on time ⇒ a temporal averaged is to be applied.

❀ This is no so simple in practice. The average must be applied during the calculation (on the fly). Transient phases must be first evacuated.

❀ Above all, the LES needs unsteady boundary conditions. The RANS computation cannot provide them directly.

425

� Generation of inlet conditions for LES. The "universal" behavior in the near-wall region 453

Equations written with internal variables

� A second possible choice consists in choosing so-called internal or wall variables, ν and u τ , to build the length and time units (note that u τ is the internal and the external velocity unit):

� Very simple formulations are also obtained for infinite Reynolds numbers:

� Now, it is seen that U and uv only depend on u τ and yu τ /ν, and not on the channel half-width h anymore.

� (97) then becomes Here:

� A transport equation for ζ i can be obtained using

527 � Denoting P θi , P k = P + G, P θ 2 the production terms of u i θ, k and θ 2 , and D θi , D k , D θ 2 their diffusion terms (molecular + turbulent), we have

� It can be noted that ❀ this is an implicit algebraic model (u i θ's appear on both sides of the relation);

❀ this model accounts for the influence of the three production terms;

❀ the scrambing term is accounted for;

❀ if only the first term is considered, the GGDH model is recovered. � Buoyancy also affects redistribution, dissipation, scrambling.

Appendix C

How to perform a good computation 537 � The acceptable blockage ratio depends on the error we can accept, and the price we are ready to pay (in terms of number of cells).

❀ The same problem is faced in experiments: wind tunnels are of limited size. In order to work with not too small scale models, it is common to accept the compromize B � 5%

❀ In CFD, it is much easier and cheaper to increase the size of the domain: cells far from the walls are large ⇒ we can choose the constraint B < 1%

✄ The minimal value for the domain height H in the present case is thus

Location of the inlet boundary

� A uniform inlet profile U = U inlet is usually applied at the inlet.

� Due to the presence of the obstacle, the streamlines are deviated upstream of the leading edge.

� Imposing a uniform inlet profile thus introduces an error.

� If the inlet is too close from the leading edge, this error can be significant.

� It can be shown for an incompressible flow of a perfect fluid that the velocity on the stagnation streamline decreases as a function of the distance d to the stagnation point

for a sphere (3D) of radius R

Cylinder Sphere

557

❀ Assuming that we can, for the purpose of this evaluation, approximate the airfoil by a cylinder of diameter e, the minimal distance between the inlet plane and the leading edge such that

❀ This is only an evaluation. This must be checked a posteriori.

❀ It is convenient to check the pressure coefficient given by the computation:

✄ Bernoulli's theorem is roughly applicable upstream of the obstacle (the flow is not turbulent, the dissipation is negligible):

such that C P varies between 0 (far upstream) and 1 (stagnation point).

✄ P � = P ∞ at the inlet due to the presence of the obstacle ⇒ C p � = 0.

✄ Since δC p � -2 δU U , the value of C p provides an estimate of the relative error on the velocity.

❀ It is wise to recompute the flow with a different inlet location.

� Example: flow through a diaphragm ❀ Here, the quantity of interest is the head loss coefficient ξ = ΔP ❀ If our criterion is to keep below 1% the error on ξ due to dicretization errors, it can be seen on the plot that the mesh is sufficiently fine.

565

� How to refine the mesh ❀ After a first computation, we want to refine the grid to investigate grid convergence. How can we proceed?

❀ The basic method consists in increasing the number of cells in each direction.

However:

✄ For Low-Reynolds number models: we have chosen the thickness of the wall-adjacent layer such that y + c � 1; it is useless to refine the mesh, except for some models that can require slightly finer meshes.

✄ For High-Reynolds number models: the wall-adjacent layer cannot be refined, since y + c > 30 is necessary.

⇒ In the wall-normal direction, it is better to keep the thickness Δ 1 of the first layer constant, and to reduce r. � Computation stops when the relative residuals are below a prescribed value.

� Be careful! Default values are often too large (� 10 -4 ). This is a major source of numerical error.

� Actually it is better to wait for the residual to reach some plateau: very often, residuals cannot be further reduced.

� Monitoring the evolution of quantities of interest (e.g., the drag or the lift for an airfoil) is necessary to avoid unpleasant surprises.