FR9710029

Production d'énergie
(hydraulique, thermique
et nucléaire)

MODELISATION ET ANALYSE NUMERIQUE DES

ECOULEMENTS DIPHASIQUES EN DESEQUILIBRE

MODELING AND NUMERICAL ANALYSIS OF

NON-EQUILIBRIUM TWO-PHASE FLOWS

97NB00078

D
I : -i
«
in,-/:;,:',», 'g
 I
FR9710029

DIRECTION DES ÉTUDES ET

RECHERCHES
SERVICE RÉACTEURS NUCLÉAIRES ET ECHANGEURS
DÉPARTEMENT PHYSIQUE DES RÉACTEURS

Avril 1997

RASCLE P.
EL AMINE K.

MODELISATION ET ANALYSE NUMERIQUE DES

ECOULEMENTS DIPHASIQUES EN
DESEQUILIBRE

MODELING AND NUMERICAL ANALYSIS OF

NON-EQUILIBRIUM TWO-PHASE FLOWS

Pages : 168 97NB00078

Diffusion : J.-M. Lecœuvre

EDF-DER
Service IPN. Département PROVAL
29- 08 ©EDF1997
1, avenue du Général-de-Gaulle
92141 Clamait Cedex ISSN 1161-0611
SYNTHESE :

On s'intéresse à l'étude numérique de modèles "bi-fluide" d'écoulements

diphasiques en déséquilibre décrits par un système à six équations. On introduit une
technique originale de décomposition du système. Cette technique permet l'utilisation
de solveurs de Riemann développés pour les écoulements monophasiques ; de plus,
elle permet une extension directe à d'autres modèles plus évolués contenant des termes
d'échanges entre phases plus précis. Les propriétés des fluides sont approchées par des
lois d'état de type gaz parfait polytropique puis étendues au cas des fluides réels.
Pour la construction des schémas numériques, l'hyperbolicité du système global
n'est pas nécessaire. En utilisant des schémas cinétiques décentrés appropriés de type
volumes finis, on montre que la méthode est capable de traiter différents régimes
d'écoulements pouvant évoluer d'un écoulement diphasique à un écoulement
monophasique et vice versa, tout en préservant les caractéristiques des grandeurs
physiques qui restent dans l'ensemble des états admissibles. Plusieurs tests numériques
modélisant des problèmes physiques sont présentés pour illustrer l'efficacité de la
méthode proposée.

97NB0W78

(HT-13/97/013/A)
EXECUTIVE SUMMARY :

We are interested in the numerical approximation of two-fluid models of

nonequilibrium two-phase flows described by six balance equations. We introduce an
original splitting technique of the system of equations. This technique is derived in a
way such that single phase Riemann solvers may be used ; moreover, it allows a
straightforward extension to various and detailed exchange source terms. The
properties of the fluids are first approached by state equations of ideal gas type and
then extended to real fluids.

For the construction of numerical schemes, the hyperbolicity of the full system
is not necessary. When based on suitable kinetic upwind schemes, the algorithm can
compute flow regimes evolving from mixture to single phase flows and vice versa. The
whole scheme preserves the physical features of all the variables which remain in the
set of physical states. Several stiff numerical tests, such that phase separation and
phase transition are displayed in order to highlight the efficiency of the proposed
method.

97NB00078

(HT-13/97/013/A)
Table des matières

1 INTRODUCTION 9

2 ECOULEMENT DIPHASIQUE, MODELISATION ET HYPERBOL-

ICITE 17

2.1 Modélisation et lois constitutives 19

2.1.1 Introduction 19

2.1.2 Lois de conservation pour le modèle bi-fluide 20

2.1.3 Formes de lois de transfert et déséquilibres 26

2.1.4 Hypothèses essentielles retenues 33

2.2 Hyperbolicité et caractère bien-posé 34

2.2.1 Caractère "bien-posé" des équations à caractéristiques complexes . . 34

2.2.2 Hyperbolicité de quelques modèles 39

3 UNE METHODE NUMERIQUE DECENTREE POUR LA RESOLU-

TION D'ECOULEMENTS DIPHASIQUES SANS TERME D'ECHANGE 57

3.1 Introduction 59

3.2 Two-fluid model of two-phase flow 60

3.3 Numerical Method 63

3.3.1 A two-step approach 63

 3.3.2 Discretization of the first step 64

3.3.3 Discretization of the second step 69

3.4 Properties of the method when based on kinetic schemes 71

3.5 Numerical Results 75

4 UN SCHEMA NUMERIQUE POUR LES ECOULEMENTS MONO-

PHASIQUES DE FLUIDES REELS 89

4.1 Introduction 91

4.2 Governing equations 92

4.2.1 Dynamical equations 92

4.2.2 Equation of state 92

4.3 Kinetic scheme 96

4.3.1 Kinetic formalism 97

4.3.2 Construction of the numerical scheme 100

4.4 Computational results 104

5 UNE METHODE NUMERIQUE DECENTREE POUR LES ECOULE-

MENTS DIPHASIQUES EN DESEQUILIBRE 115

5.1 Introduction and description 117

5.2 Two-fluid nonequilibrium model of two-phase flows 118

5.2.1 The physical model 118

5.3 Numerical method 121

5.3.1 Discretization of the first step 123

5.3.2 Discretization of the second step 126

5.4 Computational results 131

6 LE TRAITEMENT DES CONDITIONS AUX LIMITES 141

6.1 Introduction * 143

6.2 Formulation of the boundary conditions for the two-fluid model 144

6.3 Perfect fluid approach 146

6.3.1 Inflow boundary (x = 0) 149

6.3.2 Outflow boundary (x = L) 151

6.4 Real fluid approach 152

6.4.1 Inflow boundary (x = 0) 154

6.4.2 Outflow boundary (x = L) 156

A The Perthame and Roe schemes 161

A.I The Perthame scheme 162

A.2 The Roe scheme 164

 THESE DE DOCTORAT DE L'UNIVERSITE PARIS 6

Spécialité : Mathématiques Appliquées

présentée par

Khalid EL AMINE

pour obtenir le grade de DOCTEUR DE L'UNIVERSITE PARIS 6

Sujet de la thèse :

Modélisation et analyse numérique des écoulements diphasiques

en déséquilibre

Soutenue le 21 mars 1997 devant le jury composé de :

M. P.-A. RAVIART Rapporteur

M. J.-P. VILA Rapporteur
M. F. COQUEL
M. Ph. EMONOT
Mme E. GODLEWSKI
M. B. PERTHAME Directeur
M. P. RASCLE
Chapitre 1

INTRODUCTION

®ggXY PAGE(S)
!©f t BLABSK
 INTRODUCTION

Cette thèse a pour objet l'étude, d'un point de vue numérique, de modèles bi-fluide
d'écoulement diphasique liquide-vapeur. De tels écoulements gouvernent le fonctionnement
des échangeurs (coeur de réacteur, générateur de vapeur et condenseur) d'une centrale
nucléaire à eau pressurisée. Les méthodes développées peuvent être étendues à d'autres
types d'applications.

Dans un coeur de réacteur à eau pressurisée (REP) en fonctionnement normal, le

fluide caloporteur (eau) qui sert à récupérer la chaleur engendrée par la réaction nucléaire
au sein du combustible (Uranium, Plutonium) reste à l'état liquide. Cependant, en cas
d'incident (par exemple : une élévation de puissance ou une dépressurisation brutale due à
une brèche), une vaporisation localisée du liquide pourrait survenir. La présence de vapeur
peut influencer considérablement les échanges thermiques et la source de puissance. En
effet, l'efficacité du refroidissement change du fait que la conductivité thermique varie
fortement en présence de vapeur. De même, plus le taux de vapeur est important dans
un volume donné, moins il y a de molécules d'eau. Comme l'eau ralentit les neutrons et
que la réaction nucléaire qui fournit l'énergie thermique est entretenue par les neutrons
lents, le nombre de fissions nucléaires diminue et par conséquent, la production de chaleur
aussi. Les études de sûreté pour les coeurs de réacteur imposent une estimation précise
de la température du combustible, ce qui nécessite une bonne connaissance des échanges
thermiques et donc de l'écoulement.
Dans un générateur de vapeur (GV), l'énergie thermique est apportée par l'eau du circuit
primaire du réacteur. Cette eau passe à travers des faisceaux de tubes constituant le GV.
Au contact des tubes chauds, l'eau du circuit secondaire s'évapore par ebullition nucléée
et un phénomène de convection naturelle met le mélange diphasique en circulation dans le
GV. Là aussi, une bonne connaissance de l'écoulement est nécessaire à l'estimation précise
des échanges thermiques et donc du rendement du GV.

Or, la qualité des échanges de chaleur entre les fluides et les solides dépend étroite-
ment des caractéristiques et de la configuration de l'écoulement. A titre d'exemple, la
présence d'un film de vapeur à la paroi du solide dégrade fortement les échanges ther-
miques. Les variations de pression, l'énergie apportée au fluide, provoquent des change-
ments de phase (vaporisation, condensation). Pendant le fonctionnement normal ou ac-
cidentel des échangeurs, divers types d'écoulements peuvent apparaître : écoulement à
bulles, écoulement à bouchons, écoulement laminaires, écoulement à gouttelettes...

On voit donc toute l'importance d'un modèle d'écoulement qui doit être suffisament
représentatif du point de vue physique, de façon à décrire les écoulements réels apparaissant

11
dans les différents types d'échangeurs, pour les différents régimes de fonctionnement.

Un fluide diphasique est constitué par des régions monophasiques limitées par des in-
terfaces en mouvement. Localement, le comportement de chaque phase est décrit complète-
ment par des équations aux dérivées partielles traduisant les bilans de conservation de la
masse, de la quantité de mouvement et de l'énergie, et par les conditions d'interface as-
sociées. Ce type de formulation est une extension directe de la formulation des écoulements
monophasiques et s'appelle : formulation locale instantanée.
Une telle formulation nous amènerait à un problème de frontières multiples ayant des po-
sitions inconnues. La résolution de problèmes diphasiques utilisant une telle formulation,
soulève trop de difficultés mathématiques, qui restent dans la plupart des cas insurmonta-
bles. De plus, pour les études liées à la sûreté des installations et à leur dimensionnement,
une description fine des phénomènes à petite échelle n'a en général pas d'importance. Par
conséquent, la résolution de problèmes diphasiques basés sur cette formulation n'est pas
envisagée. En pratique, on moyenne ces équations locales dans le temps et dans l'espace afin
d'aboutir à un système utilisable pratiquement et obtenir l'échelle qui nous intéresse. Le
modèle obtenu est encore un système de lois de conservation et comporte trois équations
de conservation pour chaque phase et s'appelle : modèle bi-fluide. Ce modèle est donc
formulé en considérant chaque phase séparément. Ainsi, il comprend des équations de
conservation de la masse, de la quantité de mouvement et de l'énergie de chaque phase,
avec des termes d'interaction qui apparaissent dans ces équations et qui sont des termes
d'échanges à l'interface et à la paroi. Le système que nous étudions dans cette thèse s'écrit

dx{akpkuk) = Fk,

dt{akpkuk) + dx(akpkul + a;i>) = Pidxak + akpkg + Mk + Mwk, ^

2 2
dt{akpk{ek + ~ ) ) + dx{akpk{ek + h -j-)^) + pdtak = akpkguk + Ek + Ewk,

avec (A: = v,l; v = vapeur, / = liquide). ak, pk, i , ek sont respectivement le taux de
présence local, la masse volumique, la vitesse et l'énergie spécifique interne de la phase
k. rk, Mk et Ek, sont les termes de transferts interfaciaux représentant respectivement la
masse, la quantité de mouvement et l'énergie acquise par la phase k. Mwk et Ewk sont
les termes de transfert entre la phase k et la paroi, respectivement de la quantité de
mouvement et de l'énergie, p est la pression supposée la même pour les deux phases, g
est l'accélération de la pesanteur, p, est la pression à l'interface. Pour fermer le système
(1.1), nous ajoutons les lois d'états constitutives du type pk = pk(pk, hk), ek = ek(pk,pk),
la relation av + ai — 1 et nous exprimons tous les termes d'échanges et la pression à
l'interface en fonction des quantités macroscopiques introduites ci-dessus.

Le modèle bi-fluide est considéré comme le plus approprié pour une description
générale et détaillée des transitoires. Cependant, le nombre de lois constitutives qu'il est
nécessaire et suffisant de spécifier pour fermer le système est relativement élevé par rapport
aux autres modèles simplifiés tels que le modèle homogène équilibré. En l'état actuel des
choses, les travaux expérimentaux et les recherches ont fourni plus d'informations sur
les lois de fermeture pour les modèles homogènes que pour les modèles bi-fluides. La

12
fermeture des modèles bi-fluide reste donc un problème ouvert, et il existe de nombreuses
formulations des différents termes d'échange, pas toujours équivalentes. Les systèmes bi-
fluide utilisés en pratique sont souvent non hyperboliques, et, au mieux, ils ne sont que
conditionnellement hyperboliques. Pendant le fonctionnement normal ou accidentel d'un
générateur, ces conditions peuvent être violées.

Dans la première partie du chapitre 2, on fait un rappel sur la modélisation physique,

et on donne une forme générale des principaux termes d'échange ainsi que leur impor-
tance sur les phénomènes mis en jeu dans les écoulements. Dans la seconde partie du
chapitre 2, on donne tout d'abord quelques résultats concernant le lien entre la non-
hyperbolicité de certains modèles et la stabilité des schémas numériques associés. On
passe ensuite, à l'étude de Phyperbolicité d'un modèle bi-fluide à six équations incluant
un terme différentiel représentant la "masse virtuelle", puis à l'analyse de l'hyperbolicité
d'un modèle bi-fluide à quatre équations dont la "correction" porte sur l'expression de la
pression à l'interface.
Dans le cas où la pression à l'interface p,- est supposée égale à la pression p et les termes
d'échange à l'interface sont modélisés par des formules algébriques, le système (1.1) est
non hyperbolique. L'expression de p, est généralement donnée par p,; = p — (p — pi) =
p — Çpc\uv — ui\2 où f = constante et pc est la densité de la phase continue. En général,
les équations de conservation de l'énergie peuvent être découplées des autres équations de
conservation, et l'étude de l'hyperbolicité revient à celle du système à quatre équations.
En supposant le liquide incompressible (hypothèse raisonnable en général, à moins que
a v soit très petit) et que £ = 1, pc = pv, le système à quatre équations est hyperbolique
sous la condition \uv — ui\ < cv, où cv est la vitesse du son dans la vapeur. En pratique,
cette condition est largement satisfaite. Puisque le liquide est faiblement compressible,
l'expression {p — Pi) = pv\uv - u/|2 sera retenue par la suite pour la résolution du système
(1.1).

Vu les difficultés liées à la modélisation des terme d'échange et à la nature "incer-

taine" du système (1.1), nous introduisons une méthode numérique ne nécessitant pas
l'hyperbolicité du système (1.1) et permettant une extension directe à la résolution de
modèles plus évolués contenant des termes d'échange entre phases plus précis, sans changer
la méthode de base.

Dans le chapitre 3, nous considérons le système (1.1) sans terme d'échange (par ex-
emple : écoulement diphasique eau-air à basse pression et sans frottement entre phases).
Pour fermer le système, on suppose que les deux phases sont régies par des lois d'état de
type gaz parfait polytropique, i.e. p— (7* — l)PkCk, avec fv = 1.4 et 7/ = 1.0005. Cette ap-
proche nous a permis de traiter des écoulements diphasiques avec présence de phase liquide
pure localement, ce qui n'aurait pas été possible si le liquide avait été supposé incompres-
sible. Nous traitons aussi le cas où le système est non hyperbolique (p — pi) et le cas
"hyperbolique" (p ^ p,-). Pour résoudre notre système, nous introduisons une méthode de
décomposition qui réduit la résolution du système (1.1) à celle de deux systèmes découplés
ressemblant au problème classique de la tuyère et à la résolution d'un système d'équations

13
différentielles ordinaires. Ainsi, dans la première étape on considère :

dt(mk) + dx(mkuk) = 0,

dt(TnkUk) + dx(mkul + pk) = Pidxak+akpkÇ, k = v,l (1.2)

dt(Ek) + dx((Ek + p~k)uk) = QkPk9uk,

u2
avec mk = akpk, Ek = akpk{ek + ~) et p~k = akpk = (7* - l)rojte*.
Notons que ces deux systèmes évoluent indépendamment. Par conséquent, les pressions
phasiques obtenues à l'issue de ce pas peuvent être différentes.
Dans la deuxième étape, on impose l'égalité des pressions et on résout
dtmk = 0,

dt{mkuk) = 0, k = v,l (1.3)

dtEk = -pdtotk-
En ignorant les indices fc, le schéma de volumes finis associé au système (1.2) s'écrit

F?_l/t) = AtH?. (1.4)

pour calculer les flux F^,1i2, nous avons utilisé un schéma de Perthame et le schéma
de Roe. Le terme H™ comprend les termes dûs à la force de gravité (qui sont pris ex-
plicitement) et la forme discrétisée de pidxa. Ce dernier terme joue un rôle très important
dans la stabilité du schéma global. Ainsi, nous transformons dxa en a ( l - a)d r /3, avec
/3 — In(j^) et nous discrétisons suivant la nature du système. Pour pi ^ p nous util-
isons (p,# r a)" = (p,)"a!>(l - <*j)db{^ji+i ~ P?-i} e t P o u r ' e c a s Pi — P n o u s utilisons

Pour des raisons de stabilité et de conservation de l'énergie totale, le système (1.3) est
discrétisé implicitement. Les étapes de calcul pour déterminer U"+1 sont de simples
opérations algébriques. Notons que le schéma global obtenu par cette méthode n'est pas
un schéma à pas fractionnaires. Le schéma de Roe utilisé pour calculer les flux numériques
n'est pas utilisé dans le cas où l'une des deux phases peut disparaître localement, à cause
de son caractère instable dans les zones d'apparition du vide (a = 0 ou 1, avec a = av). En
utilisant le schéma de Perthame, nous avons obtenu de très bons résultats. Le schéma global
est capable de décrire des situations où l'écoulement passe du diphasique au monophasique.
Aussi, dans le cas non hyperbolique et en combinaison avec le schéma de Perthame, le
schéma global préserve, sous une condition CFL raisonnable, la positivité des masses vo-
lumiques p^, de la pression p et le contrôle des fractions volumiques, i.e. 0 < a^ < 1, ce
qui est d'importance considérable pour l'interprétation physique des résultats. Ce travail,
en collaboration avec F. Coquel. E. Godlewski, B. Perthame et P. Rascle, est soumis au
J. Comp. Phys.

Nous nous intéressons dans le chapitre 4 à l'étude d'un écoulement monophasique

d'un fluide réel, ce qui peut être considéré comme un cas particulier d'un écoulement

14
diphasique de fluides réels (a = 0 ou 1). La dynamique d'un fluide réel est gouvernée par
un système de lois de conservation dont la forme est la même que celle de la dynamique
des gaz habituels. Seule la loi d'état change. Les propriétés thermodynamiques des fluides
considérés sont données par une loi d'état tabulée de la forme p = p(p,h). Comme la
densité et l'énergie interne apparaissent dans le vecteur des variables conservatives, nous
remplaçons p par p(h — e) et nous résolvons l'équation non linéaire p = p{p(h — e), h) sur
h par une méthode de Newton. Nous définissons ensuite 7 = h/e et nous approchons la
loi d'état par p = (7 — l)pe. Nous introduisons aussi le coefficient adiabatique 7 = (c2p/p)
où c est la vitesse du son physique donnée par la table thermodynamique. En utilisant
les quantités 7 et 7, nous développons un schéma cinétique avec fonctions d'équilibre à
support compact. Ce schéma est une extension du schéma de Kaniel développé pour la
dynamique des gaz parfaits isentropiques aux cas des fluides réels non isentropiques. Nous
donnons un résultat de stabilité L1 dans les régions physiques qui nous intéressent. La
robustesse du schéma est testée sur deux tubes à choc de gaz parfait polytropique et sur
un tube à choc pour l'eau pure. Ce schéma peut donc être appliqué dans l'étude numérique
des explosions sous-marines.

Le chapitre 5 est consacré à l'étude d'écoulements diphasiques de fluides réels en

déséquilibre. Les termes sources sont pris en compte pour modéliser les éventuels
déséquilibres mécaniques et thermodynamiques entre les phases mais aussi les termes
d'échange entre chacune des phases et la paroi. On considère alors le système (1.1) complet.
Les lois d'état sont données par une table thermodynamique. On a adopté des expressions
simples pour modéliser les termes d'échange.
On résout le système obtenu par la méthode de décomposition introduite au chapitre 3.
Dans la première étape, on ajoute les termes d'échange entre chaque phase et la paroi. Les
deux systèmes obtenus sont encore découplés. Dans ce chapitre, seul le cas p, ^ p est con-
sidéré, et la discrétisation du terme pidxa est tout simplement centrée. Cette discrétisation
est suffisante pour stabiliser les calculs. Ceci est probablement dû au terme de frottement
à l'interface (traînée) que nous avons pris en compte dans (1.1). Les flux numériques sont
calculés grâce au schéma cinétique introduit au chapitre 4.
Dans la deuxième étape, on prend en compte les termes d'échange à l'interface. Ces termes
sont raides et surtout ceux associés aux transitions de phase, ils sont donc discrétisés im-
plicitement. On résout ensuite le système non linéaire obtenu par une méthode de Newton-
Raphson.
Nous testons d'abord la méthode sur un tube à choc diphasique. En utilisant le schéma
de Lax-Friedrichs dans la première étape, nous constatons que le schéma cinétique est
beaucoup moins diffusif. Nous considérons ensuite un problème d'écoulement transitoire
très violent. Ce test est caractérisé par une dépressurisation brutale dans le système qui
se traduit par une vaporisation totale de l'eau. Les résultats obtenus sont satisfaisants
et la méthode numérique utilisée est en mesure de traiter des déséquilibres mécaniques
et thermodynamiques simultanés des deux phases, ainsi que le passage d'un écoulement
monophasique liquide à un écoulement monophasique vapeur par transition de phase.

Le chapitre 6 est consacré au calcul des conditions aux limites. Le traitement des
conditions aux limites est basé sur la méthode introduite pour résoudre le modèle bi-fluide
(chapitres 3 et 5). Ainsi, il ne concerne que la première étape, puisque la deuxième ne

15
fait intervenir que des dérivées temporelles. On analyse le cas où la loi d'état est du type
gaz parfait, mais aussi le cas de fluides réels. Pour une paroi solide, on utilise la notion
d'état miroir pour calculer le flux numérique au bord. Dans le cas d'une frontière fluide,
on adapte le formalisme de Dubois qui consiste à résoudre un demi problème de Riemann
pour calculer l'état au bord. On utilise ensuite cet état avec l'état interne bordant la
frontière pour déterminer le flux numérique.

Dans cette thèse, on n'utilise que des schémas numériques "explicites" d'ordre 1. Il
serait important d'impliciter la première étape pour accroître les performances du calcul,
de considérer des schémas d'ordre supérieur et d'étendre la méthode proposée à plusieurs
dimensions d'espace.

16
Chapitre 2

ECOULEMENT DIPHASIQUE,
MODELISATION ET
HYPERBOLICITE

17

Sef t
 ECOULEMENT DIPHASIQUE, MODELISATION ET
HYPERBOLICITE

2.1 Modélisation et lois constitutives

2.1.1 Introduction

Dans des conditions données, le comportement d'un système physique donné est en général
représenté par un modèle mathématique. Les modèles mathématiques d'écoulement sont
basés sur des axiomes idéalisant certaines lois universelles (en pratique, principes de con-
servation de la masse, de la quantité de mouvement et de l'énergie) et sur des lois de
comportement, idéalisant les propriétés des fluides considérés. Les lois de comportement
servent essentiellement à fermer le système d'équations constituant le modèle. En effet,
les équations exprimant les principes de conservation contiennent les propriétés des flui-
des considérés. Ces propriétés doivent être spécifiées et s'accompagner de préoccupations
physiques. Cependant, d'un point de vue pratique, on met en évidence d'autres con-
traintes : les lois de comportement doivent être telles que le système d'équations soit
effectivement fermé, qu'il ait une solution, que cette solution soit unique en général et
qu'elle ait un sens physique. De plus, pour une étude numérique du modèle, les lois de
comportement doivent permettre de poser correctement le problème au sens d'HADAMARD
(hyperbolicité).

Les lois de comportements (ou lois constitutives), qu'il est nécessaire et suffisant
de spécifier pour fermer le système des équations moyennées, traduisent, non pas le com-
portement intrinsèque des fluides considérés, mais leurs comportements dans le cadre du
système étudié (compte-tenu des différents fluides en présence, des conditions aux limites,
etc). Il est alors important de distinguer :

• Les lois donnant les propriétés physiques (moyennes), au sein de chaque milieu con-
tinu.
• Les lois de transfert, soit aux parois, soit entre milieux continus.
• La loi topologique de taux de vide : les opérations de moyenne, qu'on fait subir
au modèle sous forme locale instantanée pour obtenir des équations moyennées,
s'accompagnent d'une perte d'information. Dans le cas des écoulements diphasiques,

19
 cette perte d'information se traduit notamment par la présence, dans les équations de
conservation moyennées, du temps de présence locale de l'une des phases moyennée
spatialement, ou de la concentration spatiale de l'une des phases moyennée dans
le temps (ou statistiquement le cas échéant). La valeur commune de ces grandeurs
appelée "taux de vide" résulte de la structure réelle de l'écoulement (géométrie des
interfaces) ; cette valeur doit être restituée au niveau du modèle par une loi qualifiée
de topologique (voir BOURÉ [2]).

2.1.2 Lois de conservation pour le modèle bi-fluide

Le modèle bi-fluide est formulé en considérant chaque phase séparément. Ainsi, le modèle
comprend des équations de conservation de la masse, de la quantité de mouvement et de
l'énergie de chaque phase, avec des termes d'interaction ou d'échange (à l'interface ou à
la paroi) qui apparaissent dans ces équations.

Les modèles bi-fluide sont obtenus en effectuant des opérations de moyenne de la

formulation locale instantanée, ces opérations de moyenne peuvent être des moyennes
en temps (ISHII [10]), des moyennes en espace (DELHAYE [3]), des moyennes spatio-
temporelles (LAHEY et DREW [13]) ou des moyennes spatio-statistiques.

Cependant, la forme résultante des équations ne diffère pas sensiblement ; ce qui

change, c'est le sens physique des quantités moyennées et l'expression des termes d'échan-
ge à l'interface, cela est un point très important qu'il ne faut donc pas perdre de vue pour
ne pas mélanger les expressions développées par différents auteurs et qui n'utilisent pas
les mêmes moyennes. Dans cette note, nous considérons un système d'équations obtenu
en effectuant des moyennes spatio-temporelles. Ainsi, en vertu de la moyenne en temps
et en espace, les deux phases sont considérées comme coexistantes en chaque point x et
en chaque instant t. Par conséquent, ces équations modélisant l'écoulement diphasique
doivent être résolues simultanément dans tout le domaine représentant l'écoulement. Elles
forment un système de 2(d+ 2) lois de conservation (où d est la dimension d'espace) :
conservation de la masse, de la quantité de mouvement et de l'énergie pour chaque phase.

Equations de conservation

+ div{akpkUk) = Tjt,
dt

uk)) + grad(akPk) = <*kPkg + Mk + Mwk, (2.1)

dt

dt
div(otkPk{ek
-\~KWK .
H .
H ~7r-)uk) = CkPkguk + Ek + Ewk-
pk 2

20
 Ces équations représentent les lois générales de conservation de la masse, de la
