Comptes Rendus Mécanique

. A stable scheme is proposed in this paper in order to obtain approximate solutions of second-momentturbulentmodelsforincompressibleflowswithorwithoutthermaltransportequation.Theanalysis of the convective terms, which includes the solution of the associated Riemann problem, enables to propose a standard projection scheme, and to get rid of spurious oscillations.


Introduction
Reynolds-averaged Navier-Stokes models are commonly used in the industrial and environmental applications in order to predict incompressible or compressible flows (see e.g.[1,2]) due to the high Reynolds numbers or Grashoff numbers encountered.These equations require closure on the Reynolds stress tensor (covariance of velocity fluctuations), but second-order closures, which consist in solving a transport equation on this symmetric tensor, are less used.Temperature or pollutant averaged transport equations also provide turbulent thermal/scalar fluxes to be modelled.All these second-order models contain first-order non-conservative terms (see e.g.[3][4][5][6]).When focusing on turbulent compressible flows, a stable numerical strategy has been proposed, which is grounded on the hyperbolic structure of the full-convective system embedded in the model (see [7][8][9][10]).More recently, this methodology has been extended to shallow-water models [6,[11][12][13].
The main purpose of the present work consists in applying similar ideas, while considering incompressible flow models including second-moment closures, in order to get rid of nonphysical oscillations which occur in many practical situations (see e.g.[14]).
We first recall the set of equations and the realizability constraint on the Reynolds stress tensor, which is mandatory to guarantee the hyperbolicity of the system.Then, we propose a basic timestepping and a special focus is given on the analysis of the full-convection part of the models.Finally, the strategy is demonstrated using a Rusanov-like solver on a wall bounded flow.This strategy is extended to anisothermal flows.

Set of equations
The governing set of equations involves mass and momentum conservations and Reynolds stress tensor R := u ′ ⊗ u ′ transport equation: with, for Newtonian fluid and incompressible flow, Σ := − p ρ 0 1 + ν 0 (∇ ∇ ∇u + ∇ ∇ ∇u T ) the mean stress tensor per mass unit, and p denotes the mean mechanical pressure.
As a covariance matrix, the Reynolds stress tensor is symmetric positive and half-definite.This refers to the so-called realizability constraint (see [1,15]).
Models for G R := u ′ ⊗ u ′ ⊗ u ′ usually rely on a classical first gradient assumption G R = −C R ν t (R, T )∇ ∇ ∇R, where C R is a positive constant, ν t is the turbulent viscosity, and T is the time scale related to turbulent dissipation.

Total kinetic energy tensor
The evolution equation on R can be combined with the kinetic energy of the mean velocity field u ⊗ u to give the total kinetic energy u ⊗ u = u ⊗ u + R evolution: Equation (3) can also be deduced directly from the mean of the instantaneous kinetic energy tensor, and emphasise that the so-called non-conservative production term −(R • ∇ ∇ ∇u + ∇ ∇ ∇u T • R) acts as energy transfer between kinetic energy of the mean velocity field and turbulent kinetic energy (note that the transfer may not be in the same direction) but is transparent with respect to the total kinetic energy.This is analogous to what happens to compressible flows with the non conservative term pdiv u term which cancels out between internal energy and kinetic energy balance.Therefore, the total kinetic energy tensor will be used in the following to derived jump conditions when dealing with non regular solutions.

Realizability constraint
To examine the realizability of the Reynolds stress tensor, we use the following Lemma: Lemma 1. Assume that the governing equation of a symmetric second-rank tensor R ∈ R 3 × R 3 is: Let δ R be the determinant of the latter tensor R. Then the governing equation of δ R is: A proof can be found in [16, Appendix 1].H is deduced from the terms R•∇ ∇ ∇u +∇ ∇ ∇u T •R, Φ r and Φ s , and the boundedness of the trace of H must be examined to ensure positivity of δ R .In case of the basic Lumley model examined in Section 3 where φ r = 0 and for incompressible flows, we have t r (H) = div u = 0 which implies that ∂ t δ R remains strictly positive for suitable initial and boundary conditions, thus all the eigenvalues of R remain strictly positive.This is also true for the same model in the comrpessible framework if div u is bounded.

A first order time scheme
Starting with initial conditions (u n , R n ) and using suitable boundary conditions, (u n+1 , R n+1 ) arises from the following first-order time scheme: A Stokes-like step has been applied in the scheme (4a) and (4b) (see, among others, [17][18][19]).The space discretization associated with right-hand-side of (4b), (4c) is derived from the analysis of the evolution system (5), which involves all convective terms:

Analysis of the evolution step
In the following, we restrict to the simple Lumley's model which does not contain so-called rapid terms Φ r = 0. We also restrict from now on to the two-dimensional framework to ease the presentation.We first recall that the latter system (5) is hyperbolic if the Reynolds stress tensor is realizable (see [20]).Eq. ( 5) is invariant under frame rotation, so we can rewrite it in the (n, τ) reference frame (where n, τ are some unit orthogonal vectors), defining x n := x • n, u n := u • n, u τ := u • τ, R nn := n T Rn, R nτ := n T Rτ and R ττ := τ T Rτ.We now define the evolution step in the n−direction, by neglecting the transverse derivatives, which yields: We set W := (u n , u τ , R nn , R ττ , R nτ ).

Hyperbolicity of the evolution system
We assume from now on that initial and boundary conditions comply with the strictly realizability constraint, which will imply that R remains strictly realizable (see [21]).As system (6) preserves the realizability, it implies that R nn > 0 whatever the unit vector n is.
Proposition 2 (Hyperbolic evolution step).System (6) is hyperbolic.It admits five real eigenvalues which are: and the associated right eigenvectors span R 5 .
We refer to [20,22] for a proof.This can be extended to the three-dimensional framework.

Riemann invariants
The Riemann invariants

Shock waves
Due to non-conservative terms in eq. ( 6), shock relations cannot be derived in a classical way.Following [3,4] approximate shock relations can nevertheless be proposed.These relations are valid in the limit of weak shocks.Assuming a linear path from left to right with respect to the variable W = (u n , u τ , R nn , R ττ , R nτ ), these relations write: where [x] := x r − x l denotes jump and x := x r +x l 2 is the arithmetic mean between left and right states of the discontinuity travelling at speed σ.
Approximate relations (10) can also be obtained from transport equations written in pseudoconservative variables R + u ⊗ u instead of R introduced in Section 1.2, as explained in [22, p. 18 and Appendix 4].

An approximate solution of the one-dimensional Riemann problem
Theorem 3 (Existence and uniqueness of the approximate solution).For the approximate shock relations (10), there exists a unique self-similar realizable solution (W(x n , t ) = w( x n t ), with W = (u n , u τ , R nn , R ττ , R nτ )) to the Riemann problem associated with system (6) supplemented with strictly realizable initial conditions W(x n , 0) = W L for x n < 0 and W(x n , 0) = W R for x n > 0, if and only if the following condition holds: As in the compressible framework [7][8][9][10], the Riemann problem structure of system (5) allows to specify wall, symmetry and inlet/outlet boundary conditions.
Proof.The proof follows the main guidelines recalled in [23].A detailed expression of intermediate states are fully given in the sequel for the anisothermal case in Section 6.
Let us consider two initial realizable states W R and W L .
Step 0. Preliminary calculation.To connect left and right states W l ,r across the 1 st and 5 t h waves we set To select the physical solution for shocks (10), we use [22,Appendix 4]).We can also write from Section 1.2: in shocks, and therefore, using (10), we get: ( u 2 n + R nn )[u n ] ≤ 0, and thus once more [u n ] ≤ 0. For the 1-shock connexion (z > 1): For the 5-shock connexion (0 < z < 1): To connect the left and right states through the 1-rarefaction wave (respectively 5-rarefaction wave) we use the Riemann invariants (9a) (respectively (9e)) still following [23].
Step 1. Solution in terms of u n and R nn variables A glance at (9b), (9c), (9d) shows that u n and R nn are constant across the three LD waves.First we can focus only on u n and R nn variables.We denote their intermediate values by u ♯ n and R ♯ nn , and link these to initial states W R,L : h 1,5 (z) are defined using eqs.(9a), (9e), ( 12) and ( 13), following [23].Note that . We set: Solving eq. ( 15) gives z 5 , hence R Step 2. Solution in terms of u τ , R nτ , and R ττ variables For these variables, the connexion from left to right can be done using z 1 and z 5 , and exploiting the Riemann-invariants across the 2 nd , 3 rd and 4 th LD waves.R nn and R nn R ττ − R 2 nτ remain constant across the 2 nd and 4 th waves (see (9b),(9d)), thus R is realizable.Moreover u τ does not vary across the 3 rd wave (see (9c)).□

