18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

Utilisation d’une méthode de suivi d’interface mixte eulérienne/lagrangienne

pour les écoulements diphasiques

Guillaume Balmigere1, Stéphane Vincent2, Jean-Paul Caltagirone2 & Erick Meillot1

1
CEA Le Ripault
BP 16, 37000 Monts
2
Université de Bordeaux 1
Laboratoire TREFLE
16, avenue Pey-Berland, 33000 Pessac
balmigere@[Link]

Résumé :

Nous présentons ici une méthode hybride de suivi d’interface s’inspirant de l’aspect particulaire des
méthodes MAC et de la représentation volumique du champ de concentration des méthodes VOF. Les
mailles du domaine sont pourvues d’un nombre de particules variable transportées de façon
Lagrangienne à partir du champ de vitesse Eulérien interpolé. Ces particules sont porteuses
d’informations de phase qui seront projetées sur le maillage Eulérien à chaque pas de temps afin
d’actualiser le champ de fraction volumique. La validation de la méthode a fait l’objet de comparaisons
sur des cas tests de tache cisaillée avec un champ de vitesse analytique, ainsi que sur des cas de rupture
de barrage 2D et 3D confirmant les performances de la méthode lors du couplage avec la résolution des
équations de Navier-Stokes.

Abstract :

We developed a hybrid Eulerian/Lagrangian method to predict the motion of free surfaces in multiphase
flows. This is inspired from the particle based MAC methods and from the Eulerian volume approach of
VOF methods. Each cell of the computational domain is provided with particles which are advected in a
Lagrangian manner by the interpolated velocity field. These particles carry the phase information which
is projected at each time step to the Eulerian grid to build the Eulerian description of the volume fraction.
The method was able to predict large interface distortions with a good accuracy, as in the case of a
prescribed velocity field corresponding to the time-reversed single vortex flow. We assessed also very
good results when coupling the method with the resolution of the Navier-Stokes equations as
demonstrated by the simulation of 2D and 3D dambreak.

Mots-clefs: Advection methods; Volume of Fluid

1 Introduction

Les représentations numériques de l’interface entre deux fluides non miscibles se

différencient principalement par la manière dont le domaine d’étude est discrétisé. La première
catégorie n’utilise pas de maillage, notamment les méthodes SPH (Smoothed Particle
Hydrodynamics) pour lesquelles la phase advectée est représentée par un ensemble de particules
dont le mouvement est assuré exclusivement par un suivi particulaire (Monaghan (1994)). Ces
méthodes ont un principe purement Lagrangien, et sont généralement utilisées pour des
problèmes d’écoulements compressibles. La deuxième famille est dite des méthodes
Lagrangiennes. Celles-ci utilisent un maillage adaptatif non-structuré permettant l’imposition
rigoureuse des conditions aux limites à l’interface (Floryan & Rasmussen (1989)). Elles donnent
un rendu très fidèle de la forme des interfaces, mais leur utilisation devient difficile dès lors
qu’elles se brisent ou se scindent et leur programmation est complexe et coûteuse en 3D. On
leur oppose généralement les méthodes dites Eulériennes qui utilisent un maillage fixe. Ces
méthodes nécessitent l’introduction d’une variable supplémentaire et d’une équation pour suivre

1
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

son évolution. Parmi elles, on distingue la méthode Front-Tracking décrite par Unverdi &
Tryggvason (1992), les méthodes de suivi de lignes de niveaux Level Set (Osher & Fedkiw
(2001)) et les méthodes Eulériennes de suivi en volume sous différentes formes telles que la
méthode MAC (Marker And Cell) présentée par Raad & Bidoae (2005), ou les variantes de la
méthode Volume of Fluid (VOF-PLIC, VOF-TVD) (Rider & Kothe (1998)). Les limites des
méthodes Eulériennes sont atteintes des lors que les échelles caractéristiques de l’interface sont
inférieures à la taille de la maille. Se produisent alors des effets de non conservation de la masse
pour la méthode Level-Set, de fractionnement artificiel pour VOF-PLIC et de diffusion
numérique importante pour VOF-TVD. Des méthodes mixtes ont récemment vu le jour (Aulisa
et al. (2003), Ishii et al. (2006)), afin de combiner les propriétés de conservation de la masse des
méthodes VOF, avec la précision des méthodes Level-Set ou particulaires pour la description
des interfaces.