quantité de mouvement et de l'énergie d'un point de vue macroscopique pour les deux
phases, vapeur (k = v) et liquide (k — /).
Afin de fermer le système (2.1), nous ajoutons les lois d'état constitutives suivantes :

k), (2.2)
k)- (2.3)
Et une loi topologique sur le taux de vide (voir section 2.1.4).
£k » Pk et pk sont reliés par l'équation d'état fondamentale suivante :

dek = Tkdsk + EjdPk. (2.4)

Pk
La température et la pression thermodynamique de chaque phase sont données par :

ds
k\Pk

Pk = Ajr • (2-6)

Les variables qui interviennent dans les lois de conservations et lois constitutives sont :

ak = taux de présence locale de la phase k,

Pk = pression pour la phase k,
Pk = masse volumique de la phase k,
Uk = vitesse de la phase k,
€k = énergie spécifique de la phase k,
hk = enthalpie de la phase k,
Sk = entropie de la phase k,
Tk = température de la phase k,
g = accélération de la pesanteur.

Fjt, Mk et Ek, sont les termes de transferts interfaciaux représentant respectivement

la masse, la quantité de mouvement et l'énergie acquise par la phase k.

Mwk et Ewk sont les termes de transfert entre la phase k et la paroi respectivement
de la quantité de mouvement et de l'énergie.

Les variables dépendantes principales choisies habituellement sont :

pg, pi pressions de la vapeur et du liquide,

ug, ui vitesses de la vapeur et du liquide,
ev, ei ou hvi hi ou sv, si comme variables thermiques,
a — av taux de présence de la vapeur.

21
Remarque 2.1 Pour une étude mathématique et numérique du système d'équations
(2.1), les variables dépendantes principales seront choisies suivant des considérations et
justifications mathématiques et physiques (par exemple, suivant que la solution du système
est continue ou non).
Selon BOURÉ [2], l'importance relative des variations de pk, Pk et hk pour les fluides
courants (eau et vapeur d'eau notamment) impose la nécessité de conserver pk parmi les
variables principales (exprimer pk en fonction de pk et de hk conduirait à des erreurs
importantes sur Pk en cas d'erreurs, même faibles, sur pk). D'ailleurs, pour celui qui veut
utiliser les lois réelles des fluides en se basant sur les tables thermodynamiques, la connais-
sance de pk est obligatoire pour calculer les autres grandeurs physiques thermodynamiques.

Cette question délicate sera discutée ultérieurement (voir chapitres suivants).

Contraintes de fermeture

Le système d'équations constituant le modèle doit être fermé. Les lois constitutives étant
les moins bien connues dans l'état actuel des choses, c'est sur elles que portent en pratique
les contraintes de fermeture. Ces contraintes dépendent en particulier :

— de la nature du problème posé (par exemple degré de régularité des conditions

initiales et aux limites envisagées) ;

— des hypothèses contenues dans le modèle (par exemple sur les coefficients de
corrélation globaux) ;

— de la façon dont le modèle est exploité (par exemple méthode numérique

utilisée).

L'étude détaillée des contraintes de fermeture est donc impossible dans le cadre
général de cette note. Elle devrait être faite pour chaque modèle particulier. Surtout que
les systèmes d'équations de conservation des modèles actuellement utilisés ne sont pas
tous rigoureusement équivalents. En général, elle est effectuée partiellement et au coup
par coup, au cours de la mise au point du modèle, à l'occasion des problèmes soulevés par
le non-respect de certaines des contraintes de fermeture.

Il existe toutefois, au moins lorsqu'on n'exclut pas le passage, par changement de

phase, d'un écoulement monophasique à un écoulement diphasique et vice-versa, des con-
traintes communes à tous les modèles :

Transfert interfacial

Lorsque o^ tend vers zéro ou l'unité, tous les termes de transfert interfaciaux doivent
tendre vers zéro. Cette contrainte, de bon sens, est respectée en général dans les modèles
existants. Cependant, la phase qui apparaît (ou disparaît) ne peut se développer (ou se

22
résorber) que s'il y a transfert de masse vers (ou en provenance de) cette phase, c'est-à-dire
si Fjt est non nul : tendant vers zéro en même temps que ak et 1 — ak, Tk doit donc ne
pas y tendre de la même façon. Autrement dit, il faut veiller à ce que :

lim > oo et lim > oo. (2.7)

arjt-M) ak a * - H 1 — QJ/fc

Selon J.BOURÉ [2], cela traduit le fait que, à l'apparition ou lors de la disparition de
la phase A:, a* et la surface spécifique S / (surface d'interface par unité de volume, au
travers de laquelle s'effectuent les transferts) ne tendent pas vers zéro de la même façon.
Dans l'exemple d'une configuration à bulles ou gouttes uniformes de rayon r, a* oc r 3 ,
S/ oc r 2 et S//ajt oc 1/r -» oo lorsque r ~» 0. On ne peut donc en particulier accepter
la formulation courante dans laquelle Fjt oc E/ oc orfc(l — ak). Son incohérence se traduit,
d'une part par la nécessité d'introduire, à la frontière d'apparition de la phase k, un a*
non nul, sous peine d'inhiber le développement de cette phase (et d'obtenir par ailleurs
une surchauffe prohibitive de l'autre phase), et d'autre part, par l'impossibilité pratique
de repasser en monophasique par vaporisation ou condensation totale (avec de même,
surchauffe ou sous-saturation prohibitive de l'autre phase). Une amélioration envisageable
de la modélisation courante serait d'adopter

Tk<xEi ce o%(l-ak)m, (2.8)

avec m < 1 (m ~ 2/3, dans l'exemple de la configuration à bulles ou à gouttes uni-

formes). Encore faudrait-il vérifier pour chaque modèle particulier, que cette forme résout
le problème ci-dessus sans induire, par ailleurs, d'autres incohérences.

Transfert pariétal :

Toujours lorsque ak tend vers zéro ou l'unité, les termes de transferts pariétaux
doivent tendre, soit vers zéro (phase apparaissant ou disparaissant), soit vers les termes
correspondants dans les écoulements monophasiques.

D'autres contraintes peuvent exister lorsque ajt est voisin de zéro ou de l'unité. De
nombreux termes devenant faibles dans les équations du modèle, le calcul numérique peut
poser des problèmes (exemple, matrices mal conditionnées). Il peut être alors souhaitable
de diviser certaines équations par ajt (en particulier les équations de quantité de mouve-
ment et d'énergie). Il faut alors veiller à ce que les limites, pour a* tendant vers zéro, des
quotients des termes constitutifs concernés par ajt soient finies.

Termes de transfert de quantité de mouvement

Termes de transferts interfaciaux

Le terme de transfert de quantité de mouvement à l'interface Mk s'écrit de la manière

suivante :
Mk = M [ + pik Vat + Mik. (2.9)

23
 est le terme de transfert de quantité de mouvement à l'interface par transfert
de masse.
Ml = Tkuik. (2.10)
Fjt est la masse acquise par la phase k, u,-fc notée aussi u£, est la vitesse moyenne phasique
à l'interface. Elle est associée aux processus d'échange de masse.
En général, on suppose que les vitesses phasiques à l'interface sont égales :

«iv = «i7 = «t- (2.11)

II y a une interprétation physique derrière cette supposition. La justification peut être

donnée en utilisant l'équation de bilan de quantité de mouvement à l'interface combinée
avec la condition de continuité de masse. Plusieurs auteurs ont développé différentes ex-
pressions pour «j :
Exemples :
Selon RANSOM [17]

{ ui si Tv > 0,

uv si r v < 0.
Selon PINET [15]
U{ = avui +
Selon STEWART and WENDROFF [26]

Etc.

•Mik est la force interfaciale généralisée pour la phase k, elle contient entre autres :
les forces de traînée, les forces transitoires, la force de Lift, la force de Faxén, la force de
Basset, la force de collisions des particules, etc. Dans ce rapport, on ne retiendra que les
forces les plus importantes, les autres seront négligées (pour plus de détails voir ISHII et ses
travaux avec d'autres auteurs comme MISHIMA, ZUBER, CHAWLA..., [8] et ses références).
La condition macroscopique de transfert de quantité de mouvement à l'interface (condition
de saut) est :

L'expression retenue pour Mik est la combinaison linéaire des forces les plus importantes
pour décrire les écoulement diphasiques dans les cœurs de Réacteurs Nucléaires.

MR + MYk + Mtk. (2.13)

Les termes apparaissant dans cette expression sont :
not
- Mil ée FD est la force stable de traînée (ou force de traînée stationnaire).

- M^. notée Mvm est la force de masse virtuelle.

24
 - Mfc noté Fi est la force de Lift.

est la force moyenne due aux tensions à l'interface, pik est la pression
moyenne phasique à l'interface. En général, on suppose que (voir section 2.1.4) :

Piv = Pu - Pi- (2.14)

Cette supposition est faite pour simplifier les calculs ; nous verrons plus loin dans quelles
conditions (ou régimes) elle est valable. On a donc finalement :

k=v,l

Termes de transferts pariétaux

•Muk est le terme de perte de charge, c'est la quantité de mouvement échangée entre
la phase k et la paroi.

Termes de transfert d'énergie

Termes de transferts interfaciaux

On décompose aussi le terme de transfert interfacial d'énergie comme suit :

Ek = ~Pik^~ + Tk(hik + ^ ) + Mtkuik + qik, (2.16)

avec :
•Pikdak/dt est le terme de transfert d'énergie dû à la pression.
U2

jfc + -•£*•) est le terme de transfert d'énergie dû au transfert de masse ; hik est
l'enthalpie moyenne de la phase k à l'interface.
•MikU{k est le terme de transfert d'énergie dû au transfert de quantité de mouvement
à l'interface.

•Ci* est la quantité (ou flux) de chaleur transférée au niveau de l'interface par con-
duction.

On ajoute la relation à l'interface suivante :

= 0, (2.17)
k=v,l

25
qui est, avec une bonne approximation, une conséquence de l'hypothèse de l'équilibre
thermodynamique local à l'interface.
Finalement les conditions de transferts interfaciaux sont les suivantes k

Tk = 0, (2-18)
k=v,l

Mk = 0, (2.19)
k=v,l

Ek = 0. (2.20)
k=v,l

Termes de transferts pariétaux :

'• il contient la puissance associée au frottement à la paroi et la quantité de

chaleur transférée au niveau de la paroi à la phase k.

2.1.3 Formes de lois de transfert et déséquilibres

Les transferts et les déséquilibres sont deux aspects inséparables de l'écoulement. Par
analogie avec de nombreux phénomènes physiques, il est raisonnable d'admettre que les
déséquilibres sont la cause des transferts. Inversement, le niveau des déséquilibres est
fonction de l'importance des transferts. On peut alors chercher à exprimer les transferts
("flux") en fonction des déséquilibres ("forces").

Lois de transfert sous forme de fonctions des variables dépendantes principales

Ces fonctions peuvent éventuellement faire intervenir des variables indépendantes. Il a été
montré (BOURÉ [l]) que le fait de représenter toutes les lois de transfert par des fonctions
conduisait à des contradictions au moins dans le cadre de l'hypothèse qui remplace clas-
siquement la loi topologique de taux de vide par l'égalité des pressions des deux phases.
Ces contradictions persistent lorsque pv ^ pi et concernent entre autres les phénomènes
de propagation (et donc de débit critique), ce qui est gênant pour l'étude des transitoires
ou, plus généralement, chaque fois qu'il apparaît des gradients (spatiaux ou temporels)
importants.
Les lois de transfert sous forme de fonctions conduisent à des vitesses de propagation de
petites perturbations indépendantes de ces lois. Ceci doit être rapproché de l'impossibilité
intrinsèque pour ces lois de transfert "fonctions"de rendre compte d'effets de proximité
(voir BOURÉ [2]).

La force stable de traînée

La force stable de traînée est due à la différence entre les vitesses phasiques. Elle est donnée

26
par [4]
FD = -~-
o
Ici le cœfficient CD est supposé dépendre de a j et Red (nombre de Reynolds), la dépendance
qui lie CD à, a<f et Red a été étudiée par ISHII et ZUBER [12] pour différents régimes
d'écoulement.
Pour des écoulements à bulles sphériques par exemple, Ai (aire interfaciale par unité de
volume) est fonction du rayon des bulles et de la fraction volumique de la phase dispersée
(vapeur).
A, = l
où Rd est le rayon moyen de courbure des éléments de la phase dispersée. Or,

cy = n-nR%,

ce qui implique :
Ai
" Rd
Ainsi, pour un écoulement à bulles,

FD = -J- ^ Ittrlttr, (2.21)

avec, ur = Ud — uc, (d : phase dispersée, c : phase continue).

Les forces de frottement à la paroi

Les forces de frottement à la paroi correspondent aux quantités de mouvement

échangées entre chaque phase k et la paroi. Les termes représentant ces échanges s'appellent
termes de perte de charge et sont notés Mwk- Pour chaque phase k, Mwk est de la forme
1
*" 8
avec

est la surface par unité de volume de la phase k en contact avec la paroi,

elle est déterminée suivant le régime de l'écoulement.

— Afc est le coefficient de frottement entre la phase k et la paroi.

Lois de transfert à termes différentiels

Les termes de transfert sous forme différentielle sont des fonctions linéaires des dérivées
partielles premières des variables dépendantes principales ; la forme de ces termes (mis à
part les fonctions ordinaires) est la plus simple possible. Elle est la seule qui ne change
ni la forme du modèle (équations aux dérivées partielles du premier ordre, en général)

27
ni le nombre d'équations aux dérivées partielles qu'il contient, ce qui a une importance
considérable (pour le calcul numérique et même analytique).
Les lois de transfert à termes différentiels rendent en général le système hyperbolique, cela
devrait donc conduire à des résultats satisfaisants quels que soient les gradients, à con-
dition toutefois que les termes du second ordre (dûs aux gradients des gradients) restent
faibles. Les instants et les régions où les gradients des gradients sont importants sont assez
limités en pratique (début d'une dépressurisation, apparition de l'ébullition, ...), mais ils
peuvent être cruciaux.
Les lois de transfert à termes différentiels ont toutefois des inconvénients sur le plan pra-
tique ; elles introduisent des coefficients inconnus subordonnés aux dérivées premières de
toutes les variables dépendantes principales. Il est impératif de limiter, ne serait-ce que de
façon arbitraire, le nombre de termes différentiels figurant en fait dans les lois de transfert,
car les coefficients inconnus jouent en pratique le rôle d'autant de paramètres ajustables.

La force de masse virtuelle et force de Lift

La force de masse virtuelle est une force transitoire qui est due à l'accélération d'une
phase par rapport à une autre.
Considérons par exemple une bulle sphérique en mouvement dans un fluide, l'accélération
de cette bulle provoque un changement d'énergie cinétique du fluide qui l'entoure, et induit
par conséquence une force résistante, cette force est calculée théoriquement dans l'article
de DREW et al. [5]. Dans un écoulement diphasique réel, la présence d'autres bulles doit
être prise en compte et l'expression générale de la force de masse virtuelle est donnée par
l'expression suivante [5] :
=Z
•*vm ^vmP

Cvm est la masse virtuelle par unité de volume de bulle,

pc est la masse volumique de la phase continue,
avm est l'accélération de la masse virtuelle.
En application du principe d'objectivité (qui traduit le fait que le comportement in-
trinsèque d'un matériau ne devrait pas dépendre de l'observateur). DREW et al. ont
développé une forme objective du terme avm, (c'est-à-dire invariante par changement de
repère galiléen) :
Dd{u Uc)
avm = ^~ + (ud - uc)[(A - 2)Vud + (1 - À)VttJ, (2.23)
qu'on peut écrire aussi :

a
™= ^ ~ ^W + (1
" x){ud ~ " c)v(Uc " Ud)' (2 24)
-
uc est la vitesse de la phase continue,
Ud est la vitesse de la phase discrète,
A est un paramètre arbitraire introduit par DREW et al. pour le développement de (2.23).
En se basant sur des formulations connues de avm pour des problèmes classiques tels
que l'accélération d'une bulle rigide sphérique unique dans un réservoir d'eau au repos, ou
d'une gouttelette unique accélérée dans un écoulement constant de vapeur et en comparant
avec (2.23) on trouve:
lim A = 2,
*0

28
 lim A = 0.
Ces deux cas asymptotiques montrent que le paramètre A doit au moins être fonction du
taux de vide. Pour des valeurs intermédiaires du taux de vide, l'expression de A (0 < A < 2)
n'est pas connue et sa valeur doit être déterminée expérimentalement.
HOUGHTON (1976) a démontré expérimentalement, pour des mélanges diphasiques dilués
que Cvm dépend de la forme des particules (Cvm = 1/2 pour une sphère seule et non
déformable), toutefois, dans le cas général et surtout pour des taux de vide plus impor-
tants, l'interaction entre les particules dispersées, qui semble diminuer le Cvm, n'est pas
bien connue. Cependant, le terme Cvm dépend du taux de vide a et doit être déterminé
expérimentalement.
La force de masse virtuelle apparaît dans les équations de quantité de mouvement sous la
forme suivante :
Mvm = adFVm = <XdPcCvm<lvm- (2.25)
Une nouvelle expression pour la force de masse virtuelle a été développée par DREW and
LAHEY [6]. Cette expression est considérée comme la plus générale :

Mvm = CvmadPc [ ( ^ + utVu^ - (j£ +UcVuc)] . (2.26)

La force s'exerçant sur une sphère en mouvement par rapport à un fluide en rotation
a été étudiée et calculée par PROUDMAN (1916) [6]. Cette force est perpendiculaire au
vecteur de rotation et au vecteur vitesse de la sphère et représente la force de Lift [6] (ou
force latérale de Lift).

FL = Mf = CLadpc (t»c - «<£) x (V x uc). (2.27)

Le cœfficient CL est appelé coefficient de Lift. Notons que la force de masse virtuelle et la
force de Lift comme elles sont présentées ci-dessus ne sont pas objectives séparément, par
contre, si nous posons Cvm = CL, nous obtenons que la somme Mvm + FL est objective.
Notons aussi que, en utilisant cette nouvelle formulation pour la force de masse virtuelle,
on supprime le paramètre A qu'on ne connaît pas, et qui apparaît dans la formulation
donnée Par DREW et al. (2.23), ce qui est de grande importance en pratique.

Puisque la force de masse virtuelle doit s'annuler dans le cas de la disparition de

l'une des deux phases, nous définissons un nouveau paramètre CM par:

<*dPcCvm = adacpCM, (2.28)

où
p = adPd + acPc- (2.29)
DlMENNA propose pour la force de masse virtuelle (ou la force de traînée dynamique),
dans le cas monodimensionnel, cette modélisation :

adacpCM\ 57 h«d ud-z~ >, (2.30)

L ot ox ox )
ce qui correspond à la forme générale développée par DREW et al. (2.23) pour A = 1.

29
Remarque 2.2 1. Pour le cas où a < 0.2 RUGLES et al. [21] ont développé une ex-
pression pour Cvm en utilisant l'expression de masse virtuelle (2.26)

2. Fvm ne peut être négligé dans le cas de très grandes accélérations spatiales (comme
dans le cas d'écoulement diphasique critique) (voir DREW et al. [5]).

Analyse de l'effet du terme de masse virtuelle pour le cas monodimensionnel

Dans ce paragraphe, on analyse dans certains cas particuliers l'influence du terme

de masse virtuelle sur la vitesse du son.

Considérons un système d'équations d'écoulement diphasique en équilibre thermique.

Le système contient quatre équations de conservation (conservation de masse et de quantité
de mouvement pour chaque phase). Ce système peut s'écrire sous la forme matricielle
suivante :

Où U = (ud, uc,p,a)T. Les termes d'échange à l'interface sous forme différentielle autres
que la force de masse virtuelle (donnée au (2.23) pour l'expression de avm) sont négligés, les
autres termes tels que la force de traînée, l'accélération gravitationnelle etc, apparaissent
dans le vecteur C à droite. Les valeurs propres associées au système (2.31) sont déterminées
en évaluant le déterminant caractéristique suivant :

(2.32)
En développant le calcul, on trouve :

( ) C v m ( T
l-Q yc yc 1 -

-~-Cvm(X - 2)(Ud - uc) - (1 + ^

Pd Vd Pd

+ { ) + r r i n I T -A) U
d-O ' (2-33)
yd 1-a (1-ar î/d (1 - <*r
avec :
yd = ud - n,
yc = uc - /x,

(1 - OL)pda\ apca\

L'équation (2.33) ne peut être résolue analytiquement à cause de sa forme complexe.

Cependant, pour certaines catégories d'écoulement, des simplifications sont faites afin de

30
la rendre abordable. Considérons par exemple le cas des écoulements à faible nombre de
Mach : si on suppose que les vitesses phasiques sont faibles devant les vitesses acoustiques,
les termes convectifs apparaissant dans (2.31) peuvent être négligés devant les termes en
dt (dérivées temporelles). L'accélération de masse virtuelle dans ce cas est réduite à :

yd et yc vérifient :
y d = y c = fi. (2.35)
L'équation (2.33) devient :

(2 36)
-
Les quatre racines de (2.36) sont réelles. Deux d'entre elles sont nulles. En utilisant la
relation (2.28) les deux autres sont liées à la propagation d'ondes acoustiques et sont
données par :
f ( l - q ) , _ot_] ~x
+ K
Pcal pda\\ '
où ac et ad sont les vitesses du son dans la phase continue et la phase discrète respec-
tivement. L'équation (2.37) définit la vitesse du son incluant l'effet du terme de masse
virtuelle sous la forme simple (2.34).

Si nous utilisons le terme de masse virtuelle (2.26), c'est-à-dire :

., n udud dud duc duc

nous obtenons les vitesses du son suivantes :

2 _ (l-Q){(1-a)Pd + apc} + PcCvm R l - a ) { a V1

(1 a)pp + PPC [ POl \

Dans le cas où l'effet de masse virtuelle est négligeable, les deux expressions (2.37)
et (2.38) coïncident et on obtient ce qu'on appelle la relation de l'onde stratifiée (ou la
vitesse de l'onde de compressibilité dans l'écoulement stratifié) développée par WALLIS
[29]:
- a)pd + ape (i-q) (2.39)
PdPc pda\

Par opposition, dans le cas où l'effet de masse virtuelle est très grand, la vitesse du
son dans un mélange homogène (obtenue en supposant que les deux fluides ont la même
vitesse) est :

" G -2 » - "lï""2. (2.40)

Pca c Pdad_

31
Pressions à l'interface

D'habitude, on posePk—pik = 0. Cela n'est pas toujours vrai ; par exemple dans les
situations où il y a expansion ou contraction des bulles. En revanche on pose pd = pu pour
la phase dispersée (vapeur). Cela est dû au fait que la phase dispersée occupe des régions
relativement petites où des grands gradients de pression ne peuvent pas être supportés.
Ainsi, on suppose que la contraction/expansion des bulles peut donner une contribution
en pk - Pik seulement pour la phase continue (liquide). De plus, la force due au terme
Pc — Pics e présente dans l'équation de quantité de mouvement sous la forme (pc — Pic)Vac.
Ainsi cette force peut avoir un effet seulement lorsque Va c ^ 0.
Dans les situations où la tension à la surface est importante, le terme pic — pu ^ 0, et les
deux pressions à l'interface sont liées par l'expression classique de Laplace :

Pic = Pu - Ha, (2.41)

avec H est la courbure moyenne de l'interface et a est la tension à la surface.

Si la phase dispersée est sphérique H = ^ , où Rj. est le rayon moyen de courbure des
éléments de la phase dispersée (bulles ou gouttelettes) (DREW [4]).
Pour un écoulement de fluide non visqueux, la pression à l'interface est donnée par
l'équation de Bernoulli :
Pic=Pc-tpcWr\2, (2.42)
où : uT — Ud — uc et £ = £ (a<j, Red), où Red est le nombre de Reynolds.
On pose en générale £ = constante pour simplifier. Notons qu'en l'absence de tension à la
surface, on a :
Pic = Pid- (2.43)

32
2.1.4 Hypothèses essentielles retenues
• Les termes visqueux et turbulents ne sont pas pris en compte dans les équations de
quantité de mouvement et d'énergie.

• L'interface est sans épaisseur (et donc sans masse F v = —F/) et sans propriété
physique.

• Vitesses moyennes au* interfaces : dans l'état actuel des choses, il y a peu de ren-
seignements sur ces vitesses moyennes. Pour l'instant, elles sont mises égales à une
même vitesse u:-, cette dernière est généralement approchée par l'expression :

U{ — avu\ + aiuv. (2.44)

Pressions dans les deux phases : on suppose que les deux phases ont la même pres-
sion :
Pv = Pl- P- (2.45)
II est devenu coutumier d'admettre pv = pi, ce qui diminue d'une unité le nombre
de variables dépendantes et remplace ainsi la donnée de la loi topologique de taux
de vide. Physiquement, le fait d'avoir pv ~ pi n'implique rien pour les valeurs en
x, t des dérivées spatiales et temporelles de pv et pi. Ces dérivées peuvent différer
sensiblement. L'hypothèse pv = pi est donc probablement trop restrictive (effets de
tension superficielle mis à part) lorsqu'on veut décrire correctement les phénomènes
de propagation et leurs conséquences, ce qui est nécessaire pour les transitoires rapi-
des et les débits proches de débits critiques ou plus généralement, chaque fois qu'il
apparaît des gradients (spatiaux ou temporels) importants.

Pressions à l'interface : on admet que les deux phases sont à la même pression à
l'interface (c'est-à-dire que la tension superficielle est négligée) :

Piv=Pii = Pi, (2.46)

comme on ne considère actuellement, dans les modèles retenus, qu'une seule pression
moyenne pour les deux phases : pv = pi = p, on est amené à définir l'écart Ap = p—pi

p-Pi = Çpc\ur\2. (2.47)

Pour les lois constitutives, on est fondé à adopter la forme la plus simple possible,
sous réserve d'y renoncer lorsqu'elle conduit à des contradictions. Ce processus de
remise en cause des hypothèses est constant et essentiel en physique.

33
2.2 Hyperbolicité et caractère bien-posé

Les systèmes d'équations établis pour décrire les écoulements diphasiques sont classés en
utilisant la méthode des caractéristiques. L'analyse montre que la plupart des modèles
proposés dans la littérature ont des caractéristiques à valeurs complexes dans les régions
physiquement intéressantes, notamment pour les écoulements diphasiques d'eau et de
vapeur. Des théorèmes bien connus (voir [19]) mènent à conclure qu'un système quasi-
linéaire d'équations aux dérivées partielles d'ordre un, qui a des caractéristiques complexes,
est mal posé en tant que problème aux valeurs initiales. Nous allons essayer d'étudier ce
phénomène dans les deux paragraphes suivants :

2.2.1 Caractère "bien-posé" des équations à caractéristiques complexes

Des caractéristiques complexes sont présumées être la cause de l'instabilité des solutions
numériques de modèles diphasiques. Ici, nous verrons que ce n'est pas toujours le cas.
Dans cette partie, nous examinons l'aspect théorique d'un mécanisme très important et qui
joue un rôle critique pour stabiliser les calculs ; c'est l'effet du coefficient de frottement
à l'interface sur les grandes longueurs d'ondes de la solution analytique. Les termes de
frottement à l'interface sont connus aussi par leurs effets stabilisant des grandes longueurs
d'ondes de la solution numérique (codes industriels tels que THYC [27], RELAP5 [20],...).
SHIEH et al. [23] ont montré analytiquement que ce terme stabilise les calculs pour toute
longueur d'onde si les mailles ne sont pas très fines. Cela résulte du fait que les petites
longueurs d'onde de la solution analytique ne peuvent être représentées sur un maillage
grossier. En conséquence, l'effet stabilisant des grandes longueurs d'onde par le terme de
frottement à l'interface devient un facteur dominant sur un maillage grossier.
L'importance du coefficient K de frottement a été constaté par SHIEH et al. en effec-
tuant trois tests différents : the Edouard's pipe problem, the Chistensen's subcooled boiling
problem, et the Oak Ridge National Laboratory (ORNL) void profile test problem.

Définition 2.1 Le problème aux valeurs initiales et aux limites "des écoulements
diphasiques" est bien posé si la solution dépend d'une façon continue de la donnée initiale
et des conditions aux bords ; autrement dit :
si u(t) et v(t) sont deux solutions correspondant respectivement aux deux données initiales
rio et VQ, alors :
\\U(t) - v(t)\\ < G\\u0 - vo\\.

Pour que, si les états initiaux sont très proches l'un de l'autre, il en soit de même pour les
états finaux pour un même temps t, il faut que le facteur G soit inférieur à 1. En effet,
pour des problèmes hyperboliques sans terme source le facteur G est en général inférieur
à 1.

34
Equation scalaire

Nous allons étudier l'effet de valeurs propres complexes sur la solution analytique d'une
équation différentielle simple :

| + f l | + K l l = 0, (2.48)

avec : p = nr + i/x,-,
fir et m sont deux constantes ne dépendant pas de x et de t.
K > 0 est le coefficient de frottement à l'interface.
Soit
u = uoexp[i(kx — wt)] (2.49)
une solution de l'équation (2.48). On a alors,

w=^tti + k^ (2.50)
i

d'où
u = uoexp[ikx — ikfxrt — (K — km)t]. (2.51)
Ainsi, l'équation différentielle (2.48) a une solution bornée quand t —> oo si : a t = to, u
est une combinaison linéaire de composantes de FOURIER exp[i(kx — wto)] avec k < j-S
Autrement dit ; si la donnée initiale du problème linéaire est suffisamment régulière, alors
l'existence de racines complexes ne cause pas de problèmes.

Ainsi, pour le cas d'une équation scalaire, le caractère "bien posé" du problème
dépend uniquement de la régularité de la donnée initiale et de la magnitude du coefficient
de frottement interfacial.

Stabilité

Nous adoptons le concept de stabilité au sens de SANZ-SERNA [22]. Ici le mot stabilité
traduit le comportement de la solution numérique pour des pas d'espace fixés, lorsque le
nombre de pas de temps augmente (ce qui est différent de la notion de stabilité au sens
de LAX-RICHTMYER [19]).
Après avoir discrétisé l'équation différentielle (2.48) par un schéma décentré, nous utilisons
le critère de stabilité suivant :
||tr+ 1 ||<||17"|| (2.52)
avec Un est la solution numérique au niveau de temps n. Généralement, on utilise la
stabilité comme elle est définie dans le livre de RICHTMYER and MORTON[19], c'est-à-dire,
si :
Un+1 = C{At)Un, (2.53)

C(&t)n doit être borné en norme pour tout n, quand At tend vers zéro et t < T .

Analyse de Stabilité pour le cas d'une équation scalaire

35
 Considérons l'équation scalaire à caractéristique complexe suivante :

O (2.54)

avec nr > 0, et des conditions aux limites périodiques.

Considérons le schéma décentré explicite de (2.54) suivant :

Un — Un TTTt TTH
j i' \ (, ir + ifii Uj
+ U
-Uj-i
+ KU?+X - 0 (2.55)
At ' Ax
qui peut se réécrire

+ KAt)U?+l = (l-o) - ib{UJ - Uf_ (2.56)

où
: At
a =

At
b =

avec a > 0.
Supposons que la solution numérique complexe U à l'instant n vérifie l'inégalité suivante :

\Vf ~ U?-i\ < L\*i ~ *;-il(|tf"l + ll^LxD/2 (2.57)

alors nous avons le résultat suivant :

Théorème 2.1 Supposons que Un satisfait le critère (2.57) et les conditions :

L\m\ <KetO< /x r (Ai/At) < 1, alors

\\Un+1\\L2<\\Un\\L2. (2.58)

Preuve :
Remarquons d'abord que toute fonction discrétisée sur un maillage de pas d'espace A i
tel que A i > j satisfait la condition (2.57), car la condition (2.57) dans ce cas est tout
simplement l'inégalité triangulaire ; cependant, les fonctions régulières peuvent satisfaire
à cette condition sur un maillage plus fin.
Nous avons :

(1 + K&t)Uf+1 = (1 - a)Uf + aC/JLi - ib{Uf - U^). (2.59)

et
En utilisant 6 = M«^ l'inégalité (2.57), on obtient :

(l + A'A02|C/;+1!2 < Kl-aWf + aUf.tf + L^A

- a)Uf + aUJ.,1 (2.60)
En utilisant l'identité : \Z\2 = ~ZZ, on a :

Wu?-i+W) ^ E I ^ I 2 + E i ^ i i2- (
3

36
Pour 0 < a < 1 , en utilisant l'inégalité (2.61) et des conditions aux limites périodiques,
on obtient :
Y,\{l-a)Vf + aU?_1\2<'£\U?\2 (2.62)
j i

de même, en utilisant l'identité ab <

(2.63)

par conséquent :

+ KAt)2\Uf+1\2 < (1 + 2L\m\At + L2\m\2At2) £ \U* (2.64)

Or, si L\pi\ < K , on obtient l'inégalité finale :

(2.65)

Système d'équations

RANSOM and HICKS [18], STEWART and W E N D R O F F [26] et STEWART [25] ont étudié le
système à quatre équations ci-dessous

d(aPv)
+ V • apvuv = 0,
dt

~ <*)PÙ+ V • (1 - a)p ui = 0,
t
Ôt
(2.66)
du
VV • uv) + aVp = K(ui — uv),

(1 _ a)pi(-~ + tt/V • m) + (1 - a ) V p = K - ti/).

Ces équations décrivent un écoulement diphasique isotherme dans une conduite droite. Les
termes à droite des équations d'impulsion sont égaux et opposés, ils représentent la force
de traînée stable sous sa forme générale comme elle est utilisée par différents auteurs. Pour
chercher les caractéristiques du système dans le cas monodimensionnel, on suppose que le
liquide est incompressible (hypothèse raisonnable, à moins que a soit très petit) et que
la loi d'état pour la vapeur est donnée par pv = pv(p). Ainsi, les racines A du polynôme
caractéristique satisfont à :

ap,{\ - m)2 + (1 - a)pv(X - uv)2 - aptcZ2(X - uv)2(X - ut)2 = 0 (2.67)

37
avec

STEWART [25], dans son article, montre que ce polynôme caractéristique admet deux
racines réelles ; une légèrement inférieure à uv — cv et l'autre légèrement supérieure à
uv + cv. La courbe du polynôme caractéristique (voir plus loin (sous-section 2.2.2)) mon-
tre que les deux autres racines réelles éventuelles doivent être comprises entre les deux
premières ; seulement, l'étude de l'équation caractéristique montre que cela est impossible
dans le cas où uv ^ uj. L'équation admet deux racines complexes conjuguées (puisque les
coefficients du polynôme caractéristique sont réels).
Si uv, ui -C cv , les deux racines complexes sont approximativement

^iTL.{Vv.ul) (2.69)

où

otpi
Le système (2.66) est discrétisé par une méthode aux différences finies semi-implicite
(voir STEWART [25], SHIEH et al. [23], HARLOW and AMSDEN [7] et LILES and REED
[14]), ces différents auteurs utilisent le même principe de discrétisation : les effets dûs au
propagations soniques sont traités d'une façon implicite, tandis que les phénomènes con-
vectifs sont traités explicitement, de façon à ce que le pas de temps At soit limité par
At < min(Ax/uv,Ax/ui). Cette discrétisation est faite sur un maillage à mailles décalées
(la pression et le taux de vide sont calculés aux centres des mailles, tandis que les vitesses
sont calculées aux interfaces j — \ ) .

Les résultats obtenus a- : des schémas similaires à celui décrit ci-dessus ont montré
que ceux-ci se comportent bit TK>ur toutes longueurs d'ondes (stables) sur un grand nom-
bre de cas pratiques. En outr si l'amplitude de K est réduite (K petit), des instabilités
apparaissent pour les grandes longueurs d'ondes. Ce phénomène est illustré par un cas test
dans STEWART [25].

Conclusions et commentaires

II est donc montré (STEWART, SHIEH et al.), qu'un calcul numérique avec un schéma
aux différences finies d'un modèle bi-fluide d'écoulement diphasique est raisonnable, s'il
y a suffisamment de transfert de quantité de mouvement entre phases et si les mailles
ne sont pas très fines, en dépit du fait que le système d'équations a des caractéristiques
complexes. De plus STEWART a donné un critère lié à la taille des mailles et l'a interprété
physiquement : pour le modèle avec termes d'échange de quantité de mouvement dans le
cas d'écoulement à bulles ou à gouttelettes, la taille critique des mailles est proportionnelle
au rayon moyen des gouttelettes ou bulles. Autrement dit, les calculs approchés peuvent
être raisonnables, si on ne s'intéresse pas à l'étude de phénomènes plus fins que l'échelle
supposée ou prise en compte dans l'expression de la loi de transfert de quantité de mou-
vement. C'est le cas où le modèle bi-fluide est destiné à un calcul approché d'un régime
d'écoulement à bulles ou à gouttelettes sans tenir compte par exemple des fluctuations
causées par une gouttelette ou une bulle individuelle.

38
 Ainsi, pour une modélisation moins fine d'écoulement diphasique, le modèle est
mathématiquement mal posé, mais les calculs approchés sont raisonnables à condition
qu'on reste dans l'échelle appropriée.

Cependant, dans le cas de modèle sans transfert de quantité de mouvement (écoule-

ment stratifié par exemple), la variation du taux de vide est liée au changement de la
surface de contact entre l'eau et la vapeur, or, d'après STEWART, toute solution approchée
autorisant le taux de vide à varier nécessite une résolution sur une échelle plus fine quand
K = 0. Dans ce cas, un problème bien posé est obtenu en prenant en compte les effets
de la tension à la surface (voir RAMSHAW and TRAPP [16]). La même chose en principe,
devrait être faite pour d'autres régimes d'écoulement, si toutefois, une résolution très fine
est souhaitée.

Dans les cas complexes du calcul d'un transitoire réel, l'instabilité des grandes
longueurs d'ondes peut avoir des effets significatifs sans que sa présence ne soit mani-
feste.

En conséquence, pour une modélisation fine, le modèle approprié est celui qui serait
mathématiquement bien posé comme problème de valeurs initiales et qui puisse décrire
l'instabilité physique de HELMHOLTZ.

2.2.2 Hyperbolicité de quelques modèles

La plupart des modèles diphasiques développés pour décrire les écoulements diphasiques
inhomogènes en déséquilibre sont représentés par un modèle bi-fluide. Le modèle bi-fluide
est un système de six équations aux dérivées partielles représentant la conservation de la
masse, de la quantité de mouvement et de l'énergie pour chaque phase.

L'objectif principal des développeurs de codes est la description d'écoulements dipha-

siques transitoires avec déséquilibres thermiques et mécaniques et éventuellement les pro-
cessus complexes d'ébullition et de condensation dûs aux transferts de chaleur. Au niveau
international, les codes les plus connus sont les deux codes américains RELAP (INEL)
et TRAC (LANL), et le code français CATHARE (CENG) ; ces trois codes utilisent des
modèles bi-fluide. Pour les processus de transferts interfaciaux, on utilise généralement des
expressions algébriques pour fermer le système d'équations : c'est le modèle de WALLIS.
La vérification et l'évaluation détaillées de ces codes, basées sur la comparaison des valeurs
expérimentales et mesurées pour un grand nombre de tests ont démontré leurs capacités à
décrire le comportement global du système dans plusieurs cas d'accidents. Cependant, cer-
tains de ces codes présentent des imperfections (ou manques) ainsi que des limitations. La
plupart de ces imperfections sont liées à la prédiction de phénomènes d'écoulement local
caractérisés par de forts gradients spatiaux ou discontinuités des grandeurs physiques ; par
exemple : propagation du front de trempe dans le cas d'un renoyage du coeur, le passage
d'une poche d'eau non borée, écoulement critique...

39
 Les raisons majeures de ces déficits peuvent être les suivantes:

1. Généralement, le caractère non-hyperbolique du système représente un problème de

valeurs initiales et aux limites mal-posé qui nécessite des termes d'amortissement
numérique pour avoir des résultats stables.

2. L'utilisation de schémas numériques d'ordre un.

3. L'utilisation de schémas numériques basés sur des discrétisations utilisant des mailles
décalées et la technique cellule donneuse "donor cell" qui est ambiguë dans le cas
d'écoulement à contre-courant.

Pour pallier ces insuffisances, divers modèles ont été développés en apportant des correc-
tions au modèle communément utilisé (modèle de WALLIS) afin de le rendre hyperbolique;
cela permettrait donc l'utilisation des méthodes numériques modernes développées pour
la résolution de systèmes hyperboliques et qui permettent une résolution plus fine.

Il existe deux approches différentes pour rendre le système d'équations bi-fluide

hyperbolique :

1. Prendre deux pressions différentes ; une pour chaque phase ; et prévoir des équations
différentielles pour décrire la dynamique de la différence de pression entre les deux
phases. Cette approche conduit à un système de 7 voire 8 équations (voir [26]) et
s'accompagne de la difficulté supplémentaire de fournir d'autres relations constitu-
tives pour fermer le système d'équations.

2. Supposer l'égalité locale des pressions pour les deux fluides, ce qui conduit à un
système à 6 équations. Et introduire d'autres termes différentiels supplémentaires
dans les équations de quantité de mouvement de chaque phase. Ces termes peuvent
être interprétés comme des forces interfaciales (force de masse virtuelle, force de
Basset, terme dû à la pression à l'interface, ...).

Avec ces deux approches, on essaye de regagner une partie de ce qui a été perdu lors des
opérations de moyenne qu'on fait subir aux équations locales instantanées.

En général, c'est la deuxième approche qui est adoptée. Les termes différentiels
habituellement retenus sont : la force de masse virtuelle et/ou le terme de pression à
l'interface.

Remarque 2.3 Dans les modèles d'écoulement diphasique avec ou sans correction (i.e.
avec terme de pression à l'interface et/ou terme de masse virtuelle), les équations d'énergie
peuvent être découplées du reste et l'étude de l'hyperbolicité se ramène à l'étude des racines
d'un polynôme caractéristique de degré quatre au lieu de six.

40
 Dans le paragraphe suivant on présente un modèle incluant un terme différen-
tiel introduit par STADTKE et al. [24]. Dans le paragraphe d'après on étudie l'effet sur
l'hyperbolicité d'un terme de pression à l'interface.

Modèle 1

Considérons le modèle diphasique à six équations suivant

( ) + -^{aipiui) = 0,

d d dv (2<71)
—(aip,ut) + -foiaiPïti) + ai-£ = F/"\

S {^- + 1 » ) + 1 {"'^' + 7,+1

R T»)+ £
Pour fermer ce système, nous ajoutons la relation

ag+ai = l, (2.72)

et les lois d'état suivantes (k = g,l) :

k), (2.73)

k), (2.74)
= defc - ~dpk, (2.75)

ainsi que

= Tk, (2.78)

41
 Pk P
—~- — —5-.
Pi Pi

Le terme F™ [k = (7, /) apparaissant dans les termes de droite du système d'équations

représente une force interfaciale non visqueuse.

-agai(agpt - aiPg)(ug - ut) l-^ - — 1

~£ (2.80)

Ce terme à été introduit par STÀDTKE et al. [24]. Les critères utilisés pour déterminer une
telle expression sont :

• Les termes non visqueux de frottements interfaciaux ne doivent pas affecter la somme
des équations de quantité de mouvement.

• Les termes non visqueux de frottements interfaciaux ne doivent pas contribuer à la

dissipation de l'énergie mécanique.

• La matrice jacobienne n'a que des valeurs propres réelles et qui admettent une in-
terprétation physique.

• La matrice jacobienne a un ensemble complet de vecteurs propres.

• le système d'équations doit traiter comme cas limites : l'écoulement monophasique

gaz/vapeur (ag —> 1), ou liquide (ag —> 0), ainsi que l'écoulement homogène
(ug = m).

• Le système d'équations doit fournir implicitement des conditions "raisonnables" d'un

écoulement critique.

L'inconnue principale choisie est le vecteur

U = (p,u s ,u;,a f l ) s 5 ,s/).

Supposons que les solutions du système soient régulières ; (2.71) est équivalent au système

42
suivant:
._2dp dag

ÔUg ÔSg
dx dx = 0

dui dsi dsi.

~
— + Ul—) = 0
(2.81)

8t

rrnv

Si nous supposons que uini — pour la phase k, nous obtenons le système (2.82)
équivalent à (2.81) suivant:

dota doia

dp dai . . 0 dp da\
+ pl + QlUl{ai)
ft -bT Yz
dp

(2.82)
dut dp
—
~ l - = Fl

9
— dx

dsi
l
dt dx"-

Remarque 2.4 On posant utnt = Uk pour la phase k, nous obtenons un système (2.82)
hyperbolique dont les valeurs propres sont données ci-après. Cependant, on n'a pas con-
servation de l'énergie totale dans le système de départ !

43
 Les valeurs propres du système (2.82) sont :

= u3 > convection
A2 = J
3
~ Ug \ "onde de taux de vide"
4 =Ui J
A5 = u + a \
= u - a >Jonde de::pression
avec
( )
u= (2.83)
C-
PgPl
u est la vitesse moyennée des deux phases :
a2 = a2 - Aa2 (2.84)

!
n2 _ OtgPl +
0
~
(2.85)
ai aï PgPl

(2.86)

Remarque 2.5 Le* mfesses de propagation correspondant aux taux de vide sont-elles
correctes du point de vue physique ?

Modèle 2

On considère ici un modèle d'écoulement diphasique simplifié ; c'est un système à qua-

tre équations constitué des équations de conservation de la masse et de la quantité de
mouvement pour chaque phase. Dans ce système, on retient le terme dû à la pression à
l'interface Ap.

0 = o,
(2.87)
d
>.«•)+£<«. o + Ap—ag = 0,
ox *
. d .
Qp iuf ) + a•t-^P-i - Ap—ai = 0.
dt^ ' ox
44
 Pour étudier Phyperbolicité de ce système, une première approche consiste à sup-
poser que le liquide est incompressible par rapport à la vapeur (i.e. p\ — constante ;
raisonnable, à moins que ag soit très petit). Pour fermer le système d'équations, on ajoute
les relations et lois d'états suivantes :

Pg = Pg(p),

^ £ 1
dp ~ c | '
Ap = (p-pi),

Le système (2.87) est de la forme :

dw
(2.88)
~W
ou
W = (otgpg, <*ipi, agpgUg,aipiui) . (2.89)

Soit W une solution régulière de (2.88) ; avec un changement de variable W = <p{U),

on a :

~ (2.90)

équivalent à

(2.91)

ou
B(U) = (U). (2.92)

Les matrices A(W) = A(<p(U)) et B(U) sont semblables, elles ont donc les mêmes valeurs
propres.
Le vecteur inconnu U est choisi de façon à avoir une forme simple de la matrice B(U) .
Le système (2.88) est formellement équivalent au système d'équations

r d_ d „ , PjK-tt/) à (1 - a) d
a
a dx
d Fi fi
l
dl
(2.93)
d Ap dd cl dd dd
Ft' a
apg, ^Z + -T^ZPg
pg + u9-*Zug = 0,

Ap d c d d
_ r ^ —-z-Pg + u/—«/ = 0.
Ul pi dx * dx
L dt (1 - a)pidx
Posons
(2.94)
le système (2.93) s'écrit donc

(2.95)
la matrice B(U) est donnée par
/
Un
a
P3 —^-P3

0 (a-1)
B(U) =
Pg <*Pg

-Ap
\ pi L - a)pi u/

