C. R. Acad. Sci. Paris, t. 331, Série I, p. 1011–1016, 2000Analyse numérique/Numerical Analysis

Comparaison de solveurs numériques pour simulerla turbulence compressible

Emmanuelle DECLERCQ a, Alain FORESTIER b, Jean-Marc HÉRARD c,d

a CEMIF, 40, rue du Pelvoux, Courcouronnes, 91020 Evry cedex, Franceb CEA Saclay, DRN/DMT/SEMT, 91191 Gif-Sur-Yvette cedex, Francec Département MFTT, division Recherche et développement, EDF, 6, quai Watier, 78400 Chatou, Franced CMI, Université de Provence, 39, rue Jolio-Curie, 13453 Marseille cedex 13, France

(Reçu le 10 décembre 1999, accepté après révision le 16 octobre 2000)

Résumé. Nous nous intéressons à la résolution numérique du système hyperbolique décrivant unécoulement turbulent, compressible, de gaz en situation isentropique. Nous construisons lesolveur de Riemann exact associé, nécessaire au schéma de Godunov, que nous comparonsau schéma Vfroe, et au schéma de Rusanov. 2000 Académie des sciences/Éditionsscientifiques et médicales Elsevier SAS

Comparison of numerical solvers for turbulent compressible flow

Abstract. We focus on the computation of the hyperbolic system describing a turbulent flow forisentropic gas. An exact Riemann solver is computed, we construct the solution of theRiemann problem and thus we derive a Godunov scheme. We provide some numericalsimulations to exhibit a comparison between Godunov scheme and two approximatesolvers: Vfroe and Rusanov scheme. 2000 Académie des sciences/Éditions scientifiqueset médicales Elsevier SAS

Abridged English version

In this Note, we focus on the non-conservative hyperbolic system arising from theK turbulentcompressible model applied to multicomponent flow and polytropic gases. Variablesρ, Y , U andKdenote respectively the mean density of the mixture, the mass fraction of one component in the mixture,the velocity and the kinetic turbulent energy of the mixtureK = 1

2 ρU′U ′. Governing equations of the

non-conservative system are given by (S). We get the entropy-entropy flux pair(E , fE) =(ρU2

2 + K +Pγ−1 ,EU + 2

3 KU + PUγ−1

)and by the entropy inequalityσ[E ] > [fE ] on shock curves ofσ discontinuity

celerity, the unique solution is exhibited. The exact solution of the Riemann problem enables to implementthe Godunov scheme. We refer to Coquel and Berthon [1,2], for a somewhat different approach to dealwith non-conservative equations, arising in systems with two distinct entropies. The aim of this study is

Note présentée par Olivier PIRONNEAU .

S0764-4442(00)01745-6/FLA 2000 Académie des sciences/Éditions scientifiques et médicales Elsevier SAS. Tous droits réservés. 1011

E. Declercq et al.

the comparison between Godunov scheme [10] and two linearized methods. The first one, Vfroe-ncv (see(1)–(5) and [3]), is very close to Godunov method, and thus gives very accurate results. It divides the CPUcost by a factor ten. However, this latter method has some drawbacks, since both pressure and turbulentkinetic energy positivity can not be ensured in all cases. This has been observed, for instance, in somestrong rarefaction waves occurring behind some bluff bodies. The Rusanov scheme (see(6)–(9) and [12])ensures positivity of the variablesρ, Y , K , but this scheme is so diffusive that it provides poorly accurateresults. We conclude that in some extreme cases, the cost of an exact solver is justified by the robustnessrequirement. Anyway, we advocate to use successively two schemes: the exact one for initialization and todeal with the boundary conditions, and a linearized one to go on with a good compromise between accuracyand computational cost.

Introduction. – On considère un écoulement isentropique turbulent de deux composants gazeux. Danscette étude on s’intéresse particulièrement au couplage pression-turbulence. La turbulence étant supposéeisotrope, le tenseur de Reynolds est décrit par l’énergie cinétique turbulente du mélangeK . La partieconvective du système étant dominante, on ne traitera que la fermeture au premier ordre du modèle (cf. [11]pour un modèle plus complet). Nous décrirons la construction du solveur de Riemann exact associé ausystème et justifierons l’utilisation du schéma de Godunov par un cas test mettant en défaut les solveursapprochés.