L’approche VOF-SM (VOF Sous-Maille) présentée ici s’inspire des caractéristiques

Lagrangiennes des méthodes MAC et de la méthodologie Eulérienne du suivi en volume VOF.
La dynamique de l’écoulement est décrite de façon Eulérienne, par la résolution du champ de
vitesse selon le modèle 1-fluide. Chaque maille de la grille contient des particules porteuses
d’informations physiques qui seront projetées sur le maillage afin d’obtenir une description
Eulérienne du champ de concentration. Nous souhaitons ainsi développer une méthode de suivi
d’interface Eulérienne, dotée d’une description sous-maille du milieu diphasique et permettant
d’appréhender statistiquement à l’échelle sous maille les phénomènes de rupture et de
reconnexion d’interfaces 3D.

2 Modèles physiques

Les équations de Navier-Stokes en formulation 1-fluide permettent la modélisation des

écoulements multiphasiques incompressibles isothermes :
∇V
. =0 (1)
⎛ ∂V
ρ⎜
∂
⎞
( )
+ (V.∇)V ⎟ = ρg − ∇p + ∇. µ ( ∇V + ∇T V ) + FTS (2)
⎝ t ⎠
∂C
+ ( V.∇ )C = 0 (3)
∂t
La méthodologie numérique utilisée est celle développée dans la bibliothèque de calcul
scientifique Aquilon. Les équations sont discrétisées sur une grille décalée cartésienne par la
méthode des volumes finis. Le couplage vitesse/pression est assuré par un algorithme de
minimisation basé sur une méthode de Lagrangien augmenté. Les équations sont discrétisées en
temps implicitement par un schéma d’Euler d’ordre 1. Le solveur est un gradient conjugué
BICG-Stab II préconditionné par une méthode de type ILU. Pour plus de détails, on se réfèrera
par exemple à Vincent & Caltagirone (2000)

3 Méthode mixte

La méthode de suivi d’interface présentée ici utilise la notion de marqueur des méthodes
MAC et la représentation Eulérienne de type 1-fluide basée sur le calcul de propriétés physiques
moyennes. Initialement, des particules dont le nombre par maille est à définir sont positionnées
dans chaque maille du domaine de calcul. Elles sont porteuses d’une information logique
relative à leur positionnement dans l’un ou l’autre des fluides en présence et d’un poids
dépendant du nombre de particules dans la maille.
La résolution des équations de Navier-Stokes fournit la connaissance du champ de vitesse
en chaque point du domaine Eulérien. Celui-ci est interpolé au niveau des particules qui sont
ensuite advectées de façon Lagrangienne. Leur nouvelle position étant connue, les informations

2
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

qu’elles portent sont projetées par des lois statistiques conservatives sur le maillage Eulérien
produisant ainsi une description Eulérienne de la fraction volumique au pas de temps suivant.
La figure 1 schématise la boucle itérative de la méthode VOF-SM.

FIG. 1 – Algorithme de la méthode VOF-SM

4 Résultats

4.1 Tache circulaire dans un champ de vitesse déformant