Hyperbolicité

Nous allons étudier l'hyperbolicité du système (2.87), ce qui est équivalent à l'étude
de l'hyperbolicité du système (2.93).
Soit
fi = {17 G JR4; pg > 0, a G]0,1[, %, t«/ G iR}. (2.96)
Soit C7 € fi ; avec le changement de variable W —> U défini par (2.94), les valeurs
propres À, de la matrice A{W) sont les mêmes que celles de la matrice B(U). Le polynôme
caractéristique P(\) de B(U) est donnée par :
P(X) det{B(U) - XI) (2.97)
Ug - À
(l-«)n
J
a a

i — A 0
(a-1)
(2.98)
us — A
<*P<? 0

-Ap
0 u/ — A
pi (1 - a)pi
Ainsi,
s2 APi
P(X) = -A)2- -A)2-

(2.99)
a pj

46
P(X) peut être écrit sous une forme plus simple en utilisant la transformation de Galilée
Au Au
2 ' 2 '
avec Au = ug — u/.
Posons donc :

(2.100)

et

(2.101)
2cg
il vient
(A - u 5 ) 2 =

P(X) est réduit à

=4
Ap
a pi
Nous allons maintenant étudier les racines de P(X). Ce qui revient à étudier les racines de

(2.104)

En posant

(2.105)

et

(2.106)

on a
(2.107)

Cas A P = 0 :
Dans ce cas, l'étude des deux fonctions f et g (voir fig. 2.1), montre qu'on a deux racines

47
 Figure 2.1: Sans correction, Ap= 0.

réelles et deux autres complexes ; autrement dit, le système d'équations (2.87) n'est jamais
hyperbolique. On peut aussi illustrer ceci par un contre-exemple. En effet, le polynôme
P(v?) peut s'écrire

(2.108)

On a un polynôme de la forme :
(2.109)
Soit
K3)
E =
2
- K2)
A =

a0 = 62{S2-2T.)

= -45A

a2 — -
peut alors s'écrire :
P(ip) = <p4 a0. (2.110)

48
Supposons A'i = K21 c'est-à-dire (1 — ot)pg = api, on a donc ai = 0. Dans ce cas

P(<p) = (p4 + 02 tp2 + a0.

P est symétrique autour de <p = 0.

OÙ O2 < 0.
Ainsi, en appliquant le Lemme 2.1 (voir ci-après), r m = 0, et la condition Q(rm) > 0 se
réduit à ao > 0,
c'est-à-dire
S2(62 - 2E) > 0, (2.111)
i.e.
6=0
ou
S2 > 2E.
5= 0 <=>• ug = ui (écoulement homogène).
Pour S ^ 0, (2.111) n'est pas satisfaite si

0 < S2 < 2E,

<*f*£ < 2E.

Or nous avons
2E = K\ + K2 = 2K2 = 2.
Il vient alors
0 < (ug - m)2 < 8c2g,

ce qui est vrai malheureusement. La relation (2.111) n'est donc pas satisfaite.

Cas AP ± 0
L'écart de pression Ap est normalement déterminé suivant le régime d'écoulement con-
sidéré. Cependant, en utilisant l'expression (2.47) avec £ = 1 et pc = pg, on a

-ui)2. (2.112)

T h é o r è m e 2.2 La relation (2.112) étant satisfaite, le système d'équations (2.87) est hy-
perbolique sur le sous-ensemble Çl° de fi, où

fi0 = {U € fi! Wg - Ufl < Cg)

49
 Figure 2.2: Avec correction, Ap ^ 0.

Preuve
II suffit d'étudier les fonctions f et g définies par (2.105) et (2.106) (voir figure 2.2) pour
montrer que le polynôme caractéristique P(v?) a toutes ses racines réelles.

Cas général
En se plaçant dans le cas général pour le système d'équations (2.87), i.e. : le liquide est
compressible et Ap n'est pas donnée sous forme explicite. Un calcul fastidieux donne le
polynôme caractéristique du système (2.87) :

P(X) =
<*sAp.
- 7 .113)

avec

Qgpiaf
En posant

27 " 7
On obtient une forme particulière du polynôme P(A), ainsi

P(A) =

50
avec

a
l "g

+ ^ (2.114)

C?(v) peut être présenté sous la forme générale suivante

Qte) = (<p2 - S2)2 - K2(<p+6)2 - Kx{v - S)2 + A3 (2.115)

Si A P = 0. le théorème énoncé ci-après suffit pour montrer que Q(<p) a deux racines
réelles et deux racines complexes.

Théorème 2.3 Considérons le polynôme suivant :

P{<p) = {y2 - S2)2 - K2{<ç + 8)2 -Ki{<p- <*)2