Numerical scheme for the evolution system
The objective of the present work is to obtain a stable numerical scheme (unlike those usually used as presented in [14], which may use partial upwinding which respect to the material mean velocity).Therefore we restrict the presentation to low order schemes such as Rusanov scheme.
Other approximate Godunov schemes and MUSL techniques might also be considered.
For the sake of brevity the scheme is presented in a one dimensional framework.The evolution step (5) can be written (using the incompressibility constraint): W n+1 i will be computed using the following Rusanov-like [24] scheme: where the time step ∆t n satisfies the CFL condition max i (|λ and the non-conservative contribution N n i are: where the discrete divergence free condition on (u n ) n i + 1 2 holds.

Wall bounded flow test case
The following test case describes a flow in the vicinity of a wall, while prescribing shear stress.The initial states W L,R and the intermediate states (I , I I , I I I , I V ) arising in the 1-dimensional Riemann problem associated with ( 6) are given in Table 1.Computations have been performed using code_saturne finite volume platform [25].The CFL parameter is set to 0.5 and the meshes contain from 250 to 32000 cells.The L1-error and the behaviour of W are shown in Figure 3, using material-upwind scheme and scheme ( 17) and ( 18) accounting for all convective effects.Table 1.Analytical solution for the tangential flow in the vicinity of a wall.

Extension to anisothermal flows with second-order turbulence models
System ( 1) is now supplemented with mean temperature θ conservation eq.(19a), turbulent heat fluxes θ ′ u ′ and temperature variance governing eqs.(19b) and (19c): Focus is still given on first-order gradient closure laws for G θ ′ 2 and G u ′ θ ′ .We consider the one dimensional Riemann problem involving eq. ( 5) with Φ r = 0, supplemented with the first-order terms of eqs.(19a) to (19c) (thus getting rid of Using the invariance under frame rotation, the counterpart of ( 6) can be derived by setting to zero the transverse derivatives, which yields: Proposition 4 (Hyperbolic evolution step).System (20) is hyperbolic.It admits nine real eigenvalues which are: and the associated right eigenvectors span R 9 .Fields associated with eigenvalues λ 2,3,4,5,6,7,8 are linearly degenerate (LD).The 1-field and 9field are genuinely non linear (GNL).
We refer to [20,22] for a proof.

Riemann invariants
The Riemann invariants I i R associated with the i th field ( f ∈ I i R , ∇ ∇ ∇ f • r i = 0) are: and consider the approximate shock relations eqs.(10) and (27).Assume strictly realizable initial conditions W(x n , 0) = W L for x n < 0 and W(x n , 0) = W R for x n > 0, then there exists a unique selfsimilar realizable solution: W(x n , t ) = w( x n t ) to the Riemann problem associated with the evolution step involved in eqs.(5) and (19) if and only if condition (11) holds.
Proof.The proof follows the main guidelines recalled in [23] and is only partially detailed in [22].We consider two initial realizable states W R , W L .We look for the four intermediate states labelled I , I I , I I I , I V from left to right separating the 5 distinct waves.
Step 0. Preliminary calculation To connect left and right states W l ,r across the 1 st and 5 t h waves we set z := R r nn R l nn > 0. To select the physical solution for shocks (10), we use [u n ] = u r n − u l n ≤ 0. For the 1-shock connection (z > 1): For the 9-shock connection (0 < z < 1): To connect the left and right states through the 1-rarefaction wave (respectively 9-rarefaction wave) we use the Riemann invariants (9a) (respectively (9e)).
Step 1. Solution in terms of u n and R nn variables A glance at (9b), (9c), (9d) shows that u n and R nn are constant across the three LD waves, as for the isothermal case.First we can focus only on u n and R nn variables.We denote their intermediate values by nn , and link these to initial states W R,L : h 1,9 (z) are defined using (23a), (23e), ( 28), (29), following [23].Note that . We set: Equation ( 32) is obtained from Riemann invariants I 2,3 R , I 7,8 R .
Step 3. Solution in terms of R ττ , θ ′2 and u ′ τ θ ′ variables.Then, we deduce the remaining two variables R ττ and θ ′2 for states I I (resp.I I I ) from Riemann invariants of the 2 − 3 wave (resp.7 − 8 wave): Both left and right states are realizable and they satisfy the condition of Theorem 5 for existence and uniqueness of the solution.
The results of the simulation where we applied the Rusanov-like scheme are presented in Figure 1.As previously we observe the absence of oscillations due to the application of the Rusanov-like scheme.Figure 1j shows a one-half convergence rate, as expected due to the contact Table 2. General solution for the tangential flow in the vicinity of a wall considering the transport of a scalar quantity.
One can also notice that the Rusanov-like scheme does not recover the exact constant state for R ττ , u ′ τ θ ′ and θ ′2 , similarly to what one gets on pressure profiles with Rusanov scheme for Euler equations with perfect gas equation of state.However this artefact disappears with space-time convergence.), R nn (c.), R R nτ (e.), θ (f.), u ′ n θ ′ (g.), u ′ τ θ ′ (h.), θ ′2 (i.) along with convergence graph (j.) obtained applying the Rusanov-like scheme.Point-source dispersion of a pollutant using a material-upwinding scheme and the Rusanov-like scheme.

Point-source dispersion of a pollutant
We present here the steady state of a point-source dispersion of a passive pollutant in a steady and homogeneous dynamic.The following hypotheses were made: • the pollutant is introduced using a Dirac mass and continuously in time at a concentration value of 1 (computationnaly, it is introduced in a single cell), • molecular diffusion is neglected in front of turbulent diffusion, One can observe in Figure 2 the results obtained by applying material-upwind scheme which present spurious spatial steady oscillations, and are suppressed by the Rusanov-like scheme.

Conclusion
A numerical strategy has been proposed to compute incompressible flows with second-order moment turbulence closures, in order to remove numerical oscillations.It is based on the analysis of the full convective part of the system.As in the compressible framework, the realizability of the Reynolds stress tensor is mandatory to obtain the hyperbolicity of the evolution system.
The numerical convergence has been demonstrated on a verification test case using a Rusanovlike scheme.It is also suitable for anisothermal flows or mean species transport equations with second-moment turbulence closures.
Of course, more accurate Riemann solvers, or relaxation solvers, and second order extension can be investigated.This methodology can be applied to projection step methods or Uzawa algorithms.
Eventually, more complex models involving non-zero contribution Φ r (∇ ∇ ∇u, R) may be considered, using the same approach (see for instance [22,Appendix 2]).

♯
nn and u ♯ n .The solution is unique since Ψ is a strictly monotonic function.Its lower and upper bounds enable to conclude that u ♯ n and R ♯ nn exist if (11) holds.

Step 2 .
Solution in terms of u τ , R nτ , and θ, u ′ n θ ′ variables As u τ , R nτ are invariants of 4-5-6 wave, we set u ⋆ τ = u I I τ = u I I I τ and R ⋆ nτ = R I I nτ = R I I I nτ , given by:

Figure 2 .
Figure 2. Point-source dispersion of a pollutant using a material-upwinding scheme and the Rusanov-like scheme.

Figure 3 .
Figure 3. Profiles of the variables u n , u τ , R nn , R ττ , and R nτ along with convergence plots obtained by applying material-upwind scheme and Rusanov-like scheme (500 cells, CFL= 0.5).