On présente une classe de  modèles

            d'écoulements diphasiques compressibles pour les simulations instationnaires. L'accent

            sera mis sur les principes de modélisation retenus pour la fermeture des grandeurs

            associées à l'interface, notamment les termes de transfert interfacial. On montre que le

            système fermé est régi par une inégalité d'entropie physique, et qu'il est hyperbolique

            au sens large. L'impact d'une fermeture particulière sera discuté, ainsi que le respect

            des conditions de positivité

HURISSE, Olivier

HÉRARD, Jean-Marc

French

I-Revues

19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
Une classe de mod̀eles diphasiques compressibles
J.M. H ÉRARD , O. HURISSE
EDF, R&D, Mécanique des Fluides, Energie et Environnement, 6 quai Watier 78400 CHATOU, FRANCE
Résuḿe :
Cet article est d́edíe à la mod́elisation deśecoulements diphasiques en approche bifluide. Les principaux éléments sont
rappeĺes, notamment les lois de fermeture pour les termes de transfert interfacial. Nous donnons ensuite l’algorithme
d’ensemble, consistant avec l’inégalit́e d’entropie, et nous précisons les techniques numériques utiliśees pour l’́etape de
relaxation en pression. Des résultats nuḿeriques compl̀etent l’ensemble de la présentation.
Abstract :
We focus in this paper on the modelling of two-phase flows whenusing the two-fluid approach. The basic ingredients are
first recalled, including the closure laws for the interfacil transfer terms. We then provide the whole algorithm that is
consistant with the entropy inequality, and examine some ways to compute the pressure relaxation step. Numerical results
complete the whole.
Mots clefs : Ecoulements diphasiques, ińegalité d’entropie, mod̀ele bifluide, relaxation en pression
1 Introduction
On examine dans cette note une classe de modèles bifluides d́edíesà la simulation deśecoulements diphasiques
compressibles. Depuis quelques annés, plusieurs auteurs s’intéressent̀a proposer un formalisme adéquatà cet
effet. On pourra notamment, sansêtre exhaustif, citer les travaux de Ransom et al ([1]), Baeret Nunziatto
([2]), Kapila et al([3]), Glimm et al ([4]), mais aussi les travaux de mod́elisation de Gavrilyuk et Saurel ([5]),
Papin et Abgrall ([6]), et Lhuillier ([7, 8]). Plus récemment, quelques travaux ayant pour objectif d’étendre
le formalisme sous-jacent aux modèles d’́ecoulement eau-vapeur ont permis de revisiter cette approche, en
consid́erant un ensemble de contraintesà satisfairea priori. Les travaux [9, 10, 11, 12] s’inscrivent dans
cette d́emarche. Cette dernière approche áegalement conduit̀a proposer des modèles de repŕesentation des
écoulements triphasiques ou multi-champs ([13, 14]), maisaussi pour la description des milieux poreux ([15,
16]), ou des milieux granulaires denses ([17]). Le principe d’ensemble est simple, et repose sur l’imposition
d’un cahier des charges. Celui-ci nécessite que le système soit tel que :
– les effets de convection soient de nature hyperbolique,
– le syst̀eme soit ŕegi par une ińegalit́e d’entropie physiquement admissible,
– les relations de saut dans la limite non visqueuse soient définies de manière unique,
– les conditions de positivité des taux de présence, densités et autres grandeurs thermodynamiquesad hoc
soient satisfaites au niveau continu.
On pŕesente tout d’abord le système d’́equations ouvert, les relations de fermeture, ainsi que lesprincipales
propríet́es du mod̀ele. On pŕesente ensuite les grandes lignes de la méthode nuḿerique, et quelques résultats.
2 Equations d’évolution du modèleà deux fluides
Les variables sont les suivantes :αk ∈ [0, 1] est la fraction volumique de phasek, sachant que :α1 + α2 = 1,
Uk, Pk, ρk, etmk = αkρk désignent respectivement la vitesse, la pression, la densité et la masse partielle de
phasek. Enfin, l’énergie totale de phasek est :Ek = ρkU2k/2 + ρkek(Pk, ρk). Le syst̀eme d’́equations ouvert,
portant sur les sept inconnuesα1 et (ρk, Uk, Pk) est alors le suivant, :



