avec 6 réel et K\ et Ki sont deux réels positifs, alors il y a deux cas où toutes les racines
de P sont réelles.
(i) K\ = 0 ou A'2 = 0 ou S = 0 , ou bien :
(ii) Ki>0,K2>0,6£0.et {K\'Z + K\'3f < AS2

Preuve
(i) Evident.
(ii) Posons X = <p — 5, il vient

= R(X)
= X2{X + 2S)2-K1X2-K2{X 2

= F(X)-G(X),
avec

F(X) = {X2-K2){X
et
G{X) = KiX2,
on a G{X) > 0, F(X) > 0 sur ] - 00, -y/K^\ U [y/KÏ, +oo[ et F(-y/K^ = F{y/KÏ) = 0.
Puisque F tend vers +oc plus vite que G quand \X\ —> 00, il existe au moins deux racines
réelles de l'équation F{X) - G(X) = 0.
P(<p) admet une racine double s'il existe <po telle que

=0 et Fivo) = 0,

51
ce qui revient à :
F{Xo)=G{Xo) et F'{Xo) = G'{Xo) (2.116)
avec Xo = <po - S.
F'{X0) = G'{X0) est équivalent à
(2.117)
Multiplions (2.117) par XQ et remplaçons le terme obtenu à gauche (c'est-à-dire G(Xo))
par F(Xo). Après un calcul simple, on obtient

Xo = V-26K2. (2.118)
Maintenant, multiplions (2.117) par (Xo + 26) il vient
o + 26) = X0(X0 + 25f + F(X0) (2.119)
et remplaçons dans (2.119) F(Xo) par G(X0). Un calcul simple donne (Xo + 26)3 = 26Ki
i.e.
X o = ^/25Ki - 26. (2.120)
De (2.118) et (2.120) on tire la condition nécessaire et suffisante pour que (2.116) soit
réalisée i.e. :
26= ^/26Ki + y/26K2
ou bien :
{I<1/3 + I<l/3f = (26)2. (2.121)
Ainsi, pour que toute les racines du polynôme P(X) soient réelles, il faut et il suffit que

(A'i /3 + K\IZf < 4<S2. (2.122)

Si AP T^ 0, la condition nécessaire et suffisante pour avoir toutes les racines réelles doit
être calculée suivant chaque expression de Ap, ou bien, il faut trouver une formule reliant
entre eux les coefficients Ki, K A3 et Ap. Ou bien, développer le polynôme Q{ip) et
le mettre sous la forme général*; aes polynômes de degré quatre et appliquer le lemme
ci-après.

Lemme 2.1 Considérons le polynôme de degré quatre et à coefficients réels suivant

Q{\) = A4 + 63A3 + 62A2 + 62A + 60.
Une condition nécessaire et suffisante pour que Q(X) ait quatre racines réelles est que :
Q ait trois racines réelles rg < rm < rj telles que les trois inégalités suivantes soient
satisfaites
Q(rg) < 0, Q(rm) > 0, Q(rd) < 0

Pour la démonstation de ce lemme, voir (D.L. HICKS [9]).

52
Bibliographie

[1] J. A. BOURÉ, Mathematical Modeling and the Two-phase Constitutive Equations,

European Two-phase Flow Group Meeting, Haifa, 1975.

[2] J. A. BoURÉ, Les lois constitutives des modèles d'écoulements diphasiques monodi-
mensionnels à deux fluides , Rapport CEA-R-4915, 1978.

[3] J. M. DELHAYE, Basic Equations for Two-phase Flow Modeling, Two-phase Flow
and Heat Transfer in the Power and Process Industries, A. E. BERGLES et ai, Eds.
Hemisphere Publishing, Washington, 1981.

[4] D. A. D R E W , Analytical Modeling of Multiphase Flows, Edited by R. T . LAHEY Jr.,

Boiling Heat Transfer, Modern Development and Advances, Elsevier Science Publish-
ers, 31-84, 1992.

[5] D. A. D R E W , I. CHENG k R. T . LAHEY, The Analysis of Virtual Masse Effect in

Two-Phase Flow, Int. J. Multiphase Flow, Vol. 5, No 4, 233-242, 1979.

[6] D. A. D R E W k R. T . LAHEY, The Virtual Masse and Lift Force on a Sphere in

Rotating and Straining Inviscid Flow, Int. J. Multiphase Flow, Vol. 13, No 1, 113-
121, 1987.

[7] F. H. HARLOW k A. A. AMSDEN, J. Computational Phys., 17, 19-52, 1975, k, 18,

440-464, 1975.

[8] G. F . H E W I T T , J . M. DELHAYE k N. ZUBER Multiphase Science and Technology,

Volume 5, Hemisphere Publishing Corporation, 1990.

[9] D. L. HICKS Hyperbolic Models for Two-Phase (or Two-Material) Flow, Sandia Na-
tional Laboratories, SAND81 0253, 1981.

[10] M. ISHII, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975.

[11] M. ISHII k K. MISHIMA, Two-Fluid Model and Hydrodynamic Constitutive Relations,

Nuclear Engineering and Design, 82, 107-126, 1984.

[12] M. ISHII k N. ZUBER, Relative Motion and Interfacial Drag Coefficient in Dispersed
Two-phase Flow of Bubbles, Drops and Particles, A. I. Ch. E. Journal, Vol. 2 5 ,
843-855, 1979.

53
[13] R. T. LAHEY & D. A. DREW, The Three Dimensional Time and Volume Aver-
aged Conservation Equations of Two-phase Flow, Advances in Nuclear Science and
Technology, Vol. 20, 1-69, 1988.
[14] D. R. LILES & H. REED, J. Computational Phys., 26, 390-407, 1978.
[15] B. PlNET, Modélisation des écoulements pour l'analyse de l'accident de perte de ré-
frigérant primaire (APRP), Entropie, 83, 1978.
[16] J. D. RAMSHAW & J. A. TRAPP, Characteristics, Stability, and Short-Wavelenth
Phenomena in Two-Phase Flow Equation System, Nuclear Science and Engineering,
66, 93-102, 1978.
[17] V. H. RANSOM, Two-Phase Flows, Ecole d'été, Analyse Numérique, 12-13 Juin, 1989.
[18] V. H. RANSOM & D. L. HICKS, Hyperbolic Two-Pressure Models for Two-Phase
Flow, J. Comput. Phys., 53, 124-151, 1984.
[19] R. D. RICHTMYER & K. W. MORTON, Difference Methods for Initial-Value Prob-
lems, 2nd éd., Interscience Publishers, John Wiley & Sons, New York, 1967.
[20] R. A. DIMENNA, J. R. LARSON, R. W. JOHNSON, T. K. LARSON, C. S. MILLER,
J. E. STREIT, R. G. HANSON & D. M. KISER, RELAP5/MOD2, models and cor-
relations, NUREG/CR-5194, EGG-2531,1988.
[21] A. E. RUGGLES, R. T. LAHEY J R . , D. A. DREW &: H. A. SCARTON, An Investi-
gation of the Propagation of Pressure Perturbations in Bubbly Air/Water Flows, J.
Heat Trans., 110, 494-499, 1988.
[22] J. M. SANZ-SERNA, Stadies in Numerical Nonlinear Instability, Why Do Leapfrog
Go Unstable, SIAM J. Sci. Stat. Comp., 6, 4 ,Oct. 1985.
[23] A. S.-L. SHIEH, R. KRISHNAMURTHY & V.H. RANSOM, Stability, Accuracy, and
Convergence of the Numerical Methods in RELAP5/MOD3, Nuclear Science and
Engineering, 116, 227-244, 1994.
[24] H. STADTKE k R. HOLTBECKER, Hyperbolic Model for Inhomogeneous Two-Phase
Flow, International Conference on Multiphase Flows, University of Tsukuba, Japan,
September 24-27, 1991.
[25] H. B. STEWART, Stability of Two-Phase Flow Calculation Using Two-Fluid Models,
J. Comput. Phys., 33, 259-270, 1979.
[26] H. B. STEWART & B. WENDROFF, Two-Phase Flow: Models and Methods, J. Com-
put. Phys., 56, 363-409, 1984.
[27] THYC, Modélisation et Méthodes Numériques, E.D.F./D.E.R., Version 3.0, 1993.
[28] A. R. D. THORLEY & D. C. WlGGERT, The Effect of Virtual Masse on the Basic
Equations for Unsteady One-Dimensional Heterogeneous Flows, Int. J. Multiphase
Flow, 11, No 2, 149-160, 1985.

54
[29] G. B. WALLIS, One-dimensional Two-Phase Flow, McGraw Hill, New York, 1969.
[30] T. WATANABE, M. HIRANO, F. TANABE & H. KAMO, The Effect of Virtual Masse
Force Term on the Numerical Stability and Efficiency of System Calculations, Nuclear
Engineering and Design, 120, 181-192, 1990.
[31] T. WATANABE & Y. KUKITA, The Effect of Virtual Masse Term on the Stability of
the Two-Fluid Model against Perturbations, Nuclear Engineering and Design, 135,
327-340, 1992.

55
Chapitre 3

UNE METHODE NUMERIQUE

DECENTREE POUR LA
RESOLUTION
D'ECOULEMENTS
DIPHASIQUES SANS TERME
D'ECHANGE

57
 A NUMERICAL METHOD USING UPWIND SCHEMES FOR THE
RESOLUTION OF TWO-PHASE FLOWS WITHOUT EXCHANGE
TERMS1

3.1 Introduction

An accurate prediction of two-phase flow phenomena is essential to safety analysis of nu-

clear reactors under off-normal or accident conditions. In general, the ability to predict
these phenomena depends on the availability of mathematical models and experimental
correlations. Among several two-phase flow models, there are two fundamentally different
formulations of the macroscopic field equations for two phase flow systems; namely the
two-fluid model and the mixture model (see [11], [22]). Here we are interested in the two-
fluid model. It is considered to be appropriate for the most general and detailed description
of transient two-phase flows. In the two-fluid model, each phase is separately described
in terms of two sets of conservation equations. The interaction terms between two phases
appear in the basic equations as transfer terms across the interphases. Because the number
of necessary closure relations are considerably higher for the two-fluid model than for the
mixture model, much more detailed experimental data are necessary to develop satisfac-
tory closure relations. However, sufficient experimental information for developing reliable
closure relations for the two-fluid model is often not available. Therefore, the present state
of the art in the two-phase flow instrumentation implies that considerable uncertainties
exist in the closure relations for the two-fluid model. In spite of these shortcomings, how-
ever, there are a number of situations where the two-fluid model is more advantageous (or
even necessary in some cases) than the mixture model (see [11]). Also, many authors agree
on the basic form of the model (equations (3.1), below), although it lacks some physical
properties. The derivation of additional terms has been widely studied.
A specific concern is that most models presently used in the large computer codes are
based on governing equations having complex eigenvalues (Wallis model [27]), but com-
plex caracteristics are in general thought to cause ill-posedness of the numerical solutions.
In order to make their models hyperbolic, some researchers introduce formulations for
interfacial coupling between the two separated momentum equations, including space and
time derivatives of phasic physical quantities (virtual mass term, interfacial pressure; see
chapter 2 and the references therein, [19], [24], . . . ) . Hence, the correct formulation of the
basic two-fluid equations and the appropriate closure relations have been discussed during
the past and, up to now, there does not exist a commonly agreed approach.
'Ce travail a été réalisé en collaboration avec F. Coquel, E. Godlewski, B. Perthame et P. Rascle.

59
Our objective in this chapter is to develop a simple numerical method capable of predicting
two-phase flow phenomena with reasonable accuracy, computational efficiency and, follow-
ing the persistent uncertainties in two-phase modeling, the implementation of new terms
may be accomplished without changing the basic solution method. That is the reason why
we only treat here the simplified system (3.1). It is solved numerically by using an upwind
numerical scheme based on a finite volume method. The unknowns are the conservative
variables and all the quantities are computed at the centers of the cells. Coupling terms
between the phases are being addressed regarding their numerical treatment within the
frame of the present method, (phase exchanges, acceleration, friction terms and also real
fluids) see chapter 5.
We propose a resolution approach for the two-fluid model which reduces the problem to
classical gas dynamics. Then we use a Roe scheme and a Boltzmann scheme involving
compactly supported equilibrium functions. The Boltzmann scheme was introduced at
first by Perthame [15] for the resolution of compressible Euler equations. This approach
ensures the positivity of the pressure and of the densities of each phase as well as the
control of the void fraction which must stay within the interval [0,1]. The conservation
of total mass, momentum and energy is also preserved by this method. Several numerical
tests assert the possibility of computing stiff phenomena with our algorithm. For example
we can compute a separation problem where fronts propagate and one of the two phases
disappears locally. For further methods on diphasic computations, we also refer to [21],
[23] and the recent extension of the work in [7].
The format of this chapter is as follows. In the second section we present the continuous
model. Then, we describe our numerical approach. In section 4, some stability properties
are proved and the last section is devoted to the illustation of the schemes on well-known
test cases.

3.2 Two-fluid model of two-phase flow

The two-fluid model is formulated by considering each phase separately in terms of two
sets of conservation equations governing the balance of mass, momentum and energy of
each phase. Since the macroscopic fields of one phase are not independent of the other
phase, the interaction terms which couple the transport of mass, momentum and energy
of each phase across the interphase appear in the field equations. The balance equations
are complemented by the state equations for the two phases and by additional correlations
for the right-hand side coupling terms. The equations contain more dependent variables
compared with the number of equations. Therefore, additionnai assumptions are necessary
in order to close the system of equations. The most common procedure (to close the system
of equations) is to postulate equal local pressure values for the two phases.

The simplest macroscopic balance equations for the equal pressure one-dimensional
two-fluid model can be written under the following form, which can be eventually supple-
mented with specific terms for practical applications,

60
 dx(avpviiv) = 0,

{ ) dx(aiptui) = 0,

dt(avpvuv) + dx(avpvul + avp) - pdxav = (pi - p)dxav + avpvg,

+ dx{aiptuf + aip) - pdxai = (p, - p)dxat + atpig, (3.1)

0 9
It D tL
dt(avpv(ev + -—-)) + dx(avpv{ev -\ (- -~)uv) + pdtav = avpvguv,
I pv A
9 1
11} D Ul
dt{aipi(ei + -—)) + dx{atpi{ei H h —)ut)+pdtai =
z pi z

(3.2)

This system of equations represents the balance equations for mass (3.1.a, 3.1.6), momen-
tum (3.1.C, 3.1.oQ and energy (3.1.e,3.1./) of the vapor and liquid phases. The variables
appearing in the above equations have the following meanings (here we set k = v, I for the
vapor and liquid phases):

ait — volume fraction of Ar-phase,

pk = density of fc-phase,
Uk = velocity of fc-phase,
et = specific internal energy of fc-phase,
p = common pressure to the two phases,
Pi — interface pressure,
g = gravity constant.

These equations have to be supplemented by two equations of state. The standard

form of the equation of state for each phase is given by a function relating the density to
the pressure and the internal energy

€k)- (3.3)

In general, there exist no explicit expressions for the relationships in equations (3.3).
In practice, they are approached by various explicit forms. But the most general relation
is simply a tabular form. From thermodynamic principles, one can use any pair of thermo-
dynamic variables, for example: {pk,hk) or (pk,Sk), where hk and s* are respectively the
specific enthalpy and the entropy of A-phase. Also, in the system of equations (3.1), the
density for each phase is one of the conservative variables. Therefore, it is computationally
easier to deal with relations under the form

€k)- (3.4)

61
Since we suppose that the two phases are in pressure equilibrium, we have

Pv = Pi = P,

i.e.
P = Pv{Pv,ev)=pi{pi,ei). (3.5)

The interfacial pressure pt- needs a closure equation, the simplest of which is p,- = p,
see (3.37) below for another example of closure relation and references. We now have
six equations for the six unknowns {c*kPk> &kPkUk, &kPk(ek + u|/2)) with k = v(vapor),
k = l(liguid) completed with the three algebraic relations (3.2) and (3.5) for a* and p.

Notice that many simplifications are encountered in practice (see [22], Wulff[ll]
pp. 85-238, ...). For instance Sainsaulieu [19], Toumi [24] assume that the liquid phase is
incompressible with constant density pi, while the vapor phase is assumed to be an ideal
gas governed by the state function

p=(f-l)pvev, (3.6)

with a constant coefficient (1 < y < 3).

For simplicity in exposition here, we assume that both phases are described by an
equation of state of ideal gas type

€k
= frv

where yk is a constant and represents the ratio of specific heat capacities of the fluid k.
Typical values for 7 are 7 = 1.4 for a diatomic gas, e.g. air, and 7 = 1.0005 for the liquid.
This is reasonable at least, in a first step, when the vapor phase is always present and
imposes its compressibility to the liquid phase. Researches are in progress to improve these
laws.

It is very difficult to solve numerically the system of equations (3.1) as it stands.

Indeed, for pi = p it is not a hyperbolic system, and even for physically relevant expres-
sions of (p — pi) it is only conditionally hyperbolic. Hence instabilities can occur in some
situations. Additional physical effects, that are not modelized here, counter balance these
instabilities. But their explicit form is not settled precisely and they can be added in each
specific case. This is what we do in specific applications. But here we restrict ourselves to
the common model (3.1). Nevertheless, the mathematical structure of this system gives
some fundamental physical properties that we would like to keep at the numerical level. As
far as stability is concerned we wish to keep nonnegative densities, temperatures, pressure,
and the void fraction (i.e. a = av) within the interval [0,1]. Eventhough we do not use it
later, notice that this system is also endowed with two entropy inequalities (one for each
phase). In the case p, = p, and for a 7-law (3.7) these entropy inequalities write,

dtakPkSk + dxakPkukSk < 0, Sk = ln^-. (3.8)

P
62
We expect that they play an essential role in the stability. But at this moment, we were
not able to recover the corresponding entropy-in-cell inequalities (see [8], [10], [14]).

Concerning accuracy, a basic property we need is to be able to recover the classical

finite-volume schemes when only one phase is present (av = 0 or a; = 0). In particular
we wish to treat the most difficult case where one of the phases disappears (see the first
numerical test in section 5 ), i.e. av = 0 or a; = 0 locally.

3.3 Numerical Method

3.3.1 A two-step approach

We develop a new method based on a decomposition of the system. This method allows to
apply advanced numerical techniques that have been recently developed and successfully
applied to single phase flows (or compressible Euler equations). Let us briefly outline the
main steps of the method. Being given the state variables at time tn, (an,pn,p%,u%), (or
the computational variables, (ro£,ra£u£,Ek)), we proceed as follows.

Step 1 (hydrodynamical step), i) we solve two classical and decoupled hydrodynam-

ical systems using usual approximate Riemann solvers and ii) include the source term
(pidxa). Each of these two problems in step 1 are close to a quasi-lD nozzle problem
where a is similar to the section of the nozzle. After the first step, we know the quantities
n+l/2 n+1/2 n+1/2 n+1/2 n+1/2 r-.n+l/2 , ., e J.\. \.
mk ' = ak ' pk ' , mk ' uk ' , Ek ' and the pressure of the two phases are
different in general.
Step 2 (restoring equality of pressures). We compute akl+1,pn+1, />£+1, e£ +1 in order
to restore the equality of the pressures between the two phases.

We now describe in more details these two steps and introduce related notations.
To simplify, we set p, = p, but later, the interfacial pressure (3.37) will be considered.

First step (hydrodynamical step)

In the first step, we solve the two following systems, one for each phase k (k = v,l):
dt{m.k) +dx{mkuk) = 0 ,

dt{mkuk) + dx{mkul+p"k) = Pkdxak + <*kPk9, (3.9)

. dt{Ek) + dx{{Ek + pk)uk)

where

63
 Pk = ockPk = (lk ~ l)mjfeCfc.

The two systems will be discretized explicitly and then, they may be solved inde-
pendently because these are uncoupled.
These two systems can be written in the following compact form

where the vector Uk is defined by Uj = [rrik, nikUk, Ek), k = v,l.

In what follows, the terms due to gravity forces will be ignored. In practice (see sectionô),
these are treated using an explicit scheme.

Second step (restoring equality of pressures).

In t h e s e c o n d s t e p , w e solve t h e t w o following c o u p l e d s y s t e m s for k = v,l

dfTTlk = 0,

dt{mkuk) = 0, (3.11)

dtEk = -pdtak.

Thus the method involves a splitting which reduces the resolution of the full system
(3.1) to classical hydrodynamics (3.9) with source terms, and a rather simple problem
(3.11).

3.3.2 Discretization of the first step

Next we explain how to solve the first step (3.10), i.e. the system of equations

1 Ji—I — H(U), (3.12)

where

U = (m, mu, E)T is the vector of conservative variables,

F(U) = {mu, mu2 + p,(E + p)u)T is the flux vector,
H(U) — (0,pdxa,0)T is the source term.

and

64
 £ = m(e + y ) ,
p = (7 - l)me.

We note that the system of equations (3.12) is composed by two physically different parts.
The first one (the left-hand side) is in conservative form, but the second (the right-hand
side) is in non-conservation form. At this level we would like to mention that contact
discontinuities can occur but the equations stay meaningful; and for the applications we
have in mind, shock waves as in classical compressible flows are not of primary interest.

We will use a classical finite volume approach and discretize (3.12) in

j (3 ^

In order to explain this, we begin by dividing the spatial domain into N cells denoted by
Kj =]xj-i/2,Xj+i/2[ called control volumes. By choosing a mesh width Ax, a time step
At and a = At/Ax we have
*;-i/2 =(j-l/2)Ax, j<EZ,
tn = nAt, n€ IN.
For simplicity we take an uniform mesh, with Ax constant, although the method we discuss
can be extended to variable meshes.

Classically the finite volume method produces approximations t/j1 of the mean value
of U on the cell [^-1/21^+1/2] a t

We assume that, at time tn, we know £/", an approximation to the average of the
solution to (3.12) on the cell ]xj-i/2,xj+i/2b a n d we wish to calculate £/"+1, the mean
value of an approximate solution at time tn+1 = tn + At. We deduce (3.13) from (3.12)
setting

; ~ - L r + 1 / 2 f/(x,i")dx, (3.14)
Ax JJX] _ll2
X]_ll2

- ^ / r F(U(xJ+1/2,t))dt, (3.15)

H (3.16)
?-7^7T / H(U(Xit))dxdt.
J
At Ax Jt» Jx,-i/2

We now detail these approximations. We begin by the fluxes F?+1 ,T As it is well

known, the numerical fluxes F?+1,2 can be obtained through approximate Riemann solvers

65
(see Harten, Lax and Van Leer [10]). In practice, we have tested Roe solver [18] and a
kinetic solver [15] (see appendix for recalling). But for tests where one of the volume
fractions can vanish, the kinetic method is prefered because of its better behavior (see also
the positivity results stated in theorem 3.1). By contrast the Roe solver is known to be
instable close to vacuum which means here a = 0 or 1.

In the case of three-point schemes, these numerical fluxes have the form

(3-17)

And for the sake of completeness, we describe the second scheme namely the kinetic
scheme.
Consider the left-hand side of the system (3.12), i.e.

In order to describe the kinetic solver, we introduce the "equilibrium functions"

associated with the state vector Un(x) which are defined by:

(3.19)

•«(*, v) = A m » ( * ) . / f - ^ X [ (v - « - ( x ) ) i / ^ ^ I , (3.20)

with

1S a
Here x nonnegative function that satisfies

f (1, w, w2)X{w)dw = (1,0,1). (3.22)

JR.

The simplest choice (we will use it in practice) is due to [15]; it reads

X{W) =
2^1 1 H<^- (3 23)
'
But other choices are possible such as the physical Maxwellian l/\/2¥exp(—w2/2). We
refer to [8], [10] for more details on the motivations for introducing this method and other
references.
We notice that the integration of the two functions defined above with respect to the v
variable yields

66
 r{x,v)dv, (3.24)
n n
m u (x)= f vfn(x,v)dv, (3.25)
JIR
En(x) = Jm {-• /"(x, «,) + gn(x, v^j dv, (3.26)
and, instead of solving the system of equations (3.18) we solve the two following linear
transport equations

[
and I (3.27)
n n
I g(x,v,t )=g (x,v)

with t > tn, (x,v) € IR2.

The following lemma is at the basis of the Boltzmann schemes (Perthame[15]).

Lemma 3.1 Let At be a small time step and t G [£n, tn + At]. The following quantities

m(x,t) = f(x,v,t)dv,
J IK
mu{x,t) - j vf{x,v,t)dv, (3.28)
JIB.
E{x,t) =

are first order in At approximations of the solutions to (3.18).

Next assume that Un(x) is piecewise constant

Un(x) = Uf = (m], m]u], Ef), forz e]xj_i/ 2 , Xj+i/2[- (3.29)

In the finite volume approach, once the numerical flux functions ^+1/2 in (3.15)
which approximate the physical fluxes at the interface between the cells j and j + 1 are
known, the formula for updating the cell averages of the conserved quantities corresponding
to (3.18) is

Assume the following CFL condition (for the choice 3.23):

At < Ax/{\un{x)\ + v / 3y / p n (x)/m"(x)), Vz e B. (3.31)

67
By using (3.28) and the condition (3.31), the integration of (3.18) gives the flux F?+l/2
(a three entries vector) which is written

(3.32)

with

mjx(w)dw. (3.33)

Here we have generically denoted

U= (m,mu,E)T. (3.34)

The flux F is obtained by replacing " > " by " < " in the above integrals. As it is
proved in [15], this scheme is consistent for any function x which satisfies (3.22).

We now describe the discretization of the term H(U) = (Q,pda/dx,Q)T. This term
plays an essential role in the stability of the global scheme. We will treat it wisely and this
following the nature of physical phenomena and the characteristics of the conservation
law system (i.e. hyperbolic (p, =p£ p) or not (p, = p)). We could have used a two fractional
step-method which consists in integrating the left part of the system (3.9) in a first step
and the right one in a second step. This method has been used for the nozzle problem by
several authors. For reasons we will see later, we have prefered to integrate the two parts
in a single step.

We propose to discretize the term dxa in the expression pdxa (or pidxa when p, ^ p)
as

dxa = a(l - a)dxln{ (3.35)

1—a
Indeed we have observed that, when a (a — av, 1 — a = a\) is very small (say 10 6 ),
all the terms in the vapor phase equations almost vanish, but not dxa and this creates
instabilities. This also appears for a close to 1 for the liquid phase equations. With the
interpretation (3.35) of pdxa (or pidxa) we get rid of this. Next, it remains to discretize
(3.35). The physical quantities appearing in (3.35) will be taken at time tn.
Let
(3.36)

68
 We will discuss here two cases.

J± Pi # P-
In this case, pdxa is replaced by pidxa = (p - (p - Pi))dxa in (3.9) (see (3.1)). The
expression of (p — pi) ([4]) is given by:

P - Pi = £Pc{uv
2
-UI) , (3.37)

£ is a constant and pc is the continuous phase density. Other expressions can be found in
the literature (see [19], [24] . . . ) . Here we have used £ = 1 and pc = pv. Notice that, in
practical tests, the resulting expression is very small numerically.

Case 2: p,- = p: the space discretization is different in the two cases.

In case 1, we use the centered formula

(3.38)

But in case 2, we use the more stable formula

{pdxa)j = Pjctj(l - a^— (3.39)

where the Minmod function is defined by

Minmod(a, b) = -(sgn(a) + sgn(b))min(\a\, \b\).

In the case 1, we could have used the Minmod instead of the centered formula (3.38).
We prefer the centered formula which gives sharper fronts. This is possible because the
formula (3.37) brings more stability (hyperbolicity) to the continuous system compared
to the case p, = p.

Finally, the quantities computed from this first step are mk and Ek, {k = v,l).
This quantities will be identified using the superscript n + 1/2.

3.3.3 Discretization of the second step

Next, it remains to discretize the two following coupled systems:

dtTnv = 0, — 0,

dtmvuv = 0, and dtmiui = 0, (3.40)

dtEv = -pdtav, dtEi — -pdtai-

69
 This step can be performed cell by cell and thus we skip the index j through out
this subsection. To ensure stability, we will use an implicit method. First, we remark that
the first two equations in each system can be omitted because they only mean that the
quantities they involve are time invariant, i.e.
Tl
mn+l _ m + 1 /2 C n+l
and I (3.41)
{mvuv)n+l ={mvuv)n+l/2, {
Thus we are reduced to studying
dtEv = -pdtav,
(3.42)
dtEi = -pdta,.
We make some further simplifications before stating the ultimate system to be discretized.
Recall that

Ek=(j?pl+ k= L

By using the equalities (3.41), the system (3.42) amounts to

dt(mvev) = -pdtav,
(3.43)
dt(miei) = -pdtan.
Adding up the two equations of (3.43) each other (in view of (3.2)) gives
dt{mvev) + dt(miei) = 0. (3.44)
To keep this fundamental conservation law, we discretize it as

(mvev)n+1 + (m,e ; ) n + 1 = (mvev)n+1/2 + (m,e,) n + 1 / 2 . (3.45)

Next, using the relationship

p'k = {jk - l)mk€k, k = v,l, (3.46)

the terms in the left hand-side of (3.45) i.e. (mtefc)""1"1 for k = v, /, are replaced by

An easy computation shows that

(mvev)n+l/2
a "-r* 1 —a
7V - 1• + 7/ - 1
where a = av and (1 — a) = aj.
The unknowns are p and a. To compute them, we need to state another equation; thus,
consider the discretized form of the first equation of (3.43)
( m v c ) n + 1 - (mvev)n+1'2 = -pn+1(an+1 - à). (3.48)

70
The value à should be chosen at time fn+1/2 i.e. à = o n + 1 / 2 and then the system (3.1)
would be solved in two fractional steps. Now, a n + 1 / 2 is unknown from the first step, since
the computed variable is m" = an+1^2p" and from this relationship we cannot
n+1 2
deduce a / . For this reason, we set â = a " and as a consequence, the system (3.1) is
solved in a single step. (To check this, add the discretized systems of (3.9) (for k = v, I)
to those of (3.11) (for k = u, /)). We now have

(mvev)n+1 ~ (mvev)n+1'2 = -pn+1(an+l - an). (3.49)

Using (3.46), we get after some algebric manipulations

= m n+l/2 e n + l/2 (2z2??L - tt«) . (3.50)

\ (fv — -U /
Finally, replacing pn+1 in (3.47) by its expression in (3.50) yields an equation with the
unique unknown an+1. Some tedious calculations give the explicit expression

m e + ^

Notice that though it does not seem to be, this expression is indeed symmetric with
respect to the two phases. Next the explicit expression of p n + 1 is obtained plugging an+1
given by (3.51) into (3.47).

We are now in a position to update all the main variables involved in the system
of equations i.e.: a, p, mv, mvuv, Ev, mi, m/uj, Ei and the variables involved in their ex-
pressions i.e.: (pk, uk, ek k = v,l).

3.4 Properties of the method when based on kinetic sche-

mes

In this section, we give several theoretical properties of the method based on the kinetic
scheme in the case 2, where the Minmod discretization (3.39) of the term pdxa was intro-
duced to make up the lack of hyperbolicity. The hyperbolic case (case 1) is not addressed
here but strong numerical evidences assessed its stability properties.
This method preserves several basic properties of the initial system. By using the finite-
volume method and the conservative variables for the computed unknowns, we have con-
servation of total mass, momentum and energy. We strengthen at first (3.31) setting

At < Ai/2(|^| + VZ^/m]). (3.52)

The weak stability result is given by the following:

71
Theorem 3.1 Suppose that at the time level tn, we have p% > 0 for k — v,l; pn > 0 and
0 < al < 1. Then, we get p%+l > 0, pn+1 > 0 and 0 < a £ + 1 < 1 under the CFL conditions
(3.52) and (3.53) (see the notations in the proof below)

< \, (3.53)

which is also implied by the condition acting only at time level n:

At :
——rnaxjk-T
3i
Ax >>

Remark 3.1 1. These estimates are crude, but give a time-step of the order of the usual
CFL condition. They guarantee that a remains in ]0,1[. The time-step could become very
small (as in usual monovhasic gas dynamics) but in practice this does not occur as we
shall see even in the sej .ration problcn (section 3.5 pb.l,).

2. To the best of our knowledge these positivity results are new and are directly
inherited from the splitting method we have proposed.

3. A0J is given by (see 3.39)

Proof of Theorem 3.1

We consider the phase associated, for instance, with the void fraction a. We know
that f(x, v, tn) and g(x, v, tn) are nonnegative functions. Since they solve transport equa-
tions we also have f(x, v, tn+l>2) > 0 and g(x,v,tn+1/2) > 0. And thus:

= XI L K*>u' t n + 1 / 2 ) d v d x * 0. (3.56)
We write

^+1/2)AlSJR.K v x
^^
n
7T L
Ax JIRxKj
9(x,v,t

72
Now, following (3.13) and (3.28) we have

U/(X tn+1/2)dXdV
•h L K 3 '"' = mT'2-r12 -|P"^, (3-58)
where p"(Aa"/Ax) is the discretization associated with pdxa. This yields

(3 59)
-

The condition (3.52) implies that \v\^ < \ and after integration of /(z, r, f n+1 / 2 ) =
fn(x - vAt,v) and g(x, v, £n+1/2) = gn(x - vAt, v) over iR x Kj we obtain after some
tedious but easy calculations

m,- > \mn3 > 0, (3.60)

and furthermore

JIRxK} *
Therefore, using (3.61) we get from (3.59):

(3.62)

Notice that pdxa has been transformed into pa(l — a)dxp (see 3.35), with /? = ln(-~;)-
Denoting A/3"/Ax the discretization associated with 5x/3, the inequality (3.62) reduces to

In practice, the velocity |u" + 1 ' 2 | is a finite quantity and the condition (3.53) is sufficient
to ensure (see (3.59))

e j + 1 / 2 > 0. (3.64)

Now, we have proved that the variables m£ + ' and e£+ ' are nonnegative.
Since 7; > 1, yv > 1 and 0 < a" < 1, we have

0< _J_ + Q" < ~^— and 0 < an < - I I - , (3.65)

7; - 1 7; - 1 7V - 1
this implies that

0 L + ^ ^ ^
a ) m e + < m e +
+1/2 l/2
Jh-m; e^ . (3.66)
Tv ~ 1

73
Using (3.51), it follows that
0 < a n + 1 < 1. (3.67)
Following (3.67) and (3.47) we get
p n + 1 > 0. (3.68)
Finally, (3.41) and (3.67) give
pnk+1>0 {k = v,l). (3.69)
Hence, we have proved the result under the condition (3.53).
Now, we will attempt to upperbound |u" + ' | in order to recover (3.54). We argue as
follows

+ / v2f,(v)dv- I «Vi+i(»)*>]• (3.70)

Jv<0 Jv<0
Following (3.58) and (3.70), we obtain

/ (i() i + (t;))dt;]-~^A^, (3.71)

v<0 iAX

and we have only to consider the case where uj and Act7} have the same sign, say
nonnegative. Then
n+l/2un+l/2 <m>n + to{f V7f._i{v)dv+ f ff{v)dv]i (3.73)
i i l Jv>0 Jv<0

and using (3.60) we have

^x/~), (3.73)

with
Mf =
m
— j U " j l •*"ej) •+• "* a a ; jV"*j-lU u j-ll "+•ej-l)i m
j + u l " j + i l "t"e j+iJM > i '' 4 J
In this formula, the "max" comes from the necessity to consider also the other possibility
when u" and A Q J are negative.
We can rewrite (3.73) as

with

M? = ^ . (3.76)

74
3.5 Numerical Results

Problem 1: Sedimentation

This problem represents a very simplified physical phenomena: separation of air and
water by gravity in a vertical tube.
The objective of this problem is to check the capability of the model and of the chosen
numerical approach to describe counter-current flow conditions and the occurence of strong
void gradients which are typical of many situations where phase separation phenomena
are dominating.
The test consists of a vertical tube of 7.5m length. In the initial state, the tube is filled
with a homogeneous two-phase mixture (air, water) of constant void fraction of a = 0.5;
the velocities are uv = u; = Oro/s. The boundary conditions at the inlet and the outlet
are nul debits for the vapor and the liquid, i.e., the velocities of the two fluids are forced
to equal zero at the two ends of the duct, since these are closed.
The phase separation is due to gravity forces only. The liquid phase falls down to the
bottom of the tube, while the gas phase comes up to the top of the tube.

All the presented results have been achieved using a constant ratio At/Ax = 6.OJ57 —
4. Figures 3.2 and 3.4 below show in the space-time plane the predicted values for the
void fraction in case 1 and case 2 (i.e. p, ^ p and p, = p). The comparison between the
steady-state solutions for the two cases is given in figure 3.5, this figure clearly indicates
that the steady-state solution in case 1 is much closer to the expected analytical steady-
state solution than in case 2.
The void fraction distribution along the vertical height at different time levels is shown
in figure 3.3, this figure indicates an upward and downward directed void front in the
transient period; the occurence of quasi-stationary conditions is achieved once the two
fronts have merged at the middle section of the tube.
Figure 3.6 shows the spatial convergence of the void fraction at the steady-state. We give
below the errors in L1 norm, choosing the finest grid solution as the reference one.

errori errori 9
0.2090 0.1069 0.9662

where errori = \\a25 — aioolli an<* error2 = H050 — Oioolli, the subscripts represent the
number of cells. Here 6 = •— * —— provides an estimate of the order of accuracy
Iog{h1/h2)
(hi = 7.5/25 and h2 = 7.5/50).
Figure 3.7 shows the predicted values for the pressure in space-time plane in case 2. One
can remark the establishement of a static pressure gradient in the gas and liquid regions.
Figure 3.8 compares the steady pressure profiles in the cases 1 and 2. Not only the case 1,
because of the centered discretization, gives an overshoot as can be seen in figure 3.8, but
also it is very sensitive to the number of computationl cells. These results can be compared
to those obtained by [21].
75
 Constant Velocity Liquid Inflow

Ambient Pressure Ambient Pressure Ambient Pressure

t t = tf

Liquid Air

Figure 3.1: Water Faucet Problem Evolution

Problem 2: Water Faucet

The water faucet problem (see [16]) consists of a vertical tube 12m in length and \m
in diameter. The top has a fixed inflow rate of water at a velocity of 10m/s, temperature
of 50°c, and a liquid volume fraction of 0.8 . The bottom of the tube is open to the
ambiant pressure and the top of the tube is closed to vapor flow. The problem is illustrated
schematically in figure 3.1.

Initially, the tube is filled with a uniform column of water at a velocity of lOm/s sur-
rounded by stagnant vapor, such that the vapor volume fraction is 0.2. The thermodynamic
properties of the system at the initial state are assumed constant at values appropriate for
air-water mixture and are 50° c for the temperature and 1.0E05(Pascal) for the pressure.
The velocity boundary conditions at the inlet are lOm/s for the liquid and 0.0m/s for the
vapor. The volume fraction of the vapor at the inlet is set constant at 0.2. The only out-
flow boundary condition at the bottom of the tube is constant pressure at 1.0.E05(Pascal).
These specific boundary conditions are implemented using the solution of a partial Rie-

76
mann problem which is described in details in chapter 6. Here, for physical considerations
(see[16]), we are led to choose p, = p. The calculation is carried out until a steady flow is
attained.
This transient problem has a particularly simple analytical solution when pressure vari-
ation in the vapor (and hence in the liquid) is ignored (see [16]), this analytical solution
was used as a code test problem in [17].
The objective of this problem is to test the stability and the convergence of the numerical
solution method. The diffusive character of the numerical method is also tested since a
discontinuity in void fraction is propagated through the solution space.

All the computations have been achieved using a constant ratio At/Ax = 5.0E — 4.
The void fraction as a function of space at different times is shown in figure 3.9. This
figure shows the dynamic propagation of the void profile down the tube until the wave has
completely passed out of the tube and the steady-state profile remains.
In order to check the convergence and the stability of the scheme, computations have been
made using a discretization with 12, 24, 48, 100 and 200 cells. The void fraction profiles
for various discretisations are compared to the above mentionned analytical solution (see
figure 3.10), which is explicitly given by:
Ofttf
1- if
(3.77)

0.2 otherwise,
where uf = lOm/s and a* = 0.8.
The L1 norm errors between the approximate and analytical solutions are displayed iin the
following table:

errori error2 errors error4 error5 e

0.4761 0.3796 0.2793 0.2141 0.1674 0.3717

with errori — ||a,-a ona /j,ti co /|| and % refers to the various discretizations (i — 1 for 12 cells,
x=
i = 2 for 24 and so on). Here 6 denotes an average value of the 0, = -*- —

We compare in the figure 3.11 the profile given by the kinetic and Roe solvers in the
first step. As announced, the Roe solver is more accurate.

Finally, the transient liquid velocity profile is given in figure 3.12.

The numerical results obtained when the Roe scheme is used in the first step are displayed
in figures 3.13, 3.14 and 3.15. One can remark that, because of the nonhyperbolicity of the
problem, when the mesh is refined, ascillations occur for the two schemes. This oscillations
appear later for the kinetic scheme because it is more diffusive.

Problem 3: Shock tube problem

In order to show that the algorithm is able to handle large differences between liquid and

77
gas (vapor) temperatures, we present a two-phase shock tube problem. The initial data is
given by the following table:

a p (bar) uv (m/s) ui (m/s) K [kj/kg) hi (kj/kg) 7v 7/

L
U 0.25 200. 0. 0. 3092.7 1338.2 1.0924 1.0182
vk 0.25 150. 0. 0. 3092.7 2338.2 1.0924 1.0182

Notice however that the mathematical meaning of this problem witn shocks is quite
imprecise (see [5]), and there is no physical interpretation for this case because the thermo-
dynamical desequilibrium requires exchange terms between phases. The computations have
been achieved using 200 cells and a constant ratio At/A.. = 5.0.E — 4. In the figures 3.16,
3.17 and 3.18, we show the pressure and enthalpies profiles at time t = 0.75E — 03.

void fraction

distance (m)

time(s)

Figure 3.2: Void fraction in space-time plane (case 1), 50 cells.

78
 1 1— -i ^i •" _ — , — _ . t

i
I /
,y
_..-••
t = 0.0 s
t - 0.5 s
0.9 • .' _ . / • ' t = 1.0s "
i _. t= 2.5S -.—-
t = 10.^-'
0.8

0.7

0.6 -

s
0.5
/
/
a
i
/
i
0.4 i
fi

/
0.3

02 - / i
i
j
0.1

••'•••••"

0
2 3 4
distance (m)

Figure 3.3: Void fraction history, 50 cells.

void fraction

distance (m)

time (s)

Figure 3.4: Void fraction in space-time plane (case 2), 50 cells.

79
 1
/ / easel
; / case 2
0.9

0.8

0.7 -

0.6

0.5

0.4

0.3

02

0.1

0
3 4
distance (m)

Figure 3.5: Comparison between steady-state solutions, t = 10 s, 50 cells.

26 cells -
50 cells --•
0.9 100 cells -

0.8

0.7

0.6

0.5

0.4

0.3

02

0.1

0
3 4
distance (m)

Figure 3.6: Spatial convergence of steady-state solution (case 1), t = 6.5 s.

80
 pressure (Pa)

distance (m)

time (s)

Figure 3.7: Pressure in space-time plane (case 2), 50 cells.

I 105000

3 4
distance (m)

Figure 3.8: Comparison between steady-state solutions, t = 10 s, 50 cells.

81
 0.6

0.55 • .-•-•i = o.3s -

W).5.s
t = 0.7s'-^~-...
t = 0.9s — -
0.5 • t = 1.1 s "
t=1.3s