La précision de la méthode VOF-SM a été éprouvée sur différents cas tests. Le premier est
celui d’une tache cisaillée par un champ analytique, dans la configuration présentée par Lopez
et al. (2005). Initialement, le centre d’un cercle de fluide de rayon 0.15 m est positionné au
point de coordonnées (0.5 ; 0.75) dans un carré de coté 1 m. Celui-ci est soumis au champ de
vitesse dérivé de la fonction courant suivante : Ψ ( x, y ) = 1 / π . sin 2 (πx) sin 2 (πy ) cos(πt / T ) où
T = 8s dans les simulations présentées ici. La déformation maximale du cercle se produit à
t = T / 2 (Fig. 3) et le cercle retrouve sa position initiale lorsque t = T (Fig. 2). La solution à
t = T est pratiquement confondue avec la solution théorique pour un maillage de 64 x 64
comme en témoigne la figure 2. La méthode VOF-SM apporte une amélioration importante par
rapport aux approches VOF-PLIC et VOF-TVD en s’affranchissant des problèmes de
dissipation numérique et de déformation importante. Les niveaux d’erreur en terme de norme L1
ont été comparés à ceux présentés par Lopez et al. (2005), basés sur une méthode VOF-PLIC
améliorée. Nous avons effectué une étude de convergence sur des maillages de 16 n, n=1..6
points par direction, avec 2, 3, 4 particules par direction (Fig. 4). La méthode proposée présente
des niveaux d’erreur et des ordres de convergence très satisfaisants selon les critères
mathématiques choisis par Lopez. L’ordre atteint est en effet compris entre 2,5 et 3.

Nous avons également considéré l’écoulement potentiel généré par la présence d’une
plaque rigide perpendiculaire à la direction principale de l’écoulement, donné par la fonction
courant Ψ ( x, y ) = 10( x − x a )( y − y a ) où ( x a , y a ) sont les coordonnées du point d’arrêt. Ce cas
illustre la nécessité d’introduire une gestion des particules des lors que le domaine possède des
frontières ouvertes. Le transport des particules entraîne en effet une raréfaction ou une
disparition de celles-ci dans certaines zones, rendant impossible la gestion de problèmes de type
injecteur notamment. Les figures 5 et 6 illustrent la mise en place de la gestion des particules
sortantes et leur redistribution dans les zones raréfiées sur ce cas de déformation de tache à
t=0.11s. La figure 7 présente l’équivalent en terme de contour C=0.5 de fraction volumique
après l’étape de projection (Fig. 1).

4.2 Rupture de barrage 2D sur fond mouillé

3
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

De nombreux travaux expérimentaux et théoriques traitent de la rupture de barrage. Ils

confèrent à ce cas tout son intérêt quand à la validation de la méthode proposée pour la
résolution des écoulements à surface libre air/eau. Stansby et al. (1998) se sont intéressés au cas
de la rupture de barrage sur fond mouillé faisant intervenir des déformations importantes de
l’interface ainsi que des phénomènes de « splashing » et de rupture. Nous avons reproduit cette
étude en deux dimensions sur un domaine de 1,2 x 0,2 m discrétisé sur une grille de 600 x 100.
Les résultats de la méthode VOF-SM sont comparés à ceux obtenus avec une méthode VOF-
PLIC classique et aux résultats expérimentaux (Fig. 8).
Initialement, nous considérons une colonne d’eau au repos de longueur 0.5m et de hauteur
0.1m. En aval, le sol est recouvert d’une épaisseur de 0.01m d’eau. Le début du calcul coïncide
donc avec la disparition artificielle de la séparation entre ces deux zones, l’eau est alors mise en
mouvement par la force de gravité. Une onde de choc se propage en aval et donne naissance à
un jet déferlant. La surface libre obtenue est comparée aux données expérimentales aux instants
t=0.24s et t=0.4s. Les résultats sont en bon accord avec les données expérimentales tant au
niveau de la hauteur d’eau que de la vitesse de propagation du jet. La méthode donne également
un rendu fidèle des structures décrites par Stansby et al. (1998) liées à l’emprisonnement d’air.

4.3 Rupture de barrage 3D sur fond sec avec obstacle

Nous avons poursuivi le cas de la rupture de barrage dans une configuration 3D de fond
sec en aval, et en présence d’un obstacle. Initialement, une colonne d’eau au repos de
dimensions 1,22 x 0,55 x 1 m se situe en amont d’une zone sèche où se trouve un obstacle
solide. Le tout est inclus dans un domaine fermé. Le domaine de calcul de dimensions 3,22 x 1
x 1 m est discrétisé sur un maillage peu raffiné de 100 x 36 x 36. La figure 9 présente une vue
3D de la surface libre à t=0.7s. La figure 10 confronte les données expérimentales issues de
Kleefsman et al. (2006) et l’évolution de la hauteur d’eau au cours du temps aux points de
coordonnées (1;0) et (2.7;0). L’évolution macroscopique de la surface libre est correctement
décrite par la méthode VOF-SM.