1. Solveur de Riemann exact et schéma de Godunov

On s’intéresse à la résolution du système portant sur les variables moyennées par le processus de Favre :ρ est la densité du mélange,U la vitesse du mélange,Y la fraction massique d’un des deux composantsetK la turbulence sommée des deux fluides. La pressionP est connue par la loi d’étatP = P (ρY ). OnnoteC = (ρ, ρY, ρU)t,W =W (x, t) = (C,K)t, avecx ∈Ω (⊂R2), t> 0.

C =

ρρYρU

, F (C,K) =

ρUρY U

ρU ⊗U +(

23 K + P

)I

,

(S)

∂tC +∇F (C,K) = 0,

∂tK +∇(KU) + 23 K∇U = 0.

Le système différentiel (S) vérifié parW s’écrit aussiW,t+∑i=1,2AiW,i = 0. Par invariance par rotation,

on se ramène au système (Sn) monodimensionnel correspondant à l’écriture de (S) dans un repère localreposant sur la normale~n. Les composantes normale et tangentielle deU sont désignées parun et ut.SoitM la matrice de passage dans le repère normal,MW = Wn = (ρ, ρY, ρun, ρut,K)t. On définitAnparAn =

∑i=1,2(MAiM

−1)ni et on note∂n =∇ · ~n. On s’intéresse au problème de Riemann suivant :

(PR)

∂tWn +An∂nWn = 0 (Sn),Wn(x, t= 0) =Wg si x · n < 0,Wn(x, t= 0) =Wd si x · n > 0.

La recherche de toutes les entropies mathématiques du système(C,K/ρ2/3,E) conduit à ne garder parmicelles-ci que l’entropie physique du systèmeE de fluxfE (cf. [6]). E entropie convexe (linéaire par rapportàK), correspond à l’énergie totale du système :

E =ρU2

2+K + ρ

∫P (ρY )

ρ2dρ, fE =

ρU2

2U +

5

3KU + ρU

∫P (ρY )

ρ2dρ+PU.

La partie admissible des ondes de choc est sélectionnée en écrivant la croissance de l’entropie aux chocs.

1012

Solveurs numériques en turbulence compressible

La vitesse du son dans l’écoulement bifluide turbulent est notéec′ =√Y c2 + 10

9Kρ et c2 = P ′(ρY ),

(Sn) admet pour valeurs propresλ1 = un − c′, λ2,3,4 = un, λ5 = un + c′.Le système (Sn) est hyperbolique (non strictement). Les champs caractéristiques sont soit linéairement

dégénérés, soit vraiment non linéaires à la condition suffisante queP (ρY ) soit une fonction convexe duvolume spécifique1/ρ. On en déduit d’après [9] que la solution est constituée d’au plus quatre étatsconstants séparés par des ondes de choc ou détente et une discontinuité de contact correspondant à l’ondetriple. Au passage de la discontinuité de contact séparant les états1 et 2, les sauts suivants sont nuls (avecla convention[ϕ]ba = ϕ(b)− ϕ(a)) : [un]21 = 0 et [2K/3 + P ]21 = 0. La construction des ondes de détenterepose sur le calcul des invariants de Riemann.

On donne l’expression des ondes de détente [14] :Dg entre l’état gaucheg et l’état1 etDd entre l’état

droit d et l’état2, avec pourk ∈ g, d, c′k(x) =√Ykc2(Ykx) + 10

9Kkρ

5/3

k

· x2/3 :

Dk(Uk) =

(ρ,Y, uτ ,K,un); ρ > 0, Y = Yk, uτ = uτk, K =

Kkρ5/3

ρ5/3k

, un = unk +

∫ ρk

ρ

ck′(x)

xdx

.

Quant aux ondes de choc, elles sont obtenues en écrivant les relations de saut de Rankine–Hugoniot aufranchissement d’une discontinuité. Dans [5], on définit le produit non conservatif, noté[ ]S , associé auchemin des lignes droitesS comme suit :