0.45 •

S 0.4 •

ï
| 0.35

0.3 - • • - • , .

. . . . . .

0.25 • • • - - . . . . - \ ^

0.2

6 10 12
distance (m)

Figure 3.9: Transient void fraction profile, 24 cells.

0.45
analytic
12 cells
24 cells
48 cells —
0.4 100 cells — -

À
200 cells — -

0.35

0.3 •

0.25 • ./•

0.2 •

0.15
6 10 12
distance (m)

Figure 3.10: Convergence void fraction profile, t = 0.4 s.

82
 0.45
analytic
Roe scheme
Kinetic scheme

0.4 -

0.35

0.3 -
• / / \

\
0.25 • \
X
02 '•

0.15
6 10 12
distance (m)

Figure 3.11: Convergence void fraction profile, t = 0.4 s, 48 cella.

O.Os
t = 0.1 s
t = 0.3s —.,-•:
t = 0.5.3-'
t
-1=1.3 S

4 6
distance (m)

Figure 3.12: Transient liquid velocity profile, 96 cells.

83
 0.65
0.0 s —
0.1 s —
0.6 t 0.3.S-"'-'-
,t*-0'.5s —

0.55 t :
t=1.1 s —
t=1.3s —
0.5

0.45

S 0.4

0.35

0.3

0.25

0.15
10 12
distance (m)

Figure 3.13: Roe scheme. Transient void fraction profile, 24 cells.

0.45
12 cells
24 cells
48 cells
i 100 cells
0.4 200 cens
/ /
\ \
j--y-
••- \ \
0.35 / / / -

0.3

1
025

*• \ " ' . .

\ ^
.
02

0.15
6 8 10 12
distance (m)

Figure 3.14: Roe scheme. Convergence void fraction profile, t = 0.4 s.

84
 19
0.0s
0.1 s
18 0.3 5..- —
o:s-
0.7 s
17 1.3s — - -

16

15

14

13

12

11

10

9
6 10 12
distance (m)

Figure 3.15: Roe scheme. Transient liquid velocity profile, 96 cells.

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 3.16: Pressure profile.

85
3.13e+06

3.126+06 -

3.05e+06 -

3.046+06
0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9
distance (m)

Figure 3.17: Vapor enthalpy profile.

2.46+06

2.26+06

1-2e+O6
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 3.18: Liquid enthalpy profile.

86
Bibliographie

[1] F. COQUEL, K. EL AMINÉ, E. GODLEWSKI, B. PERTHAME k P. RASCLE, Numerical

method using upwind schemes for the resolution of two-phase flows, to apear in J.
Comput. Phys.
[2] F. CORON & B. PERTHAME, Numerical passage from kinetic to fluid equations, SIAM
J. Numer. Anal., 28, pp. 26-42, 1991.
[3] F. D E VUYST, Schémas non-conservatifs et schémas cinétiques pour la simulation
numérique d'écoulements hypersoniques non visqueux en déséquilibre thermochimique,
Thèse de l'Université Paris VI, france, 1994.
[4] K. E L AMINÉ & P. RASCLE, Ecoulement Diphasique, Modélisation et Hyperbolicité,
Rapport interne EDF HT-13/96/004, 1996.
[5] G. DAL MASO, P. LE FLOCH k F. MURÂT, Definition and weak stability of non-
conservative products, J. Math. Pures Appl., 74, 483-458, 1995.
[6] T. GALLOUËT, An introduction to finite volume methods, cours CEA-EDF-INRIA,
(France) (1992), to appear in finite volume methods by R. EYMARD, T. GALLOUËT
k R. HERBIN, Handbook of Numerical Analysis, P. G. CIARLET k J.-L. LIONS
(Eds.), North-Holland, Amsterdam.
[7] J.M. GHIDAGLIA, A. KUMBARO k G. LE COQ, Une méthode volumes finis à flux
caractéristiques pour la résolution numérique des systèmes hyperboliques de lois de
conservation, C. R. Acad. Sci. Paris, t. 322, série I, p. 981-988, 1996.
[8] E. GODLEWSKI k P.-A. RAVIART, Hyperbolic systems of conservation laws, Ellipses,
Mathématiques k applications, vol. 3/4, 1990.
[9] E. GODLEWSKI k P.-A. RAVIART, Numerical approximation of hyperbolic systems
of conservation laws, AMS 118, Springer-Verlag, New-York, 1996.
[10] A. HARTEN, P.D. LAX k B. VAN LEER, On Upstream Differencing and Godunov-
type Schemes for Hyperbolic Conservation Laws, SIAM Review, Vol 25, n° 1, pp.
35-61, 1983.
[11] G.F. HEWITT, J.M. DELHAYE k N. ZUBER, Multiphase Science and Technology,
Volume 5 Hemisphere Publishing Corporation, 1990.

87
[12] M. ISHII, Tkermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975.
[13] R.T. LAHBY & D.A. DREW, The Three Dimensional Time-and Volume Averaged
Conservation Equations of Two-phase Flow, Advances in Nuclear Science and Tech-
nology, Vol.20, 1-69, 1988.
[14] S. OsHER & E. TADMOR, On the convergence of difference approximations to scalar
conservation laws , Math. Comp., Vol 181, pp. 19-51, 1988.
[15] B. PERTHAME, Second-order Boltzmann schemes for compressible Euler equations
in one and two space dimenions , SIAM J. Numer. Anal., Vol.29, n° 1, pp. 1-19,
February 1992.
[16] V.H. RANSOM, Numerical Benchmark Tests, Multiphase Science and Technology,
Volume 3, Edited by G.F. HEWITT, J.M. DELHAY k N. ZUBER, Hemisphere, 1987.
[17] V.H. RANSOM & V. MOUSSEAU, convergence and Accuracy of the Relap5 Two-phase
Flow Model, Proceeding of the ANS International Topical Meeting on Advances in
Mathematics,Computations and Reactor Physics, Pittsburg, Pennsylvania , 1991.
[18] P.L. ROE, Approximate Riemann Solvers, Parameter Vectors, and Difference
Schemes, J. Comput. Phys., 43, 357-372, 1981.
[19] L. SAINSAULIEU, Contribution à la modélisation mathématique et numérique des
écoulements diphasiques constitués d'un nuage de particules dans un écoulement de
gaz, Thèse d'Habilitation à diriger des recherches, Université Paris VI , 1995.
[20] A.S.-L. SHIEH, R. KRISHNAMURTHY $ V.H. RANSOM, Stability, Accuracy, and Con-
vergence of the Numerical Methods in RELAP5/MOD3, Nuclear Science and Engi-
neering, 116, 227-244, 1994.
[21] H. STÂDTKE, G. FRANCHELLO & B. WORTH, Numerical Simulation of Multi-
Dimensional Two-Phase Flow Based on Hyperbolic Flow Equations, preprint, 30th
Meeting of the European Two-phase Flow Group Piacenza, June 6-8, 1994, Italy.
[22] H.B. STEWART & B. WENDROFF, Two-Phase Flow: Models and Methods, J. Comput.
Phys., 56, 363-409, 1984.
[23] I. ToUMI, A Weak Formulation of Roe's Approximate Riemann Solver, J. Comput.
Phys., 102, 360-373, 1992.
[24] I. TOUMI, An Upwind Numerical Method for a six Equation Two-Phase Flow Model,
Rapport DMT/94/168, SERMA/LETR/94/1619, C.E.A., 1994.
[25] J.A. TRAPP &c R.A. RlEMKE, A Nearly-Implicit Hydrodynamic Numerical Scheme
for Two-phase Flows, J. Comput. Phys., 66, 62-82, 1986.
[26] P. VlLLEDIEU, Approximation de type cinétique du système hyperbolique de la dy-
namique des gaz hors équilibre thermochimique, Thèse de l'Université Paul Sabatier
(Toulouse), france, 1994.
[27] G.B. WALLIS, One-dimensional Two-Phase Flow, McGraw Hill, New York, 1969.

88
Chapitre 4

UN SCHEMA NUMERIQUE
POUR LES ECOULEMENTS
MONO-PHASIQUES DE
FLUIDES REELS

89

PAÛE(S)
!®f t BLASitC
 A NUMERICAL SCHEME FOR ONE-PHASE FLOW OF REAL
FLUIDS

4.1 Introduction

The numerical schemes proposed for the computation of inviscid compressible gas dynam-
ics in which the properties of the fluid are represented by the ideal equation of state, have
been extended to the real gas case by several authors. The main difference with regard
to the special case of an ideal gas is that the ideal gas law is replaced by an equation of
state of the form p = p(p, e) or p — p(p,pe). Colella and Glaz [1] extended the numerical
procedure for obtaining exact Riemann solution to the real gas case, by introducing an
iterative solution procedure to handle shock waves. Dubois [4] developed an algorithm for
the computation of rarefaction-compression waves using tables of Srinivasan et al. [21].

The two approaches listed above are "exact" Riemann solvers. For the extension
of other approximate Riemann solvers and flux-vector splitting techniques, we refer the
reader to [5], [13] and [15].

Here, we will be particularly concerned with the development of an approximate

solution algorithm for the Riemann problem for the medium pure water with a general
equation of state. The algorithm is obviously viable for the case of a real gas (and also
for vapor, which may be considered as a real gas). Notice that, a one-phase flow of steam
or water is a particular case of two-phase flow when the void fraction is zero or one.
The first practical application of this algorithm is in the field of reactor nuclear safety,
for example, in the normal operation of a Pressurized Water Reactor (PWR), the coolant
fluid remains in general in the liquid state. Another important application is the numerical
study of underwater shock phenomena. Using the so-called Tait equation (see [2] or [10]),
which is a simplified equation characterizing the properties of the pure water, Holl [10] has
established that the Riemann problems for the water and the ideal gas can be solved using
identical methods. Based on this result, Cooke and Chen [2] adapted the Roe method to the
numerical calculation of shock wave propagation in liquid media with a general equation
of state, called the UNESCO equation of state. However, since a relation specifying the
dependence of internal energy upon the other state variables is not available, the governing
equations were formulated in terms of enthalpy. This formulation leads to the appearance
in the energy equation of a term involving the time derivative of pressure, which is treated
as a source term. Unfortunately, this leads to a Riemann problem for a system where the
energy equation is not a conservative law.

91
 In this chapter, we develop an algorithm for approximating the solution of the sys-
tem (4.1) for real fluid in its original conservation form. This algorithm is based on the
construction of a Boltzmann scheme involving compactly supported velocity distribution
functions. The proposed algorithm is positively conservative. This means that, the density
and internal energy remain positive throughout the computational process.

4.2 Governing equations

4.2.1 Dynamical equations

In fluid dynamics, the evolution of the fluid is governed by the equations of conservation
of mass, momentum, and energy. For an inviscid non-heat-conducting flow, the balance
equations in one-dimension can be written in the form

| | 0 , (4.1)

where the vectors U and F(U) are given by

U = (p,pu,E), (4.2)

F(U) = {pu,p + pu2,{E + p)u), (4.3)

where t is the time, x the distance, p the density, u the velocity, p the pressure and
E = p(t + u 2 /2) is the total energy (e is the specific internal energy).
In order to close the system of equations (4.1), the pressure p must be specified through
a constitutive equation, the EOS, which characterizes the fluid material. Notice that, the
dynamical behavior of a fluid is significantly affected by the properties of p (see [14]).

4.2.2 Equation of state

Real fluid

The thermodynamic properties of a fluid are given by an equilibrium equation of state

e = e(v,s), (4.4)

called the complete equation of state, expressing the internal energy as a function of the
specific volume v — \/p (p is the density) and the specific entropy s. Physical principles of
thtrmodynamics require the equation of state (4.4) to satisfy certain convexity constraints
ant; to have appropriate asymptotic properties. The constraints are expressed in terms of
various derivatives of e. For the fluids considered here (pure phases), the thermodynamic

92
constraints on the domain of definition and the asymptotic conditions are satisfied (see
[14] for more details).

In the following, we introduce a particularly important second derivative of the

internal energy which determine the speed of propagation of waves, namely, the adiabatic
exponent

The pressure and the temperature are defined by

and T(v,s) = £ (4.6)

OV\s OS\

7 can also be expressed as

Thermodynamic stability requires only that 7 > 0, but the strict hyperbolicity of (4.1)
requires strict positivity of 7 i.e.: 7 > 0 (see below). For many materials and in particular
for gas/steam and water, 7 > 1 (see figure (4.2)), away from phase transitions. Here, we
restrict ourselves to pure phases (gas/steam or water), for phase transitions, see chapter
5.

The model for fluid flow considered here, writen in the form (4.1) neglects viscosity,
heat conduction, and radiation. As a result, the dynamics require only partial specification
of the thermodynamics of the material, as contained in an incomplete equation of state.
In other words, only a relation of the form
p = p(v,e), (4.8)
is needed. From equation (4.8), one can also compute the adiabatic exponent. Indeed
d ., . d ,, . de d ., .