∂α1
∂t
+ VI
∂α1
∂x
= φ1(W ) ;
∂mk
∂t
+
∂mkUk
∂x
= 0 ;
∂mkUk
∂t
+
∂mkU
2
k
∂x
+ αk
∂Pk
∂x
+ (Pk − PI)
∂αk
∂x
= Ik(W ) ;
∂αkEk
∂t
+
∂αkUk(Ek + Pk)
∂x
+ PI
∂αk
∂t
= VIIk(W ) .
(1)
1
19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
Les termes de transfert interfacialφ1(W ) et I1(W ) doivent faire l’objet de lois de fermeture, ainsi que la
vitesseVI et la pressionPI , sachant queI1(W ) + I2(W ) = 0.
2.1 Lois de fermeture
Afin d’obtenir une ińegalit́e d’entropie physiquement admissible, il est nécessaire que les termes de trans-
fert interfacial soient choisis judicieusement. La forme suivante ”classique” pour le transfert de quantité de
mouvement est en particulier admissible :
Ik = (−1)
k m1m2
τU (m1 + m2)
(U1 − U2) (2)
de m̂eme que le terme suivant afférentà l’équation gouvernant la fraction volumique :
φ1 = α1α2(P1 − P2)/(τP Π0) (3)
Leséchelles de tempsτP , τU sont fournies par la litt́erature, le termeΠ0 désignant une constante ayant la di-
mension d’une pression.
L’in égalit́e d’entropie fournit́egalement une information essentielle. En effet, si l’on postule la forme consis-
tante :
VI = βU1 + (1 − β)U2 (4)
où la fonction scalaireβ resteà d́eterminer, il s’ensuit que la contrainte suivante doitêtre v́erifiée :
PI = µP1 + (1 − µ)P2 (5)
où µ = a1(1−β)
a1(1−β)+a2β
, avecak = 1/Tk, Tk désignant la temṕerature de phasek.
Le troisìeme point du cahier des charges se traduit par une contraintesur la nature du champ associé à une
des valeurs propres du système convectif associé à (1), qui estλ = VI . Il importe effectivement que celui-ci
soit linéairement d́eǵeńeŕe, afin de garantir que les relations de saut par champ isolé soient d́efinies de manière
unique. Les fermetures du type (4) satisfaisant ce crit` re, qui peuvent̂etre exhib́ees ([9, 11]), correspondent de
fait à : β = 0, β = 1, β = m1/(m1 + m2). Ces ŕesultats sont intéressants̀a plusieurs titres. On remarque
tout d’abord que le casβ(1 − β) = 0 correspond au modèle utiliśe par Baer et Nunziatto ([2]), Kapila et al
([3]), Glimm et al ([5]) et d’autres auteurs, la formeétant un peu distincte dans [5]. La proposition hybride
β = m1/(m1 + m2), examińee dans [11], a un comportement dynamique. Par contre, asymptotiquement,
elle conduità privilégier pourVI la phase la plus présente, ce qui peut surprendrea priori. Nous pŕesentons
maintenant les principales propriét́es du syst̀eme ferḿe.
2.2 Principales propriétés
Propriét́e 1 :
Le mod̀ele (1), muni des fermetures(VI , PI) préćedentes, est hyperbolique sans condition sur l’espace des
états physiques, dès lors que les vitesses relatives ne sont pas soniques, i.e.|U1−U2| 6= ck. Les valeurs propres
sont :λ1 = VI , λ2 = U1, λ3 = U2, λ4,5 = U1 ± c1, λ6,7 = U2 ± c2.
Propriét́e 2 :
Les solutions ŕegulìeres du mod̀ele (1) muni des fermetures (2), (3), (4), (5), satisfont l’inégalit́e d’entropie
suivante :
∂
∑
mksk
∂t
+
∂
∑
mkuksk
∂x
≥ 0. (6)
où sk désigne l’entropie phasique.
Propriét́e 3 :
Les ondes associées aux valeurs propresλ1, λ2, λ3 sont lińeairement d́eǵeńeŕees. Dans chaque onde simple,
les relations de saut sont uniques.
Remarque 1 : L’ensemble de ces remarques sera utilisé pour la recherche d’approximations numériques
convergentes. Il faut noter que la dernière propríet́e est essentielle, et n’a pas d’équivalent dans le cadre
des mod̀eles bifluides classiquesà relaxation de pression locale instantanée. Sans cette dernière propríet́e, des
sch́emas d’approximation distincts vont converger vers des solutions distinctes.
2
19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
3 Algorithme global
3.1 Le sch́ema temporel
Nous ne d́etaillons pas ici l’ensemble de l’algorithme, et renvoyonse partieà la littérature aff́erente en ce
domaine. Une techniquèa pas fractionnaires est utilisée, qui consistèa traiter successivement :
– uneétape d’́evolution du syst̀eme convectif hyperbolique associé à (1)
– l’ensemble des termes sources de relaxation en pression eten vitessse associésàφ1 et Ik.
Le sch́ema de Rusanov et le schéma de Godunov approché ([18, 11]), mis en oeuvre pour obtenir des approxi-
mations stables et consistantes des solutions du système convectif (́etape d’́evolution), sont d́ecrits dans [11].
D’autres sch́emas peuvent̂etre consid́eŕesà ce stade, en particulier des schémas de relaxation tels que celui
propośe dans [19], ou d’autres solveurs de Riemann ([20]).
Au sein de la secondétape, la prise en compte des effets de trainée statique ne pose pas de difficulté parti-
culière. Par contre, l’étape de relaxation en pression n’est pas triviale, celle-ci couplant de manière non lińeaire
trois EDO distinctes.
3.2 Deux sch́emas de simulation de l’́etape de relaxation en pression
L’ étape isoĺee de relaxation en pression contribueà augmenter le bilan entropique global. Elle s’´ crit :
τP
∂α1
∂t
= α1α2(P1 − P2)/Π0
∂mkek
∂t
+ PI
∂αk
∂t
= 0
∂mk
∂t
=
∂mkUk
∂t
= 0
(7)
Sachant quePI = P1, il est clair que l’́ecart de pression(P1 − P2) suit au sein de l’́etape (7) l’́evolution :
1
2
∂(P1 − P2)
2
∂t
= −(α2ρ1c
2
1 + α1(γ2 − 1))P1 + α1P2)(P1 − P2)
2/τP /Π0 ≤ 0 (8)
dès lors que la phase vapeur2 suit une loi de gaz parfait :ρ2e2 = (γ2 − 1)P2. Ce ŕesultat est totalement
conformeà celui de [11]. Nous d́ecrivons ci-dessous deux algorithmes distincts d’intégration du système (7).
Ces deux algorithmes sont tels que :
– Le principe du maximum pour le taux de présenceα1 est assuŕe,
– On respecte le principe de conservation local de l’énergie totale.
3.2.1 Algorithme semi-discret
L’algorithme semi-discret́evalue, pour une donnéeW̃ , l’ état impliciteW suivant :
∆P = P1 − P2 = (P̃1 − P̃2)e
−C̃∆t
α1/(1 − α1) = α̃1/(1 − α̃1)e
(P1−P2)∆t
Π0τP
(9)
sachant que :C = (α2ρ1c21 + α1ρ2c
2
2 + α1(P1 −P2)(ρ2
∂e2
∂P2
)−1)/(τP Π0). La pressionP1 est ensuite calculée
en ŕesolvant l’́equation scalaire implicite non linéaire suivante :
H(P1) = m̃1e1(P1,
m̃1
α1
) + m̃2e2(P1 − ∆P,
m̃2
1 − α1
) = m̃1ẽ1 + m̃2ẽ2 (10)
Pour des thermodynamiques classiques, la fonctionH(P ) est croissante, garantissant l’unicité de la solution ;
l’existenceétant assuŕee pour des pas de temps∆t suffisamment grands. Le principe du maximum discret est
clairement assuré pourα1 (voir (9)), et l’algorithme est encore globalement conservatif enénergie totale, par
construction.
3
19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
3.2.2 Algorithme discret
L’algorithme discret calcule, pour une donnéeW̃ , l’ état impliciteW satisfaisant :
τP (α2 − α̃2) + ∆tα1α2(P1 − P2)/Π0 = 0
m1e1(P1,
m1
1−α2
) + P1(α̃2 − α2) = m̃1ẽ1
m2e2(P2,
m2
α2
) − P1(α̃2 − α2) = m̃2ẽ2
(11)
sachant que :mk = m̃k, et :Uk = Ũk. Consid́erant pour simplifier la présentation deux lois phasiques du type
suivant :ρkek = (γk − 1)Pk, on peut alors montrer que le système discret (11) admet une unique solutionα2
dans l’intervalle[0, αmax], sachant queαmax =
1+(γ1−1)α̃2
γ1
, quelque soit la valeur de∆t. Les pressionsP1, P2
sont alors positives. L’algorithme est ici encore globalement conservatif eńenergie totale. Ce résultat est en
parfait accord avec celuiénonće dans [17].
Remarque 2 : On noteráegalement que lorsqueτP = 0, l’algorithme discret s’identifie avec l’algorithme de
relaxation instantańee d́ecrit dans [21]. On rappelle cependant que l’utilisation d’une telle valeur conduirà
une divergence du calcul, lorsque le maillage est raffiné, i d́ependamment du cas test considéré.
4 Résultats nuḿeriques
On pŕesente ci-dessous deux résultats simples. Le premier correspondà l’évaluation de la précision obtenue
lors de l’approximation des solutions de l’étape de relaxation en pression, le secondà la simulation d’un
probl̀eme de Riemann diphasique.
4.1 Un cas test analytique
Ce cas test analytique est extrait de [22]. L’échelle de temps considéŕee pour ce test est :τP = 10−6. Le temps
d’intégration estT = 10−5, etΠ0 = 21 × 105. On montre l’erreur de convergence en normeL1 obtenue sur
un cas analytique avec les deux algorithmes d´ crits ci-dessus (voir les figures (1,2)).
-3 -2 -1 0 1
5
10
FIG. 1 – Convergence des deux schémas de relaxation en pression : algorithme discret (cercles) et semi-discret
(carŕes), pourPl etPv (pointillés), avecτP = 10−6. On trace iciLog(||erreur||1) en fonction deLog( δtT ).
4.2 Un problème de Riemann
On consid̀ere maintenant un problème de Riemann, pour lequel les conditions initiales sont les suivantes :
W (x < 500, t = 0) = WL et :W (x > 500, t = 0) = WR, avec :
αl ρl Ul Pl ρv Uv Pv
Etat L 0.95 800 10 150 · 105 80 15 145 · 105
Etat R 0.90 700 10 145 · 105 70 15 140 · 105
4
19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
-3 -2 -1 0 1
-8
-6
-4
-2
0
FIG. 2 – Convergence des deux schémas de relaxation en pression : algorithme discret (cercles) et semi-discret
(carŕes), pour la fraction volumiqueαl, avecτP = 10−6.
La figure3 illustre à l’instantT = 0.012 le comportement des vitesses et pressions de phase, sachantque
τU = τP = 10
−3, la simulation ayant́et́e ŕealiśee sur un maillage comportant1000 mailles.
0 200 400 600 800 1000
10
11
12
13
14
0 200 400 600 800 1000
1,44e+07
1,46e+07
1,48e+07
1,5e+07
FIG. 3 – Profils de vitesse (à gauche) et de pression (à droite) pour la simulation du cas2, pour un maillage
comportant1000 cellules, avec :CFL = 0.5, sachant queτU = τP = 10−3.
5 Conclusion
La prise en compte des termes de transfert de masse et de chaleur interfaciaux rentre parfaitement dans le
cadre du pŕesent formalisme ([13]). La classe de modèles consid́eŕes permet en outre d’envisager le couplage
interfacial instationnaire avec d’autres modèles diphasiques, notamment avec les modèles homog̀enes. Les po-
tentialit́es d’utilisation sont assez vastes, et le contrˆ le entropique est un outil puissant pour la modélisation,
mais aussi pour la recherche d’approximations stables et consistantes. Signalons enfin qu’une possibilité d’ex-
tension dynamique de l’approche d´ crite dans cette note est faite dans [23], l’objectifétant cette fois de pouvoir
gérer la transition d’́ecoulement. La formulation pourrait alors devenir une alternativeà la formulation multi-
champs utiliśee dans le contexte nucléaire civil.
Remerciements
Ce travail áet́e ŕealiśe dans le cadre du projet NEPTUNE, qui est financé par le CEA (Commissariat` l’Énergie
Atomique), EDF, l’IRSN (Institut de Radioprotection et de Sûret́e Nucĺeaire) et AREVA-NP.
5
19 èmeCongr̀es Français de Ḿecanique Marseille, 24-28 aôut 2009
Références
[1] Ransom V. and Hicks D. Hyperbolic two-pressure model fortwo-phase flow. J. Comp. Phys., 53, 124–
151, 1984.
[2] Baer M. and Nunziato J. A two-phase mixture theory for thedeflagration to detonation transition (DDT)
in reactive granular materials. Int. J. Multiphase Flow, 12(6), 861–889, 1986.
[3] Kapila A., Son S., Bdzil J., Menikoff R., and Stewart D. Two-phase modeling of a DDT : structure of the
velocity relaxation zone. Phys. of Fluids, 9(12), 3885–3897,1997.
[4] Glimm J., Saltz D., and Sharp D. Two-phase flow modelling of a fluid mixing layer. J. Fluid Mech., 378,
119–143, 1999.
[5] Gavrilyuk S. and Saurel R. Mathematical and numerical modelling of two-phase compressible flows with
micro-inertia. J. Comp. Phys., 175, 326–360, 2002.
[6] Papin M. and Abgrall R. Fermetures entropiques pour les syst̀emes bifluides̀a septéquations. C. R.
Mécanique, 333, 838–842, 2005.
[7] Lhuillier D. A mean field description of two-phase flows with pase changes. Int. J. of Multiphase Flow,
29, 511–525, 2003.
[8] Lhuillier D. Evolution of the volumetric interfacial areain two-phase mixtures. C. R. Ḿecanique, 332,
103–108, 2004.
[9] Coquel F., Galloüet T., H́erard J., and Seguin N. Closure laws for a two-fluid two-pressure model. C. R.
Acad. Sci. Paris, I-332, 927–932, 2002.
[10] Hérard J. Numerical modelling of turbulent two phase flows using the two-fluid approach. AIAA paper
2003-4113, http ://www.aiaa.org/, 2003.
[11] Galloüet T., H́erard J., and Seguin N. Numerical modelling of two phase flows using the two-fluid two-
pressure approach. Math. Mod. Meth. in Appl. Sci., 14(5), 663–700, 2004.
[12] Hérard J. and Hurisse O. Some recent numerical advances for two-phase flow modelling in NEPTUNE
project. Int. J. Num. Meth. Fluids, 59, 285–307, 2009.
[13] Hérard J. A three-phase flow model. Mathematical and Computer Modelling, 45, 732–755, 2007.
[14] Hérard J. A relaxation scheme to compute three-phase flow models. AIAA paper 2007-4455,
http ://www.aiaa.org/, 2007.
[15] Hérard J. Un mod̀ele hyperbolique diphasique bifluide en milieu poreux. C. R. Mécanique, 336, 650–655,
2008.
[16] Girault L. and H́erard J. A two-fluid hyperbolic model in a porous medium. soumis, 2008.
[17] Galloüet T., Helluy P., H́erard J., and Nussbaum J. Hyperbolic relaxation model for granular flows.
soumis, 2008.
[18] Buffard T., Gallouet T., and H́erard J. A sequel to a rough Godunov scheme. Application to real gases.
Computers and Fluids, 29(7), 813–847, 2000.
[19] Ambroso A., Chalons C., Coquel F., and Galié T. Relaxation and numerical approximation of a two-fluid
two-pressure diphasic model. Math. Mod. Numer. Anal., 2008.
[20] Schwendeman D., Wahle C., and Kapila A. The Riemann problem and high resolution Godunov method
for a model of compressible two-phase flow. J. Comp. Phys., 212,490–526, 2006.
[21] Hérard J. and Hurisse O. A simple method to compute two-fluid models. Int. J. of Comp. Fluid Dyn.,
19(7), 475–482, 2005.
[22] Hérard J. and Hurisse O. Schémas d’approximation pour une classe de modèles bifluides. Rapport interne
EDF, à paraitre.
[23] Hérard J. Une classe dynamique de modèles diphasiques bifluides. Rapport interne EDF,à paraitre.
6


Une classe de modèles diphasiques compressibles

http://documents.irevues.inist.fr/bitstream/2042/36725/1/658.pdf

Une classe de modèles diphasiques compressibles

Abstract

Similar works

Full text

Available Versions

I-Revues