[K∂nun]S =

∫ 1

0

(Kg + s(Kd −Kg)

)(ud − ug) ds=

Kg +Kd

2[un].

Les ondes de choc sont données avec la notationCk(Uk) pourk ∈ g, d avecCg entre l’étatg et l’état1etCd entre l’état2 et l’étatd :

Ck(Uk) =

(ρ,Y, uτ ,K,un); ρ > 0, Y = Yk, uτ = uτk, K =

4ρ− ρk4ρk − ρ

Kk,

un = unk −

√(ρ− ρk)

[23 (K −Kk) + P −Pk

]ρρk

.

La partie admissible des ondes de choc a été sélectionnée en écrivant la croissance de l’entropie auxchocs. La démonstration de l’existence de la solution du problème de Riemann repose sur la constructionparamètrée de la solution, du même type que celle de la méthode de Smoller [13] (p. 353). On évitecependant de paramétrer le passage des états gauche et droit aux états intermédiaires par des fonctionsexponentielles, trop coûteuses dans les implémentations informatiques. La résolution du problème deRiemann de façon exacte se ramène alors à un problème de recherche du zéro d’une fonction non linéairedeR2, qui grâce à ses propriétés de monotonie par rapport à ses deux arguments, permet une résolutionnumérique par une méthode de type Newton [7,6].

On cherche à résoudre ce système sur un maillageΩi=1,...,N (Ω =⋃iΩi, Ωi ∩ Ωj = ∅ si i 6= j,

i ∈ 1, . . . ,N). SoitWni la solution discrète sur la mailleΩi, Wn

i =def

1|Ωi|

∫ΩiW (x,n∆t) dx, où ∆t est

le pas de temps soumis à une condition CFL. La construction de la solution exacteWR du problème deRiemann, est nécessaire au calcul des flux du schéma de Godunov [10]. On noteW ∗gd = (C∗gd,K

∗gd) =

WR(Wg,Wd, x/t = 0), `ij = |∂Ωi ∩ ∂Ωj| et V (i) est l’ensemble des cellules voisines de la mailleΩi.nij est la normale à l’interface des maillesΩi, Ωj , orientée de la mailleΩi versΩj .

Ki =

( ∑j∈V (i)

K∗ij

)/( ∑j∈V (i)

1

).

1013

E. Declercq et al.

Le schéma de Godunov est alors donné par la formule :Cn+1i =Cni −

∆t

|Ωi|∑

j∈V (i)

F (W ∗ij) · nij `ij ,

Kn+1i =Kn

i −∆t

|Ωi|∑

j∈V (i)

(K∗iju

∗nij +

2

3Kiu

∗nij

)`ij .

Cette méthode numérique s’applique sur tout type de maillage, structuré ou non [7].

2. Le schéma Vfroe-ncv

Une confrontation du schéma de Godunov à des solveurs approchés permet de rendre compte de l’intérêtd’utiliser un solveur exact même plus coûteux. On s’est intéressé tout d’abord au solveur linéarisé duschéma Vfroe-ncv. Le schéma Vfroe a été introduit en 1996 dans [8]. Il est basé sur une résolution localedu problème de Riemann linéarisé. Une extension de ce schéma aux variables non conservatives a été faitedans [3]. Ici, pour préserver les invariants de Riemann aux franchissements des discontinuités de contact,nous prendrons la variableZ = (Y,un, ut,K,P )t, τ = 1/ρ. On résout le problème de Riemann linéarisé(PRL) associé à la matrice constanteAn(Z) avecZ = (Zg +Zd)/2, moyenne des étatsg etd.

(PRL)

∂tZ +An(Z)∂nZ = 0,Z = Zg si x · n < 0,Z = Zd si x · n > 0.

(1)

Avec c ′2 = (γP + 109 K)τ , les valeurs propres deAn(Z) sontλ1(Z) = un− c ′, λ2,3,4 = un, λ5 = un + c ′

Zd = Zg +∑i=1,5

αi(Z)ri(Z). (2)

On noteZ∗gd(Zg,Zd) la solution à l’interfacex/t= 0 du (PRL), avecZ1 etZ2 les deux états intermédiaires,Z∗gd = Zg si un − c ′ > 0, Z∗gd = Z1 si 0< un < c ′,

Z∗gd = Z2 si−c ′ < un < 0, Z∗gd = Zd si un + c ′ < 0,(3)

Z1 = Zg +

1

2c ′

[un]dg −

τ

c ′

[2

3K + P

]dg

r1(Z),

Z2 = Zd −1

2c ′

[un]dg +

τ

c ′

[2

3K + P

]dg

r5(Z).

(4)

L’extension du schéma Vfroe au système non conservatif est donnée par :Cn+1i =Cni −

∆t

|Ωi|∑

j∈V (i)

F (Z∗ij)nij`ij ,