P
dv |e de
Therefore we get:
. _ v dp _ v (dp dp
7
~ ~"p"d~o\s ~~p \lh\e ~ P fc|,

Consider now the system of conservation laws (4.1), where the pressure is given by
the incomplete equation of state (4.8). If we set W = (p,u,e), the system of equations
(4.1) can be written in the following form
d

93
Expressed in nonconservation form, the system (4.10) reads

h'(W)$-W + f'(W)~W = 0. (4.11)

01 OX
The systems (4.10) and (4.11) are equivalent for smooth flows. Setting B = h'(W) and
A = f'{W), the characteristic polynomial resulting from det(B~1A — XI) = 0 is given by
(« - A)3 - (u - X)(v2{ppe - pv)) = 0, (4.12)
where

By definition, the sound speed c is given by

c2 = v2(ppe-pv). (4.14)
By setting v = l/p and from the relation settled in (4.9), we can derive the relation
c2 = _v2dj> = dp
dv\s dp

we find again the classical definition of the sound speed which says that: c2 is the derivative
of pressure with respect to density at constant entropy. From (4.7) and (4.15) we obtain

c2 = _J_|p =jp
p1 OV\s p
The adiabatic exponent 7 is now seen to be a dimensionless measure of the sound speed c.

The sound speed is real because 7 > 0 as required by thermodynamic stability.

Moreover, the sound speed is strictly positive, except where 7 = 0, which can happen
only in three-phase regions (i.e., triple points) and at special points in two-phase regions.
Away from such points, the system (4.1) is strictly hyperbolic and its three eigenvalues
are u — c,u,u + c.

Ideal gas

For an ideal gas, pv and e are independent of v at fixed T, according to the experiments
of Boyle and of Joule-Thompson. The internal energy is a function of temperature alone,
e = e(T), and pv is related to T by the ideal gas law,

pv = RT, (4.17)
where R is a constant. A polytropic gas is an ideal gas for which the specific heats at
constant volume and pressure, i.e., cv and cp are constant. In this case, we have
e = cvT, (4.18)

94
and

h = cpT, (4.19)

where h — e + pv is the specific enthalpy.

Note that, by the ideal gas law (4.17), we have

cp - c = R, (4.20)

and the equation of state takes the form

p=(7-l)ep, (4.21)

where 7 is the ratio of specific heats, 7 = cp/cv.

For a polytropic gas, the entropy s is defined up to an additive constant by

s = cJogiTv"1'1)) + const
= cvlog{—)1 + const.
P"
Thus, the pressure can be expressed as a function of entropy s and density p,

p=Kealc*p">, (4.22)

where K is some constant. In addition, the sound speed in the gas (polytropic gas) is given
by

| = 7AV /cV-! = IP. (4.23)

°P\s P

Remark 4.1 It should be emphasized, however, that the adiabatic exponent 7 is not the
ratio of specific heats 7 for a general equation of state.

Let us return to the real fluid case. In analogy with the special case of a polytropic
gas (4.21), we define a dimensionless quantity 7 by

7(p,e)=l + ^. (4.24)

Recall that, for polytropic gases, the ratio of specific heats 7 is equal to the adiabatic
exponent 7.

The quantities 7 defined in (4.24) and 7 defined in (4.16) were used by several authors
([1], [13], [15]) for the construction of approximate solution algorithms for the Riemann
problem for real gases. Here, these two quantities will be used also for the Riemann problem
for pure water.
In numerical calculations of fluid flow, the EOS can be specified by an analytical formula,
an empirical fit, or a table. As far as we are concerned, the thermodynamic properties of

95
the fluid are given by a table [19]. Notice that, in this table, the input thermodynamic
variables can be choosen only among the following three pairs (p, h), {p,T) and (p,s).
Before presenting our Riemann solver, we indicate that these two quantities satisfy

1<7<3 and (4.25)

on the set of regions in pressure and temperature which covers the physical range of
interest. The temperature and the pressure ranges are respectively 5 to 526°C and 0.00871
to 400 bar. Figures 4.1 and 4.2 show the distribution of the variables 7 and 7 in the p-T
plane (pressure-temperature) for the liquid region. We see that 7 can be very large.

4e+07
3.5e+07
3e+07
2.5e+07
100 1J (Pa)
200 1e+07
300 5e+06
Temperature (Celsius) 400
500

Figure 4.1: Variations of the quantity 7 in the liquid region (water).

4.3 Kinetic scheme

The kinetic scheme we will attempt to constuct here, is an extension of the Kaniel scheme
[11] (derived for isentropic ideal gas) to non-isentropic real fluids (for another extention,
see [24]).
In order to approximate the system of equations (4.1), we consider the linear kinetic equa-
tions (Boltzmann equations) with appropriate initial data, which modelize the transport
of equilibrium functions whose the moments are first order approximations in time of the
solutions of (4.1). This formalism is then used to build a first-order explicit scheme of Flux
Vector Splitting type (FVS).

96
 4e+07
3.5e+07
3e+07
2.56+07
2e+07
100 LSe+OPressure (Pa)
200
300 1e+07
5e+06
Temperature (Celsius) 400 500

Figure 4.2: Variations of the adiabatic exponent 7 in the liquid region (water).

4.3.1 Kinetic formalism

From a given state vector at time tn

Un(x) = (pn(x),pnun(x),En(x)),
and the nondimensional thermodynamic variables 7 and 7 such that

and
p"(x) = (7n(x)-l)pn(x)en(x),
we define the associated "equilibrium functions" as:

pn(x) v - un(x)
P(x,v) =
7n(a;) - 1 cn(x) cn{x)
(4.26)

gn{x,v) =

where
(4.27)

97
Note on the one hand that thanks to the bounds (4.25) // defined in (4.27) remains noneg-
ative. On the other hand, since / n (.,.) is compactly supported, for 7 > 1, /"(.,.) and its
moments are integrable and so is gn(.,.).
The support of these two functions is built in such a way that, the velocity v is bounded
by the propagation velocities of the physical sound waves, u — c < v < u + c. It is easy to
check that, the integration with respect to v of the functions defined above gives

Pn(x) = jmfn(x,v)dv, (4.28)

pnun(x) = fvr{x,v)dv, (4.29)

JIn

En(x) = JmhÇfn(x,v)+gn(x,v)\dv. (4.30)

To the initial data (4.26) we associate the two following linear transport equations

and I (4.31)
{ g(x,v,tn)=gn(x,v)

with t > tn, (x,v) e m2.

The exact solutions of the two initial value problems (4.31) are respectively:
f(x,v,t) = fn(x -vt, v) and g(x,v,t) = gn(x - vt,v). Clearly f(x, .,*), g(x, .,t) and their
moments are integrable. The principle of the method is based on the following results.

Lemma 4.1 Let At be a small time step and t € [tn,tn + At]. The following quantities

p{x,i) = f{x,v,t)dv,
J In
pu{x,t) = /' vf{x,v,t)dv, (4.32)
J In

E(x,t) -

are first order in At approximations of the solutions to (4-1)•

Proof

Let Û{x,t) = {p{x,t),p~u{x,t),Ë(x,t)) be the exact solution of the Cauchy problem

associated to the system of conservation laws (4.1):
o ~ o
Ft +dx
~ (4.33)

98
with t >tn and x 6 JR. Un(x) is assumed to be smooth. The first order accurate Taylor
expansions in time for p{x,t), pu(x,t) and E(x,t)) in the neighborhood of {x,tn) are:

p(x,t) = p{x,tn)-At—pu(x,tn)+O(At2),

pu(x,t) = ^

Ë(x,t) = E(x,tn)~ &t-j^

To prove lemma 4.1, it suffies to show that the time expansions of the approximate
solutions (4.32) differ from those of the exact solution by O(At2). Indeed we write

POM) = Jmf(x,v,t)dv

m Q

= p(x, tn) - At v-z-fix, v, tn)dv + O(At2)

JIR
im ox

and

pu(x,t) = / vf(x,v,t)dv
JIR

= pu(x, n - At J v2^f(x, v, tn)dv + O(At2)

And finally,

E(x,t) =

O{At2)

= E(x, tn) - At j m [Y^/(«, v, tn) + v^g(x, v, tn)] dv + O(At2),

99
after some tedious calculations we obtain

E(x,t) = E(x,tn)-At+
Ox
This completes the proof.

4.3.2 Construction of the numerical scheme

By using lemma 4.1 and the finite-volume method, we will build a first-order flux splitting
scheme for the system of conservation laws (4.1). The spacial domain is divided into cells
denoted by Kj =]z,-_i/2.#j+i/2[> where Xj+1/2 = {j + 1/2)Ax with j e 2Z, and Ax is the
mesh width. Assume that the initial state Un(x) — U(x,tn) is a piecewise constant over
the cells Kj i.e.,

lT(x) = Uf = (roj,m]u],EJ) for* €]Xj_1/2,xj+l/2[ (4.34)

We wish to compute the average denoted by C/"+1 of U(x,tn+l) over the cells Kj at time
tn+1 =tn + At,

/"+1 = ~ T J + I / 2 U(x, tn+1)dx. (4.35)

The expression for updating cell averages of the conserved quantities can be written as,

At
'~~ -J7-1/2). (4-36)

Lemma 4.2 Under the CFL condition

~ ( | < | + c?)<l, Vj€Z (4.37)

the finite volume scheme (4-36) is given by

U*+l) - F{U*.X,U?)), (4.38)

where the numerical flux function F(Uf t^j+i) is split as

F(U?,U?+i) = F+(U?) + F"(t/; + 1 ). (4.39)

The expressions of F + and F~ are given in (4-45)-(4-48)-

Proof

By integrating over ïïtx Kj x [tn,tn+1] we have

100
this gives
r
t-1/2

+ (0,0, 1)T(<7(*, v, tn+1) - g{x, v, tn))] dvdx

K"2. y)
+ (0,0, v)T{g(xj+1/2, v, t) - g(xj_1/2, v, t))]
Using lemma 4.1 and the exact solutions of (4.31) we get

Cj-i/2 - v(t ~ f "), v)) + (0,0, v)T{gn

-9n{xj.1/2-v(t-tn),v))]dvdt. (4.41)

Since the functions /"(.,.) and gn(.,.) given in (4.26) have a bounded support, the
u-integration is reduced therefore to those such that \v — u\ < c, i.e. |u|(t-t n ) <
and by the C.F.L. condition we have |u|(t— tn) < Ax i.e.
x
j±i/2 - Ax < xj±l/2 - v(t - tn) < xj±1/2 + Ax.
Thus the expressions of /"(.,.) and gn(.,.) involved in (4.41) reduce to

fn(xj+1/2-v(t-tn),v) =
fn(xj+1,v) ifv
and

gn{xj>v) if v > 0

gn{xj+1,v) ifv<0
where
3-i
V- UT
1 Pi (4.42)
T" — 1 c" v-un I

101
and

(4.43)

The functions / n (x ; _i/2 — v(t — tn),v) and 5n(Xj_i/2 - v(t — tn),v) are obtained in the
same way.

By setting

*(£?, v) = («, v2, ~)Tr{xh v) + (0,0, v)Tgn(Xj, v), (4.44)

the expression (4.41) takes the form

-I h(U^v)dv- f h(U?,v)dv]
J
Jv>0 Jv<0
At

For the numerical computations, we develop the expressions of F + and F by using

the change of variable w = (v — u)/c (we drop the indexes j and n for simplicity), and
after some computations we obtain:

F+(U) = , (4.45)

{ f(c3F+ + Z<?uF+ + 3cu2F+ + u + + uF+)

where
U=(p,pu,E)T. (4.46)
The functions F / , 0 < d < 3 are denned by

(4.47)

where B = moi(mtn(-u/c, 1), —1).

The flux F~ is obtained by replacing " + " by " — " in the expression (4.45), and Fj" are
defined by

0 < d < 3. (4.48)

dj-d+2
102
We will now show how to update all the variables involved in the construction of the
equilibium functions. At time level tn+1, the elementary variables p, u,e are known from
the state vector U of the conservative variables, for simplicity in exposition, we skip the
indexes j and n + 1. The computation of these variables is achieved using an iterative
method. Indeed: The table from reference [19] gives the density as a function of one of
the following three pairs of thermodynamic variables (p, h), (p, T) or (p, s). Since for real
fluids, an explicit relation specifying the dependence of temperature T or entropy s upon
the other thermodynamic variables of state is not available, we choose the first pair, i.e.
(p, h), where h is given by

h = e+p/p, (4.49)

and the EOS is now of the form

p = f(p,h). (4.50)

The problem now consists in finding the pair (p, h). By using the relation (4.49), the
problem is reduced to finding the enthalpy h satisfying

p = f(p(h-e),h). (4.51)

This leads to a nonlinear equation on the enthalpy. It is solved by using the iterative
algorithm of Newton-Raphson.
We have choosen the enthalpy h rather than pressure p to solve the equation (4.50),
because the density is virtually insensitive to pressure variations in one-phase regions.
The computation of the other variables is now straightforward: the pressure is obtained
from (4.49), the sound speed is given by [19], and the nondimensional quantities 7 and 7
follow thanks to (4.24) and (4.16).

The stability results are established by the next theorem.

Theorem 4.1 Suppose that, 1 < 7 < 3 and 1 < 7. If p] > 0 and e] > 0 for all j <E Z,
then p]+1 > 0 and e" +1 > 0 for all j € Z, and the scheme (4.38) is L1 stable.

Proof. To prove this theorem, we follow the ideas of Perthame developed in [17].
The hypothesis of theorem 4.1 implies that fn{x,v) > 0 and gn(x,v) > 0 (see (4.26)). As
an outcome, the solutions of (4.31) are positive, i.e. /(x,u,£) > 0 and g(x,v,t) > 0 for
t > tn.
Using lemma 4.1, we have

[ f(x,vX+1)dvdx > 0,

with h = Ax. We also have

n+1 2

En+l — i(u Lï
103
 * ïr +1/a L5((«'-«r l ) a +(«r l ) 2 + 2 «r l (»-«r i
and hence e" +1 > 0.

Our next concern will be the L1 stability of the scheme. Since the scheme is conser-
vative, we have
h
= £ \PT\ =

\\E(x,tn+1)\\LHlR) = h £ |E; +1 | = ft £ £; + 1 = ft £ £7 =
J6^ jtZ jtZ

According to the Cauchy-Schwarz inequality, it follows that,

and the proof is complete.

4.4 Computational results

For quality assurance, the scheme is first tested on two classical problems for the ideal gas.
In the first two cases, we set 7 = 1.4 and 7 = 3 which correspond to the maximum speed
of sound in the gas. Next, a Riemann problem for pure water is investigated. In this later
case , the values of the state equation are obtained using a thermodynamic table given in
reference [19].
All the numerical computations were done with a uniform grid, and the number of dis-
cretization cells is 200.

Test case 1: Sod shock tube. The initial data is given by the following table:

104
 e P u
2.5 1. 0.
UR 2. 0.125 0.

The length of the tube is 1. The initial discontinuity is located at the point x — 0.45.
Figures 4.3-4.6 show the results for the solution at time t = 0.25.

Test case 2: Lax shock tube. The initial data is given by the following table:

V P u
L
U 3.528 0.445 0.698
UR 0.571 0.5 0.

For this test also, the length of the tube is 1, but the initial discontinuity is now located
at i = 0.5. The final time is t = 0.16. The figures 4.7-4.10 give the corresponding results.
We refer the reader to [17] for a comparaison of the results obtained for these two test
cases.

Test case 3: we consider now a shock tube problem for pure water. The left and
right states are defined by the following table:

p(bar) p(kg/m>) u (m/s)

400. 1014.34 0.
k 0.955 997.162 0.
u
The initial discontinuity is located in the midle of the tube, i.e. at x = 0.5. Notice that
the values tabulated above correspond to the same temperature for each side of the dis-
continuity, T = 25°c. This test is similar to the one presented in [2]. The difference is that
the pressure in the left state is limited to 400 bar (rather than 1000 bar as in [2]), and this
because, the table from reference [19] is specifically designed to fit the state equation on
a set of regions in pressure which covers the physical range of interest for the analysis of
flow problems in nuclear reactors. Pressures much above 400 bar do not appear in such
problems, and so, these are not available.
Figures 4.11-4.16 show the results for the solution at time 2.08 10~4 obtained by our nu-
merical scheme. In figure 4.12, density variations across the contact discontinuity are weak,
although its position is apparent in Figure 4.14. The backward moving rarefaction wave
and the forward moving shock appear to be moving at approximately the same speed,
because the fluid velocity is so small in comparison to the speed of sound in water. Figures
4.15 and 4.16 show the important variations of 7 and 7 and the difference between them,
(in the perfect gas case, 7 = 7). Indeed we notice that 7 — 1 can be very small while 7 is
always greater than 50. However this does not affect the numerical results. This underline
the capabilities of the proposed method.

105
0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.3: Pressure, Sod shock tube.

0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.4: Velocity, Sod shock tube.

106
0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.5: Density, Sod shock tube.

0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.6: Internal energy, Sod shock tube.

107
s
a.

0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.7: Pressure, Lax shock tube.

•02
0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.8: Velocity, Lax shock tube.

108
 0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.9: Density, Lax shock tube.

2>
tu

0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.10: Internal energy, Lax shock tube.

109
4.56+07

0.4 0.5 0.6 0.7 0.8 0.9

distance (m)

Figure 4.11: Pressure of the liquid at time 2.08 10"4.

0.1 0.3 0.4 0.5 0.6 0.7 0.8 0.9

distance (m)

Figure 4.12: Density of the liquid at time 2.08 10~4.

110
 0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.13: Velocity of the liquid at time 2.08 10"4.

145000

140000 -

135000

130000

125000

120000

115000

110000 •

105000
0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9

Figure 4.14: Enthalpy of the liquid at time 2.08 10"4.

Ill
1.35

1.25 •

1.15 -

1.05

0.1 02 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.15: 7 of the liquid at time 2.08 10~4.

7000

6000 -

5000

4000 •

3000

2000

1000

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 4.16: 7 of the liquid at time 2.08 10~4.

112
Bibliographie

[1] P H . COLELLA& H. M. GLAZ, Efficient Solution Algorithms for the Riemann Problem
for Real Gases, J. Comput. Phys., 59, 264-289, 1985.
[2] C. H. CoOKE & T.-J. CHEN, On shock capturing for pure water with general equation
of state, Communications in Applied Numerical Methods, Vol 8, p. 219-233, 1992.

[3] F. COQUEL, K. EL AMINE, E. GODLEWSKI, B. PERTHAME & P. RASCLE, Numerical

method using upwind schemes for the resolution of two-phase flows, to apear in J.
Comput. Phys.
[4] F. DUBOIS, Evaluation du flux d'Osher pour l'air à l'équilibre chimique, Rapport
Aérospatiale, Programme Tebaldi, Code Euler 3D en hypersonique et supersonique
élevé, S/DT/MI n° 41 747, avril 1989.

[5] P. GLAISTER, An Approximate linearised Riemann Solver for the Euler Equations
for Real Gases, J. Comput. Phys., 74, 382-408, 1988.
[6] E. GODLEWSKI & P.-A. RAVIART, Hyperbolic systems of conservation laws, Ellipses,
Mathématiques & applications, vol. 3/4, 1990.
[7] E. GODLEWSKI & P.-A. RAVIART, Numerical approximation of hyperbolic systems
of conservation laws, AMS 118, Springer-Verlag, New-York, 1996.

[8] A. HARTEN, P.D. LAX & B. VAN LEER, On Upstream Differencing and Godunov-
type Schemes for Hyperbolic Conservation Laws, SIAM Review, Vol 25, n° 1, pp.
35-61, 1983.
[9] D.L. HICKS, hyperbolic models for Two Phase (or Two Material) Flow., Sandia Na-
tional Laboratories. SAND81 0253. 1981.
[10] R. HOLL, Wellenfokussierung in Fluiden, Ph. D. Dissertation, Technishen Hochschule
Aachen, West Germany, 1982.

[11] S. KANIEL, A Kinetic Model for the Compressible Flow equations, Indiana University
Mathematics Journal, Vol. 37, N° 3, p. 537-563, 1988.
[12] R.J. LEVEQUE, Numerical methods for conservation laws, Birkhàuser, Bale, 1990,
Lecture Notes in Mathematics.

113
[13] M. S. Liou, B. VAN LEER & J. S. SHUEN, Splitting of Inviscid Fluxes for Real
Gases, J. Comput. Phys., 87, 1-24, 1988.
[14] R. MENIKOFF & B. J. PLOHR The Riemann problem for fluid flow of real materials,
Reviews of Modern Physics, Vol. 61, TV" 1, p. 75-130, January 1989.
[15] J.-L. MONTAGNE, H. C. YEE k M. VINOKUR Comparaison study of high-resolution
shock-capturing schemes for real gas, NASA Technical Memorandum 100 004, July
1987, Ames Research center, Moffett Field, California.
[16] S. OsHER h E. TADMOR, On the convergence of difference approximations to scalar
conservation laws , Math. Comp., Vol 181, pp. 19-51, 1988.
[17] B. PERTHAME, Second-order Boltzmann schemes for compressible Euler equations
in one and two space dimensions, SIAM J. Numer. Anal., Vol.29, n° 1, pp. 1-19,
February 1992.
[18] V.H. RANSOM & D.L. HICKS, Hyperbolic Two-Pressure Models for Two-Phase Flow,
J. Comput. Phys., 53, 124-151, 1984.
[19] P. RASCLE & O. MORVANT, THETIS : Calcul par interpolation des fonctions ther-
modynamiques de fluides diphasiques, EDF HT-13/94/014, 1994.
[20] P.L. ROE, Approximate Riemann Solvers, Parameter Vectors, and Difference
Schemes, J. Comput. Phys., 43, 357-372, 1981.
[21] S. SRINIVASAN, J. C. TANNEHILL & K. J. WEILMUENSTER, Simplified Curve Fits
for the Thermodynamic Properties of Equilibrium Air, NASA R. P. n° 1181, août,
1987.
[22] H.B. STEWART & B. WENDROFF, Two-Phase Flow: Models and Methods, J. Comput.
Phys., 56, 363-409, 1984.
[23] I. TOUMI, A Weak Formulation of Roe's Approximate Riemann Solver, J. Comput.
Phys., 102, 360-373, 1992.
[24] A. ZELMANSE, Formulation cinétique et schémas de Boltzmann pour le calcul
numérique en mécanique des fluides, Thèse de Doctorat de l'Université Paris XIII,
1995.

114
Chapitre 5

UNE METHODE NUMERIQUE

DECENTREE POUR LES
ECOULEMENTS
DIPHASIQUES EN
DESEQUILIBRE

115

!®f t
 AN UPWIND NUMERICAL METHOD FOR NON-EQUILIBRIUM
TWO-PHASE FLOWS

5.1 Introduction and description

In most practical situations, the flow field and the topographical distribution of the phases
must be described using average properties. Indeed, the presence of the interfaces separat-
ing the phases and the interfacial exchange of mass, momentum, and energy, complicate
considerably the analysis of two-phase flow in comparison to single-phase flows. An exact
formulation of two-phase flow problem would have required the description of the evolution
in time and space of the fields of each phase, together with a prediction of the geometry
of the interfaces. Such an approach is totaly impractical, except in some relatively simple
situations, such as the case of horizontal stratified flow.

In the general case of two-fluid nonequilibrium models, the flow fields of the two
phases and the exchanges occuring at the interfaces are explicitly taken into consideration.
Proper averaging of the local instantaneous conservation equations ([6], [19], . . . ) yields
two sets of phase conservation equations (one for each phase) in terms of phase-average
properties and balance equations that describe what happens at the interface (interfacial
exchange of mass, momentum and energy). The dynamic of interactions between the two
phases are fully described by laws governing the interfacial mass, momentum and energy
exchanges. This approach results in the so-called "six-equation models". No assumptions
are made in this case regarding, for example, thermal equilibrium or the velocity ratio.
Various models dealt with fewer equations and the choice between the various alternatives
depends primarily on the nature of the problem to be solved. For example, by forcing one
of the phases or both to remain saturated, or by assuming equal average phase velocities,
the number of conservation equations is correspondingly reduced. In this case, we end up
with four or five equation models. In the case of the mixture models, the existence of
interfacial exchages is not explicitly considered; the evolution of the mixture is a priori
specified using relationships linking the behaviour of variables such as the void fraction
and the departures from thermal equilibrium of the phases to local variables such as the
pressure. In this case, only three mixture conservation equations are needed. The simplest
possible mixture model is the homogeneous-equilibrium model (in which the phases are
intimately mixed), in this later, equal phase velocities and thermal equilibrium between
the phases are assumed. The full six-equation model might be needed, for example, for
establishing the critical flow condition or for calculating the evolution of the mixture during
a fast transient during which strong departures from equilibrium are expected to occur.

117
Detailed descriptions of two-phase models were published by several authors: Delhaye [6],
Ishii [16] etc. The working set of one-dimensional two-fluid conservation equations that
is given below, is similar to the sets used in the modern codes for the analysis of Light
Water Reactors (LWR) under off-normal and accident conditions and constitutes a rather
complex set of six coupled partial differential equations. The solution of this set requires
external specification of several laws governing the exchanges between the phases and
between the phases and the wall.

5.2 Two-fluid nonequilibrium model of two-phase flows

5.2.1 The physical model

The two-fluid nonequilibrium model consists of individual mass, momentum and energy
conservation equations for both the liquid and vapor/gas phase. The interfacial and wall
transfer terms appear in the field equations as source terms. In one-dimension, this model
can be written under the following form,

dt(avpv) + dx(avpvuv) = r v ,

dt{<*ipi) + dx{aipiui) = Ti,

dt(avpvuv) + dx(avpvul + avp) - Pidxav = avpvg + Tvuiv +

dx{aipiuf + atp) - Pidxai = aiptg + Tiuit + M$

uv2 p uv2
dt(avpv(ev + -=-)) + dx(avpv(ev -\ 1- -~-)uv) + pdtav =
2 pv 2

avpvguv + rv(hiv + ^f-) -H M,Çulv + </,-„ + ç wt ,,

2 2
• T ) ) + «r(a/Pi(e/ + - + - j

The volume fractions are related by,

<*„ + <*/ = !. (5.2)

In these equations, the v and / subscripts refer to the vapor and liquid phases, re-
spectively. The variables appearing in the above conservation equations have the following
meanings:

118
 Qfe = volume fraction of k-phase,
p = the common pressure,
Pi = interfacial pressure,
Pk = density of k-phase,
Uk = velocity of k-phase,
ejt = specific internal energy of k-phase,
hik = enthalpy of k-phase at the interface,
Uik = velocity of k-phase at the interface,
g = force of gravity.

The source terms Tk, Mj^, Mwk, qik and qwk are the mass transfer rate at the interface,
interfacial friction (drag force), wall friction, interfacial heat transfer and the wall heat
transfer into phase k respectively. Here we have assumed gravity to be the only external
body force acting. In order that mass, momentum and energy be conserved at the interface
between phases, the following relations must hold:

r
* = o. (5-3)
k=v,l
MPk=0, (5.4)

£ Tkihik + ^) + MPkuik + qik = 0. (5.5)

Jfc=v,/

In order to close the system of equations (5.1), additional relationships are needed. We
refer to these relationships as constitutive equations. These equations generally specify
how the materials (i.e., phases) interact with themselves and with each other. A general
form for each of these equations is given below; specific forms and assumptions used in
our computations are also presented.

Constitutive equations

The first of these constitutive equations are the equations of state: for the moment, we
assume that a relationship giving the vapor (respectively the liquid) density as a function
of pressure and the vapor (respectively the liquid) specific enthalpy is available:

Pk = Pk{p,hk). (5.6)

These relations are given by thermodynamic table [26]. In practice, these are usually
approached by explicit equations (see section 5.3).

In addition to these thermodynamic relations, specifications for the interfacial and

wall transfer terms are required for closure.

119
 In this paper, thermal and mechanical non-equilibrium between the phases and the
processes involving interfacial mass and energy exchanges and interfacial friction will be of
interest. The interphase mass transfer term F describes the net transfer of mass from one
phase to the other. The most important basic physical processes by which this can occur
are condensation and boiling, (while flashing, i.e., vapor generation solely as the result of a
reduction in system pressure, can be considered as a special type of boiling). The processes
of condensation and boiling require that heat be transferred between phases; hence the
mass transfer rate can be computed by dividing the energy transfer that goes into the mass
exchange by the heat of vaporization. The coupling between the mass and heat transfer
processes suggest that the driving potential for both mass and energy transfer must be
related to the degree of thermal nonequilibrium. Thus, we set the following algebraic
relationships:

i- Tv),

qu = Ha(Ti - 7}),
Hiv and Hu are the heat transfer coefficients from the interface. These coefficients can
be written under the form Hik — <Xv<*iPkCpk/tik where cpk is the specific heat at constant
pressure of phase k and tj£ can be considered to be a rate constant for the interfacial
heat transfer process. Models of any desired degree of complexity can be constructed by
appropriate calculations of Uk. In our case, we have assumed *,-* = constant. Replace
cpk(Ti — Tk) by hik — hk, the resulting form for qik is
qik = avaipk(hik - hk)/tik. (5.7)
This form was chosen not because of any expected physical validity, although it should be
adequate as a rough approximation, but because of its simplicity.
The interfacial heat fluxes and the mass exchange rates Tk are usually linked through the
following thermal energy jump condition

J2 (ft* +r **«) = °> (5-8)

k=v,l
thus, the interphase vaporisation (or condensation) rate is given by,

r, = -f^±fL, (5.9)
«i "il
i.e., the vapor source term in the vapor continuity equation (Tv) is equal to the heat
transfer out of the liquid divided by the latent heat of vaporisation.

To derive a practical set of constitutive equations, it is necessary to make some

further assumptions:
the interfacial temperature is assumed to be the saturation temperature,
Ti = r,,

120
the interface enthalpies are taken to have the saturation values of the phase in question,

hik = h$k(p),

and the phasic velocities at the interface are equal,

uiv = uu = Ui, (5.10)

where u, is some average interfacial velocity which we assume at the form (see chapter 2)