FIG. 2 – Tache cisaillée,

contour C=0.5, t=T

FIG. 3 – Tache FIG. 4 – Courbe de convergence basée sur le calcul de l’erreur L1

cisaillée, contour C=0.5,
t=T/2

4
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

FIG. 5 – Champ de particules FIG. 6 – Champ de particules FIG. 7 – Contour C=0,5 de

à t=0,11s sans redistribution à t=0,11s avec redistribution fraction volumique à t=0,11s

(a)

(b)

FIG. 8 – Rupture de barrage 2D sur fond mouillé. t=0.24s (a), t=0.4s (b)

FIG. 9 – Rupture de barrage 3D sur fond mouillé avec obstacle. Représentation de la surface
libre au temps t=0.7s.

5
18ème Congrès Français de Mécanique Grenoble, 27-31 août 2007

(a) (b)

FIG. 10 – Rupture de barrage 3D sur fond mouillé avec obstacle. Evolution de la hauteur d’eau
au cours du temps, dans le plan y=0 m aux positions x= 1 m (a), x=2,7 m (b)

5 Conclusions

Nous avons développé une méthode de suivi d’interface hybride Eulérienne/Lagrangienne

d’ordre élevé et plus précise que les méthodes VOF rencontrées, donnant non seulement une
description macroscopique de l’interface mais aussi une description sous-maille. Celle-ci est
utilisable en 3D et une optimisation de l’algorithme de redistribution des particules devrait
permettre de diminuer le coût lié à la description particulaire du milieu diphasique. Un effort
important doit maintenant être fait sur la modélisation des phénomènes sous-maille.

Référence

Monaghan J.J. 1994 Simulating free surface flows with SPH, J. Comput. Phys. 110, 399-406
Floryan, J.M. & Rasmussen, H. 1989 Numerical method for viscous flows with moving
boundaries. Appl. Mech, Rev. 42(12), 323-341
Unverdi, S.O. & Tryggvason, G. 1992 A front tracking method for viscous, incompressible,
multi-fluid flows J. Comput. Phys. 100, 25-37
Osher, S. & Fedkiw, R. 2001 Level set methods: an overview and some recent results. J.
Comput. Phys. 169, 463-502
Raad, P.E. & Bidoae, R. 2005 The three-dimensional Eulerian-Lagrangian marker and micro
cell method for the simulation of free surface flows. J. Comput. Phys. 203, 668-669
Rider, W.J. & Kothe, D.B. 1998 Reconstructing volume tracking. J. Comput. Phys. 141, 112-
152
Aulisa, E., Manservisi, S. & Scardovelli, R. 2003 A mixed markers and volume-of-fluid method
for the reconstruction and advection of interfaces in two-phase and free-boundary flows, J.
Comput. Phys. 188, 611-639
Ishii, E., Ishikawa, T. & Tanabe Y. 2006 Hybrid Particle/Grid Method for predicting motion of
micro-and macrofree surfaces, Journal of Fluids Engineering 128, 921-930
Vincent, S. & Caltagirone, J.-P. 2000, One cell local multigrid method for solving unsteady
incompressible multi-phase flows, J. Comput . Phys. 163, 172-215
Lopez, J., Hernandez, J., Gomez, P. & Faura, F. 2005, An improved PLIC-VOF method for
tracking thin fluid structures in incompressible two-phase flows, J. Comput . Phys. 208, 51-74
Stansby, P.K., Chegini, A. & Barnes, T.C. 1998 The initial stages of dam-break flow, J. Fluid
Mech. 374, 407-424
Kleefsman, [Link]., G. Fekken, G., Veldman, A.E.P., Iwanowski, B. & Buchner, B. 2006 A
volume of fluid bades simulation method for wave impact problems, J. Comput. Phys. 206, 363-
393