Kn+1i =Kn

i −∆t

|Ωi|∑

j∈V (i)

(K∗iju

∗nij +

2

3Kiu

∗nij

)`ij .

(5)

Le solveur linéarisé peut générer des états intermédiaires qui ne respectent pas la positivité des variablesphysiquesρ, P ouK . Par exemple, à l’état miroir d’une double raréfaction symétrique (correspondant àune condition limite de paroi en détente) pour les conditions initiales suivantes :(Y,un, ut,K,P ) pourx< 0 et (Y,−un, ut,K,P ) pourx> 0 avecun < 0, ~n étant la normale sortante, on obtient :

K1 =K

(1 +

5

3

unc ′

), P1 = P

(1 + γ

unc ′

), doncK1 < 0 pour

unc ′

<−5

3etP1 < 0 pour

unc ′

<− 1

γ.

De plus, on ne peut pas garantir la positivité de maille des valeursρ, Y , P , K .

1014

Solveurs numériques en turbulence compressible

3. Le schéma de Rusanov

On considère le système hyperbolique non conservatif suivant :

∂tW +∑i=1,2

(Hi(W )

),i

+∑i=1,2

BiW,i = 0, (6)

dont la matrice associée estAn :

An =∑i=1,2

∂Hi(W )

∂W· ni +

∑i=1,2

niBi. (7)

On notera les valeurs propresλk et on introduitSi,j = maxi,j(maxk(|λk(Wi)|, |λk(Wj)|)

)et Wij =

(Wi +Wj)/2.Le flux numériqueFR du schéma de Rusanov [12] a alors l’expression suivante :

FR(Wi,Wj , nij) =1

2

(F (Wi)nij + F (Wj)nij − Si,j(Cj −Ci)

). (8)

Le schéma de Rusanov, appliqué à notre système (S) s’écrit :Cn+1i =Cni −

∆t

|Ωi|∑

j∈V (i)

FR(Wni ,W

nj , nij

)`ij ,

Kn+1i =Kn

i −∆t

|Ωi|∑

j∈V (i)

((Kun)ij −

1

2Sij(Kj −Ki) +

2

3Kiuijnij

)n`ij .

(9)

Le schéma de Rusanov a été choisi pour notre système parce qu’il respecte la positivité de maille [14] deρ, ρY etK et parce qu’il assure également à la fraction massique de rester comprise dans le bon ensemble0< Y < 1, sous condition CFL :∆t

∑j∈V (i)

(Sij + 2

3 ujnij)`ij < 2|Ωi|.

4. Résultats numériques

Grâce au schéma linéarisé Vfroe, le temps de calcul se trouve divisé par10, tout en donnant unebonne précision des résultats. Si l’on compare la normeL2 d’une solution2D avec celle calculée par laméthode de Godunov on constate un pourcentage d’erreur inférieur à1.5. Avec le schéma de Rusanov, ondivise encore le temps de calcul d’un facteur3. De plus, le principe du maximum pour la concentrationest respecté. On traite le cas d’une tuyère à rétrécissements (cf. figure 1) pour des conditions initiales :(ρ,Y, u, v,K,P ) = (1.1,1,600,0,1000,100000) sur maillage régulier de600 quadrangles. Ce cas testrévèle les avantages et inconvénients des différentes méthodes. Avec le schéma de Rusanov aucun problèmede positivité n’apparaît, mais l’évaluation deK se révèle assez approximative. Quant au schéma Vfroe-ncv, il donne de bons résultats, mais pour le traitement des fortes détentes, notamment lors de simulationd’arrière corps où l’écoulement se raréfie, il est préférable d’utiliser le schéma de Godunov pour garantir lapositivité de la densité et de la turbulence. On donne ici une comparaison de la turbulence calculée par laméthode de Godunov avec celle obtenue par le schéma Vfroe-ncv. On donne en (a) la turbulence calculéepar le schéma de Godunov et en (b) le pourcentage d’écart, calculé sur la turbulence, du schéma Vfroe-ncvpar rapport au schéma de Godunov.