U{ = avv,[ + cniuv.

All these assumptions are open to question, but are simplest to make at this time.

In two-fluid modeling each phase may have an indépendant velocity and, in fact,
may move in opposite directions from each other (countercurrent flow). The interphase
momentum transfer during this process is the result of the forces exerted on the phases
by each other. The forces can be considered as two basic components: i) the friction force
between phases M $ which is similar to the friction between the individual phases and
the wall and ii) the interference between the two phases. The second (transient) category
accounts for the fact that both phases are in contact with each other and when one phase
moves the other phase must be displaced; it is commonly called the "virtual mass force"
M?k (see chapter 2) and will not be considered in the next. Both components of the
interphase momentum transfer tend to reduce the slip between the phases, for most real
and transient applications. The general expression for M-% is given in chapter 2, here we
have used the following expression

Now, it remains to give an explicit form for the wall transfer terms. For qwk we use
Qwk — <*kQw where qw is given as a function of x, this simplified model neglects the heat
content of the wall. For the frictional terms Mwk, we use Mwk = —FwkakpkUk, Fwk is a
frictional coefficient at the wall.

5.3 Numerical method

To solve the system of equations (5.1), we will use the splitting method introduced in
chapter 3. We now give a brief outline of the method and the advantages it offers for
the numerical approximation. First, we split the system of equations into two "fractional
steps" or two stages. In the first step, we solve simultaneously two systems of equations

121
which look like those of the quasi-lD nozzle problem with source terms and in the second
step, we solve one system of ordinary differential equations (O.D.E.) in order to take into
account the coupling of the two phases.

First step (hydrodynamical step)

In this step, we solve the two following systems, one for each phase k (k = v, I):

dt(mk) + dx{mkuk) = 0,

dt{mkuk) + dx{mkul + pk) = pidxak + akpkg + Mwk, (5.11)

dt (Ek) + dx({Ek + pk)uk) = akpkguk + qwk,

where

mk = akpk,

u2
Ek ~ akpk{ek + •—),

Pk = <xkpk.

The expression of the interfacial pressure p, is given later (see (5.28)).

As in chapter 3, we will use an explicit method to discretize these two systems. This
lead to two uncoupled problems which may be solved separately. Recall also that, since
the two systems evolve independently, the phasic pressures obtained from this step may
have different values. In the next, the algebraic terms in the right-hand side (gravity forces
and exchange terms between the phases and the wall) will be ignored. In practice, these
terms are treated explicitly. The remaining systems to be solved are now

where Uk = {mk, mkuk, Ek)T, F{Uk) = (mkuk,mkul + pk, (Ek + pk)uk)T and H(Uk) -
(0,pidxak,0)T, k = v,l.

Second step (restoring equality of pressures)

In this step, we tie the two sets of equations of the first step using the fact that
the pressure of the two phases are the same (pv = pi), and we solve the following coupled

122
system with interfacial exchange terms.

dt{mvuv) = rvUi

dtEv = -pdtav +
(5.13)
dtmt = Ti,

= -pdta, + r,{hsi + f ) + Mffw + qu.

This is a splitting which reduces the resolution of the full system (5.1) to classical
hydrodynamics (5.11) with source terms, and a problem (5.13) of ordinary differential
equations.

5.3.1 Discretization of the first step

The problems (5.12) will be solved using the same methods, so we skip the index k. The
system to be considered writes now
ÔU dF(U)
(5.14)

where

U = (m, mu,E)T is the vector of conservative variables,

F(U) = (mu, mti2 + p, (E + p)u)T is the flux vector,
H(U) — (0,pidxCK,0)T is the source term,

and

Y
p = ap.

The system (5.14) must be supplemented by an explicit constitutive equation relating

the pressure p to the other thermodynamic variables i.e.: p and e. This equation is called
the EOS (equation of state) and characterizes the fluid material.
For an ideal gas, the EOS reduces to p^= (7 — l)pe where 7 is the ratio of enthalpy to
internal energy, and is a constant. Several authors [1], [20] and [21], have advocated the

123
use of an equivalent 7 for real gases. In our case, we will use the same strategy for the two
phases i.e. vapor and liquid.
As in the ideal gas, we define

7=^, (5.15)

h is the specific enthalpy. Combining this with the definition h — e + p/p gives

p=( 7 -l)pe. (5.16)

This is identical in form to the EOS of the ideal gas, but now we have 7 = 7(p,e).
In terms of the variables involved in our system we have

p=(7-l)me. (5.17)

And the general thermodynamic equation given in (5.6) is approached by

(7 can be expressed as a function of (p, h)).

It will be convenient to introduce another dimensionless variable called the adiabatic
exponent

y = -V-(9/) = P
-c\ (5.19)
p \dvj3 p
where v = l/p, is the specific volume, 5 is the specific entropy and c = (dp/dp)l is the
speed of sound. The variable 7 is seen to be a dimensionless measure of the sound speed,
and it will be employed for the contruction of the numerical scheme. For most materials,

To dicretize the system (5.14), we use an explicit upwind scheme based on finite-
volume method. Let Ax denote the mesh width and At the time step, we set tn = nAt,
n € W; Xj— jAx and ij_i/2 = {j - l/2)Ax, j e 7Z\ and a = At/Ax.

Assume that, at time tn, we know f/j1, an approximation to the average of the
solution U to (5.14) at time tn across the cell ]ij_i/2i x j+i/2[- Hence f/"+1, defined in the
same way as U? but at time tn+1 = tn + At is given by

with

[/;-—/ U{x,t*)dx. (5.21)

Ax JXj_lf2
In (5.20), H" représente an approximation of the source term H(U), (see below).
The numerical flux function F?+1,2 plays the role of an average flux through Sj+1/2; it is
usually noted by

124
This flux can be obtained through the approximate Riemann solvers developped for the
numerical solution of real gas. Here, we have used the kinetic solver introduced in chapter
4. To compute the flux F?, t ,2, it suffices to replace the density p by the partial density m
and the pressure p by the partial pressure p. The flux writes then
(5.22)

where

F+(U) = 2cuF? + u2F+) , (5.23)

f ? + u3F0+) + fi(7)p(cF+ + UF0+)

with
U = (m, mu, E)T. (5.24)

The functions FJ1" are defined by

F+ = (5.25)
d dy-d +
where B — max(min(—u/c, 1), —1).
The flux F~ is obtained by replacing " + " by " - " in the expression (5.23). FJ are defined
by

F~
d = (5.26)
dj-d +2

In this chapter, we consider only the case p,- ^ p and write Pidxa in (5.11) (see (5.1))
as
Pidxa = {p-(p- Pi))dxa. (5.27)
The expression of (p - p.) (see chapters 2 and 3) is given by:
p-Pi = Pv{uv-ui)2, (5.28)

Next, it remains to discretize (5.27). The quantities appearing in (5.27) will be taken
at time tn and its discretized form is given by the following centered formula
1
( \

This discretization is sufficient to stabilize the computations (see section 5.4).

The quantities computed from this step will be denoted by m£ +1 , {m.kUk)n+1^2 and

125
5.3.2 Discretization of the second step

Now, we are left with the task of discretizing the following system

dt(mvuv) = Mvi,

dtEv — —pdtav +
(5.30)

dt{miui) = Mu,

dtEi = -pdtat + En,

where, for k = v, I,

Mki= %
Eki = Tk(hsk + u?

For simplicity, we drop the index j through out this subsection since the computations
are performed cell by cell. Since the correlations describing interfacial transport processes
have often very small time constants, the source terms require a strictly implicit treatment.
The discretised system associated to (5.30) is then

m7vlT1 - " + 1 = 0,

n+1 1
{mvuv) - = 0,

1
^ = 0,
(5.31)
= 0,

= 0,

The motivation for the discretization of the terms pdtak in the energy equations is given
in chapter 3.

126
 In the system of equations (5.31), replacing rrik and Ek by their respective expres-
sions (defined above), ejt by hk — p/pk in the energy equations, and adding the relation
av + Q/ = 1 we obtain:

- mnv+1/2 = 0,

{avPvuvy+l - {mvuv)n+ll* - = 0,

= 0,

(5.32)

= 0,

= 0,

= 0.

The interfacial exchange terms involved in (5.32) are expressed as

where fc^
From the principle of conservation at the interface between phases we have:
pn+l _ _pn+l

+ =

The relaxation time constants t,jt (k — v,l) are assumed to have the same value
Uv = Ul = 0.005 sec. When written in these form, the interfacial exchange terms are
implicit enough to be stable.

Since we deal with nonanalytic equations of state, the conservative variables are not
chosen as fundamental unknowns to solve the nonlinear system of equations (5.32). For this
reason, the vector of fundamental unknowns is composed with elementary variables. The
appropriate fundamental unknowns are (a v ,u v ,/i v ,p, aj, u/,/i/). Notice that, the solution
method remain conservative.

127
 The variables pk appearing in the conservative variables and the expression of the
interfacial terms (see above) are eliminated by using the constitutive equations, these are
obtained from a thermodynamic table [26] as functions of {hk,p). Densities were not chosen
as fundamental unknowns because these are not input variables in thermodynamic table
and this because the pressure equation of state is a very sensitive function of the density
for an all-liquid system (see figure 5.1). Thus, if density were a fundamental unknown
small errors in the density would produce large errors in pressure. This would lead to a
less stable algorithm.

Notice that, the thermodynamic table from reference [26] is specifically designed
for efficiently computing one or two-phase flows in thermodynamic equilibrium. Moreover,
in two-phase region the equation of state (see figure 5.1) gives only the density of the
mixture as a function of the pressure and the specific enthaly of the mixture. Now, six-
equation two-fluid models require an equation of state for each phase. However, in a
thermodynamic nonequilibrium flow (superheated liquid or/and supercooled vapor) the
enthalpies may lies outside the pure phase regions (see figure (5.2)). In this case, the
equations of state corresponding to each pure phase are simply extended by extrapolation
across the saturation boundary (physically corresponding to a supersaturated metastable
phase), and so the densities are computed. The saturated enthalpies as function of the
pressure are given by [26].

Now, it remains to solve the nonlinear system (5.32).

In the case where interfacial exchange terms are neglected, the system (5.32) reads

(<*vPvK)n+1 + (a v p v = 0.
+1
a; + a?+1 - 1 = 0, (5.33)

(a//>/u,) n+1 =

(a,M<)n+1 + = 0.

The solution of this nonlinear system reduce to the one of a nonlinear equation as we shall
see below. From (5.33.1, 5.33.2) and (5.33.5, 5.33.6) we deduce the velocities,

{mvuv) n+l/2
(5.34)
l

128
 500000
1e+06
1.5e+06

ilpy(J/Kg)
4e+0
3.5e+07

1e+07
Pressure (Pa) 5e+O6

Figure 5.1: Equilibrium equation of state p = p(p, h).

IT
Press

critical point
CL
1
saturated vapor boundary
saturated liquid boundary

y
Two-phase \ Vapor region
Liquid region /
region

h (Enthalpy)

Figure 5.2: Phase diagram in the p-h plane.

129
and

l =

Next, using the energy equations gives the enthalpies as functions of the pressure,

(5-36)

and

' ~mr i/2 2

mr i/2 ' l
'
Recall that, the densities are obtained from the thermodynamic table as functions of
enthalpies and pressure. Now, according to (5.36) and (5.37) it follows that the densities
are only functions of the pressure

l=ff(p n + 1 ), (5.38)
and

PÏ+1 = p?+1{pn+1,h?+1) = p? + 1 (p n + 1 ,/i? + 1 (j> n + 1 )) = h{pn+1). (5.39)

Since we will use an iterative method to solve the problem, it is convenient to work with
nondimensional equations. For this reason, we consider the equation (5.33.4). The volume
fractions are obtained from the mass equations

(5.40)
g{pn+1Y
and

Substituting (5.40) and (5.41) into (5.33.4) gives

This nonlinear equation is solved by a Newton method. The algorithm writes

where the subscript i indicate successive iterative approximations. The derivative /'(pf
at the point p " + 1 is given by

1 + (5 44)
« W ) ) ] V^PbK^PIpr^j ' '
130
with

-1

and
n+l/2 -1

The partial derivatives above are computed at the points (p" +1 , ^£ +1 (p" +1 )) (f° rfc= u,/)
and are obtained from the table in reference [26]. The starting point is simply chosen
among p", p"_ 1 or p" + 1 where j and n refers respectively to the cell j and time tn.

We now return to (5.32). This nonlinear system can be written in the form:

F(X) = 0, (5.45)

with
X = (Q V , UV, hv, p, ai,ui, hi), (5.46)

and

(5.47)

Since we deal with a stiff problem, the components of X and F(X) are first scaled by
their typical values. The resulting system is then solved using a globally convergent
method based on an iterative Newton-Raphson algorithm (see [7]). The Jacobian matrix
dF(X)/dX is difficult to compute analytically, therefore we have used finite differences to
compute the partial derivatives.

From the components of X (5.46), the computation of all the other variables is
straightforward. Recall that, the speed of sound for each phase is given by c* = Ck{p,hk),
for k = v,l, and is obtained from thermodynamic table.

5.4 Computational results

Test case 1: This test was proposed by Toumi in [34]. It consists of a two-phase shock
tube problem. The initial state is given by the following table:

a p (bar) uv (m/s) ui {m/s) K (kj/kg) hi (kj/kg)

LO.25 200. 0. 0. 3092.7 1338.2
UR 0.25 150. 0. 0. 3092.7 1338.2

131
 The purpose is to compare the results obtained by two different schemes: the Lax-
Friedrichs and the proposed Kinetic schemes when used in the first step of our method.
As we can see in the Figures 5.3 and 5.4 the Lax-Friedrichs scheme has much more nu-
merical diffusion than the Kinetic one but the results obtained with the two schemes are
in reasonable agreement. These results can also be compared to those obtained by Toumi
[34] which use an approximate Riemann solver (the Roe one) on the full system.

Test case 2: Edward's Pipe Blowdown.

This test contains many features of the depressurization and blowdown that occur in pos-
tulated loss of coolant accident in a reactor system . It tests the two-phase model through
a wide range of void fraction variation and involves the models for interphase friction,
interphase heat transfer and mass transfer. The test consists of a horizontal pipe of 4 m
length. In the initial conditions, The pipe is filled with a saturated liquid at 7 x 106 Pa.

At time zero, the pipe was blown down to atmospheric pressure. The fast depres-
surization results in a strong evaporation of the liquid phase "flashing". The void fraction
distribution along the pipe length at different times is shown in Figures 5.5 and 5.6. The
initial phase of the depressurization is characterized by the propagation of a rarefaction
wave from the opening into the pipe and its reflection at the closed end (see Figure 5.7).
When the steady-state is established (see Figure 5.8), the pressure in the pipe is equal
to the atmospheric pressure and all the thermodynamic variables have saturated values
(see for example the liquid density displayed in figure 5.11). Figures 5.9 and 5.10 show
the phasic velocities, we can remark that the liquid velocity is always below the velocity
of the vapor phase.

L-F scheme »
Kinetic scheme —

0.2 0.4 0.8 1 1.2

distance (m)

Figure 5.3: Pressure, comparaison between Lax-Friedrichs and Kinetic schemes.

132
 70
L-F scheme o
Kinetic scheme
60

SO

I .

-10
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
distance (m)

Figure 5.4: Vapor velocity, comparaison between Lax-Fried richs and Kinetic schemes.

I = 0.004 s
I = 0.008 s
0.9 <=-OOT6 s
t = 0.062 s

0.8 t=jfc4s — - -

0.7

0.6
S
u 0.5

0.4

0.3

0.2

0.1

0
0.5 1 1.5 2 2.5 3 3.5
distance (m)

Figure 5.5: Void fraction history.

133
 void traction

1.5 distance (m)

«me(s) 0.5

Figure 5.6: Void fraction in space-time plane.

t«0.000s
t-0.002 s
1 = 0.004 s
7e+06 -••.. T^.U.UlbS

t»o!toos —
^-__ \ t-0.190s
6e+06 /" ~ —-^ -. t = 0.400 s "
\
\

5e+06 L

I 4e+06
\ \

38+06
--- _ \\v
2e+06 -
"A..^
1e+06 -

0
0.5 1.5 2 2.5 3.5
distance (m)

Figure 5.7: Pressure history.

134
 pressure
8e+06-

1.5 distance (m)

time (s) 0.5

Figure 5.8: Pressure in space-time plane.

250
= 0.000 s
= 0.002 s - —
= 0.010 S
= 0.040 s
200 = 0.100 S /-
= 0.200 s — •/
= 0.300 S •••/•
t = 0.400 S -/•- /

150

I
100

2
S
50

-50
0.5 1.5 2 2.5 3.5
distance (m)

Figure 5.9: Transient liquid velocity profile.

135
250
t = 0.000 s /
t = 0.002 S '
t = 0.010 s • ••-•/
\ = 0.040 s —*•
200 t = 0.100 s -•-'- -
t = 0.200 s ?—!
t = 0.300 S/'--••/
t « 0.400,4
y
-•••*
/
S /
150 .'•
y
J
!

y L
100
y'
/—-j.
50

-50
0.5 1.5 2 2.5 3.5
distance (m)

Figure 5.10: Transient vapor velocity profile.

0.05 0.1 1.5 distance (m)

0.15 -,
Ci O*\
«me (s) 0.3 0.35
0.4

Figure 5.11: Liquid density in space-time plane.

136
Bibliographie

[1] P H . COLELLA k H. M. GLAZ Efficient Solution Algorithms for the Riemann Problem
for Real Gases, J. Comput. Phys., 59, 264-289, 1985.

[2] F. COQUEL, K. EL AMINÉ, E. GODLEWSKI, B. PERTHAME k P . RASCLE, Numerical

method using upwind schemes for the resolution of two-phase flows, to apear in J.
Comput. Phys.
[3] K. EL AMINÉ k P. RASCLE, Ecoulement Diphasique, Modélisation et Hyperbolicité,
Rapport interne EDF HT-13/96/004, 1996.
[4] J.A. BOURÉ, Mathematical Modeling and the Two-phase Constitutive Equations, Eu-
ropean Two-phase Flow Group Meeting, Haifa, 1975.

[5] J.A. BOURÉ, Les lois constitutives des modèles d'écoulements diphasiques monodi-
mensionnels, à deux fluides , Rapport CEA-R-4915, 1978.
[6] J.M. DELHAYE, Basic Equations for Two-phase Flow Modeling, Two-phase Flow
and Heat Transfer in the Power and Process Industries, A.E. BERGLES et ai, Eds.
Hemisphere Publishing, Washington, 1981.

[7] J. E. DENNIS k R. B. SCHNABEL, Numerical Methods for Unconstrained Optimiza-

tion and Nonlinear Equations, (Prentice-Hall, Engelwood Cliffs, NJ), 1983.
[8] D.A. DREW, Analytical Modeling of Multiphase Flows, Edited by R.T. LAHEY Jr.,
Boiling Heat Transfer, Modern Development and Advances, Elsevier Science Publish-
ers, 31-84, 1992.
[9] R. L. FERCH, Method of Characteristics Solutions for Non-Equilibrium Transient
Flow-Boiling, Int. J. Multiphase Flow, Vol. 5, 265-279, 1979.

[10] P. GLAISTER, An Approximate linearised Riemann Solver for the Euler Equations
for Real Gases, J. Comput. Phys., 74, 382-408, 1988.
[11] E. GODLEWSKI k P.-A. RAVIART, Hyperbolic systems of conservation laws, Ellipses,
Mathématiques k applications, vol. 3/4, 1990.
[12] E. GODLEWSKI k P.-A. RAVIART, Numerical approximation of hyperbolic systems
of conservation laws, AMS 118, Springer-Verlag, New-York, 1996.

137
[13] A. HARTEN, P.D. LAX & B. VAN LEER, On Upstream Differencing and Godunov-
type Schemes for Hyperbolic Conservation Laws, SIAM Review, Vol 25, n° 1, pp.
35-61, 1983.
[14] G.F. HEWITT, J.M. DELHAYE k N. ZuBER Multiphase Science and Technology,
Volume 5 Hemisphere Publishing Corporation, 1990.
[15] D.L. HICKS hyperbolic models for Two Phase (or Two Material) Flow., Sandia Na-
tional Laboratories. SAND81 0253. 1981.
[16] M. ISHII, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975.
[17] M. ISHII k K. MISHIMA, Two-Fluid Model and Hydrodynamic Constitutive Relations,
Nuclear Engineering and Design, 82, 107-126, 1984.
[18] M. ISHII k N. ZUBER, Relative Motion and Interfacial Drag Coefficient in Dispersed
Two-phase Flow of Bubbles, Drops and Particles, A. I. Ch. E. Journal,, vol 25,
843-855, 1979.
[19] R.T. LAHEY k D.A. DREW, The Three Dimensional Time-and Volume Averaged
Conservation Equations of Two-phase Flow, Advances in Nuclear Science and Tech-
nology, Vol.20, 1-69, 1988.
[20] M. S. Liou, B. VAN LEER k J. S. SHUEN Splitting of Inviscid Fluxes for Real Gases.
J. Comput. Phys., 87, 1-24, 1990.
[21] J.-L. MONTAGNE, H. C. YEE k M. VINOKUR Comparaison study of high-resolution
shock-capturing schemes for real gas, NASA Technical Memorandum 100 004, July
1987, Ames Research center, Moffett Field, California.
[22] S. OSHER k E. TADMOR, On the convergence of difference approximations to scalar
conservation laws , Math. Comp., Vol 181, pp. 19-51, 1988.
[23] V.H. RANSOM, Numerical Benchmark Tests, Multiphase Science and Technology,
Volume 3, Edited by G.F. HEWITT, J.M. DELHAY and N. ZUBER, Hemisphere,
1987.
[24] V H. RANSOM, Two-Phase Flows, Ecole d'été, Analyse Numérique, 12-13 Juin, 1989.
[25] V.H. RANSOM k D.L. HICKS, Hyperbolic Two-Pressure Models for Two-Phase Flow,
J. Comput. Phys., 53, 124-151, 1984.
[26] P. RASCLE k 0 . MORVANT, THETIS : Calcul par interpolation des fonctions ther-
modynamiques de fluides diphasiques, EDF HT-13/94/014, 1994.
[27] P.L. ROE, Approximate Riemann Solvers, Parameter Vectors, and Difference
Schemes, J. Comput. Phys., 43, 357-372, 1981.
[28] L. SAINSAULIEU, Contribution à la modélisation mathématique et numérique des
écoulements diphasiques constitués d'un nuage de particules dans un écoulement de
gaz, Thèse d'Habilitation à diriger des recherches, Université Paris VI , 1995.

138
[29] H. STADTKE k R. HOLTBECKER, Hyperbolic Model for Inhomogeneous Two-Phase
Flow, International Conference on Multiphase Flows, University of Tsukuba, Japan,
September 24-27, 1991.
[30] H. STADTKE k R. HOLTBECKER, Numerical Simulation of Two-phase Flow Based on
Hyperbolic Flow Equations, Sixth International Topical Meeting on Nuclear Reactor
Thermal Hydraulics, NURETH-6, October 5-8, 1993, Grenoble France.
[31] H.B. STEWART, Stability of Two-Phase Flow Calculation Using Two-Fluid Models,
J. Comput. Phys., 33, 259-270, 1979.
[32] H.B. STEWART k B. WENDROFF, Two-Phase Flow: Models and Methods, J. Comput.
Phys., 56, 363-409, 1984.
[33] I. TOUMI, A Weak Formulation of Roe's Approximate Riemann Solver, J. Comput.
Phys., 102, 360-373, 1992.
[34] I. TOUMI, An Upwind Numerical Method for a six Equation Two-Phase Flow Model,
Rapport DMT/94/168, SERMA/LETR/94/1619, C.E.A., 1994.
[35] J.A. TRAPP k R.A. RlEMKE, A Nearly-Implicit Hydrodynamic Numerical Scheme
for Two-phase Flows, J. Comput. Phys., 66, 62-82, 1986.
[36] G.B. WALLIS, One-dimensional Two-Phase Flow, McGraw Hill, New York, 1969.

139
SSIXT §»Aûë(S)
!®f t BLAHK
Chapitre 6

LE TRAITEMENT DES
CONDITIONS AUX LIMITES

141
M&WÏ g»A6E(S)
l®f t 1LA&8K
 THE TREATMENT OF BOUNDARY CONDITIONS

6.1 Introduction

In preceding chapters, the basic equations associated with the two-fluid model were dis-
cussed. The solution of this set of equations has been attempted by different authors
applying a variety of numerical methods. One common problem hindring all this efforts
has been the existence of complex characteristics associated with the equation set. Indeed,
as it was seen previously, the two-fluid model is not hyperbolic for unequal phase veloc-
ities and, even for physically relevant expressions of the interfacial coupling between the
two phases including space and time derivatives of phasic physical quantities, it is only
conditionally hyperbolic. Theoretically, for quasi linear partial differential equations, the
existence of complex characteristics is supposed to imply that the basic equation set is not
well posed. This property has been recently associated with inconsistent treatment of the
interfacial exchange terms (the so-called, jump conditions) and the assumption that both
phases are at equal pressure.
Another numerical problem associated with the solution of this set of equations is the
proper treatment of the boundary conditions, i.e.: the required boundary conditions of
physical origin to be imposed (and their number), the specification of the remaining vari-
ables and, the formulation and discretization of these conditions. Due to the "wave-like"
nature of two-phase flows, if care is not taken, any boundary error propagates throughout
the computational domain and, in most cases, may result in a severe instability problem,
eventhough the integration scheme for the interior points is stable. Notice that the number
of boundary conditions may change during a transient.
Since we deal with problems characterized by propagation phenomena, the analysis of the
number and type of conditions at the boundaries is theoretically linked to the properties
of the characteristics. For complex two-fluid applications, the usual mathematical method
of defining the correct number of boundary conditions is by analysing the direction of the
emanating characteristics.

In this chapter, we will give a detailed discussion of the boundary treatment for
the two-fluid model, when it is solved by the splitting method described in the previous
chapters (3 and 5). This method in some extent bypasses the possible lack of real charac-
teristics and in general allows for a convenient and rather simple handling of the problem.
Indeed, the method proposed in preceding chapters and upon which the analysis of the
boundary conditions is developed below, does not require the hyperbolicity of the global
system. One of its advantages is the ease with which boundary conditions can be derived

143
and implemented, based on the physics of the problem, without being concerned with the
complicated structure of the original system.

6.2 Formulation of the boundary conditions for the two-

fluid model

Let us recall the splitting method introduced for approximating the two-fluid model (see
chapters 3 and 5 for the notations).
In the first stage, we solve two uncoupled systems, one for each phase k (k = u,/):

{ dt{mk) + dx{mkuk) = 0,

dt{mkuk) + dx{mku\ + pk) = pkdxak, (6.1)

where

mk = akpk,
u2
Ek =QkPk(

Pk = <xkpk,

and in the second stage, we solve one system of ordinary differential equations (ODE).
The second step is not concerned with the treatment of the boundary conditions, since it
involves only time derivatives.

From now on, we ignore the right-hand side. These two systems reduce to

where the vector Uk is defined by [/J = {nik,mjtwt,Ek), k = v,l.

If the index k is skipped, the system (6.2) writes

where

U = (m, mu, E)T is the vector of conservative variables,

F(U) = (mu, mu2 + p, (E + p)u)T is the flux vector,

144
and

u\

p = ap.

The pressure is given by an equation of state of the form p = p(p, e).

The set of equations (6.3) is solved in a finite domain 0 < x < L, where we shall
always assume that x = 0 is the inflow boundary and x — L the outflow boundary. The
problem is an initial boundary value problem, because both initial data in the region
0 < x < L and the boundary conditions at x = 0, L are needed. Notice that, for most hy-
drodynamics problems, it is well known that a state U(t) can not be completely prescribed
at a given boundary as we shall see below.

We will consider that the space interval [0, L] is divided into JV volumes of length
Ax. The finite-volume approach of the system (6.3) is given by (see previous chapters)

i_ (C /n + i _ u?) + J-(F?+l/2 - Ff_l/2) = 0, j = 1,.. .N. (6.4)

For the internal cells we have classically
U?+1), j = l,...N-l, (6.5)
where F(.,.) is some numerical flux function and, Uf is a constant vector over the volume
Kj =] z j-i/2) a; j+i/2[- The value of the function F(f/", U™+1) can be obtained by an exact
solution of the Riemann problem (Godunov) or an approximate Riemann Solver (Roe [14],
Osher [10], Colella-Glaz [1]). Another way to compute the numerical flux is the use of the
Flux Vector Spliting method (Perthame [11], Van Leer [18]), and the flux in this case is
expressed as

Since we use a finite-volume method, an evaluation of the boundary conditions in

terms of a numerical flux at the boundaries is needed. This fact was first recognized by
Godunov et al. [7].

i) For the case of a solid wall boundary, we use the well known notion of a "mirror
state" of the states defined at the extreme cells of the computational domain and hence,
the flux at the boundary is defined from the state at the cell edging the boundary and its
mirror state. For example, if the inflow boundary is a solid wall, we shall define below a
mirror state UQ of the state in the first cell £7" and, the flux across the boundary is given
by:

F?/2 = F(US,U?). (6.6)

ii) In the case of a free field boundary, we use the approach developped by Dubois-Le
Floch [2], [3] and which consists in the computation of a boundary state Ub obtained by

145
the solution of a partial Riemann problem (see below) and the flux across the boundary
(we consider only the inflow boundary) is given by:

F?/2 = F(Ub,U?). (6.7)

In order to compute the numerical boundary state Ub, recall that in usual fluid dynamics,
four cases should be analyzed from a physical point of view (see [2], [3], [4]):

Supersonic inflow. Ub is entirely given.

Subsonic inflow. The incoming flow is characterized by its static thermodynamical

properties, i.e.: total enthalpy and entropy, or by its temperature and its momentum
(or velocity). The state Ub is the unique state that can be connected to the interior
state by at most one contact discontinuity and a 3-wave.

Subsonic outflow. The static pressure is given. Ub is the unique state which can be
connected to the interior state by a 1-wave.

Supersonic outflow. No data is required and, Ub is identical to the interior state.

Here, for numerical purposes, the above mentionned 1- and 3-waves are always approxi-
mated using rarefaction-compression waves. This allows for the use of associated Riemann
invariants.

We will now discuss the different situations encountered in two-phase flow for a given
boundary (inflow or outflow), and the problems connected with the implementation of the
boundary conditions. Notice on the one hand that, in two-phase flow problems, it appears
that for physical reasons, the void fraction can be imposed only at inflow boundaries,
but not at outflow boundaries (see Wulff [8] pp. 85-238). This is an important point to
keep in mind. On the other hand, since the systems of equations considered here are
solved in conservation form, the boundary conditions have to be expressed as a function
of conservative variables, whose expressions involve the volume fractions.
The physical boundary conditions are derived in many cases from experimental set-ups
and are given in terms of phasic measurable quantities such as the primitive variables pk,
Uk, Pk- Hence, additional data or assumptions on the void fraction of physical origin have
to be added in order to define completely the numerical problem.

In two-phase flow applications, supersonic inflow (respectively outflow) is unusual;

hence, only subsonic inflow (respectively outflow) will be analysed here.

6.3 Perfect fluid approach

Consider the system of equations (6.2), i.e.:

146
 H (6.8)
dt
assume now that both phases are described by an equation of state of ideal gas type

Pk = (7* - l)/>*efc, (6.9)

where -yjt is a constant; pk is given in this case by

Pk = ilk - l)mkek. (6.10)

If we let qk represent the momentum qk = mkuk, then the flux is expressed by

U ^ ^ ^ ^ \ , (6.11)
and the Jacobian matrix A(Uk) = dF(Uk) of the system (6.8) is given by

0 1 0

Çk
A(Uk) = (3-7*)
mi mk

-(7*-1)4
where

A direct calculation gives the eigenvalues of A(Uk)

Ajt,2 = «it»

t,3 = «it + Cjt,

where ck = \Z{ykp~k)/mk is some numerical local speed of sound. One can remark that
ck = ck = \/{jkPk)/Pk, where ck is the phasic speed of sound associated to the 7-law.
The three right eigenvectors of A(Uk) are given by

t.i = (1, «it - cfc, Hk - ukck)T,

rk,3 = (1, «fc + cjt, ^

The three pairs of Riemann invariants are given by the table 6.1.

147
 Eigenvalues

2c fc
Riemann uk Uk uk-
invariants
Pk Pfc
Pk
m?

Tableau 6.1: Riemann invariants.

U k.b

x=0
Inflow boundary

Figure 6.1: Inflow boundary conditions.

148
6.3.1 Inflow boundary (z = 0)

Three cases must be discussed.

Case 1: the inflow is closed to the two phases.
As it is well known, the boundary conditions for the two systems of equations requires
zero velocities at the wall i.e. ujt = 0. Numerically (see above) we introduce the "mirror
states", (we consider only the case x = 0) which are defined by:

Pk,o = Pfc.ii

«fc,o= -Ufc.i, for k = v,l

Pkfi = Pfc.l •

Recall that by convention we have used the index 0 for the mirror state and 1 for the inner
state.
Now, at time tn, we have pv — Pi = P this leads in the equations above to only five
conditions on the boundary. In order to complete the number of boundary conditions,
data on the void fraction is required. If we make the hypothesis that <*o = <*i, the mirror
states are completely prescribed, and in term of conservative variables we have

qk,o =-Çk,i, for k =

the last equatity can be replaced by pk,o = pk,i •

Case 2: the inflow is closed to one of the two phases.
In this case, the void fraction a must be given. The boundary conditions of the phase for
which the inflow is closed are specified by a null inflow velocity. For the other phase, the
computation of the boundary state is given below.

To show more clearly the computation of the boundary conditions in this case, we
consider the problem described in chapter 3, namely the Faucet flow. In this problem (see
chapter 3), the physical variables to be imposed at the inflow boundary are, uv£ = 0 (the
inflow boundary is closed to the vapor), a/,6, iufi and /i/^. As noted above, in order to
derive the boundary conditions for the vapor phase, it suffices to determine the mirror
state UVtQ. This is easily achieved by setting:

UVlo = [mVfl, gv,o, EVto)

where pvJ0 = pvA, uv,0 - -uVti, pv,0 = pv,i and avfi = 1 - Qi,b-
For the liquid phase, the boundary state Uij, is computed as follows:

149
since the velocity of the vapor is null at the inflow, the numerical boundary flux FVtl/2 for
the vapor has the form
= F{UVf(hUv,i) = (0,Pv,6 (6.12)
where pv$ is the partial pressure of the vapor on the boundary. If we assume that p/^ = pVtb
at the boundary, the partial pressure of the liquid is given by

Pi,b - (6.13)

The remaining conservative variables are

Pi,b
7/ - 1 Kb '

Çl,b —

Case 3: the inflow is open to the two phases.

In this case also, the void fraction a at the boundary noted aj must be given. For simplicity,
we restrict ourselves to the case where aj, = 0 or 1, i.e., only one phase enters the
computational domain (say the vapor phase). In this case, we ignore the void fraction,
and only the dynamics of the pure-phase in question is considered. Assume the vapor
velocity uVjb and enthalpy hVyb are given at the boundary. Thus, the computation of the
boundary state Uv<b is obtained by the solution of a partial Riemann problem, which
reduces to the solution of the following system:

Pv,b

(6.14)
-V — 1
7v-l
lv —A

Pi Pv,l
By eliminating the intermediate variables and, after some calculations one obtains

Iv
Pv,b =

150
 "îv Pv,b
Pv,b =

Çv,b =

v,b— Pv,b\ 1
7v

Remark 6.1 Notice that, we have assumed that Uk,b are given and Uk,b > 0, for k =
v, / at the inflow boundary. For two-phase flow, other boundary conditions such that the
pressure and the enthalpies may be imposed instead of the velocities. In this case, during
a transient the sign of the velocities may change. Therefore, the number and the analysis
of the boundary conditions must be adapted following the problem to be solved.

6.3.2 Outflow boundary (x = L)

k,b

U k,N

x=L
Outflow boundary

Figure 6.2: Outflow boundary conditions.

Case 1; the outflow is closed to the two phases.

This case is analogous to the one described previously for the inflow boundary.
Case 2; the outflow is open to the two phases.
In this case, the void fraction cannot be imposed, because the void fraction information
propagates from the interior of the computational domain to its exterior. We therefore set
ftf, = O^r. (O.lO)

The pressure pf, at the boundary is given. Since the outflow is open to the two phases, the
phasic pressures are equal, and we have

Pv,b = Pl,b — Pb- (6.16)

151
The partial Riemann problems are solved with only 1-waves (one for each phase) issued
from the interior states Uk,N- The computation of the boundary states Uk,b is achieved by
solving:

Uk,N i Cfc ^ = Ukb

7Ck,N «fc,6 ^rcfc,<"
(6.17)
Pk,N _ P
Pk,b

From the second equation in (6.17), we have

~ - mk,N, (6.18)
Pk,N I
by using (6.15) we obtain
/ \ I/T*
mk,b = — rnk,N. (6.19)
\Pk,Nj
The velocities are deduced from the first equation in (6.17),

(6.20)
1
"*fc,/v y "**:,& j

and the specific internal energy is given by

This completes the computation of the states of the two phases.

6.4 Real fluid approach

The thermodynamic properties of a real fluid may be embodied in the relation

P = p(p,e), (6.22)

expressing the pressure p as a function of the density p and the specific internal energy.
In some cases, it is convenient to work with the internal energy per unit volume £ — pe
[9] rather than specific internal energy. In this case, the equation of state takes the form

P = p{p,e). (6.23)

Assume now that the partial pressure p\ in the system of equations (6.2) can be expressed
as a function of the partial density m* and the partial internal energy per unit volume
i = mkCki viz. as

Pk — Pk{rnkiêk)i (6.24)

152
this holds indeed if pk, the fc-phase pressure depends linearly on pk (case of thermally
perfect fluids); also, it is possible if ak is considered as a parameter.

In order to avoid errors in the computation of the Jacobian matrix, we express the
equation (6.24) in terms of all variables Uk, viz. asp* = pk(mk{Uk), £k{Uk)) • Furthermore,
we shall assume that the first derivatives of the partial pressure are available.

In terms of the partial internal energy per unit volume ëk and m*, the derivatives
of the partial pressure pk in (6.24) are defined as
£., I UPk \ /f. rtr\
\ O£k )

For an ideal fluid, Xk — 0, and ëk = {"/k — 1)-

We now construct the Jacobian matrix A(Uk) = dF(Uk)/dUk and find its eigenvalues and
right eigenvectors. Defining the partial momentum qk = rnkUk, the flux function can be
written as

F(Uk) = («,mk- f + + %- (6.26)

where

Pk = Pk{mk, Ek - (6.27)
2mk'
The expression of the Jacobian is given by
/ 0

A{Uk) = ~(2-«*)-T {2-kk)uk

t ~ - Hk)uk
where the total enthalpy, Hk, is defined by Hk = hk + u\/2, and hk = + Pk/pk is the
specific enthalpy. The three eigenvalues of A(Uk) are found to be
Ajt,l = Uk~ Cfc,

where c\ = \k + «Jfe/ifc is some numerical local speed of sound.

The corresponding right eigenvectors are given by
T
rk,i = (1, uk - ck, Hk -

nt,3= (l,uk + ck,Hk +

153
6.4.1 Inflow boundary (x = 0)

Only the last case will be discussed here. The treatment of the cases where the inflow
is closed to the two phases or only for one phase, is achieved in the same way as in the
perfect fluid case.
Case 3;
Let us put ourselves under the conditions stated in case 3 for the perfect fluid. Consider
here the liquid phase. To compute the boundary state £//,&, we use the Riemann invariants
for the liquid phase. Therefore we have,

ui,b = given,

hi tb — given,

U* = Ulfi,

(6.28)
PÎ ~ Pi,b,

dp

By eliminating the intermediate variables, the system (6.28) reduce to solving

dp
(6.29)
uifi = "U +

where the uknown is the liquid-phase pressure p/,i of the boundary state. This is a nonlinear
equation which can be writen under the form
f(xi) = 0 (6.30)

where, x/ = p;^ and

dx
f(xi) = ulA - m,b (6.31)

The nonlinear equation (6.30) is solved by an iterative method of Newton-Raphson type,

combined with the trapezoidale rule for the integration. For each pair (x,s), the corre-
sponding density and sound speed are obtained from the thermodynamic table [13].

The integral in (6.29) can be expressed in terms of the adiabatic exponent. Indeed,
the speed of sound of a fluid is defined by (see chapter 4)

c2=dp 21
P'

154
Thus we have
dp = c2dp along an isentrope.
Now, the pressure p decreases along a 3-rarefaction curve, so p also.
By adapting the method of Smoller [15], we introduce a parameter x by setting

x > 0 since p < Pi, where p\ is the pressure in the cell 1. Hence
dp = —pie~xdx = —pdx,
and
dp _ dp _ —pdx _ dx
p c2p jp 7 "
Integrating this last quantity gives

Pi
We turn next t o the integral in (6.29), we have

FrPi,b

'pi,\
dp _ /•*«.» dp

-dx
y/li

For sufficiently weak nonlinear boundary conditions, i.e., C//,i not too far from C//,b,
the variations of the adiabatic exponent can be neglected, and one can assumes that
ll = 7/,6 = 7/,i- So it comes

I — — — — ^ — ^ — ™ — / 1 / ' ^ 2^^ dx
Pl.i P'(Pi s)cl{Pi s
) \s=s,A ]j Pl,lîl,l JO

7/4 ~ 1
27/.1 Pl,b
' — 1
(7/,i - 1) V PM7/.1

and p^b, which is the uknown of the equation (6.29) is given by

27/.1
(7/4 - 1)
iiL + i
.(7U-1)0''1

155
6.4.2 Outflow boundary ( i = L)

case 1; the outflow is closed to the two phases.

This case is analogous to the one described previously for the inflow boundary.
case 2; the outflow is open to the two phases.
If a is constant (see case 2 of the outflow boundary in the case of perfect fluid), the Riemann
invariants are the same as those of the one-phase flow. Thus, under the assumption (6.15),
the entropies Sk and the quantities
k
rk-Uk jf °^k Apk

are constant along 1-rarefaction curves. By using (6.16), it follows that

&k,b = s
k,Ni

dp (6-32)

Using (pb, Sfc.b) as input thermodynamic variables in table [13] we get pk,b and hk,b- The
specific internal energy is then

ek,b = hk,b - — • (6.33)

Pk,b
The integral involved in the expression of the particle velocity is calculated numerically
by the trapezoidal rule. The densities pk and the speed of sound Ck are obtained form the
table in reference [13] for each pair of thermodynamic variables (pk, Sk), k = v, I.
Finally, the boundary states are

Uk,b = (mfc>6, çfc,6, Ek,b), for k = v, /, (6.34)

where

The integral in (6.32) can also be expressed in terms of the adiabatic exponent.
Hence, the pressure p and the density p decreases along a 1-rarefaction curve. To use the
method of Smoller [15], we introduce a parameter x by setting

— = e~ x < 1, therefore x > 0,

PN

where p^ is the pressure in the cell N. Thus we have

dp = —pwe~xdx = —pdx,

156
and
dp __ dp dx
P ~~ ÎP 7 '
Integrating this quantity gives

7
PN

Now, the integral in (6.32) is expressed as

rPb dp fP*> dp
Jpk.N

" Jo
-dx

For computational efficiency, if the states Uk,N and f4,(, are sufficiently close, the
integral in (6.32) can be calculated analytically by assuming jk,b = îk,N- We then obtain

r dp
,N Pk(p, s)ck(p, s) {s=3kN
Pk,N
,N7k,N
pk,N7k Jo
r

-1
~ 7k,N) V Pk,N7k,N

and finally we get

- 1
Pb
- 1

which complete the computation of the boundary states.

157
Bibliographie

[1] P. COLELLA & H. M. GLAZ, Efficient Solution Algorithms for the Riemann Problem
for Real Gases, J. Comput. Phys., 59, 264-289, 1985.
[2] F. DUBOIS, Conditions aux limites fortement non linéaires pour les équations d'Euler,
Cours CEA-EDF-INRIA, "Méthode de différences finies et équations hyperboliques",
organisé par P . L. LIONS, novembre 1988.
[3] F. DUBOIS & P . LE FLOCH, Boundary Conditions for Nonlinear Hyperbolic Systems
of Conservation Laws, Journal of Differential Equations, 71, 93-122, 1988.
[4] F. DUBOIS, Résolution numérique des équation d'Euler Multidimensionnelles, Cours
de DEA d'Analyse Numérique, Ecole Polytechnique, 1993.
[5] E. GODLEWSKI & P.-A. RAVIART, Hyperbolic systems of conservation laws, Ellipses,
Mathématiques & applications, vol. 3/4, 1990.
[6] E. GODLEWSKI & P.-A. RAVIART, Numerical approximation of hyperbolic systems
of conservation laws, AMS 118, Springer-Verlag, New-York, 1996.
[7] S. K. GODUNOV & COLL, Résolution Numériquedes Probème Multidimensionnels de
la Dynamique des Gaz, Edition MIR, Moscou, 1979, 1976.
[8] G.F. HEWITT, J.M. DELHAYE & N. ZUBER, Multiphase Science and Technology,
Volume 5 Hemisphere Publishing Corporation, 1990.

[9] J.-L. MONTAGNE, H. C. YEE & M. VINOKUR Comparaison study of high-resolution

shock-capturing schemes for real gas, NASA Technical Memorandum 100 004, July
1987, Ames Research center, Moffett Field, California.
[10] S. OSHER, Riemann Solvers, the Entropy Condition, and Differece Approximations,,
SIAM J. Numer. Anal., 21, n°2, 217-235, 1984.
[11] B. PERTHAME, Second-order Boltzmann schemes for compressible Euler equations
in one and two space dimenions , SIAM J. Numer. Anal., Vol.29, n° 1, pp. 1-19,
February 1992.
[12] V.H. RANSOM, D.L. HICKS, Hyperbolic Two-Pressure Models for Two-Phase Flow,
J. Comput. Phys., 53, 124-151, 1984.

158
[13] P. RASCLE & 0 . MORVANT, THETIS : Calcul par interpolation des fonctions ther-
modynamiques de fluides diphasiques, EDF HT-13/94/014, 1994.
[14] P.L. ROE, Approximate Riemann Solvers, Parameter Vectors, and Difference
Schemes, 3. Comput. Phys., 43, 357-372, 1981.
[15] J. SMOLLER, Shock Waves and Reaction-Diffusion Equations, Springer Verlag, New
York, 1983.

[16] H. STADTKE & R. HOLTBECKER, Numerical Simulation of Two-phase Flow Based on

Hyperbolic Flow Equations, Sixth International Topical Meeting on Nuclear Reactor
Thermal Hydraulics, NURETH-6, October 5-8, 1993, Grenoble France.
[17] H.B. STEWART k B. WENDROFF, Two-Phase Flow: Models and Methods, J. Comput.
Phys., 56, 363-409, 1984.
[18] B. VAN LEER, Flux-vector Splitting for the Euler Equations, Lecture Notes in Physics,
n°170, (Krause Ed), Springer, Berlin, 507-512, 1982.

159
NEXT PAGKS)
!®fft &LANK
Annexe A

The Perthame and Roe schemes

Consider the following Cauchy problem for the unknown U = U(x, t) with t > tn and
x eR
du | dF(U) = Q
dt di (A1)

where:
[/ = (m, mu, £ ) r is the vector of conservative variables,
F(U) = (mu, mu2 + p,(E + p)u) is the flux vector,
t/ n (x) is a given state at time tn,
p and £" are given by: p = (7 — l)me and £" = m(e + u2/2). 7 is a constant 1 < 7 < 3.
All the variables appearing above are defined in sections 2 and 3 of chapter 3. In this
appendix, p and p appearing in the Euler equations are replaced by m and p in order to
take into account the volume fraction of the considered phase.
In order to built first-order shemes, we assume that Un(x) is a piecewise constant over
cells Kj=]xj_1/2,xj+1/3[
Un(x) - U* = (m],m"u",E") forx € Kj. (A.2)
Uf is defined by

We wish to calculate U"+1, the approximate solution at time £n+1 = tn -\- At in cell Kj.

The generic formula for updating cell averages of the conserved quantities corre-
sponding to (A.I) is

161
In the next, we compute the fluxes •fT+i/2 associated to the Perthame and Roe schemes.

A.I The Perthame scheme

The Perthame scheme is one of the several numerical methods based on kinetic theory
which have been proposed for the computation of inviscid compressible flows. These meth-
ods use the Boltzmann or the Boltzmann-like equation of the kinetic theory of gases. They
are based on the connection between the Boltzmann equation and the Euler equations gov-
erning the dynamics of inviscid compressible flows.

In order to compute the fluxes F±(Uj) given in (3.33), we follow a similar course of
reasoning as that used by Perthame [15] in the ideal gas case.
By integrating over the time interval [tn,tn+1] we have

U(x,tn+1) ~ U(x,tn) = f dtU(x,t)dt

Jt"
The implementation of the upwind principle for flux splitting is performed at the Boltz-
mann level and then mapped to the Euler level. Thus, by using the lemma 3.1 which gives
the connection between the vector U and the moments of the Boltzmann function, we
have

f
Jt"
at(7(i,t)dt = f f
Jtn f |(l,t;,^)Tat/(x,t;,
JR I 2

= -jfB J\(v,v2,j)Tdxf(x
After that, we integrate over each cell Kj in order to find the approximation to the average
of the solution across the cell Kj

?+1-U?)Ax = fX3^/2(U(x,tn+1)-U(x,tn))dx

+(0,0, v)Tdxg{x, v, t)]dxdvdt

= ~f f[{v,v2,^-)T(f(xJ+1/2,v,t)-f(xj.1/2,v,t))

+(0,0, v)T{g(xj+1/2, v, t) - s(ij-i/2, v, t))]dvdt

162
The exact solutions to (3.27) for t >tn are given by

f(x,v,t) = f»{x-v(t-tn),v) (A.5)

n n
g(x,v,t) = g (x-v{t-t ),v) (A.6)

so it comes

(U»+1-U?)Ax = -

T n
-fn(*j-i/2 - «(* - tn),v)) + (0,0, v)T{gn(x
( j+1/2 - v{t - tn),v)

-9n(xj.1/2-v(t-tn),v))]dvdt (A.7)

fn and gn are proportional to x- Now, x has a bounded support (3.23), the v-integration is
"reduced therefore to v such that \(v—u)\/m/p\ < \ / 3 , i.e. M(£ — *n) < At(\u\ + y/p/m y/3),
and by the C.F.L. condition (3.31) we have \v\(t - tn) < Ax i.e.

< Xj±1/2 - V(t - tn) < Xj±if2 + AX

Thus the expressions of / " and gn reduces to

if u > 0

if v < 0

and

if v > 0
9n(xj+i/2-v{t-tn),v)=*

n n n
/ n (xj_!/ 2 - u( - < ), u) and <7 (xj_i/2 — u(t — t ), v) are expressed in the same way.
By setting

h(U?,v) =
A(0,0, v ] x((v - (A.8)

163
The expression (A.7) takes the form

U +1 U
?
3
= f-^l ÙLX
J
JV>O
h(U?,v)dv+ JV<O
J
f h(U?+1,v)dv J

- I h(U?_ltv)dv- f h{Uf,v)dv]
Jv>0 Jv<0

= UJ - ~ [ ( F " + (Uj) + Fn~ (Uj+1)) - (F" + ([/,_!) + F""(Uj))}

= UJ - —[(^(Uj, Uj+1) - FiUj-uUj)]

The numerical flux function Fn(Uj, Uj+i) plays the role of an average flux through £j + 1 / 2 ;
it is usullay noted by

and the final expression for updating celle averages of the conserved quantities is

At
(A.9)

for the numerical computations, we develop the expressions of F+ and F by using the
change of variable w = (v — u) y/m/p then, we find the expression given in (3.33)

m,jx{w)dw. (A.10)

A.2 The Roe scheme

The Roe scheme [18] is one of the most popular Riemann solvers for the gas dynamics.
The use of the Roe scheme was extended to other hyperbolic systems of conservation laws
by several authors, and also to hyperbolic nonconservative systems (see for example [23]
and [24]). Here, we will adapt the Roe method to solve our system of equations given in
(A.I).

164
The Jacobian matrix of this system is given by
/ 0 1 0 \

dF(U) ( 3 _ 7 )JL
A(U) = 2 ro2
dU
7 - 1 <73 <
2 m3 m w
m
where ç = mu, and the total enthalpy, H is defined by

In terms of the primitive variables, the matrix A(U) writes:

0 1 0 ^\

A(U) = 7-1

yu
\ 2
The eigenvalues of A(U) are given by

Ai = u - c,

A2 = u,
A3 = u + c,
where c = \/{jp)/m is the numerical speed of sound.
The right eigenvectors of A(U) are

r3=

The idea of Roe is, instead of solving the nonlinear problem (A.I) with the state
vector U™ at time tn given in (A.2), he introduce a local linearization of the system at
each interface Xj+1/21 and solve Riemann problem associated to the linear system defined
by

~U = 0 (A.11)
where the matrix A{Uf, U"+ï), known as the Roe averaged matrix, is some Jacobian matrix
depending on the states t/J1 and U™+1. For convenience, we denote Uf by ULand C/"+1 by
UR. The Roe averaged matrix is contructed to have the following properties:

165
 1. A(UL, UR).{UR - UL) = F{UR) - F{UL)

2. A(UL,UR) is diagonalizable with real eigenvalues

3. A(UL, UR) —-» A{U) smoothly as UL, UR —» U.

To derive an averaged Jacobian satisfying the conditions listed above, we introduce a

parameter vector W — {w\,wi,W3)T, and we construct two matrices B(WL,WR) and
C(WL, WR) such that
UR - UL = B(WL, WR).(WR - WL), (A.12)
and
F(UR) - F{UL) = C(WL, WR).(WR - WL). (A.13)
The Roe parameter vector in our case is defined by

u VÛ2 (A.14)

\ « / \
In terms of the component of W, U is given by
/ to? \

U = (A.15)

- p(W)
where
-1,
(A.16)

The matrices B and C satisfying the desired properties (A.12) and (A.13) are found to be
0 0 \

B(WL, WR) = 0
(A.17)
7 - 11^_
\ 7 7 7 /
and

7-1 _ 7
C(WL,WR) = (A.18)
7 7 7

0 w3 w2 )

166
where the overbars denote arithmetic means
1 .
V = 2 \VL + VR) • (A.19)

Finaly, the matrix defined by

A(WL,WR) =C(WL,WR).(B(WL,WR))-1 (A.20)

satisfies the conditions of Roe . Its expression is

0 1 0

A(UL,UR) = 7-1 (A.21)

2__lû 3 - ÛH H - (7 - l)û 2

where the averaged quantities are defined by

(A.22)
+ y/mR

The eigenvalues and right eigenvectors of A(UL, UR) are

Ai = û - c,

A2 = û,

X3-û + c,
and

r"i = (1, û — c, H — ûc) ,

r3=

where ë = ( 7 - f
The expression of the flux function for Roe's solver is now given by

?) + F(Uf+l)) - ^|

Since the system (A.ll) is linear U"+i — U™ may be decomposed as

(A.23)

167
and hence

(A-24)

= J2$i\Xi\r'i (A.25)

where the coeficients /?,• are given by

2c

2c
= Am- -rf.
c
(.) denote the jump in a quantity, and m is defined by

m =

168