Conclusion. –Nous préconisons, en général, une utilisation du solveur Vfroe corrélée à la résolutionpar solveur exact pour le traitement des frontières du domaine. Pour certains cas physiques plus« pathologiques » (fortes détentes), nous préférerons le schéma de Godunov, plus robuste. Ce travail sepoursuit par la prise en compte des termes visqueux (laminaires et turbulents) tout en gardant un traitementconvectif couplé (pression-turbulence). Un travail similaire, sur un modèle de turbulence au second ordre,a montré la supériorité d’une telle approche couplée, par rapport à une approche numérique basée sur undécouplage des équations, de fortes instabilités au niveau des conditions limites de paroi apparaissant avecle modèle découplé [4].

1015

E. Declercq et al.

Figure 1. – (a)K par le schéma de Godunov ; (b) Écart surK entre Vfroe-ncv et Godunov.

Figure 1. –(a)K by Godunov scheme; (b) Difference betweenK by Vfroe-ncv andK by Godunov scheme.

Références bibliographiques

[1] Berthon C., Contribution à l’analyse numérique des équations de Navier–Stokes compressibles à deux entropiesspécifiques. Application à la turbulence compressible, Thèse de doctorat, Université Paris-VI, 1999.

[2] Berthon C., Coquel F., Convection diffusion system with first and second order in non-conservation form, (enpréparation).

[3] Buffard T., Gallouët T., Hérard J.-M., A sequel to a rough Godunov scheme: application to real gases, Computersand Fluids 29 (7) (2000) 813–847.

[4] Brun G., Hérard J.-M., Jeandel D., Uhlmann M., An approximate Roe-type Riemann solver for a class of realizablesecond order closures, Int. J. of Comput. Fluid Dynam. 13 (3) (2000) 213–249.

[5] Dal Maso G., Le Floch P., Murat F., Definition and weak stability of nonconservative products, J. Math. Pures etAppl. 74 (6) (1995) 483–54.

[6] Declercq E., Forestier A., Hérard J.-M., Part 1: An exact Riemann solver for a multicomponent turbulent flow,Part 2: Comparison of numerical solvers for a multicomponent turbulent flow, soumis à Int. J. of Comput. FluidDynam.

[7] Forestier A., Hérard J.-M., Louis X., Solveur de type Godunov pour simuler les écoulements turbulentscompressibles, C. R. Acad. Sci. Paris, Série I 324 (1997) 919–926.

[8] Gallouët T., Masella J.-M., Un schéma de Godunov approché, C. R. Acad. Sci. Paris, Série I 323 (1996) 77–84.[9] Godlewski E., Raviart P.-A., Numerical Approximation of Hyperbolic Systems of Conservation Laws, Applied

Mathematical Sciences 118, Springer-Verlag, 1996.[10] Godunov S.K., A difference method for numerical calculation of discontinuous equations of hydrodynamics,

Math. Sb. 47 (1959) 217–300.[11] Mohammadi B., Pironneau O., Analysis of theK-Epsilon Turbulence Model, Masson, Paris, 1994.[12] Rusanov V., Calculation of interaction of non-steady shock waves with obstacles, J. of Comput. Math. and Phys.

USSR 1 (1961) 267–279.[13] Smoller J., Shock Waves and Reaction–Diffusion Equations, 2nd ed., Springer-Verlag, 1994.[14] Xeuxet-Declercq E., Comparaison de solveurs numériques pour le traitement de la turbulence bi-fluide, Thèse de

doctorat, Université d’Evry, 1999.

1016