Grenier These
Grenier These
Grenier These
École Doctorale
Sciences Pour l’Ingénieur, Géosciences & Architecture
Thèse de Doctorat
Jury
Étrange coutume que d’ouvrir ce manuscrit de thèse sur la partie qui vient clôturer
effectivement ces 3 années de travail. Mais tellement indispensable pour remettre en pers-
pective l’état d’esprit qui a abouti à ce résultat.
Je remercie aussi parmi les personnes de Saipem (l’entreprise qui m’a financé au cours
de cette thèse Cifre) tout particulièrement Jacques Ruer. Il a su me faire confiance, lire
et écouter patiemment mes rapports et exposés un peu techniques sur la résolution de la
partie hyperbolique par solveur de Riemann du modèle à fraction de volume. . . Il m’a aussi
laissé une grande liberté de recherche en me permettant de rester au sein du laboratoire.
Sont à remercier aussi Jean-Marc Bonnissel, puis Nicolas Butin, qui ont suivi les dé-
veloppements techniques et qui ont eux aussi ingurgité de la SPH et toutes ses subtilités.
D’une base initialement purement Nantaise, les travaux m’ont amenés à explorer
d’autres variantes de la SPH développées ailleurs. C’est ainsi qu’après avoir éprouvé
quelques difficultés sur la méthode à fraction de volume (quoi ? quatre mois sur ce MUSCL ?
broutille ! !), je suis allé rencontrer Jean-Paul Vila à l’Insa Toulouse. Je tiens à le remercier
de nous avoir accueilli, David et moi, ces quelques jours et d’avoir pris le temps de regarder
en détail, les mains dans le cambouis, les raisons de mon blocage sur cette méthode.
Parallèlement à cela, mon exploration m’a aussi permis de collaborer avec Andrea
Colagrossi, de l’Insean à Rome. Bien que je n’ai jamais foulé le sol de la cité aux 7 col-
lines (enfin celle transalpine ; la provencale Nı̂mes ne m’est pas inconnue !), les nombreux
échanges avec lui ont été très fructueux. L’implémentation de sa méthode, mathématique-
ment moins complexe que celle qui me donnait un peu plus de fil à retordre, m’a permis
d’avoir assez rapidement des résultats et de dissiper, à moins d’un an de la fin de thèse,
les nuages qui en assombrissaient le dénouement. Je suis particulièrement reconnaissant
d’avoir aussi travaillé avec lui.
Je remercie aussi Françis Leboeuf d’avoir bien voulu être rapporteur de mon manuscrit
de thèse et d’avoir fait d’excellentes remarques en tant que spécialiste de la mécanique
des fluides compressibles. Mes remerciements vont aussi à Loı̈c Boudet et Yves Lecointe
d’avoir voulu, respectivement, participer et présider à mon jury de thèse.
Enfin, je tiens à remercier beaucoup d’autres personnes, qui sont intervenues moins
formellement dans le cadre de ma thèse (ils ne figurent pas dans la bibliographie, ni sur
la couverture du manuscrit) mais n’en demeurent pas moins indispensables à la réussite
de celle-ci.
Il y a tout d’abord ma famille qui m’a vu émigrer dans le nord (vous savez au dessus
d’une ligne Bordeaux-Lyon. . .), ne comprenant pas pourquoi j’échangeais le soleil nı̂mois
contre la pluie bretonne. Pendant trois années, ils ont patiemment essayé de comprendre,
à travers mes récits un peu obscurs sur la SPH, mes histoires de sorties voiles ou vtt ou
nocturnes. . . Puis un beau matin de février, arrivant sous le soleil dans la cité des ducs,
écoutant ma soutenance pieusement et s’émerveillant devant la houle magique de Maı̂tre
Fébulle, ils ont vu cette fabuleuse vie nantaise.
Cette thèse s’est aussi bien déroulée grâce à la très bonne ambiance de travail au
labo, à ses collègues et sa dynamique équipe de doctorants. Il y a tout d’abord Lionel
Gentaz qui m’accueilli en stage de Master (avec Alain Clément), toujours partant pour
le badminton et le vtt. Pierre Ferrant qui a suivi de loin, avec un regard non SPH, ma
thèse. Anne et Thuy pour leur gentillesse.
Et puis tous les jeunes, qui sont partis, revenus, restés, arrivés, titularisés. Qui n’ont
jamais été croisés qu’au deuxième étage du bâtiment D, à la cafét’, ou revus jusque
sur la terrasse du Remorqueur. Pèle-mêle : Mathieu Do. (le père de la SPH au LMF),
Yannoux (ça bouine ?), Micky, Jsex, Romulus Hélène et Velcro (t’as faim ?), Babar et
Elise (Vélocampus Powa !), Alban, PM (le pécheur de loup), JB, Guigs, Aurélien le DAH !,
Adrien Hellfest en kayak, Joss, OG, Ally, Matt Du. sailor master, Elich’, Hakim et Vonvon
LMF surf crew, Matt DL (on se disait avec Charles), Charles donc chess master, Clément,
Chadi, Bruno, Auré D., Guillaume F., Ilyes, Lionel H., Wawann.
Et surtout les 3 aı̂nés qui m’ont recueilli comme pauvre stagiaire de Master, perdu
seul en D01 et dans Nantes. Avec qui j’ai découvert Nantes (de nuit surtout...) : Duke
Ducro, Pierrot a.k.a. Rdr et Ouss.
Une mention toute spéciale à mon cobureau, Fébo, le normalien le plus cool jamais
vu, le roi de la vague (en bassin) et du vent en voile légère.
A Saipem, je tiens aussi à remercier Xavier Riou, docteur pleinement heureux dans
une boite d’ingénierie. Merci à ses conseils de début de thèse. Je décerne aussi une mention
spéciale à la TEDE young team, avec ses nombreux jeunes ingénieurs et stagiaires, qui ont
vu débarquer le temps d’une présentation un extraterrestre les inondant de numérique et
de SPH pour aussitôt repartir dans son labo.
Je remercie aussi Jérome et Emilie pour m’avoir permis de m’initier à l’enseignement
durant ma thèse.
Et finalement, je fais un clin d’oeil à ceux qui n’ont pas travaillé directement avec moi,
mais qui ont été plus ou moins au courant de ma thèse ces années nantaises : Lolotte,
Pablo et Lilie, Oliv’, Claire, Céd J80, Ingrid, Justine. Et les meilleurs : Palpak, Guigui et
Zaza, Tibo et Eve.
Table des matières
Introduction 1
i
TABLE DES MATIÈRES
ii
TABLE DES MATIÈRES
Conclusion 153
Bibliographie 156
iii
Introduction
1
INTRODUCTION
foré, et grâce, entre autres, à la pression élevée régnant dans la roche, le réservoir expulse
naturellement l’huile par le puits. Cette récupération primaire est cependant limitée dans
le temps et par conséquent ne permet pas d’extraire toute la quantité contenue dans le
réservoir. Pour augmenter celle-ci on a donc recours à une récupération assistée : par
l’usage de pompes sur terre et/ou par injection de fluides dans le réservoir sous la poche
d’huile (pour « pousser » celle-ci). C’est cette dernière technique qui est employée en
production en mer grâce à la disponibilité de l’eau. Cependant un réservoir géologique
n’est pas une entité parfaitement étanche, et la chute de pression d’une part et l’injection
d’eau d’autre part, provoquent la présence d’eau en sortie, dans le puits de production.
Or, la proportion de cette eau mélangée à l’huile ne cesse de croı̂tre au cours de la vie de
production d’un puits : de traces à son début, la part d’eau peut atteindre 99% du volume
de mélange extrait en fin de vie (proportion atteinte suivant les conditions économiques
du lieu de production). Ce mélange n’étant pas valorisable directement, il est nécessaire
de séparer l’eau de l’huile : on peut alors récupérer d’un côté l’huile pour la revendre
et de l’autre l’eau pour la ré-injecter dans le réservoir pour en améliorer le rendement
(sous réserve de pureté vis-à-vis de la porosité de la roche). Ce processus de séparation
est fait dans des séparateurs gravitaires (par décantation) ou dans des hydrocyclones (par
centrifugation) suivant les conditions de mélange. Ces systèmes sont actuellement disposés
en surface sur les plate-formes pétrolières.
Or, il a été vu plus haut que les exploitations en mer doivent se faire dans des profon-
deurs d’eau toujours plus grandes (> 1000m). Parmi toutes les contraintes supplémentaires
que cet environnement extrême apporte, certaines ont un impact direct sur le mélange
eau-huile, et d’autres sur les infrastructures, par exemple les séparateurs. En effet, comme
cela a été mentionné précédemment, la quantité d’eau présente dans le mélange peut deve-
nir prépondérante : se pose alors le problème de la pertinence économique et énergétique
de remonter à la surface un fluide non valorisable sur de grandes hauteurs pour l’éva-
cuer ensuite (ou le redescendre pour le ré-injecter). Cela a des répercutions évidentes en
terme de dimensionnement du réseau de transport. De plus, à grande profondeur règnent
de très faibles températures (de l’ordre du degré Celsius) qui contrastent avec celles de
l’huile sortant du réservoir (de l’ordre de quelques dizaines de degrés). De par sa com-
position chimique complexe et ses propriétés, le mélange eau-huile-gaz peut réagir à ce
choc thermique par la formation de composés solides (les hydrates), qui peuvent le cas
échéant boucher les conduites. Bien qu’il existe des procédés pour prévenir (par isolation
ou chimiquement) ces dépôts, leur déploiement dans ces conditions peut devenir pénalisant
économiquement.
Les problématiques amenées par l’extraction d’un mélange eau-huile à grande pro-
fondeur sont donc posées. La séparation apparaı̂t comme pouvant être un point clé d’une
exploitation viable si elle peut être effectuée au plus près de la tête de puits du gisement. Si
théoriquement, rien n’empêche de faire de la décantation avec des conditions thermiques
et mécaniques assez sévères, cela impose pratiquement une conception différente prenant
en compte toutes ces contraintes. Vient s’ajouter aussi un facteur contraignant supplé-
mentaire qui est la difficulté des opérations à ces grandes profondeurs : toute intervention
de l’homme est alors relativement complexe et sans commune mesure avec ce qui peut
se passer sur une plate-forme, que ce soit pour la mise en œuvre initiale du séparateur,
son exploitation ou sa maintenance. Toutes ces contraintes font qu’un séparateur doit être
conçu de façon la plus optimale possible tout en respectant une efficacité donnée (c’est-à-
dire la capacité à séparer correctement les deux fluides ; celle-ci est quantifiée en mesurant
2
INTRODUCTION
Une fois déterminés les modèles physiques à étudier, il est nécessaire de les valider
sur le cas d’étude choisi ou dans un premier temps sur des configurations simplifiées (par
exemple l’écoulement d’une bulle isolée). Pour effectuer cette démarche de validation des
modèles s’offrent au physicien deux choix : soit la comparaison directe de ses modèles à
la réalité, soit une étude virtuelle avec des outils de simulation informatique. La première
option, l’expérimentation, est la plus évidente : la confrontation au réel a longtemps été la
seule façon de confirmer la viabilité des modèles. Mais la lourdeur de mise en place d’essais
dans certaines configurations ainsi que l’émergence de nouveaux moyens technologiques
ont donné plus d’importance à un nouveau moyen de validation : la simulation par des
méthodes numériques. L’avènement de l’informatique a ainsi permis de faire résoudre par
des machines les équations régissant certains phénomènes physiques (faisant intervenir
des notions mathématiques compliquées à résoudre manuellement) et de confronter leur
solution aux hypothèses physiques. Et l’essor de ces méthodes de résolution a été portée
par la rapidité grandissantes des ordinateurs à fournir les simulations et par le faible
coût matériel de celles-ci. Cependant, il ne faut pas perdre de vue que l’expérimentation
reste indispensable à la validation de la réponse numérique ou pour valider les concepts
physiques hors de portée de la simulation actuellement.
Ayant fait le choix de la modélisation numérique, et parmi les outils développés dans
l’équipe Hydrodynamique et Génie Océanique du Laboratoire, deux d’entre eux sont po-
tentiellement susceptibles de nous permettre d’étudier notre problématique. Le premier,
VOFNI, est basé sur la méthode des Volumes Finis et une méthode de capture d’interface
(Volume Of Fluid). Il s’agit d’une méthode reposant sur des maillages eulériens, dont les
limites de la manipulation sont difficiles à repousser lorsqu’on étudie des grands mouve-
ments de corps au sein du fluide (tels que le déplacement d’un navire dans de la houle,
sujet au cœur de la thématique de l’équipe). C’est pour cette raison qu’un deuxième outil
a été développé : la méthode SPH dont l’une des principales caractéristiques est l’absence
d’utilisation de maillage. Cette méthode a été initialement développée à la fin des années
70 dans le domaine de l’astrophysique, avant d’être appliquée aux problèmes d’hydrody-
namique à surface libre au milieu des années 90. Quatre doctorants (Doring (2005), Oger
(2006), Deuff (2007) et Guilcher (2008)) et plusieurs stagiaires de Master se sont succédés
dans l’équipe pour développer cette méthode et prolonger son application à travers dif-
3
INTRODUCTION
Dans un premier temps, il a consisté à étudier les formulations utilisées dans un en-
vironnement monofluide pour évaluer leurs capacités à étudier des écoulements bifluides.
Sur cette base, des formulations multifluides ont été proposées. Elles résultent de discus-
sions et de travaux menés en collaboration avec Andrea Colagrossi (INSEAN Rome) d’une
part, et avec Jean-Paul Vila (INSA Toulouse) d’autre part. Elles sont présentées dans les
deuxième et troisième chapitres de la première partie.
Puis, ayant l’objectif d’ajouter les phénomènes physiques d’effets de viscosité et de
tension superficielle aux formulations choisies, il en a été étudié différentes modélisations.
Celles-ci ont été ensuite validées sur des cas simples ne mettant en jeu que ces effets. Ces
aspects sont détaillés dans la deuxième partie.
Dans la troisième partie (chapitres 1 à 4) sont exposés des cas-tests de validation plus
complets, mettant en jeu de façon progressive tous les éléments étudiés précédemment.
Les résultats des formulations SPH sont confrontés, suivant les cas, à des résultats expéri-
mentaux ou analytiques, ou lorsque ceux-ci ne sont pas disponible à cause de la complexité
de l’écoulement, à des résultats issus d’autres méthodes numériques dont la validation a
été effectuée par ailleurs.
Dans le dernier chapitre, une configuration très simplifiée d’un séparateur gravitaire
est présentée afin de démontrer la capacité de la méthode SPH proposée à traiter ce genre
de problématique.
Enfin, on réalise un bilan de ces travaux en esquissant des pistes de développements
futurs.
4
Première partie
5
Chapitre I.1
Les écoulements envisagés font intervenir des fluides incompressibles avec des effets
visqueux et de tension superficielle. Nous allons rappeler les équations instationnaires
régissant ces mouvements de fluides que nous considérons compressibles : les équations de
Navier-Stokes.
∂ΩS
∂ΩSL
∂ΩI xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Y
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
X
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
7
CHAPITRE I.1. LES ÉQUATIONS PHYSIQUES
où μ est la viscosité dynamique et λ est la viscosité de volume du fluide qui vérifient
l’hypothèse de Stokes (2μ+3λ = 0). Dans les cas où μ est constant et le fluide incompressible
⃗ ⋅ u⃗) = 0), la divergence du tenseur prend la forme :
(i.e. (∇
pK = p(ρK ).
L’équation d’état couramment choisie en SPH pour modéliser des fluides quasi-incompressibles
est l’équation de Tait :
c20 ρ0K ρK
γK
pK = K [( ) − 1] + p0K
γK ρ0K
où γK est le coefficient polytropique du fluide considéré, ρ0K et p0K ses masse volumique et
pression de référence et c0K est la vitesse du son lorsque ρK = ρ0K . Cette dernière vitesse
permet de régler le degré de compressibilité du fluide.
8
I.1.4. VITESSE DU SON ET PSEUDO-COMPRESSIBILITÉ
Dρ∗ u⃗∗ ∗ 1 1 ⃗∗ 1
∗ = −div(p I) + divTv∗ + S et p∗ = (ρ∗ γ − 1).
Dt Re F r2 γM a2
∂pK
c2K = .
∂ρK
9
CHAPITRE I.1. LES ÉQUATIONS PHYSIQUES
ce qui se traduit par une limitation du nombre de Mach : M a < 0, 1. Pour nous assurer
de cette contrainte, tout en ayant une dynamique de l’écoulement imposée (caractérisée
par sa vitesse V0 ), le paramètre sur lequel il nous faut intervenir est la vitesse du son
caractéristique c0 . La vitesse du son utilisée numériquement pour les calculs est donc
totalement arbitraire (tout en respectant la condition énoncées ci-dessus).
Pour résumer, pour modéliser des écoulements incompressibles nous effectuons la ré-
solution d’un système d’équations compressibles avec une vitesse du son distincte de la
vitesse du son réelle au sein du fluide compressible.
1 ∂ρ
χT = )
ρ ∂p T
et qui vaut, pour une loi d’état linéaire, χ−1 T = ρc . L’utilisation d’une vitesse du son
2
artificielle (telle que justifiée ci-dessus) introduit donc une compressibilité artificielle que
nous contrôlons via le nombre de Mach.
Pour les écoulements bifluides, les vitesses du son dans chacun des fluides peuvent
être choisies de manière indépendantes pourvu qu’on vérifie M a < 0, 1 pour chaque fluide.
Or, d’après Colagrossi et Landrini (2003), les coefficients de compressibilité ne peuvent
être indépendants pour pouvoir opérer dans la même gamme de pressions. Il en déduit
alors : √
c20X ρ0X c20Y ρ0Y γX ρ0Y
= soit c0X = c0Y .
γX γY γY ρ0X
Bien que d’un point de vue physique deux fluides peuvent coexister avec des compressibi-
lités différentes, notre expérience numérique tend à considérer l’hypothèse de Colagrossi
importante pour des simulations stables pour les deux fluides.
10
I.1.6. CONDITIONS AUX LIMITES
Cette relation traduit le fait que pour une interface courbe en équilibre, les forces de
tension sont compensées par une différence de pression de part et d’autre de cette interface.
Le lien de cette relation avec l’équation de quantité de mouvement (I.1.1) sera fait
ultérieurement.
Dans cette étude les effets de capillarité (liés aux effets de tension superficielle en pré-
sence d’une paroi solide) seront éludés. A priori on suppose que, dans les écoulements étu-
diés, soit les interfaces ne rejoignent pas une paroi (avec présence d’un angle de contact),
soit ce contact existe mais les effets sont négligeables. Ainsi près d’une paroi solide les
effets capillaires se traduisent par un ménisque de mouillage dont la longueur capillaire lc
peut être évaluée par la relation suivante :
√
σ
lc =
ρg
si les fluides sont soumis à l’accélération de la gravité g. Cette longueur (de l’ordre du
millimètre pour un domaine air-eau sous pesanteur terrestre) est faible devant les dimen-
sions caractéristiques des écoulements étudiés où une interface rejoint une paroi solide (de
l’ordre du mètre).
⃗ = u⃗Y ⋅ n
u⃗X ⋅ n ⃗.
11
CHAPITRE I.1. LES ÉQUATIONS PHYSIQUES
Cette condition peut aussi se voir comme une condition de non-pénétration à l’interface :
chaque fluide ne peut se mélanger à l’autre (cette interprétation s’étend aussi au cas des
parois ∂ΩS que le fluide ne peut pénétrer). Le seul degré de liberté restant sur la vitesse
est celui sur la vitesse tangentielle à l’interface.
Sur l’interface ∂ΩI (respectivement sur la surface libre ∂ΩSL ), cette condition cinéma-
tique peut se reformuler de la manière suivante : tout point matériellement sur l’interface
(resp. la surface libre) doit y rester.
u⃗X = u⃗Y .
Cette condition permet le développement d’une couche limite au niveau de l’interface ∂ΩI
entre les deux fluides ou au voisinage d’une paroi ∂ΩS .
12
Chapitre I.2
Dans ce chapitre on présente d’abord la SPH telle qu’elle a été formulée par Monaghan
et telle qu’elle a été utilisée dans le contexte du laboratoire Doring (2005), Oger (2006).
On s’attardera sur la première tentative d’extension aux écoulements bifluides faite par
Oger. Il sera mis en évidence les défauts de cette formulation. Ceux-ci sont à l’origine
d’une réflexion sur la modélisation des écoulements multifluides qui s’est traduite par une
reformulation de la méthode SPH (travail mené avec la collaboration Andrea Colagrossi
de l’INSEAN). On présentera cette reformulation dans un deuxième temps.
On a donc égalité entre la valeur du champ f au point r⃗ et sa valeur interpolée < f >.
Cependant, pour une question de dérivabilité de la fonction test (de Dirac), on remplace
celle-ci par une fonction de forme (couramment appelée noyau) plus régulière. On introduit
alors la longueur de lissage h qui est une mesure de l’étalement du noyau autour du point
d’intérêt. Cette grandeur doit tendre vers 0 de manière à faire tendre la fonction de forme
vers un Dirac. Nous perdons donc l’égalité entre le valeur du champ et son interpolation :
< f (⃗
r) >h = ∫Ω f (⃗ r − x⃗, h)dV.
x)W (⃗
13
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
r − x⃗, h)d⃗
∫Ω W (⃗ x=1 (I.2.1)
⃗ (⃗
< ∇f ⃗ (⃗
r) >h = ∫ ∇f r − x⃗, h)dV
x)W (⃗
Ω (I.2.2)
= ∫ f (⃗ r − x⃗, h)⃗
x)W (⃗ ndS + ∫ f (⃗ ⃗ (⃗
x)∇W r − x⃗, h)dV.
δΩ Ω
Nous voyons ainsi l’intérêt de ce schéma d’interpolation : il est possible d’évaluer le gra-
⃗ > à partir des valeurs connues (du champ f , du noyau W et de son gradient)
dient < ∇f
sans avoir à connaı̂tre la valeur du champ ∇f .
De la même manière que précédemment, en imposant la condition suivante :
⃗ (⃗
∫Ω ∇W r − x⃗, h)d⃗
x=0 (I.2.3)
on peut montrer Mas-Gallic et Raviart (1987) que l’approximation du gradient est aussi
d’ordre 2 :
⃗ (⃗
< ∇f r) >h = ∇f ⃗ (⃗
r) + O(h2 ).
Nous ne discuterons pas ici d’une variation du schéma d’interpolation où la longueur
de lissage h peut varier en temps et en espace. Se reporter à Doring (2005), Oger (2006)
pour plus d’informations.
14
I.2.1. SPH CLASSIQUE
Pi
Voisinage de Pi D
Une fois introduit cet espace de discrétisation, nous pouvons définir l’approximation
de l’intégrale d’une fonction f sur le domaine Ω par la formule de quadrature suivante (de
type trapèze) :
∫Ω f (⃗
x)d⃗
x ≈ ∑ f (⃗
xi )ωi . (I.2.4)
i∈P
Cette formule de quadrature est valable pour les intégrales volumiques, c’est-à-dire que
l’on ne peut pas approximer avec cette formulation les intégrales surfaciques apparais-
sant dans la dérivation du produit de convolution (I.2.2). Grâce au caractère compact
15
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
xi ) >h
< f (⃗ = ∑ f (⃗ xi − x⃗j )ωj
xj )W (⃗
j∈D
et sa dérivée :
⃗ (⃗
< ∇f xi ) >h = ∑ f (⃗ ⃗ (⃗
xj )∇W xi − x⃗j )ωj . (I.2.5)
j∈D
∑ Wij ωj ≠ 1 et ⃗ ij ωj ≠ ⃗0.
∑ ∇W (I.2.6)
j∈D j∈D
16
I.2.1. SPH CLASSIQUE
1.03 1D 1.03 2D
<1>
1.02 1.02
< ∇x >
1.01 1.01
1 1
0.99 0.99
0.98 0.98
0.97 0.97
1 1.5 2 2.5 3 1 1.5 2 2.5 3
h/Δx
1.03 3D 600
1D
550
2D
Nombre de voisins
1.02 500 3D
450
1.01 400
350
1 300
250
0.99 200
150
0.98 100
50
0.97 0
1 1.5 2 2.5 3 1 1.5 2 2.5 3
h/Δx h/Δx
Tel qu’il a été mentionné, ces tests ont été effectués sur une distribution régulière
des.points de discrétisation. Des techniques de correction seront présentées ultérieure-
ment pour augmenter la précision lorsque cette distribution se déstructure (à cause du
mouvement lagrangien) .
17
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
Étant dans le formalisme Lagrangien (les volumes de contrôle ωi sont suivis au cours du
temps et restent constants) d’une part, et ayant conservation de la masse localement (et
donc globalement) d’autre part, on peut remarquer que les masses mi sont constantes
au cours du temps. Les points Pi ne sont plus seulement des points nécessaires à la
définition des volumes ωi pour la formule de quadrature (I.2.4) mais ce sont aussi des
points matériels qui sont identifiés comme étant des particules. D’où la qualification de
méthode particulaire pour la méthode SPH.
Dρω
= 0. (I.2.8)
Dt
Ainsi, de cette dernière équation on en déduit une équation d’évolution des volumes au
cours du temps :
Dω
= ωdiv(⃗ u).
Dt
Une discrétisation intuitive de l’équation de quantité de mouvement serait, avec ce
qui a été mentionné précédemment sur la dérivée de l’opérateur d’interpolation discret :
Dρ⃗
u 1
= ⃗ >i +⃗
− < ∇P g = ⃗ ij ωj + g⃗.
− ∑ Pj ∇W (I.2.9)
Dt i ρ j
En pratique, ces discrétisations sont peu précises et mènent à des instabilités numériques.
La principale raison en est la violation du principe d’action-réaction entre deux particules
(où normalement l’effort F⃗j→i = −mi ρ1i Pj ∇W
⃗ ij ωj de la particule j sur i est opposé à l’effort
réciproque F⃗i→j ). Avec ce schéma nous avons :
1 1
− mi ⃗ ij ωj
Pj ∇W ≠ −mj ⃗ ji ωi
Pi ∇W
ρi ρj
ce qui ne garantit pas la conservation locale de quantité de mouvement et de moment
angulaire.
Pour pallier cela on peut réarranger les équations de manière à les symétriser et
retrouver la réciprocité des interactions. Nous ne détaillerons pas ici tous les travaux
concernant ces aspects (e.g. Bonet et Lok (1999), Vila (2000), Monaghan (1992)) car tous
ne respectent pas la consistance du schéma complet (i.e. que le système discrétisé converge
vers le système continu ; un moyen de la démontrer est de savoir si la discrétisation choisie
pour l’équation de continuité génère de manière formelle une discrétisation donnée pour
l’équation de quantité de mouvement). Nous allons seulement développer la méthode
variationnelle présentée par Bonet et Lok (1999) et reprise dans Colagrossi (2005), qui
garantit la consistance du schéma. Cette approche est basée sur une étude variationnelle
et une approche thermodynamique.
La conservation du moment pour un système de particules non soumis à des forces
extérieures peut s’écrire de la façon suivante :
D⃗
ui
∑ mi = ∑ T⃗i (I.2.10)
i Dt i
18
I.2.1. SPH CLASSIQUE
E = ∑ mi ei . (I.2.12)
i
Le travail des forces internes T⃗i dû à un déplacement arbitraire δ w⃗i est directement relié
à la variation incrémentale DE de l’énergie potentielle élastique totale E :
La variation élémentaire de densité Dρi est discrétisée selon la manière suivante (ce choix
met en évidence la symétrie du terme de la divergence de v⃗) :
⃗ ij ωj .
Dρi = −ρi ∑(δ w⃗j − δ w⃗i ) ⋅ ∇W
j
Dρi = ⃗ ij ωj − ρi ∑ δ w⃗j ⋅ ∇W
ρi ∑ δ w⃗i ⋅ ∇W ⃗ ij ωj . (I.2.15)
j j
⃗ ij = −∇W
Grâce à la propriété d’antisymétrie du gradient du noyau (∇W ⃗ ji ), la variation
incrémentale de l’énergie DE s’écrit finalement :
mi mj
⃗i ]
DE[δ w = ∑∑ ⃗ ij ⋅ δ w⃗i .
(pi + pj )∇W
i j ρi ρj
Ainsi grâce au principe des travaux virtuels (I.2.13), on en déduit la force interne T⃗i
s’exerçant sur chaque particule i :
mi mj
T⃗i = −∑ ⃗ ij .
(pi + pj )∇W
j ρi ρj
19
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
In fine, grâce à (I.2.10), nous pouvons en déduire une discrétisation consistante de l’équa-
tion de quantité de mouvement (sans terme source). Le système complet discrétisé s’écrit
donc :
⎧
⎪
⎪
⎪
⎪
D⃗ u
⃗ ij ωj + ρi g⃗
⎪
⎪ ρi = − ∑(pi + pj )∇W
⎪
⎪
⎪
Dt i j
⎪
⎪
⎪
⎨ Dρ
= −ρi ∑(⃗ ⃗ ij ωj
uj − u⃗i ) ⋅ ∇W (I.2.16)
⎪
⎪
⎪ Dt i
⎪
⎪
⎪
j
⎪
⎪
⎪ D⃗
x
⎪
⎪
⎪ = u⃗i
⎩ Dt i
Le traitement des conditions aux limites est un point délicat à effectuer pour les
méthodes particulaires en général et la méthode SPH ne déroge pas à la règle. Les dif-
ficultés proviennent de deux points distincts : l’un concerne l’application à proprement
parler des conditions aux limites sur les points mobiles que sont les particules. L’autre
difficulté vient de l’incomplétude du support compact en limite de domaine. En effet, en
s’approchant d’une des frontières du domaine (où s’applique une condition aux limites), le
point de discrétisation voit son support d’interpolation associé intersecter cette frontière.
L’interpolation des gradients, qui est basée sur une quadrature volumique et surfacique
(§I.2.1.1.2), ne peut alors plus s’appliquer car on n’approxime pas cette dernière com-
posante. Pour contourner cela on considère les cas où les conditions aux limites doivent
s’appliquer :
– s’il s’agit d’une paroi, les conditions à imposer sont celles de non-pénétration et de
glissement ;
– s’il s’agit de la surface libre, les conditions dynamiques et cinématiques sont à
imposer.
20
I.2.1. SPH CLASSIQUE
⃗
n
Fi i Fi i
u⃗Fi u⃗i Fi Fi
Fig. I.2.3 – Particules (●) et leurs particules (●) fantômes. Paroi plane et coin concave.
cule génératrice i (⃗
n est la normale à la paroi P ) :
⎧
⎪ x⃗Fi = x⃗i + 2(⃗
xP ⋅ n⃗ − x⃗i ⋅ n
⃗ )⃗
n
⎪
⎪
⎪
⎪ u⃗Fi = u⃗i − 2(⃗ ⃗ )⃗
ui ⋅ n n
⎨
⎪
⎪ ρFi = ρi
⎪
⎪
⎪
⎩ ωFi = ωi
Cette technique permet de vérifier la condition de glissement pourvu que les parti-
cules fantômes soient répliquées à partir des particules i situées à une distance inférieure
au rayon du support compact (2h ou 3h) afin de compléter entièrement le support d’in-
terpolation des particules dans le domaine fluide. Pour les points P de la paroi, on a :
⃗> =
< u⃗P ⋅ n ∑ (⃗ ⃗ )WiP ωi +
ui ⋅ n ∑ (⃗ ⃗ )WjP ωj
uj ⋅ n
i∈particules j∈f antomes
= ∑ (⃗ ⃗ )WiP ωi +
ui ⋅ n ∑ (−⃗ ⃗ )WjP ωj
uj ⋅ n
i∈particules j∈particules
= 0.
⃗ = 0)pour les équations
u⋅n
Cependant, en dérivant la condition cinématique à la paroi (⃗
d’Euler, on obtient une accélération non nulle pour les points P sur la paroi :
D⃗
uP 1
< ⃗> =
⋅n ∑ (pi + pP )(∇W⃗ iP ⋅ n
⃗ )ωi
Dt ρP i∈particules
1
+ ⃗ jP ⋅ n
∑ (pj + pP )(∇W ⃗ )ωj + g⃗ ⋅ n
⃗
ρP j∈f antomes
1
= ∑ (pi + pP )(∇W⃗ iP ⋅ n
⃗ )ωi
ρP i∈particules
1
+ ⃗ jP ⋅ n
∑ (pj + pP )(−∇W ⃗ )ωj + g⃗ ⋅ n
⃗
ρP j∈particules
= g⃗ ⋅ n
⃗.
Afin de corriger ce défaut (qui se manifeste dans les simulations par des particules qui
ont tendance à traverser les parois perpendiculaires à la gravité) , il est donc nécessaire de
prendre en compte le terme source g⃗ dans l’extrapolation de la densité pour les particules
fantômes :
D⃗
uP ∂p ρ0 c20 ρ γ
< ⃗ >= 0
⋅n ⇐⇒ = ρ⃗ ⃗ avec p(ρ) =
g⋅n [( ) − 1]
Dt ∂n γ ρ0
21
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
Cette dernière équation se résout à une constante près, déterminée par la valeur de la
densité de la particule (cf. paragraphe sur les conditions initiales pour la méthode de
résolution) :
1
ρ0γ−1 (γ − 1) γ−1
La condition dynamique de surface libre (une pression patm est imposée sur la surface
libre ΩSL ) est plus difficile à satisfaire étant donné qu’on ne dispose par de points de
discrétisation sur cette surface, ni même de localisation explicite des particules dans son
voisinage sauf à utiliser des algorithmes spécifiques élaborés. Une manière de procéder
serait alors, de la même façon que pour les parois, de compléter le support avec des
particules fantômes sur lesquelles on appliquerait la condition pFi = patm .
Cette technique serait extrêmement complexe à mettre en œuvre dans la mesure où il
faudrait être capable de localiser la surface libre, d’en déterminer une normale locale, de
générer des fantômes de manière assez précise et consistante. Sauf pour le cas particulier où
la pression atmosphérique est nulle (à une constante près) : chaque particule fantôme ayant
une pression nulle il n’est plus nécessaire de les créer. Pour se retrouver donc dans cette
configuration, on opère un changement de référentiel de pression (via l’équation d’état, en
imposant ρ = ρ0 sur la surface libre). On est alors ramené à un problème d’interface un
peu spécifique où le fluide supérieur est considéré toujours à pression et vitesse nulles. On
est alors confronté aux difficultés numériques spécifiques à un problème d’interface (cf. §
22
I.2.1. SPH CLASSIQUE
suivant).
Il est à noter que cette manière d’envisager les choses n’est valable que si l’on ne
vient pas corriger le noyau d’interpolation pour interpoler dans le fluide, en le considérant
incomplet. Dans le raisonnement précédent, il est complet avec une partie de son support
pour laquelle les champs interpolés sont nuls.
Dans ce cas particulier d’une interface où le fluide supérieur est à pression et vitesse
nulles, on peut s’intéresser plus finement à la prise en compte de la présence de la surface
libre à travers la consistance et la précision des schémas de discrétisation des opérateurs
différentiels (la divergence de la vitesse dans l’équation de continuité et le gradient de
pression dans l’équation de quantité de mouvement) et du schéma discret pour les points
de discrétisation proche de la surface libre. Cette démarche, où l’on étudie notamment
l’influence de l’approximation de l’intégrale surfacique au niveau continu puis discret, est
relativement complexe et n’a pas été abordé dans cette étude (cf. travaux de Colagrossi
et alii (2008)).
D r⃗
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
x
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxx
r⃗
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
+ xx r⃗
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Y
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxxxxxxx
23
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
du champ de densité n’est lui, pas désirable et on verra au §I.2.3 qu’on peut modifier le
schéma pour éviter son lissage.
h
Δtac ≤ CCF Lac
c
où CCF Lac est un nombre adimensionnel permettant de pondérer cette restriction suivant
l’ordre du schéma d’avance en temps.
Une autre condition à vérifier est basée sur les termes sources Monaghan (1992) :
√
h
Δt ≤ 0.4
s
f
où f est le module de la force s’exerçant sur une particule (par exemple la gravité). Dans
la pratique, ce pas de temps Δts est largement supérieur (de plusieurs ordres) au pas de
temps acoustique Δtac et n’est jamais calculé.
Ces conditions sont à vérifier sur chacune des particules à chaque pas de temps. Le
calcul global est donc guidé par le pas de temps minimal sur l’ensemble des particules du
domaine de discrétisation :
Δt = min(Δtac
i , Δti ).
s
i
24
I.2.1. SPH CLASSIQUE
D⃗
u 1
= − ⃗ ij ωj + g⃗
∑(pi + pj + Πij )∇W
Dt i ρi j
où les quantités f¯ij et fij sont définies par (notations couramment utilisées par la suite) :
fi + fj
f¯ij = et fij = fi − fj .
2
Le coefficient α est déterminé en fonction du problème à résoudre : sa valeur est classi-
quement comprise entre 0,01 et 1 (si l’on est précautionneux, on peut utiliser des valeurs
comprises entre 0,01 et 0,1).
Cette approche est relativement simple et peu coûteuse pour améliorer la stabilité des
simulations. Cependant, elle atteint rapidement des limites pour plusieurs raisons. Tout
d’abord, le choix du coefficient α se fait dans l’optique d’éviter la formation d’instabilités
tout en minimisant la perte d’énergie (le système doit rester le plus conservatif possible).
Ensuite, la détermination du coefficient dépend tout aussi bien du problème à traiter
(rupture de barrage, ballottement, instabilités, . . .), que de la modélisation du fluide (sa
vitesse du son). Enfin, pour une configuration donnée, la dissipation à apporter au schéma
ne sera pas nécessairement la même dans toutes les zones (un zone d’impact nécessite un
coefficient plus grand qu’une zone de propagation de houle par exemple) et au cours du
temps : bien que basé sur la vitesse de l’écoulement, ce terme de viscosité artificielle ne
prend en compte que partiellement la partie acoustique de l’écoulement (via la vitesse du
son). Ainsi les discontinuités physiques (au sens hyperbolique) ne sont pas correctement
modélisées.
Pour rendre la tâche encore un peu plus complexe, la simulation d’écoulements bi-
fluides nécessite la détermination d’un coefficient de viscosité artificielle pour chaque fluide
(puisque chacun possède ses propriétés {densité et vitesse du son}).
25
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
Bien qu’il existe d’autres formes du terme de viscosité artificielle Balsara (1995), un des
moyens de décentrer le schéma de manière satisfaisante (du point de vue des instabilités et
de manière cohérente par rapport à l’hyperbolicité du système) tout en restant conservatif
consiste à introduire des solveurs 1 de Riemann. Ceci sera abordé au chapitre suivant.
⎪
⎪ ρI (γY −1) (γY − 1)g(HI − y) Y γ −1
⎪
⎪
⎪ ρ(y) = [( ) + ] pour 0 < y < HI
⎩
ρ 0Y c20Y
A titre d’exemple, sur la figure (I.2.5) sont représentés les profils de pression hydrosta-
tique (air et eau) incompressible et compressible dans une cuve de 2 × 1m. Afin d’exagérer
les différences, les vitesses du son sont prises égales à 1m.s−1 pour les deux fluides (ce
qui est loin de √l’hypothèse de pseudo-compressibilité qui serait basée ici sur une vitesse
caractéristique gHcuve ).
1
Barbarisme provenant de l’anglais solver mais largement adopté par les praticiens français à la place
de résolveur.
26
I.2.2. LES INSTABILITÉS À L’INTERFACE
0 5
pX 10 15
2
1.5
y
0.5
0
15 5000 9985
pY
Fig. I.2.5 – Profils de pression incompressible (−−) et compressible (2) dans la cuve dans
chaque fluide.
27
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
0 -0.6 0
(Pi+Pj) - Interface
Pj - Interface
-50
-2000
-0.8 -100
-4000
ρg ∇y P
ρg ∇y P
∇y P
(Pi+Pj) - SL -150
Pj - SL
-6000 (Pi+Pj) - Interface
1
Pj - Interface
-1 -200
-8000
-250
-10000
-1.2 -300
-4 -2 0 2 4 -4 -2 0 2 4
y/h y/h
(a) Composante verticale du gradient de pression (b) Comparaison entre le terme moteur de l’écou-
lement et l’accélération de pesanteur
par la particule la plus proche de l’interface côté eau est égale à 25% de la gravité contre
270 fois l’accélération de la pesanteur pour celle se trouvant du côté de l’air (voir figure
I.2.6(b))Oger (2006).
L’approximation des conditions dynamiques de surface libre se traduit donc concrè-
tement par de fortes accélérations des particules « légères » situées à proximité (dans un
rayon d’interaction 2h) de l’interface.
Malgré le fait que les champs de gradient de la pression (et de la vitesse) soient dis-
continus à l’interface, le traitement par l’opérateur d’interpolation des gradients spatiaux
des variables continues (pression et vitesse normale) permet d’effectuer un lissage de ces
champs gradients discontinus. Cette perte de précision est alors compensée par une sta-
bilité du schéma global (les mécanismes exacts de ce transfert d’erreur n’ont pas été mis
en évidence à ce jour).
28
I.2.3. FORMULATION SPH-MULTIFLUIDE
conduit ainsi à utiliser une renormalisation du gradient de ce noyau qui diffère de ce qui
est présenté classiquement avec ce type de noyau.
Enfin, se basant sur la méthode SPH dite classique, on n’évoquerons dans cette section
que les modifications apportées. Les points n’ayant pas été abordés sont considérés comme
traités de la même manière que dans la méthode classique.
∂ x⃗
J = det(F) avec F =
∂X ⃗
⃗
où x⃗ est la position spatiale d’un point matériel quelconque X.
L’équation de continuité s’écrit donc finalement :
D log J ⃗
ρ0 (X) v(X)⃗
= u) avec J =
div(⃗ = . (I.2.18)
Dt ρ(X)⃗ ⃗
v0 (X)
Le Jacobien peut donc être vu comme le rapport entre la masse volumique initiale ρ0 (X) ⃗
d’un point X ⃗ et sa masse volumique ρ(X) ⃗ à l’instant t, ou comme le rapport entre son
volume spécifique v(X) ⃗ et son volume initial v0 (X).⃗
Cette écriture de l’équation de continuité en terme du Jacobien J au lieu de la masse
volumique ρ est réalisée pour des raisons numériques qui seront expliquées ultérieurement.
⃗ ∈ K) :
avec la loi d’état fermant le système (pour X
29
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
où le nouveau noyau W S est connu sous le nom de noyau Shepard dans la littérature
Belytschko et alii (1998).
À partir du schéma d’interpolation discret classique (I.2.5) et en effectuant une ma-
nipulation sur le gradient de la fonction à évaluer (on utilise la propriété < ∇(f.1) >= f <
∇1 > +1 < ∇f >), on peut obtenir un nouveau schéma d’interpolation du gradient d’une
fonction en introduisant le noyau Shepard :
xi ) ≻ = ∑ fj ∇W
≺ ∇f (⃗ ⃗ ijS ωj − ⃗ ijS ωj
fi ∑ ∇W
j j
(I.2.21)
1 fi
= ⃗ ij ωj
∑ fj ∇W − ⃗ ij ωj
∑ ∇W
ci j ci j
où fi = f (⃗
xi ), ci = c(⃗ ⃗ ij = ∇W
xi ) et ∇W ⃗ (⃗
xi − x⃗j , h). Il est à noter que cette formulation est
légèrement différente de celle connue comme étant le gradient du noyau Shepard ∇W ⃗ ijS
Belytschko et alii (1998) :
1 ≺ fi ≻
⃗ ijS ωj
∑ fj ∇W = ⃗ ij ωj
∑ fj ∇W − ⃗ ij ωj
∑ ∇W
j ci j ci j
où l’on remarque que le second terme fait intervenir l’interpolation du champ f au point x⃗
avec le noyau Shepard. Cette dernière est obtenue par différenciation de (I.2.20) alors que
(I.2.21) est une correction de (I.2.5) pour en améliorer la précision, en particulier proche
de la surface libre, ce qui est notre objectif.
Ainsi le schéma d’interpolation utilisé est :
⎧
⎪ Wij
⎪
⎪
⎪ ≺ f (⃗
x) ≻ = ∑ fj ωj ;
⎪
⎪ j ci
⎨ ⃗ ij
∇W
⎪
⎪
⎪ ⃗ (⃗
⎪
⎪ ≺ ∇f x ) ≻ = ∑ (fj − fi ) ωj .
⎪
⎩ j ci
30
I.2.3. FORMULATION SPH-MULTIFLUIDE
si l’on définit la masse mj des particules comme étant le produit de leur masse volumique
ρj par leur volume élémentaire ωj . Cette approche, utilisée dans la formulation originale
de la SPH, ne permet pas de prendre en compte de manière précise les discontinuités de
densité entre deux fluides. Pour pallier cela, Hu & Adams se sont appuyés sur une autre
évaluation de la masse volumique :
< ρ(⃗
x) > = mi ∑ Wij ou < ρ(⃗
x ) > = m i ni avec ni = ∑ Wij . (I.2.23)
j j
Cette formulation est bien capable de prendre en compte une discontinuité de densité
car l’évaluation de la densité au point x⃗ ne dépend que de sa propre masse mi et d’une
contribution géométrique ni de la part des particules j voisines.
Avec cette formulation, une variation incrémentale de densité Dρi peut se mettre sous
la forme :
Dρi = mi ∑ ∇W ⃗ ij ⋅ (δ w⃗i − δ w⃗j ). (I.2.24)
j
pi
⃗ i ] = ∑ mi
DE[δ w Dρi
i ρ2i
pi
= ∑ ⃗ ij ⋅ (δ w⃗i − δ w⃗j )
∇W
2 ∑
i ni j
pi pj D⃗
ui
= ∑ ∑ ( 2 + 2 ) ∇W ⃗ ij ⋅ δ w⃗i = − ∑ T⃗i ⋅ δ w⃗i = − ∑ mi ⋅ δ w⃗i
i j ni nj i i Dt
Wij
≺ ρ(⃗
x) ≻ = ∑ mj WijS avec WijS = et ci = ∑ Wil ωl pour tout x⃗ ∈ K . (I.2.25)
j∈K ci l∈K
31
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
Le calcul du terme ci = c(⃗x) est basé sur la sommation de quantités pour des particules
appartenant au même fluide K que le point d’intérêt x⃗ (d’où la restriction ∀⃗ x ∈ K). De
cette manière les discontinuités de densité sont traitées explicitement.
Il est intéressant de remarquer que cette évaluation de la densité (I.2.25) est de concep-
tion différente de la formule (I.2.23) car, étant basée sur le noyau Shepard, elle est fonction
des volumes ω qui doivent être calculés via l’équation de continuité (I.2.18). Ainsi, l’utili-
sation de l’équation (I.2.25) induit une relaxation dans le lien entre la masse mi , le volume
ωi et la masse volumique ρi d’un point i quelconque. Il faudra donc, par la suite, manipuler
la relation mi = ρi ωi avec précaution. En particulier, le champ de densité calculé direc-
tement à partir de la masse et du volume (i.e. mi /ωi ) ne coı̈ncide plus avec l’évaluation
donnée par (I.2.22) ; l’évaluation (I.2.25) donne seulement une représentation lissée ≺ ρ ≻
du champ de densité (i.e. ≺ ρ ≻ < mi /ωi ).
Ainsi cette nouvelle formulation nécessite l’intégration temporelle de l’équation de
continuité (I.2.18) pour l’évaluation des volumes tandis que l’évaluation de la densité est
faite avec le noyau Shepard (I.2.25). La divergence de la vitesse est évaluée par :
1
ui ) =
div(⃗ ∑(⃗ ⃗ ij ωj
uj − u⃗i ) ⋅ ∇W avec di = ∑ Wil ωl (I.2.26)
di j l
où cette fois-ci les sommations sont faites sur tout le voisinage, sans tenir compte du type
de fluide. De cette manière, les discontinuités de vitesse à l’interface sont régularisées par
(I.2.26), tandis que les discontinuités de densité sont explicitement prises en compte par
(I.2.25).
Avec cette nouvelle formulation, une variation incrémentale de densité Dρi peut se
mettre sous la forme :
ρi
Dρi = −ρi div(δ w ⃗i ) = − ∑(δ w⃗j − δ w⃗i ) ⋅ ∇W ⃗ ij ωj (I.2.27)
di j
pi
⃗ i ] = ∑ mi
DE[δ w Dρi
i ρ2i
pi
= − ∑ mi ⃗ ij ⋅ (δ w⃗j − δ w⃗i ) ωj
∑ ∇W
i ρ i di j
mi mj pi pj D⃗
ui
= ∑∑ ( + ) ∇W ⃗ ij ⋅ δ w⃗i = − ∑ T⃗i ⋅ δ w⃗i = − ∑ mi ⋅ δ w⃗i
i j ρi ρj di dj i i Dt
Termes géométriques Ces termes n(⃗ x) (cf. équation I.2.23) peuvent être introduits
de la même manière que par Hu & Adams Hu et Adams (2006). La principale différence
provient de ce qu’on raisonne sur une approximation de la distribution des volumes au
lieu des masses des particules. Ainsi, au début d’une simulation on choisit d’aligner carté-
siennement les particules, leurs volumes sont donc identiques. D’autre part, travaillant sur
32
I.2.3. FORMULATION SPH-MULTIFLUIDE
des fluides faiblement compressibles, nous pouvons toujours supposer que les volumes ωj
des particules proches de la particule d’intérêt i sont presque identiques à ωi (i.e. que le
gradient spatial de la distribution des volumes est négligeable devant l’échelle de longueur
assimilée à la taille du support compact). Avec cette hypothèse, il est possible d’écrire :
Wj (⃗
xi ) ωj Wj (⃗
xi ) ωj Wj (⃗
xi ) Wj (⃗
xi )
= ≃ = .
d(⃗xi ) ∑k Wk (⃗
xi ) ωk ∑k Wk (⃗
xi ) n(⃗xi )
⃗ ij
∇W pi pj
uj − u⃗i ) ⋅
ui ) = ∑(⃗
div(⃗ ⃗ i = ∑(
et ∇p ⃗ ij .
+ ) ∇W
j ni j ni nj
Cette écriture met en évidence ce qui a déjà été mentionné plus haut, à savoir que l’éva-
luation de la masse volumique ρ(⃗ xi ) et le volume ωi sont reliés grâce à une distribution
lissée Mi des masses (qui ne coı̈ncide pas avec la masse ponctuelle mi ).
Ainsi, les termes géométriques l(⃗x) permettent d’évaluer la masse lissée M tandis que
les volumes ω sont évalués à partir de l’intégration temporelle de l’équation de continuité
(I.2.18). Puis, grâce à ces deux champs M et ω, on peut évaluer le champ de masse
volumique ρ avec l’équation (I.2.29), et ensuite le champ de pression (grâce à l’équation
d’état).
Remarque :
Cette approche avec termes géométriques a été choisie dans ce travail pour rester proche
des développements effectués par Hu & Adams et ainsi mener des comparaisons dans un
cadre restreint. Toutefois, à titre complémentaire après tous les cas de validations étudiés,
quelques calculs ont été menés sans l’approximation faite sur la faible variation spatiale
des volumes, et en comparaison avec ces mêmes calculs opérés avec cette approximation,
les différences n’ont pas été décelables qualitativement. Pour des travaux ultérieurs, il
serait donc envisageable de ne pas utiliser cette approximation.
33
CHAPITRE I.2. À BASE DE SPH CLASSIQUE
Puis une fois la densité et la pression connues, on peut évaluer les dérivées temporelles :
⎧
⎪ D⃗
x
⎪
⎪
⎪ = u⃗i ; ni = ∑ Wik
⎪
⎪
⎪ Dt i
⎪
⎪
⎪
k
⎪
⎪
⎪
⎪
⎪
⎪ D log J ⃗ ij
∇W
⎨ uj − u⃗i ) ⋅
= ∑(⃗ ;
⎪
⎪
⎪
Dt i j ni
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪ D⃗
u pi pj
⎪
⎪ ρ = − ∑ ( + ) ∇W ⃗ ij + S⃗ .
⎪
⎪ Dt i n n
⎩ j i j
Ces dérivées permettent alors de mettre à jour les volumes, vitesses et positions des par-
ticules.
On peut remarquer que si l’on ne considère qu’un seul fluide, le modèle reste inchangé
et l’on a alors li = ni .
34
Chapitre I.3
La SPH hybride
Dans ce chapitre nous allons reconsidérer la méthode SPH par une approche plus
mathématique, menée principalement dans les travaux de Vila (2000), Lanson et Vila
(2001) (repris par Deuff (2007) et Guilcher (2008)). Le système d’équations présenté plus
haut (les équations d’Euler pour un fluide compressible) présente un caractère hyper-
bolique, ce qui change radicalement l’approche que l’on peut avoir pour l’étude d’un
écoulement incompressible (que l’on cherche toujours à étudier dans notre cadre de tra-
vail, la pseudo-compressibilité n’étant qu’un moyen d’approcher la solution asymptotique
incompressible). En effet dans un système d’équations pour un fluide incompressible, le
caractère elliptique des équations induit que l’information se propage instantanément dans
tout le domaine. Ce qui n’est plus cas avec un système hyperbolique, où cette information
va avoir une vitesse et une direction de déplacement propres dans l’espace. De plus, le
caractère compressible du fluide permet aussi l’existence de discontinuités dans les champs
même si l’écoulement étudié n’en présente pas initialement. Ceci remet en cause le cadre
d’une approche différentielle du système d’équation, telle que nous l’avons précédemment
présentée.
Ainsi ces problématiques mathématiques et numériques naissant de l’étude et de la
modélisation d’un fluide compressible vont être abordées dans ce chapitre. Des solutions
vont être proposées et étudiées, provenant en partie du domaine des Volumes Finis. L’on
parlera ainsi de SPH hybride, méthode à mi-chemin entre la conception classique de la
SPH (à partir d’une réflexion sur les méthodes non-maillées) et la méthode des Volumes
Finis.
35
CHAPITRE I.3. LA SPH HYBRIDE
S = (0, ρ⃗
g ) le terme source. L’opérateur de transport Lu est définit par :
∂φ ∂
Lu (φ) = + ∑ (ui φ).
∂t i=1,d ∂xi
Ce système (I.3.1) voit donc sa solution exacte satisfaire la propriété de conservativité
(i.e. conservativité de la masse et de la quantité de mouvement) :
D
( φ dx) = ∫ S dx.
Dt ∫Rd
(I.3.2)
Rd
De plus, contrairement au schéma SPH classique où l’on raisonne en terme de masses
ponctuelles (cf. §I.2.8), on observe ici l’évolution temporelle des volumes ω transportés
par le champ u⃗ :
Dω
= ω div(⃗u).
Dt
Enfin, on adjoint l’équation de transport des volumes matériels, qui traduit le caractère
lagrangien du schéma :
D⃗x
= u⃗.
Dt
I.3.1.1.1 Formulation ALE
Le formalisme ALE (pour Arbitrary Lagrangian-Eulerian) consiste à modifier l’opé-
rateur de transport Lu afin que celui-ci soit lié à un champ de vitesse arbitraire u⃗0 :
∂φ ∂
Lu0 (φ) = + ∑ i
(u0 i φ).
∂t i=1,d ∂x
Par conséquent le flux intervenant dans l’équation (I.3.1) est modifié de la sorte :
F = F E − u0 φ
où F E est le flux eulérien. Les équations de transport des points matériels et d’évolution
des volumes sont aussi modifiées :
D⃗x Dω
= u⃗0 et = ω div(⃗
u0 ).
Dt Dt
Ce choix de formulation permet donc de « dégénérer » le système d’équations en une
approche purement eulérienne (⃗ u0 = ⃗0), ou en une approche strictement lagrangienne
u0 = u⃗), ou encore en une approche mixte (où la vitesse u⃗0 peut être une vitesse lissée
(⃗
localement).
36
I.3.1. LA FORMULATION MONOFLUIDE
Dh f (⃗
x) ⃗ (⃗
= < ∇f x) >h − f (⃗ ⃗ 1 (⃗
x) < ∇ x) >h
Dh f (⃗
xi ) = ∑ (f (⃗
xj ) − f (⃗ ⃗ ij ωj
xi ))∇W (I.3.3)
j∈D
x) = ⃗0.
et qui garantit finalement que Dh 1(⃗
Dωi φi Dωi φi
+ ωi Dh Fi = ωi Si soit ⃗ ij ωj = ωi Si
+ ωi ∑(Fj − Fi )∇W
Dt Dt j
Cette discrétisation a été étudiée par Lanson et Vila (2001) et Oger (2006) et pré-
sente une bonne précision au sein des fluides. Cependant elle présente deux inconvénients
majeurs. Le premier est que ce schéma n’est pas conservatif (cf. éq. I.3.2) dans le sens où :
D
(∑ ωi φi ) ≠ ∑ ωi Si . (I.3.4)
Dt j j
Le deuxième inconvénient concerne les opérateurs Lu et divergence qui ne sont pas dé-
finis au niveau des discontinuités (qui existent au sein du fluide compressible ; cf. §I.3.1.1 ).
L’idée est donc d’utiliser une formulation variationnelle et de reporter ces opérateurs
de dérivation sur les fonctions-tests associées pour obtenir une solution faible.
37
CHAPITRE I.3. LA SPH HYBRIDE
D∗h f (⃗
xi ) = ∑ (f (⃗ ⃗ ij − f (⃗
xi )∇W ⃗ ji ) ωj .
xj )∇W
j∈D
On retrouve ainsi la discrétisation consistante du schéma SPH classique sous une forme
légèrement différente. Ce schéma vérifie donc la conservativité discrète (i.e. le contraire
de l’équation (I.3.4)).
Cependant, la principale différence avec la même démonstration faite plus tôt (outre
le formalisme mathématique qui permet par la suite d’étudier certaines propriétés de la
méthode, entre autres démontrer la convergence sous certaines hypothèses Vila (2000))
va concerner l’étude de la stabilité de ce schéma.
38
I.3.1. LA FORMULATION MONOFLUIDE
Nous avons en effet déjà évoqué ce problème de stabilité pour le schéma SPH classique,
dû à l’utilisation d’une avance en temps explicite avec la discrétisation spatiale centrée.
La solution retenue pour le schéma classique avait été de décentrer le schéma en espace
au moyen d’une pseudo-viscosité rajoutée dans l’équation de quantité de mouvement.
Mais avec le cheminement effectué avec cette formulation variationnelle, nous remar-
quons que les équations sont fortement liées entre elles par l’évaluation de la divergence
des flux. La dissipation à appliquer pour décentrer le schéma devra donc s’appliquer, a
priori, sur toutes les composantes du flux.
Dh x⃗(⃗
xi ) = ∑(⃗ ⃗ ij ωj .
xj − x⃗i ) ⊗ ∇W
j
Cette estimation doit, par ailleurs, être égale à la matrice identité I. On peut donc intro-
h
duire un nouvel opérateur d’interpolation renormalisé DR définit par :
h
DR f (⃗
xi ) = xi ) ⋅ ∑(f (⃗
B(⃗ xj ) − f (⃗ ⃗ ij ωj
xi ))∇W
j
39
CHAPITRE I.3. LA SPH HYBRIDE
convergence est assurée sous la seule condition de raffinement du pas spatial (Δx → 0).
Remarques :
∂φv ∂φv ⎛ u⃗ ⋅ n
⃗ ρnx ρny ⎞
+ An = 0 avec An = ⎜c2 nx /ρ u⃗ ⋅ n
⃗ 0 ⎟
∂t ∂n ⎝c2 ny /ρ 0 u⃗ ⋅ n⃗⎠
40
I.3.1. LA FORMULATION MONOFLUIDE
des trajectoires particulières qui sont les courbes caractéristiques. On peut les déterminer
par l’étude du système diagonalisé :
où ∂W est le vecteur accroissement des variables caractéristiques (obtenu par multiplica-
tion du vecteur des variables convectives avec la matrice de passage du système diagonal)
et Λ la matrice diagonalisée de An . Ce système, correspondant à la diagonalisation des
équations d’Euler, peut se réécrire sous forme scalaire :
∂wi ∂wi
+ λi = 0. (I.3.7)
∂t ∂n
Et on en déduit, par intégration, pour chaque composante wi de W , l’invariant de Riemann
dn
(qui reste constant le long de la courbe caractéristique = λi ) :
dt
⎛ u⃗ ⋅ t⃗ ⎞
⎜⃗ ⃗ 2c ⎟
⎜u ⋅ n − ⎟
W =⎜ γ − 1⎟ .
⎜ 2c ⎟
⎜ ⎟
⎝u⃗ ⋅ n
⃗+
γ − 1⎠
t
dt = u⃗ ⋅ n
⃗
dn
dt =
u⃗
dn
⋅ n⃗
−
c
w1 constant
⃗+
dn = u⃗ ⋅ n
c
w2 constant w3 constant dt
n
0
Fig. I.3.1 – Représentation des ondes caractéristiques dans le plan (⃗
n, t).
Sur ce diagramme, on peut avoir une représentation des régions d’influence des ondes
dans le plan espace-temps. Chaque droite correspond à une onde simple (liée à un invariant
de Riemann), telle que la quantité ∂wi reste constante au travers de cette droite. Et chaque
droite sépare deux régions où la solution est uniforme.
41
CHAPITRE I.3. LA SPH HYBRIDE
Ainsi on est capable de quantifier, en tout point de l’espace et à chaque instant, l’état
du fluide affecté par la perturbation causée par une discontinuité initiale.
D D.C. D.C.
C
D
C
Fig. I.3.2 – Problème de Riemann initial et certaines configurations parmi toutes celles
possibles dans le plan (x, t) avec D signifiant détente, C choc et D.C. discontinuité de
contact.
42
I.3.1. LA FORMULATION MONOFLUIDE
u
u
−
c
φ1/3 φ2/3
φ1/6 R.H.
Fig. I.3.3 – Cas particulier d’un problème de Riemann avec les différents états φ et
les relations de passage (« R.H. » pour Rankine-Hugoniot et « inv. Riemann » pour les
invariants).
Nous allons donner le principe de résolution du problème de Riemann pour les équa-
tions d’Euler isentropiques monodimensionnelles sous forme conservative, sans rentrer
dans les détails (se reporter à Ivings et alii (1998)). La fermeture du système se fait grâce
à l’équation de Tait, avec une évaluation de la pression absolue grâce à la masse volumique
(i.e. p = Bργ ) :
⎧
⎪ ∂ρ ∂ρu
⎪
⎪
⎪ ∂t + ∂x
⎪ = 0
⎨
⎪
⎪
⎪ ∂ρu ∂(ρu2 + Bργ )
⎪
⎪ + = 0
⎩ ∂t ∂x
2c0 2c1/3
u0 + = u1/3 +
γ−1 γ −1
43
CHAPITRE I.3. LA SPH HYBRIDE
Onde de valeur propre λ = u : Sur cette onde, qui est une discontinuité de contact,
les relations de Rankine-Hugoniot nous donnent :
u1/3 = u2/3
{
ρ1/3 = ρ2/3
On appellera cette zone entre les deux ondes non-linéaires région étoilée *.
Solution complète : On peut réécrire cette dernière relation sur les vitesses de la façon
suivante :
f (ρ∗ ) = (u1 + f1 (φ1 , ρ∗ )) − (u0 − f0 (φ0 , ρ∗ )) = 0.
Ainsi le problème de Riemann peut être totalement résolu en trouvant la solution ρ∗ de
cette équation linéaire. Ceci peut se faire grâce à la méthode de Newton-Raphson, qui en
quelques itérations (3 ou 4 typiquement), trouve une solution à 10−6 près. On initialise le
calcul avec la solution ρ∗ d’un solveur de Riemann approché que nous détaillerons plus
loin. De plus, au cours des itérations il faut remarquer que les solutions intermédiaires
influent sur la considération des ondes à prendre en compte (il faut actualiser l’état étoilé
à chaque itération suivant le type d’onde rencontré).
La vitesse normale u∗ est obtenue grâce à la demi-somme des relations valables pour
un choc à gauche et à droite (ou grâce à celle pour une détente à gauche et à droite) :
1 1
u∗ = (u0 + u1 ) + (f1 (φ1 , ρ∗ ) − f0 (φ0 , ρ∗ )).
2 2
La vitesse tangentielle, quant à elle, ne peut être déterminée avec l’approche choi-
sie (monodimensionnelle). On prolonge donc les vitesses tangentielles des états gauche
(φ1/3 = φ0 ) et droit (φ2/3 = φ1 ) jusqu’à la discontinuité de contact qui est une condition de
glissement (⃗ u1/3 ⋅ t⃗ ≠ u⃗2/3 ⋅ t⃗).
Pour compléter la détermination de tous les états possibles, il reste à préciser les états
à l’intérieur des détentes (φ1/6 si la détente se trouve à gauche ou φ5/6 si elle est à droite) :
⎧
⎪ 2 γ−1 x
⎪
⎪ u1/6 = (c0 + u0 + )
⎪
⎪
⎪ γ+1 2 t
⎨ 2
⎛ c1/6 ⎞
1/(γ−1)
⎪
⎪
⎪
⎪
⎪ ρ 1/6 = où c1/6 = u1/6 − x/t
⎪
⎩ ⎝ γB ⎠
⎧
⎪ 2 γ−1 x
⎪
⎪
⎪ u5/6 = (−c1 + u1 + )
⎪
⎪ γ+1 2 t
et ⎨ 2 1/(γ−1)
⎪
⎪
⎪
c
⎪
⎪ ρ5/6 = ( 1 ) où c5/6 = x/t − u5/6
⎪
⎩ γB
44
I.3.1. LA FORMULATION MONOFLUIDE
Les détentes, qui sont des ondes qui s’étalent dans le temps, sont comprises entre :
dx dx
u0 − c 0 < < u∗ − c ∗ et u∗ + c∗ < < u1 + c 1 .
dt dt
Les chocs sont exactement positionnés aux vitesses déterminées par les relations de Rankine-
Hugoniot :
ρ ∗ u∗ − ρ0 u0 ρ ∗ u∗ − ρ1 u1
s0 = et s1 =
ρ ∗ − ρ0 ρ ∗ − ρ1
On peut remarquer que le problème de Riemann est autosimilaire, c’est-à-dire que sa
solution, notée (ρe , u⃗e ), est valable quel que soit le rapport x/t.
t
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
, φj )
F (φixxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx F (φ , φ )
xxxxxxxxxxx
xxxxxxxxxxx
ixxxxxxxxxxx
xxxxxxxxxxx j
i F (φi , φj )
xxxxxxxxxxxxxxxxxx xxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx xxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
j xxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxx φ
xxxxxxxxxxxxxxxxxx
j
xxxxxxxxxxxxxxxxxx xxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
⃗ ij
n φ
xxxxxxxxxxxxxxxxx
i
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
φ xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxx
ixxxxxxxxxxxxxxxx
xxxxxxxxxxx
φ
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
j
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxx nij xxxxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxx nij
Les équations d’Euler peuvent se mettre sous la forme discrète suivante (en utilisant
le formalisme Volumes Finis) :
Dωi φi
⃗ ij lj
+ ∑ Fij ⋅ n = 0 (I.3.10)
Dt j
où Fij = F (φi , φj ) est le flux sur la face commune à i et j, lj est la longueur de cette face,
⃗ ij son vecteur normal (orienté vers l’extérieur de i) et ωi le volume de la maille i.
n
Ainsi pour calculer le flux numérique Fij on se ramène à un problème local (cf. figure
I.3.4) où la face commune devient une interface entre deux états discontinus φi et φj . Avec
la méthode de Godunov, on détermine donc ce flux en résolvant exactement le problème
de Riemann.
Remarque :
45
CHAPITRE I.3. LA SPH HYBRIDE
Solveur de Riemann à 2 détentes : Pour ce type de solveur, nous allons faire une
approximation sur le calcul de l’état étoilé. Nous allons supposer que le problème de
Riemann est posé avec 2 ondes correspondant à des détentes (appelé aussi TRRS pour
Two-Rarefaction Riemann Solver). Ainsi à partir des équations (I.3.8) et (I.3.9) issues des
relations de Rankine-Hugoniot, on obtient :
⎧
⎪
⎪
⎪
⎪
1 1
⎪
⎪ u∗ = (u 0 + u 1 ) + (c0 − c1 )
⎨ 2 γ−1
⎪
⎪
⎪ 1 γ−1
⎪
⎪
⎪ c∗ = (c0 + c1 ) + (u0 − u1 )
⎩ 2 4
L’équation d’état et la définition de la vitesse du son permettent d’avoir :
c2∗ 1/(γ−1)
ρ∗ = ( ) .
γB
La vitesse tangentielle dans la région étoilée est déterminée de la même façon que dans le
solveur exact, à savoir par prolongation des états gauche et droit.
C’est l’état étoilé déterminé grâce à cette approximation qui est utilisé comme valeur
initiale dans la méthode itérative de Newton-Raphson du solveur exact de Godunov.
46
I.3.1. LA FORMULATION MONOFLUIDE
Solveur de Riemann HLL : Ce type de solveur n’est plus basé sur le calcul approché
de l’état étoilé mais sur une estimation directe des flux. Pour cela, Harten, Lax and Van
Leer, Harten et alii (1983), ont supposé que la configuration des ondes est telle que seules
deux ondes séparent les trois régions du plan (x, t) correspondant à des états constants,
cf. figure I.3.5.
S0 t S1
B C
A D x
De plus, en supposant que les vitesses de ces ondes sont données par un autre algo-
rithme, et en appliquant la formulation intégrale des lois de conservation, on obtient une
expression approchée du flux numérique.
Le flux Fij est estimé directement par la relation suivante :
⎧
⎪ x
⎪
⎪
⎪ F0 si S0 ≧
⎪
⎪
⎪ t
⎪
⎪
⎪ S1 F0 − S0 F1 + S0 S1 (φ1 − φ0 )
F HLL = ⎨ x
⎪ si S0 < < S1 (I.3.11)
⎪
⎪ S1 − S0 t
⎪
⎪
⎪
⎪
⎪
⎪
x
⎪ F1 si S1 ≦
⎩ t
où F0 = F (φ0 ) et F1 = F (φ1 ). On peut aussi obtenir une estimation de l’état étoilé
équivalent (ce qui sera nécessaire pour estimer le flux lagrangien) :
⎧
⎪ x
⎪
⎪
⎪ φ0 si S0 ≧
⎪
⎪
⎪ t
⎪
⎪
⎪ S1 φ1 − S0 φ0 + F0 − F1
φHLL = ⎨ x
⎪ si S0 < < S1 (I.3.12)
⎪
⎪ S1 − S0 t
⎪
⎪
⎪
⎪
⎪
⎪
x
⎪ φ1 si S1 ≦
⎩ t
Une estimation simple des vitesses S0 et S1 des ondes est donnée par Davis (1988) :
S0 = min(u0 − c0 , u1 − c1 ) et S1 = min(u0 + c0 , u1 + c1 ).
Cette méthode HLL, qui est très efficace et robuste, est très populaire dans la commu-
nauté de la CFD3 des écoulements compressibles. Cependant elle ne prend pas en compte
la présence d’une discontinuité de contact, ce qui est incorrect dans le cas des équations
d’Euler.
3
Computational Fluid Dynamics
47
CHAPITRE I.3. LA SPH HYBRIDE
Solveur de Riemann HLLC : Ce solveur, Toro (1997), est basé sur le même raison-
nement que précédemment sauf qu’on suppose ici que 3 ondes sont présentes, en ajoutant
la discontinuité de contact en plus des deux ondes existantes (ce qui ajoute un C pour
Contact à l’acronyme HLL) .
Les développements donnent des résultats semblables au solveur HLL, sauf en ce qui
va concerner la composante tangentielle de la vitesse et du flux. Ainsi les deux premières
composantes du vecteur des variables conservatives de l’état étoilé équivalent φHLLC sont
donnés par la relation (I.3.12), la densité et la vitesse normale se conservant de part et
d’autre de la discontinuité de contact. Pour la vitesse tangentielle on a donc :
⎧
⎪
⎪
⎪
⎪ ρ0 uy0 si S0 ≧
x
⎪
⎪
⎪ t
⎪
⎪
⎪
⎪
⎪
⎪
x
⎪
⎪ ρ
HLL u
y0 si S0 < < S∗
(ρuy )HLLC = ⎨ t
⎪
⎪
⎪ x
⎪
⎪
⎪ ρHLL uy1 si S∗ < < S1
⎪
⎪
⎪ t
⎪
⎪
⎪ x
⎪
⎪
⎪ ρ1 uy1 si S1 ≦
⎩ t
où ρHLL est évalué grâce à la première composante de (I.3.12).
De la même manière, les deux premières composantes du vecteur flux F HLLC sont
données par la relation (I.3.11) tandis que le flux du moment tangent est donné par :
⎧
⎪
⎪
⎪
⎪ ρ0 ux0 uy0 si S0 ≧
x
⎪
⎪
⎪ t
⎪
⎪
⎪
⎪
⎪ x
⎪
⎪ (ρux )
⎪ HLL u
y0 si S0 < < S∗
(ρux uy )HLLC = ⎨ t
⎪
⎪
⎪ x
⎪
⎪
⎪ (ρux )HLL uy1 si S∗ < < S1
⎪
⎪
⎪ t
⎪
⎪
⎪ x
⎪
⎪
⎪ ρ1 ux1 uy1 si S1 ≦
⎩ t
où (ρux )HLL est évalué grâce à la deuxième composante de (I.3.12).
Les vitesses des ondes sont basées sur les estimations suivantes :
⎧ ⎧
⎪ si p∗ ≦ pK
⎪
⎪ S = u0 − q 0 c 0 ⎪
⎪
1
⎪ 0 ⎪ 1/2
⎨ S ∗ = u∗ avec qK =⎨ (ρ∗ /ρK ) − 1
γ .
⎪
⎪
⎪ ⎪
⎪
⎪ [1 + ] si p∗ > pK
⎩ S 1 = u1 + q 1 c 1 ⎪
⎩ γ(1 − ρK /ρ∗ )
Les quantités étoilées u∗ et ρ∗ (pour obtenir p∗ ) sont issues d’un solveur approximatif
simple tel qu’un TRRS ou un TSRS.
48
I.3.1. LA FORMULATION MONOFLUIDE
⃗ ij
∇W
⃗ ij =
où n ⃗ ij ∥ .
∥∇W
Tandis que le formalisme Volumes Finis discrétisé s’écrit :
Dωi φi
⃗ ij lj
+ ∑ Fij ⋅ n = 0.
Dt j
On trouve donc une correspondance entre le flux numérique Volumes Finis Fij ⋅ n ⃗ ij et
la demi-somme 12 (Fi + Fj ) n ⃗ ij du schéma SPH, ainsi qu’entre la mesure lj de la face de
la maille et la quantité ∥∇W ⃗ ij ∥ n⃗ ij ωi ωj (cette dernière correspondance est valable quelle
que soit la dimension de l’espace). Cependant il faut noter la différence notable entre les
sommations effectuées dans les deux schémas. En effet, dans un cas (Volumes Finis) la
sommation s’étend aux voisins adjacents à la maille i tandis que dans l’autre (SPH) la
sommation s’étend à la notion de voisinage du point i dans le support compact.
Une autre grande différence entre les deux approches est que la formulation SPH est
conçue pour être employée de manière lagrangienne alors que les Volumes Finis sont clas-
siquement employés sur des maillages fixes. Les idées développées par Harten et Hyman
(1983) et généralisées à la SPH par Vila (2000) vont être brièvement exposées.
On raisonne sur la formulation ALE (cf. §I.3.1.1.1) du schéma SPH qui permet de
mettre en évidence le flux eulérien :
qui est un problème de Riemann en x⃗ij dans le repère mobile lié au champ de transport
v⃗0 . Or on sait résoudre, en eulérien, le problème suivant :
⎧
⎪
⎪
⎪
⎪ ∂φ ∂
⎪
⎪
⎪ + ⃗ ij )
(F E (xij , t, φ) ⋅ n = 0
⎪ ∂t ∂x
⎨
⎪
⎪
⎪
⎪
⎪ φ si x < xij
⎪
⎪ φ(x, 0) = { i
⎪
⎩ φj si x > xij
x + X 0 (t) t
φ = φE ( , φi , φj ) où X 0 (t) = ∫ u⃗0 (xij , τ ) ⋅ n
⃗ ij dτ
t 0
49
CHAPITRE I.3. LA SPH HYBRIDE
u t u0
u
−
c
F E (xij , φi , φj ) GE (xij , φi , φj )
u+c
φi φj
x
xij
0
qui traduit le fait que l’on ne résout plus le problème de Riemann en xt mais en x+Xt (t) .
La résolution du problème de Riemann sur la raie caractéristique λ0ij = u⃗0 (xij , t) ⋅ n ⃗ ij
0 0
donne un vecteur solution φE (λij , φi , φj ) et un flux eulérien égal à FE (φE (λij , φi , φj )). Le
flux lagrangien s’écrit finalement :
⃗ ij , on le
Le problème de Riemann étant résolu monodimensionnellement sur la direction n
reprojette dans l’espace :
gE (⃗
nij , φi , φj ) = ⃗ ij .
GE (φi , φj ) ⋅ n
50
I.3.1. LA FORMULATION MONOFLUIDE
On remplace alors les valeurs utilisées par le solveur de Riemann à l’interface (φi et
φj ) par les vecteurs φij et φji de variables conservatives corrigés par les gradients :
φij = φi + DR
h
xij − x⃗i )
φi ⋅ (⃗
{
xij − x⃗i )
φji = φj + DR φj ⋅ (⃗
h
Cependant, cette correction ne peut être utilisée telle quelle car le calcul du gradient
peut être inexact localement (entre autres, au voisinage des discontinuités physique où le
gradient lui-même n’est pas défini). Ceci peut créer des extrema locaux qui vont générer
des oscillations dans la solution, oscillations qui peuvent exister numériquement mais qui
ne correspondent pas à la solution physique de problème qui est monotone (à cause de la
viscosité réelle du fluide). L’idée est donc de limiter les flux numériques dans les zones de
forts gradients tout en respectant une condition moins restrictive que la monotonicité, la
Variation Totale Diminuante (Total Variation Diminishing, Harten (1983)) .
Cette condition impose que la variation totale de φ :
∂φ
T V (φ) = ∫ ∣ ∣dx
Ω ∂x
ne doit pas augmenter dans le temps. Il en découle les propriétés suivantes : le nombre
d’extrema de la solution n’augmente pas ; un extremum local n’augmente pas non plus ; un
schéma monotone vérifie cette condition. Un type de limiteur satisfaisant cette condition
TVD est le Minmod (figure I.3.7). Il corrige chaque composante des gradients indépen-
damment avec un facteur β l tel que :
φlij − φli = βijl (φlj − φli )
{
0 ≤ βijl ≤ 1
51
CHAPITRE I.3. LA SPH HYBRIDE
où βijl est en pratique borné par le limiteur le Minmod défini par :
⎧
⎪ DRh l
xij − x⃗i )
φi ⋅ (⃗
⎪
⎪
⎪ β l
= minmod (1, )
⎪
⎪ ij
φj − φli
l
⎨
⎪
⎪
⎪
⎪
⎪
⎪ minmod(a, b) = min(0, max(a, b)) + max(0, min(a, b))
⎩
φlji φlj
φli φlij h l
DR φj
h l
DR φi
x
φlj
φlij
φli φlji
Il existe deux manières d’éviter cet écueil. La première consiste effectivement à com-
parer les états φlij et φlji de manière à détecter cette inversion du problème de Riemann
et par la suite la corriger Marongiu (2007) :
φlij = φli
si (φli − φlj )(φlij − φlji ) < 0 alors {
φlji = φlj
52
I.3.1. LA FORMULATION MONOFLUIDE
Ceci a pour corollaire que le limiteur peut devenir moins efficace si la discontinuité n’est
pas inversée (il permet de réduire, au mieux, celle-ci de moitié).
φl φl
φlj φlj
φli φli
x x
Fig. I.3.9 – Région d’influence des deux variantes du limiteur Minmod : à gauche version
classique, à droite version « papillon ».
p 1
Etot = ei + ec = u∥2 .
+ ∥⃗
ρ(γ − 1) 2
53
CHAPITRE I.3. LA SPH HYBRIDE
La comparaison sera donc effectuée entre la solution exacte et celle de chaque solveur :
1 N ∣Etot
exact
− Etot
solveur
∣
err = ∑ exact
N i=1 Etot
Cas à 2 chocs : Les deux états gauche et droit sont définis par :
⎛1100⎞ ⎛1000⎞
φG = ⎜ 200 ⎟ ; φD = ⎜ 0 ⎟
⎝ 0 ⎠ ⎝ 0 ⎠
La masse volumique de référence ρ0 du fluide est prise égale à 1000kg.m−3 (avec γ = 7).
La vitesse du son est égale à 1466m.s−1 . La solution est formée de deux ondes de choc
voyageant vers la gauche et la droite aux vitesses respectives 1788m.s−1 et 1827m.s−1 .
3.5E+08
3E+08
2.5E+08
2E+08
P
1.5E+08
1E+08 Godunov
HLL
HLLC
5E+07 TRRS
TSRS
Solution exacte
0
-0.02 0 0.02
x
Fig. I.3.10 – Répartition de pression dans le tube à choc no 1 à l’instant t = 10−5 s.
Cas à 2 détentes : Les deux états gauche et droit sont définis par :
⎛ 1100 ⎞ ⎛1000⎞
φG = ⎜−200⎟ ; φD = ⎜ 0 ⎟
⎝ 0 ⎠ ⎝ 0 ⎠
Les propriétés du fluide restent identiques mais cette fois-ci deux détentes vont être
créées, se propageant dans des directions opposées aux vitesses respectives 2151m.s−1
et 1466m.s−1 .
54
I.3.1. LA FORMULATION MONOFLUIDE
3E+08
2.5E+08 Godunov
HLL
HLLC
2E+08 TRRS
TSRS
Solution exacte
1.5E+08
P
1E+08
5E+07
-5E+07
-0.02 0 0.02
x
Fig. I.3.11 – Répartition de pression dans le tube à choc no 2 à l’instant t = 10−5 s.
Cas à un choc et une détente : Les deux états gauche et droit sont définis par :
⎛1100⎞ ⎛1000⎞
φG = ⎜ 0 ⎟ ; φD = ⎜ 0 ⎟
⎝ 0 ⎠ ⎝ 0 ⎠
Les deux états sont donc au repos. La discontinuité de pression existante va générer une
onde de choc qui va se déplacer vers la droite à une vitesse de 1627m.s−1 et une onde de
détente qui va vers la gauche à 1951.s−1 .
3E+08
Godunov
2.5E+08
HLL
HLLC
TRRS
TSRS
2E+08 Solution exacte
P
1.5E+08
1E+08
5E+07
-0.02 0 0.02
x
Fig. I.3.12 – Répartition de pression dans le tube à choc no 3 à l’instant t = 10−5 s.
55
CHAPITRE I.3. LA SPH HYBRIDE
Les différents solveurs donnent des solutions très proches. Les différences entre eux sont
négligeables devant l’erreur qu’ils commentent tous par rapport à la solution exacte. Cette
erreur relative sur l’énergie entre la solution exacte et les solutions données par les solveurs
est néanmoins relativement faible : de l’ordre de 0.9%, 9.9% et 8.1% respectivement pour
les trois cas tests.
3E+08
1er ordre
2.5E+08 Limiteur strict
Limiteur papillon
Avec inversion
2E+08
Solution exacte
P
1.5E+08
1E+08
5E+07
-0.02 0 0.02
x
Fig. I.3.13 – Répartition de pression dans le tube à choc no 3 à l’instant t = 10−5 s. Influence
du schéma M.U.S.C.L. et des différents limiteurs.
Sont confrontés à la solution exacte, un calcul avec un solveur de Riemann exact 1er
ordre (c’est-à-dire sans utilisation du schéma M.U.S.C.L) puis des calculs avec différents
limiteurs du schéma M.U.S.C.L. : un limiteur Minmod strict, le même en version « pa-
pillon » et un dernier en version sans traitement particulier des inversions des problèmes
de Riemann. On constate cette fois des différences notables entre les solutions obtenues.
Δerr × 10−6
Avec inversion réf.
1er ordre -278008
Limiteur strict -165571
Limiteur papillon -21088
Tab. I.3.1 – Différence d’erreur relative de l’énergie totale Etot de chaque type de limi-
teur par rapport à celle du solveur de Godunov avec un limiteur Minmod autorisant les
inversions.
56
I.3.1. LA FORMULATION MONOFLUIDE
D’un point de vue énergétique (cf. table (I.3.1)), il ressort clairement que l’absence de
schéma M.U.S.C.L. engendre beaucoup de dissipation dans le calcul de la solution numé-
rique. L’ajout de ce schéma est donc indispensable. La façon de le limiter via le Minmod,
et plus particulièrement le traitement des inversions, est aussi une source de dissipation :
adopter un limiteur strict paraı̂t être un peu trop restrictif (et par conséquent dissipatif).
En comparaison le limiteur « papillon » est environ 8 fois moins dissipatif. Quant à ignorer
complètement les inversions, cela permet de minimiser l’erreur sur l’énergie du système.
D’un point de vue qualitatif (figure I.3.13), l’utilisation d’un limiteur Minmod sans
traitement des inversions paraı̂t être la meilleure façon de capturer les discontinuités de
pression, tant pour les chocs que dans les détentes.
57
CHAPITRE I.3. LA SPH HYBRIDE
P
0.8 Godunov 6000
1.81
5333.33
Y
0.6 H2 H1 1.8
4666.67
1.645
1.65
1.655
X 4000
Y
0.4
3333.33
2666.67
2000
0.2
1333.33
666.667
0
2
X
2.5 3 0
Y
0.4 0.4
0.2 0.2
0 0
2 2.5 3 2 2.5 3
X X
0.4 0.4
0.2 0.2
0 0
2 2.5 3 2 2.5 3
X X
Fig. I.3.14 – Champ de pression pour les différents solveurs au temps t(g/H)1/2 = 6, 75
Temps CPU
Godunov ref.
HLL -13,1%
HLLC -0,4%
TRRS -4,6%
TSRS -5,4%
Tab. I.3.2 – Écart relatif par rapport au solveur de Godunov sur le temps CPU global de
la simulation.
Les gains restent toutefois assez faibles (de l’ordre de 5%), ce qui n’est pas suffisant
pour être décisif dans le choix du solveur utilisé par la suite. On privilégie la précision en
58
I.3.2. EXTENSION À DEUX FLUIDES ET DIFFICULTÉS
Godunov
HLL
HLLC
TRRS
TSRS
Godunov
HLL
HLLC
TRRS
TSRS
Fig. I.3.15 – Hauteur d’eau totale en (x/H)1 = 0, 825 et (x/H)2 = 1, 653 (distances indi-
quées à partir du bord droit). Comparaison aux calculs de Colagrossi et Landrini (2003) :
bifluide ligne continue noire sans symboles, monofluide ligne pointillée sans symboles et
aux expériences de Zhou et alii (1999) : (△).
59
CHAPITRE I.3. LA SPH HYBRIDE
0
Godunov
-0.05 HLL
HLLC
-0.1 TRRS
(Etot − Etot0 )/ΔE
TSRS
-0.15
-0.2
-0.25
-0.3
-0.35
-0.4
0 2 4 6 8
t(g/H)1/2
Fig. I.3.16 – Évolution de l’énergie au cours du temps pour les différents solveurs.
travail effectué plus haut sur la construction du solveur de Riemann, en prenant soin de
ne pas faire de simplifications (à cause des lois d’état différentes pour chaque fluide). Il
en résulte alors de nouveaux flux, asymétriques, entre les particules de part et d’autre de
l’interface.
On autorise donc, par cette approche de solveur de Riemann interfacial, des flux entre
les particules et tout particulièrement des flux de masse à travers l’interface (comme
ailleurs). On pourrait penser, à priori, que cela a une faible incidence dans les configura-
tions choisies. Or quand on fait un simple test de maintien hydrostatique d’une colonne
eau/air, on obtient une solution divergente, l’effort appliqué d’un fluide sur l’autre crois-
sant rapidement.
Pour comprendre ce phénomène, on revient à un cas monofluide. On peut en effet
observer ce phénomène en étudiant le cas simple d’un fluide soumis à la pesanteur dans
un domaine (1m × 1m) avec une surface libre et en observant la répartition de la masse,
cf. figure I.3.17.
Sur ce cas particulier, le schéma M.U.S.C.L. n’a pas été employé pour accentuer le
phénomène. L’absence de ce schéma conduit à résoudre les problèmes de Riemann entre
des états discontinus alors que le champ moteur de l’écoulement (le champ de pression
hydrostatique) est en première approximation linéaire (aux effets de compressibilité près).
On introduit alors une dissipation excessive qui force le transfert de masse.
Enfin sur ce cas, la présence d’une condition aux limites de type mur avec glissement
est traitée par la technique des particules fantômes. Comme cela a, par ailleurs (§I.2.1.3.1),
été mentionné, ces particules fantômes voient leurs masses corrigées par la variation de
densité hydrostatique. Ceci crée donc un flux de masse entrant à la paroi qui va traverser
le domaine. Les parois verticales ne faisant pas intervenir de correction particulière (la
pesanteur est orthogonale à la paroi), le flux ne peut sortir par ces frontières. La dernière
frontière est la surface libre, où la condition de non-pénétration impose un flux de masse
nul à travers celle-ci. Il n’y a donc pas de condition de flux sortant sur ce domaine : la
masse va s’y accumuler. Les particules de la surface libre n’étant pas contraintes (condi-
tion de surface libre cinématique), celles-ci vont pouvoir se déplacer en hauteur et, par
60
I.3.3. LA FORMULATION À FRACTION DE VOLUME
Y
0 0.5
1.51
1.2
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
X X
61
CHAPITRE I.3. LA SPH HYBRIDE
α1 + α2 = 1.
Cette condition de saturation traduit le fait que le domaine en volume est complètement
occupé par les deux fluides (figure I.3.18).
α2
α1
Fig. I.3.18 – Schéma de principe de la description de l’espace par deux fractions de volume.
Nous nous intéressons à un cas particulier de cette condition (celle qui concerne les
écoulements à interface libre) où α1 est nul dans toute la partie occupée par le fluide 2
(et inversement pour α2 dans le fluide 1). Cette condition est vérifiée à l’approximation
numérique près (la fraction volumique est proche de zéro). On peut ainsi supposer qu’il
n’existe qu’un seul champ de vitesse u⃗ commun aux deux fluides.
∂ ρ̃K
⃗ ⋅ (ρ̃K u⃗)
+∇ = 0. (I.3.14)
∂t
62
I.3.3. LA FORMULATION À FRACTION DE VOLUME
On peut aussi introduire la densité volumique de masse de mélange ρ (bien qu’au niveau
de la description continue les fluides ne se mélangent pas) telle que ρ = ρ̃1 + ρ̃2 . Elle vérifie
l’équation :
∂ρ
+∇ ⃗ ⋅ (ρ⃗
u) = 0.
∂t
Cependant, avec l’utilisation des lois d’état qui seront introduites ultérieurement, cette
équation seule n’est pas suffisante pour décrire correctement le problème : les deux équa-
tions (I.3.14) sont nécessaires.
où la grandeur mélangée Xm est le barycentre des grandeurs des espèces 1 et 2, affectées
respectivement de leur fraction volumique (notées plus légèrement avec α = α1 et α2 =
1 − α).
Cette loi de mélange introduit donc dans le système une nouvelle inconnue α qui peut
suivre une relation algébrique (reliant cette fraction aux masses volumiques partielles) ou
une équation aux dérivées partielles. Pour ce genre de relation, exprimant la conservation
de la fraction volumique (ou densité volumique de volume) du fluide 1, les modèles VOF,
Hirt et Nichols (1981), adoptent une équation de transport :
∂α
⃗
+ u⃗ ⋅ ∇α = 0.
∂t
63
CHAPITRE I.3. LA SPH HYBRIDE
Pour les fluides incompressibles, celle-ci peut se mettre sous la forme d’une équation de
conservation :
∂α
+∇⃗ ⋅ (α⃗
u) = 0.
∂t
Mais pour les fluides pseudo-compressibles, cette équation de conservation n’est pas for-
cément justifiée. Une autre relation pour la fraction de volume sera introduite ultérieure-
ment.
Pression de mélange : Avec la loi de mélange introduite éq. (I.3.17), nous pouvons
maintenant définir la pression de mélange p∗ :
p∗ = αp1 + (1 − α)p2 .
64
I.3.3. LA FORMULATION À FRACTION DE VOLUME
Viscosité de mélange : De la même façon qu’une pression de mélange a été définie pour
l’équation de quantité de mouvement, il est nécessaire de préciser le tenseur des contraintes
visqueuses Tv . On suppose qu’il est défini comme celui d’un seul fluide Newtonien de
mélange (on néglige la partie antisymétrique du tenseur parce qu’on ne s’intéresse qu’à la
solution incompressible des écoulements) :
Tv = ⃗ ⊗ u⃗)
μ(∇
avec une viscosité dynamique de mélange vérifiant la loi de fermeture :
μ = α∗ μ1 + (1 − α∗ )μ2 . (I.3.18)
65
CHAPITRE I.3. LA SPH HYBRIDE
Dω D⃗
x
= ω div(⃗
u) et = u⃗.
Dt Dt
⎧
⎪
⎪
⎪ ∂α p1 − p2
⎪
⎪ ⃗
+ u⃗ ⋅ ∇(α) =
⎪ ∂t
(R) ⎨
⎪
⎪
⎪ ∂φ ∂F ∂G
⎪
⎪ + + = S
⎪
⎩ ∂t ∂x ∂y
⎧
⎪
⎪
⎪
⎪
∂α
⃗
⎪
⎪ + u⃗ ⋅ ∇(α) = 0
⎨ ∂t
⎪
⎪
⎪ ∂φ ∂F ∂G
⎪
⎪ + + = S
⎪
⎩ ∂t ∂x ∂y
La solution issue de la première équation est notée α̃. Nous verrons à l’étape suivante qu’il
n’est pas nécessaire de l’expliciter (car elle ne correspond pas à l’équilibre des pressions) et
par conséquent de définir un schéma numérique pour sa résolution. On ne s’occupe donc
que de la résolution de la deuxième équation.
Cette équation hyperbolique étant la même que celle sur laquelle a été basée la dé-
monstration du schéma SPH hybride monofluide, nous pouvons en reprendre les résultats
66
I.3.3. LA FORMULATION À FRACTION DE VOLUME
en terme de discrétisation. Ainsi le schéma discret associé à cette étape de transport sera :
⎧
⎪
⎪
⎪
⎪
D⃗
xi
⎪
⎪ = u⃗i
⎪
⎪
⎪
Dt
⎪
⎪
⎪ Dωi
⎨ = ωi ∑(⃗ ⃗ ij ωj
uj − u⃗i ) ⋅ ∇W
⎪
⎪
⎪
Dt j
⎪
⎪
⎪
⎪
⎪
⎪
Dωφi
⃗ ij ∥ ωj + ωi S
⎪
⎪ = −ωi ∑ 2 gE (⃗ nij , φi , φj ) ∥∇W
⎪
⎩ Dt j
Ce schéma discret est résolu numériquement avec un schéma d’avance en temps de type
Runge-Kutta d’ordre 4 et un schéma de décentrement en espace de type solveur de Rie-
mann exact (donnant le flux gE (⃗ nij , φi , φj )). Ce dernier sera détaillé spécifiquement pour
les équations étudiées au prochain paragraphe.
La solution issue de ce schéma est notée φ̃ et va servir au calcul de la deuxième étape.
√
⎧
⎪ q − q̃(ρ̃1 , ρ̃2 ) + (q − q̃(ρ̃1 , ρ̃2 ))2 + 4ρ̃1 ρ̃2 c21 c22
⎪
⎪
⎪ γ + (ρ̃1 , ρ̃2 ) =
⎪
⎪
⎪ ρ̃2 c22
γ + (ρ̃1 , ρ̃2 ) ⎪
⎪
⎪
α∗ (ρ̃1 , ρ̃2 ) = avec ⎨ q̃(ρ̃ , ρ̃ ) = ρ̃2 c22 − ρ̃1 c21
1 + γ + (ρ̃1 , ρ̃2 ) ⎪
⎪
⎪
1 2
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩ q = ρ02 c22 − ρ01 c21
Il est important de noter que le schéma discret de la méthode proposée ne fait alors
appel à aucun schéma numérique d’évolution de la fraction de volume.
67
CHAPITRE I.3. LA SPH HYBRIDE
Onde de valeur propre λ = u − c : Pour une onde de détente (soit p∗g ≤ pg ), les
invariants de Riemann sont :
⎧
⎪ cg ρ1∗g ρ2∗g cg ρ1g ρ2g
⎪
⎪
⎪ u∗g + log ( ) = ug + log ( )
⎪
⎪ 2 ρ 01 ρ02 2 ρ 01 ρ02
⎨
⎪
⎪
⎪ ρ1∗g ρ1g
⎪
⎪
⎪ =
⎩ ρ2∗g ρ2g
On en déduit les densités partielles de l’état étoilé :
⎧
⎪ ug − u∗g
⎪
⎪
⎪ ρ1∗g = ρ1g exp ( )
⎪
⎪ cg
⎨
⎪
⎪ ug − u∗g
⎪
⎪
⎪ ρ 2∗g = ρ 2g exp ( )
⎪
⎩ cg
où la vitesse de l’état étoilé est donnée par :
pg − p˜0g
u∗g = ug + fg (φ′g , p∗g ) où fg (φ′g , p∗g ) = cg log ( ).
p∗g − p˜0g
Si, cette fois-ci, une onde de choc (p∗g > pg ) est présente à gauche, les invariants de
Riemann sont :
⎧
⎪ ρ1∗g (u∗g − s) = ρ1g (ug − s)
⎪
⎪
⎨ ρ2∗g (u∗g − s) = ρ2g (ug − s)
⎪
⎪
⎪
⎩ ρg (ug − s)(ug − u∗g ) = −(pg − p∗g )
68
I.3.3. LA FORMULATION À FRACTION DE VOLUME
Avec les deux premières relations et après quelques calculs, nous pouvons en déduire les
densités partielles :
⎧
⎪ ug − s
⎪
⎪
⎪ ρ1∗g = ρ1g
⎪
⎪ u∗g − s
⎨
⎪
⎪
⎪ ug − s
⎪
⎪
⎪ ρ2∗g = ρ2g
⎩ u∗g − s
où la vitesse s de l’onde de choc est déduite via la pression de mélange :
u∗g − ug
s = u∗g + ρg c2g .
p∗g − pg
En pratique, celle-ci n’est pas calculée numériquement avec cette relation pour éviter les
problèmes dans les discontinuités de contact (où u∗g = ug et p∗g = pg ). On utilise plutôt la
relation suivante : √
2 ρg
s = u∗g − cg .
p∗g − pg
Enfin, le troisième invariant de Riemann permet d’obtenir la vitesse de l’état étoilé :
p∗g − pg
u∗g = ug + fg (φ′g , p∗g ) où fg (φ′g , p∗g ) = − √ .
ρg (p∗g − p˜0g )
69
CHAPITRE I.3. LA SPH HYBRIDE
De même, on peut montrer que le saut de fraction volumique est aussi quelconque à
travers la discontinuité de contact. Par contre, ce saut devient nul à travers un choc ou
une détente. Pour résumer, si l’on se trouve à gauche de la discontinuité on a α∗g = αg
(puisque α est constante à travers la 1-onde), et inversement on a α∗d = αd à droite de la
discontinuité.
Solution complète : On peut réécrire la relation sur les vitesses au niveau de la dis-
continuité de contact de la façon suivante :
1 1
u∗ = (ug + ud ) + (fd (φ′d , p∗ ) − fg (φ′g , p∗ )).
2 2
À l’intérieur des détentes, nous avons les relations suivantes, pour la détente gauche
(état « 1/6 ») :
⎧ ⎧
⎪ ug − u
⎪
⎪
⎪
x ⎪
⎪
⎪ ρ = ρ exp ( )
⎪ u = t + cg
⎪ ⎪
⎪ 1 1g
cg
⎨ avec ⎨
⎪
⎪
⎪ ⎪
⎪
⎪ ug − u
⎪
⎪ p = αg p1 (ρ1 ) + (1 − αg )p2 (ρ2 ) ⎪
⎪
⎪ ρ 2 = ρ 2g exp ( )
⎩ ⎩ cg
⎧ ⎧
⎪
⎪
⎪ x ⎪
⎪ ud − u
⎪
⎪ u = t − cd
⎪ ⎪
⎪ ρ1 = ρ1d exp (− cd )
⎪
⎨ avec ⎨
⎪
⎪
⎪ ⎪
⎪
⎪ ud − u
⎪
⎪ p = αd p1 (ρ1 ) + (1 − αd )p2 (ρ2 ) ⎪
⎪
⎪ ρ 2 = ρ 2d exp (− )
⎩ ⎩ cd
Les détentes, qui sont des ondes qui s’étalent dans le temps, sont comprises entre :
dx dx
ug − cg < < u∗ − c g et u∗ + cd < < ud + c d
dt dt
toutes les vitesses ayant été définies auparavant.
70
I.3.3. LA FORMULATION À FRACTION DE VOLUME
⃗ i
∇φ
∇ ⃗ i∥ , n
⃗ l φi = min(∥∇φ ⃗ j )⃗
⃗ i ⋅ ∇φ ni ⃗i =
où n
jD ⃗ i∥
∥∇φ
Ce nouveau gradient limité ∇ ⃗ l φi est par la suite utilisé dans le schéma M.U.S.C.L..
Cette correction a été pensée pour être utilisable en deux dimensions. Elle reste coû-
teuse (elle oblige à effectuer une sommation supplémentaire sur le voisinage, entre le calcul
des gradients et le calcul final des dérivées temporelles) mais s’avère indispensable pour
mener à bien les simulations.
71
CHAPITRE I.3. LA SPH HYBRIDE
i
)n⃗
φj
⋅ ∇⃗
i
i
(n⃗
φi φj
72
I.3.3. LA FORMULATION À FRACTION DE VOLUME
5E+08
Tait avec γ = 7
Tait avec γ = 1, 1
4E+08
3E+08
Linéarisée
2E+08
p 1E+08
-1E+08
950 1000 1050 1100 1150
ρ
Fig. I.3.20 – Variation de la pression pour les différentes lois d’état dans la gamme de
masse volumique étudiée. La vitesse du son vaut 1466m.s−1 et la masse volumique de
référence 1000kg.m−3 .
discontinuité physique de manière à pouvoir utiliser la formulation dans une version où
le gradient nécessaire au schéma M.U.S.C.L. n’est pas limité (correction présentée plus
haut). Néanmoins, la fraction de volume ne peut pas être prise égale à 0 ou 1 pour des
raisons numériques : elle est donc initialisée à 10−7 uniformément sur tout le domaine.
Les résultats sont présentés sur les figures I.3.21, I.3.22 et I.3.23.
2.5E+08
2E+08
1.5E+08
P
1E+08
SPH-VF MUSCL non limité
SPH-VF MUSCL limité
Riemann Tait MUSCL
SPH-VF sans MUSCL
5E+07 Riemann Tait sans MUSCL
Solution exacte
73
CHAPITRE I.3. LA SPH HYBRIDE
1E+08
P
5E+07
-5E+07
1E+08
5E+07
-0.02 0 0.02
x
lors de l’étape de relaxation (α ∉ [0, 1] causé par des densités partielles négatives . . .). Cette
formulation est clairement fausse mais néanmoins elle fournit des résultats intéressants
(symboles respectifs ◯ et ◇) : ceux-ci sont proches de la formulation équivalente mono-
fluide (avec solveur de Riemann et équation de Tait). Ceci est à analyser prudemment
mais laisserait à penser que la résolution de la partie hyperbolique du schéma SPH-VF
soit en un sens correcte.
Enfin, nous pouvons remarquer la forte dégradation de la précision apportée par la
74
I.3.3. LA FORMULATION À FRACTION DE VOLUME
correction qui vient limiter les gradients avant la procédure M.U.S.C.L. (symboles ◻). On
retrouve alors une précision légèrement supérieure à celle d’une formulation sans schéma
M.U.S.C.L. mais encore fortement éloignée des résultats utilisant le schéma M.U.S.C.L.
non limité. Mais, comme signalé au moment de la présentation de cette correction, celle-ci
reste indispensable pour que le schéma numérique SPH-VF fonctionne correctement.
60
50
40
ρu
30
20 h/Δx = 1
h/Δx = 1, 4286
10 h/Δx = 2
Solution exacte
0
-0.2 0 0.2 0.4 0.6
x
Fig. I.3.24 – Influence du rapport h/Δx (avec pour un ratio croissant, respectivement 2,
4 et 6 voisins) sur la précision de la résolution du problème de Riemann du tube à choc
no 5.
Les propriétés des deux fluides sont les suivantes : le fluide 1 possède une masse
volumique au repos de 1kg.m−3 et sa vitesse du son est de 3m.s−1 . La masse volumique
du fluide 2 vaut 1000kg.m−3 et sa vitesse du son est de 15m.s−1 . Les états initiaux sont
donnés sous forme du vecteur de variables convectives φv = (α, ρ1 , ρ2 , u)t .
75
CHAPITRE I.3. LA SPH HYBRIDE
Cas d’une discontinuité de contact convectée : Les deux états gauche et droit sont
définis par :
⎛1 − 10 ⎞
−7 −7
⎛ 10 ⎞
⎜ 1 ⎟ ⎜ 1 ⎟
φvG = ⎜ ⎟ ; φvD = ⎜ ⎟
⎜ 1000 ⎟ ⎜1000⎟
⎝ 0, 15 ⎠ ⎝ 0, 15 ⎠
La discontinuité est initialement positionnée en x = 0, 25m. Cet état initial crée une
discontinuité de contact qui va être simplement advectée dans le domaine à une vitesse
de 0, 15m.s−1 . Cela permet de tester la diffusion numérique du schéma.
Sans le schéma M.U.S.C.L, on remarque sur la figure I.3.25 que la formulation eulé-
rienne du schéma SPH-VF est proche de celle Volumes Finis. Avec le schéma M.U.S.C.L.,
elle devient beaucoup diffusive qu’en Volumes Finis, du fait de la correction.
Par contre, dans la formulation SPH-VF purement lagrangienne, on observe comme
attendu que la discontinuité est positionnée exactement quelle que soit l’utilisation du
M.U.S.C.L., et ce sans diffusion même sur la fraction volumique (cf. figure I.3.26). On
voit ici l’intérêt du caractère lagrangien couplé à cette formulation.
Cas à un choc et une détente : Les deux états gauche et droit sont définis par :
⎛1 − 10 ⎞
−7 −7
⎛ 10 ⎞
⎜ 100 ⎟ ⎜ 1 ⎟
φvG = ⎜ ⎟ ; φvD = ⎜ ⎟
⎜ 104 ⎟ ⎜1000⎟
⎝ 0 ⎠ ⎝ 0 ⎠
76
I.3.3. LA FORMULATION À FRACTION DE VOLUME
x
Fig. I.3.26 – Répartition de la fraction volumique dans le tube à choc no 4 à l’instant
t = 3, 333s.
La discontinuité est initialement positionnée en x = 0, 3m. Cet état initial va créer une
détente qui va se déplacer vers la gauche et un choc qui va aller à droite. La discontinuité
de contact se déplace peu et reste à sa position initiale (x ≃ 0, 3m).
Pour la discontinuité de contact, on peut remarquer (voir figure I.3.28) le même phé-
nomène que pour le cas-test précédent : le caractère lagrangien du schéma SPH-VF permet
de bien la capturer même si elle commence à diffuser légèrement du fait de sa faible vitesse.
Quant au choc et à la détente, on observe (figure I.3.27) néanmoins une plus grande
diffusion qu’avec le schéma Volumes Finis.
On notera que les décalages en espace qui peuvent être observés entre nos résultats et
ceux des Volumes Finis proviennent certainement d’un recalage incorrect de ces dernières
données car, sur notre formulation SPH-VF, les discontinuités de contact sont positionnées
exactement (respectivement pour chaque tube à x = 0, 15 × 3, 333m et x = 0, 3m).
77
CHAPITRE I.3. LA SPH HYBRIDE
ρu
x
Fig. I.3.27 – Répartition de la quantité de mouvement dans le tube à choc no 5 à l’instant
t = 0, 03s. Légende des symboles : cf. figure ci-dessous.
_ Solution exacte
78
Deuxième partie
79
Dans ce chapitre nous allons aborder la modélisation de deux phénomènes physiques
distincts : la viscosité (réelle cette fois-ci) et les effets de tension superficielle. Dans l’étude
des équations de Navier-Stokes, on peut séparer l’étude de la partie hyperbolique de ces
équations (Euler) et les parties diffusives et capillaires. Les modélisations de ces effets
diffusifs et capillaires proposées dans cette partie pourront s’appliquer aux deux approches
distinctes précédentes pour modéliser la partie hyperbolique.
81
82
Chapitre II.1
La viscosité
1
Tvij = μij 2
xi − x⃗j ) ⊗ (⃗
[(⃗ ui − u⃗j ) + (⃗
ui − u⃗j ) ⊗ (⃗
xi − x⃗j )]
rij
où μij est la viscosité dynamique entre les deux particules i et j, tandis que rij est la
distance entre elles.
Puis on obtient l’effort visqueux volumique (pour rappel :ρ Du⃗ = divTv ) agissant sur
Dt
une particule en appliquant simplement l’opérateur dérivée d’interpolation discret au ten-
seur (et en remarquant que celui-ci est symétrique, i.e. Tvji = Tvij ) :
divTvij = Dh (Tv ).
1 1
divTvij = ∑ μij ( ⃗ ij
+ ) Tvij ⋅ ∇W
j ni nj
divTvij = ⃗ ij ωj .
∑ μij 2Tvij ⋅ ∇W
j
83
CHAPITRE II.1. LA VISCOSITÉ
Cette constante peut être calculée analytiquement et vaut : C th = Δ(f (⃗ x)) = 4. La formu-
lation de Morris P.S.E.R. s’écrit finalement :
4 (⃗ ⃗ ij ωj
xi − x⃗j ) ⋅ ∇W
divTvij = ∑ 2μij 2
ui − u⃗j ).
(⃗
Ci j rij
84
II.1.2. VISCOSITÉ DYNAMIQUE INTER-PARTICULAIRE
On pourra se référer à Cueille (2005) pour une étude de la précision des 3 derniers
schémas (Monaghan, Morris et Morris P.S.E.R.) sur des fonctions scalaires sur différents
types de maillages.
On peut noter que pour un fluide donné (masse volumique, viscosité dynamique, vitesse
du son), le raffinement de la taille h des éléments de discrétisation (via h/Δx) peut mener
très rapidement à restreindre le pas de temps général à cause de la condition « visqueuse »
(terme en h2 ). Pour éviter cela, Cueille (2005) propose d’impliciter la contrainte visqueuse
et obtient une réduction notable du temps de calcul. Cependant, la prise en compte des
effets de tension superficielle va imposer une condition sur les pas de temps sensiblement
aussi restrictive (terme en h3/2 ). Dans cette première approche on n’a donc pas cherché à
spécifiquement reformuler le problème comme Cueille pour s’en affranchir.
85
CHAPITRE II.1. LA VISCOSITÉ
II.1.5 Validation
Nous allons ici valider ces différents schémas et estimer leur précision sur un cas
académique d’écoulement visqueux : l’écoulement de Poiseuille plan.
L’écoulement s’effectue dans un canal périodique (pour que le profil stationnaire s’éta-
blisse tout en limitant le domaine de calcul) à un régime laminaire. L’élément moteur est
un gradient de pression assimilable à un champ de force volumique parallèle à l’axe du ca-
nal (celui-ci peut être vu comme un champ de pesanteur, l’écoulement étant alors vertical
descendant).
Le bilan de quantité de mouvement projeté sur l’axe e⃗y s’écrit donc :
∂P ∂ ∂v
= 0 soit 0 = (μ(x) ) − ρ(x)g.
∂y ∂x ∂x
Cette formulation reste généraliste (en terme de répartition de viscosité et densité dans
l’espace) de manière à pouvoir en dériver la solution aussi bien dans le cas monofluide
(avec le même fluide dans tout le canal) que bifluide (avec deux fluides X et Y occupant
des régions distinctes, l’interface étant évidemment parallèle à l’axe d’écoulement).
Les différentes conditions cinématiques en découlant vont servir de condition aux
limites pour l’intégration de l’équation précédente. Sur les parois (en x = 0 et x = L), nous
avons une condition de vitesse nulle (v0 = 0 et vL = 0), tandis qu’à l’interface (en x = xI ),
si présente, on a continuité de la vitesse tangentielle (vIX = vIY ).
Dans le cas monofluide (fluide X ), le profil de vitesse analytique du régime stationnaire
s’écrit, cf. e.g. Cueille (2005) :
ρX g
v(x) = − x(L − x).
2μX
Dans le cas bifluide (avec le fluide X entre 0 et xI , et entre xI et L le fluide Y), le profil
devient, cf. e.g. Chanteperdrix (2004) :
ρ(x)g τI
v(x) = (x − xI )2 − (x − xI ) − vI
2μ(x) μ(x)
gxI ρX μY − ρY μX ρY g x2I τI
τI = et vI = + xI .
2 μX + μY μY 2 μY
Une fois calculé ce profil analytique de vitesse, nous pouvons en déduire le débit à
travers la section du canal :
L
Q = ∫0 ρ(x)v(x)dx.
On peut aussi en déduire la vitesse maximale dans le canal ainsi que le nombre de
Reynolds correspondant. Pour le cas monofluide, la vitesse maximale vm est atteinte en
x = L/2 :
2
ρX g L 2 ρX L3
vm = X et Re = ( X ) g .
2μ 4 μ 8
86
II.1.5. VALIDATION
Pour le cas bifluide, la recherche d’extremum de la vitesse analytique donne comme posi-
τI
tion : xm = xI + . On remarque alors que, si τI est positif (i.e. ν Y > ν X ), le maximum
ρ(x)g
global se situe alors pour x > xI , dans le fluide Y :
Y τI2 ρY L Y
vm = −( + vI ) et ReY = ( )v .
2μY ρY g μY 2 m
Pour l’autre fluide X , le maximum local se situe donc à l’interface, en x = xI , puisque la
solution est monotone.
X ρX L X
vm = vI et ReX = ( X ) vm .
μ 2
Et inversement si ν Y < ν X .
87
CHAPITRE II.1. LA VISCOSITÉ
inférieure ; mais à faible résolution spatiale cela peut poser problème lors de la création
des particules fantômes). Les points sont initialement disposés cartésiennement. Le noyau
est la gaussienne (avec un rapport h/Δx = 1, 3298). La formulation choisie ici par défaut
est la formulation SPH-Multifluide.
Remarque :
Exceptionnellement, nous avons choisi ici la loi d’état linéarisée pour la formulation SPH-
Multifluide. Sur les autres cas présentés par la suite, il a été vérifié que le choix de la loi
d’état n’influe pas sur les résultats. La loi d’état sera donc celle de Tait pour la SPH-
Multifluide, tandis que sera celle linéarisée pour la SPH-VF (sauf mention contraire).
Nous allons étudier l’influence de divers paramètres numériques sur la précision (celle-
ci sera calculée relativement à la solution analytique détaillée plus haut). L’écoulement
s’effectuera à un nombre de Reynolds de 1 (soit μ = 1000kg.m−1 .s−1 et vm = 1m.s−1 ).
où errv est l’erreur relative sur la vitesse par rapport à l’expression analytique sur l’en-
semble des N particules échantillonnées dans la largeur du canal. On peut alors étu-
dier la convergence spatiale et on observe (voir figure II.1.2) que la formulation SPH-
Multifluide converge à un ordre proche de 2. L’utilisation de la renormalisation pour le
noyau (§I.3.1.2.4) permettrait d’atteindre ce taux de convergence. On observera son effet
avec l’autre formulation SPH-VF.
Remarque :
La convergence dont il est question ici (et par la suite) est une convergence vers la solution
analytique. L’ordre de convergence en espace devrait reposer sur une convergence vers
une solution numérique « saturée » (où l’augmentation de la résolution spatiale n’aurait
88
II.1.5. VALIDATION
0 0.1
0.08
20 points
-0.02
40 points
0.06 80 points
160 points
errQ
errV
0.04
-0.04
20 points
40 points 0.02
80 points
160 points
Solution exacte 0
-0.06
-0.02
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3
y t
(a) Erreur relative sur la vitesse à t = 3s (b) Erreur relative sur le débit
-2
10
-3
errRM S
10
SPH-Multifluide
10
-4 1er ordre
2e ordre
3e ordre
20 40 60 80 100 120 140 160
N /L
aucune incidence sur la précision de la solution). Ceci n’a été fait qu’incomplètement
ici puisqu’on n’observe pas encore de saturation de la courbe de convergence aux plus
hautes résolutions étudiés. Il serait donc nécessaire d’effectuer des calculs plus fins (et
plus coûteux). Néanmoins, même si l’on n’obtient pas l’ordre de convergence absolu du
schéma, cela permet d’avoir une idée de la vitesse de convergence aux résolutions utilisées
en pratique.
89
CHAPITRE II.1. LA VISCOSITÉ
0.15
0.4
Flekkøy
Monaghan
Flekkøy Morris
0.1 Monaghan Morris P.S.E.R.
Morris
Morris P.S.E.R.
errQ
errV
Solution exacte
0.05 0.2
0
-0.05
-0.1
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3
y t
(a) Erreur relative sur la vitesse à t = 3s (b) Erreur relative sur le débit
errRM S
Flekkoy 2, 1.10−2
Monaghan 5, 4.10−3
Morris 4, 3.10−4
Morris P.S.E.R. 2, 9.10−4
Tab. II.1.1 – Erreur relative quadratique moyenne.
Les schémas les plus précis sont ceux basés sur la formulation de Morris : ils sont 1 à 2
ordre de magnitude plus précis que les schémas de Monaghan et Flekkøy (respectivement
précis à 0, 5% et 2% contre au mieux 0, 03% pour Morris). La renormalisation P.S.E.
apporte effectivement une amélioration de la précision.
90
II.1.5. VALIDATION
0.005 0.1
0 0.08
Gaussien
-0.005 0.06
Cubique
errQ
errV
-0.01 0.04
Gaussien
-0.015
Cubique 0.02
Solution exacte
-0.02 0
-0.025 -0.02
0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3
y t
(a) Erreur relative sur la vitesse à t = 3s (b) Erreur relative sur le débit
91
CHAPITRE II.1. LA VISCOSITÉ
10-2
errRM S
10-3
SPH-Multifluide
SPH-VF
SPH-VF renormalisé
10-4
pente 2e ordre
pente 2e ordre
20 40 60 80 100 120 140 160
N /L
On observe (voir figure II.1.5) que la formulation SPH-VF converge avec un ordre
légèrement supérieur à 2 (à noyau identique). Néanmoins, on peut remarquer qu’à faible
résolution, cette formulation est moins précise que la formulation SPH-Multifluide. La
92
II.1.5. VALIDATION
dissipation excessive induite par le limiteur de gradient (cf. §I.3.3.2.3) pourrait expliquer
ce manque de précision.
Quant à la renormalisation du gradient du noyau, son utilisation semble mener au
contraire du résultat attendu (augmentation de la précision spatiale qui permettrait de se
rapprocher du 2e ordre théoriquement).
Remarque :
Par ailleurs nous avons étudié aussi l’influence des coefficients CCF La et CCF Lv pour les
deux formulations. Avec le schéma d’avance en temps employé (Runge-Kutta 4e ordre),
le coefficient de la condition restrictive sur l’acoustique CCF La peut être égal à 2,5 pour
la formulation SPH-Multifluide (cf. Colagrossi et Landrini (2003)) et égal à 4,6 pour la
formulation SPH-VF (voir l’analyse faite par Guilcher (2008)). Ceci a été vérifié sur ce cas.
Il a aussi été testé l’influence du coefficient de restriction sur les effets visqueux CCF Lv :
la différence relative entre un calcul mené avec ce coefficient égal à 0,1 et un autre calcul
avec CCF Lv = 1 est inférieure à 2.10−4 .
Fluide ρ0 μ vm Re
X 1000 5 2, 7.10−3 0,769
Y 100 5.10−2 7, 7.10−3 0,027
La vitesse du son est prise égale à 0.1m.s−1 dans les deux fluides. La discrétisation spa-
tiale est de 40 points dans la largeur du canal. La formulation utilisée est SPH-Multifluide
avec un noyau gaussien.
93
CHAPITRE II.1. LA VISCOSITÉ
0.4 0.4
Moyenne harmonique Moyenne harmonique
Moyenne arithmétique Moyenne arithmétique
0.3 0.3
0.2 0.2
errV
errV
0.1 0.1
0 0
Fig. II.1.6 – Erreurs relatives sur le canal de Poiseuille bifluide suivant la viscosité dyna-
mique inter-particulaire.
vont être aussi légèrement influencées par l’autre viscosité dans le cas de la moyenne har-
monique. Mais globalement cette moyenne permet d’avoir une meilleure précision globale
(l’erreur relative quadratique moyenne vaut 0,1% contre 1,6% pour la moyenne arithmé-
tique).
Pour la formulation SPH-VF (figure II.1.6(b)), on remarque que l’erreur sur la vi-
tesse est, quelle que soit la moyenne utilisée, proche de la moyenne arithmétique. Afin
d’expliquer une partie de l’origine de cette erreur, il est nécessaire de rappeler que, dans
cette formulation, la viscosité dynamique de chaque particule est déjà estimée comme un
barycentre de la viscosité dynamique de chaque fluide (I.3.18), pondéré par la fraction vo-
lumique. Or celle-ci va avoir tendance à diffuser à l’interface de telle sorte que la viscosité
dynamique de chaque particule va être déjà influencée par la viscosité de l’autre fluide
(via la fraction volumique).
L’observation de cette fraction volumique en fin de simulation (figure II.1.7(a)) montre
que celle-ci aura plus diffusé à l’interface avec l’utilisation d’une moyenne harmonique pour
la viscosité dynamique inter-particulaire. On peut supposer que ce schéma étant plus précis
sur la vitesse que celui de la moyenne arithmétique (cf. figure II.1.6(b)), il va engendrer une
diffusion plus forte à l’interface. Ce qui aura pour conséquence de lisser encore davantage le
saut de viscosité dynamique (figure II.1.7(b)) et de minimiser la contrainte de cisaillement
à l’interface. D’où un schéma légèrement plus précis. Néanmoins, l’erreur reste grande dans
le fluide Y et diminuera donc la précision de cette formulation.
94
II.1.5. VALIDATION
1.2 6
Moyenne harmonique Moyenne harmonique
Moyenne arithmétique Moyenne arithmétique
1 5
0.8 4
0.6 3
μ
α
0.4 2
0.2 1
0 0
-0.2 -1
0.035 0.04 0.045 0.05 0.055 0.06 0.065 0.035 0.04 0.045 0.05 0.055 0.06 0.065
y y
(a) Fraction volumique à t = 10s (b) Viscosité dynamique particulaire à t = 10s
-2
10
errRM S
10-3
SPH-Multifluide
SPH-VF
SPH-VF renormalisé
pente 2e ordre
-4
pente 2e ordre
10
20 40 60 80 100 120 140 160
N
Figure (II.1.8) on retrouve les tendances observées lors des simulations monofluides
pour la formulation SPH-Multifluide. Celle-ci converge avec un ordre proche de 2. Par
contre, comme mentionné plus haut pour la formulation SPH-VF, l’erreur générée à l’in-
terface se propage dans le fluide Y et dégrade fortement la précision du calcul.
95
CHAPITRE II.1. LA VISCOSITÉ
Remarque :
Conclusion
Concernant la validation de la modélisation des effets visqueux, nous avons donc mené
celle-ci sur les deux formulations et dans deux configurations, monofluide et bifluide. Si, sur
le cas monofluide, les différences entre les deux formulations ne sont pas très importantes,
elles le deviennent par contre sur le cas bifluide, et sont toujours à l’avantage de SPH-
Multifluide. Comme cela a été souligné, la dissipation excessive induite par la correction
faite dans la formulation SPH-VF impacte fortement ses résultats en bifluide.
Il serait aussi intéressant d’effectuer d’autres cas de validation des effets purement
visqueux sur des problèmes possédant une solution analytique mais présentant une moins
forte directionnalité de l’écoulement. Un cas couramment employé est celui de la cavité
entrainé. Cependant il ne peut pas être décliné en version bifluide. Nous verrons plus loin
(au chapitre III.2) le cas des instabilités de Rayleigh-Taylor qui peut être utilisé en tant
que validation des effets visqueux bifluides.
96
Chapitre II.2
97
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
⃗
n
x⃗
Tts ⋅ x⃗
t⃗
∂ΩI
Fig. II.2.1 – Propriété de projection sur la tangente t⃗ du tenseur Tts .
Comme il a été déjà mentionné, cette approche a été aussi adoptée par la commu-
nauté SPH. Il a été alors remarqué que cette force volumique de tension de surface pouvait
engendrer localement des pressions négatives Morris (2000). Or la méthode SPH est for-
tement sensible à la présence locale de pressions négatives qui conduit à des appariements
indésirés de particules Monaghan (2000). Il a été proposé par Hu et Adams (2006) de
modifier le tenseur Tts pour éviter ces pressions négatives de la façon suivante :
1
⃗ ( I−n
Tts = σ∥∇C∥ ⃗I ⊗ n
⃗I )
d
où d est la dimension de l’espace. Hors des bénéfices numériques attendus (examinés plus
loin lors de la validation), il est important de remarquer que l’introduction du terme 1/d
change la propriété de projection sur le plan tangent à l’interface du tenseur des effets de
tension de surface : celle-ci n’est plus respectée et l’on ne peut donc plus parler stricto
sensu d’efforts de tension de surface.
98
II.2.2. MODÉLISATION ALTERNATIVE
Dans les deux formulations, le tenseur discret Ttsi est calculé à partir du gradient de
⃗ i avec la relation (II.2.1).
la fonction indicatrice ∇C
99
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
Suivant Brackbill et alii (1992), Morris (2000) et Hu et Adams (2006) proposent une
condition restrictive sur les pas de temps basée sur la relation suivante :
1/2
min(ρX , ρY ) h3
Δtts ≤ CCF Lts ( )
2πσ
où ces auteurs prennent la constante CCF Lts égale à 0, 25.
D’autre part, les lois d’état choisies imposent une condition restrictive sur le choix de
la vitesse du son numérique dans l’écoulement Chanteperdrix (2004). Celle-ci doit assurer
que les variations de pressions engendrées par les effets capillaires restent négligeables :
σ
≪ 1. (II.2.3)
ρc2 R
II.2.5 Validation
Nous allons aborder ici la validation du schéma C.S.S.. Nous ne présenterons ici que
des résultats issus de l’implémentation du schéma avec la formulation SPH-Multiphase,
n’ayant pas du tout réussi à coupler ce schéma avec la formulation SPH-VF (les investi-
gations sur ce problème n’ont pas été poussées plus loin faute de temps).
100
II.2.5. VALIDATION
Δx N∞ N2 N1
1/96 4, 96.10−2 7, 74.10−3 2, 28.10−3 avec 1/d
1/128 5, 06.10−2 8, 20.10−3 2, 42.10−3
1/160 5, 71.10−2 8, 31.10−3 2, 43.10−3
1/192 5, 32.10−2 8, 42.10−3 2, 46.10−3
101
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
0.15 0.15
0.025 0.025
0.1 0.1
0.05 0.05
y
y
0 0
-0.05 -0.05
-0.1 -0.1
-0.15 -0.15
-0.1 0 0.1 -0.1 0 0.1
x x
(a) Résolution Δx = 1/96 (b) Résolution Δx = 1/192
Fig. II.2.2 – Courants parasites à t = 2.10−3 s (sans 1/d). Le trait plein matérialise l’inter-
face pour faciliter l’appréciation du phénomène.
restent donc présents dans les simulations qui seront présentées avec ce schéma par la
suite.
102
II.2.5. VALIDATION
1.5
0.2
Y
-0.2
-0.2 0 0.2
X
103
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
vitesse
0.08 0
0.075
-0.2
0.07
-0.4
0.065
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
t t
(a) Convergence sur la position (b) Convergence sur la vitesse
Remarque :
104
II.2.5. VALIDATION
Comme noté dans Morris (2000), à haute résolution, le schéma C.S.S. implémenté en SPH
peut générer localement des pressions négatives. Pour éviter cela, il peut être utilisé la
procédure de contrôle de l’instabilité de tension Monaghan (2000). Au lieu d’implémenter
cette correction supplémentaire purement numérique, nous avons préféré tirer parti du fait
que ces simulations ne possédant pas de surface libre, on peut librement ajouter une pres-
sion de fond 2 . Cela ne change pas la solution théoriquement mais son introduction dans
l’équation d’état permet d’éviter d’avoir des pressions relatives négatives et les instabilités
numériques qui s’ensuivent.
Cependant, il a été remarqué que l’utilisation d’une pression de fond de valeur trop
grande (typiquement la pression atmosphérique) ne permettait pas de poursuivre correc-
tement les simulations. En effet, lors de l’adimensionnalisation des équations (cf. §I.1.3),
la longueur de référence choisie est, pour des raisons numériques (liées à l’utilisation du
ratio r/h, cf. §I.2.1.1.1), la longueur de lissage h du schéma d’interpolation. Celle-ci est
forcément très petite devant la longueur physique de référence du problème (dans le sens
où la résolution numérique du problème impose des éléments de discrétisation petits de-
vant la taille du domaine de calcul). La pression de fond peut donc être très différente de
la pression de référence découlant de l’adimensionnalisation, et ainsi poser problème d’un
point de vue purement numérique.
Ainsi, la détermination de cette pression de fond n’est pas aisée. Dans ce cas de
validation en apesanteur, la seule pression de référence du problème (à savoir le saut de
pression à travers l’interface, ΔP = σ/R) ne suffit pas à avoir des simulations stables, et
la pression de fond choisie est prise égale à 50P a.
105
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
0.35
Analytique
Numérique
0.3
Erreur relative
0.25
τ [s], errτ
0.2
0.15
0.1
0.05
0
0.6 0.8 1 1.2 1.4
σ
Fig. II.2.5 – Comparaison entre la période d’oscillation donnée par la relation analytique
(II.2.5) et les résultats de la simulation numérique et erreur relative. Résolution spatiale :
60 × 60 points.
Les résultats numériques (cf. figure II.2.5) obtenus sont satisfaisants. La tendance
issue de la relation (II.2.5) est correctement respectée, à 3% près.
L’ordre de convergence en espace de cette formulation a été étudié : voir figure (II.2.6).
Celui-ci est compris entre le 1er et le 2e ordre.
0.1
0.08
0.06
0.04
errτ
0.02
SPH-Multifluide
pente 1er ordre
pente 2e ordre
106
II.2.5. VALIDATION
x - avec x - avec
y - avec y - avec
x - sans x - sans
y - sans y - sans
0.09 0.09
position
position
0.08 0.08
0.07 0.07
0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 0.3 0.4 0.5
t t
(a) 900 points (b) 14400 points
Fig. II.2.7 – Influence du terme 1/d sur la période d’oscillation suivant la résolution
spatiale. Cas à faible ratio de densité.
107
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
0.35
Analytique
Numérique (avec 1/d)
Numérique (sans 1/d)
0.3 Erreur relative (avec 1/d)
Erreur relative (sans 1/d)
0.25
τ [s], errτ
0.2
0.15
0.1
0.05
0
0.6 0.8 1 1.2 1.4
σ
Fig. II.2.8 – Comparaison entre la période d’oscillation donnée par la relation analytique
(II.2.5) et les résultats de la simulation numérique suivant la présence du terme 1/d et
erreur relative. V0 = 2m.s−1 .
D’autre part, il n’a été observé aucune différence sur la forme de l’interface suivant la
présence de la correction.
Par contre, à plus forte résolution (sur le premier cas à faible ratio de densité, des
simulations à 56700 et 230400 points ont été conduites), des instabilités apparaissent effec-
tivement à l’interface (figure II.2.10), ce qui dégrade fortement la précision du modèle de
tension superficielle. Mais le fait de rajouter le terme correctif de Hu & Adams n’améliore
pas la qualité de la pression à l’interface.
Conclusion
Nous avons mené ici la validation des effets de tension superficielle sur des cas tests
de goutte oscillante en apesanteur avec deux ratio de densité différents. Il en ressort
que le schéma C.S.S. couplé à la formulation SPH-Multifluide donne des résultats très
satisfaisants en terme de prédiction de la pulsation d’oscillation de la goutte, et ce malgré
le fait que des courants parasites (artefact numérique) existent et n’aient pas été corrigés.
Toutefois, il convient d’être prudent dans l’utilisation du schéma à haute résolution
où les corrections proposées n’apportent pas de réponses satisfaisantes. Il est nécessaire
d’être aussi attentif dans la détermination de la pression de fond.
Dans la dernière partie sera présenté un cas de validation plus complexe faisant inter-
venir ce schéma C.S.S. avec succès (cf. §III.4.2).
108
II.2.5. VALIDATION
Fig. II.2.9 – Influence du terme 1/d sur le champ de pression de la goutte oscillante à
différents instants dimensionnels. Faible ratio de densité. Résolution de 14400 points. Les
deux colonnes de droite (avec 1/d) sont aux mêmes instants que les deux colonnes de
gauche (sans 1/d).
Fig. II.2.10 – Instabilités à haute résolution (230400 points). Influence du terme 1/d sur
le champ de pression de la goutte oscillante. Faible ratio de densité.
109
CHAPITRE II.2. LES EFFETS DE TENSION SUPERFICIELLE
110
Troisième partie
111
Chapitre III.1
Ballottement linéaire
y
h1 ⃗
a
ρX
g⃗
0
ρY
−h2 x
0 L
113
CHAPITRE III.1. BALLOTTEMENT LINÉAIRE
a0 L 4
ξ= (x − + ∑ 2
cos(ω2n+1 t) cos(k2n+1 x))
g 2 n≥0 Lk2n+1
où
⎧
⎪ nπ
⎪
⎪ kn =
⎪ L
⎨ gkn Δρ
⎪
⎪
⎪ ωn =
⎪
⎩ ρ1 coth(kn h1 ) + ρ2 coth(kn h2 )
Cette solution reste valable tant que l’hypothèse de régime linéaire des oscillations est
vérifiée, c’est-à-dire que l’amplitude du mode le plus énergétique (pour n = 0) reste petite
devant la largeur L de la cuve et devant la hauteur h1 et inférieure à h2 , soit :
a0 a0 L a0 L
≪1 , ≪ 1 et < h2
g gh1 g
La fraction volumique est initialisée telle que dans la partie supérieure du domaine on
ait α = 1 − et dans la partie inférieure α = avec = 10−6 .
Le nombre de modes utilisés pour comparer les résultats numériques à la solution
analytique mentionnée ci-dessus est pris égal à n = 100. On s’intéresse aux élévations de
surface libre sur les parois, soit en x = 0 et en x = L.
La discrétisation spatiale est choisie de façon à avoir un certain nombre d’éléments par
longueur d’onde du premier mode (soit deux fois la largeur de la cuve). Deux résolutions
ont été étudiées de manière à effectuer les comparaisons avec Chanteperdrix (2004) (la
convergence en espace n’a donc pas été le critère déterminant pour justifier ces résolu-
tions) : la première avec 40 points dans la largeur de la cuve et la deuxième avec deux fois
plus de points. Le pas spatial est identique pour chaque direction de l’espace.
Remarque :
Il est nécessaire de noter que le pas de résolution choisi par Chanteperdrix est différent
dans les deux directions de l’espace (il choisit un nombre d’éléments égal à 40 × 100 et
80 × 150 alors que le domaine a pour dimensions 1 × 2, 25m). De manière à assurer la plus
grande précision possible de l’interpolation SPH, nous avons choisi d’avoir une répartition
114
III.1.3. RÉSULTATS
III.1.3 Résultats
III.1.3.1 Avec interface
Dans un premier temps, nous nous plaçons dans la configuration décrite ci-dessus afin
de pouvoir comparer les résultats de la formulation SPH-VF avec ceux issus du formalisme
Volumes Finis (cf. Chanteperdrix (2004)).
t t
Fig. III.1.2 – Comparaison de l’élévation de surface libre sur les parois gauche (x =
0 courbes du bas) et droite (x = L courbes du haut) entre la solution analytique, la
formulation Volumes Finis à fraction de volume Chanteperdrix (2004) et SPH-VF
On peut observer (cf. figure III.1.2) que sur un maillage grossier l’amplitude des os-
cillations est très fortement atténuée dès la deuxième oscillation. Sur un maillage plus
fin, l’amortissement est plus réduit pour la première oscillation mais devient non négli-
geable pour les suivantes. On retrouve le comportement de la formulation tel qu’il avait
été observé lors de l’étude monodimensionnelle des tubes à choc.
Pour compléter cette observation et pour écarter le biais dans la validation que peut
constituer la diffusion de l’interface (phénomène qui n’est pas quantifié à proprement
115
CHAPITRE III.1. BALLOTTEMENT LINÉAIRE
parler avec le choix de suivi des particules proches de l’interface), le même cas va être
conduit sans interface mais avec une surface libre.
Remarque :
Il convient de noter que, pour les trois formulations SPH-VF utilisées (avec et sans
limiteur, et (pi + pj )) , la partie du domaine modélisant l’eau (soit le bas de la cuve)
contient aussi l’autre fluide (l’air) de manière négligeable (exactement du point de vue
de la fraction volumique).
Tous les schémas ont été utilisés avec un noyau renormalisé. Comme cela a déjà été
évoqué (cf. §I.3.1.2.4), le critère de passage de la matrice de renormalisation à la matrice
identité est basé sur une des valeurs propres de cette matrice (limite prise ici égale à 0, 6).
Les résultats sont présentés sur les figures III.1.3 et III.1.4. Ils sont quantifiés dans le
tableau III.1.1 et divers enseignements peuvent en être tirés.
Tout d’abord, la formulation « SPH-VF (pi +pj ) », qui correspond au schéma SPH-VF
non décentré (sans viscosité artificielle ni solveur de Riemann), possède un comportement
extrêmement proche de la formulation SPH classique sans viscosité artificielle, comme
attendu. La formulation avec fraction de volume est donc aussi stable (voire plus précise)
que la formulation SPH classique pour cette configuration particulière. Cette dernière est
précise sur un maillage grossier mais le raffinement ne lui fait pas gagner beaucoup plus
de précision.
116
III.1.3. RÉSULTATS
0 0.5 1 1.5 2
SPH hybride monofluide 0
0.005
0
ξ
-0.005
0 0.5 1 1.5 2
t
Fig. III.1.3 – Comparaison de l’élévation de surface libre sur les parois gauche (x = 0
courbes du bas) et droite (x = L courbes du haut) entre la solution analytique et les
diverses formulations. Résolution grossière.
errξ
grossier fin
0 L 0 L
SPH classique 0,13 0,12 0,13 0,11
SPH hybride monofluide 0,19 0,21 0,12 0,12
SPH-VF (pi + pj ) 0,13 0,12 0,10 0,10
SPH-VF sans limiteur 0,26 0,38 0,14 0,16
SPH-VF + limiteur 0,68 1,03 0,30 0,40
Tab. III.1.1 – Erreur quadratique moyenne en temps de l’erreur relative d’élévation de
surface libre sur chaque paroi 0 et L.
Ensuite, lorsque pour ces deux formulations (SPH-VF et SPH classique) on introduit
un solveur de Riemann pour décentrer le schéma en espace (et deviennent donc « SPH-VF
sans limiteur » et « SPH hybride monofluide »), on observe une légère dissipation numé-
rique. L’écart aux résultats des formulations non décentrées se réduit avec l’augmentation
de la résolution spatiale. La formulation SPH-VF sans limiteur est toutefois légèrement
plus dissipative que la formulation SPH hybride monofluide.
Enfin, lorsque le limiteur de gradient est introduit dans la formulation SPH-VF, la
dissipation devient non négligeable (on a entre 30% et 100% d’erreur). Néanmoins il faut
noter que cet ajout est nécessaire pour pouvoir simuler des écoulements bifluides (cf.
remarques lors de son introduction §I.3.3.2.3).
117
CHAPITRE III.1. BALLOTTEMENT LINÉAIRE
0 0.5 1 1.5 2
SPH hybride monofluide 0
0.005
0
ξ
-0.005
0 0.5 1 1.5 2
t
Fig. III.1.4 – Comparaison de l’élévation de surface libre sur les parois gauche (x = 0
courbes du bas) et droite (x = L courbes du haut) entre la solution analytique et les
diverses formulations. Résolution fine.
Remarque :
À titre d’illustration, le temps CPU pour ce cas est donné tableau III.1.2.
Temps CPU (s)
SPH classique 188
SPH hybride monofluide 377
SPH-VF (pi + pj ) 1126
SPH-VF sans limiteur 2481
SPH-VF + limiteur 3939
Tab. III.1.2 – Temps de calcul CPU pour la résolution grossière. Le coefficient CFL
acoustique CCF Lac est égal à 1 pour toutes les formulations.
Il est important de noter que le coefficient CFL a été pris identique pour tous les
calculs, alors qu’a priori il peut être beaucoup plus important pour les formulations avec
solveur de Riemann (facteur proche de 2 en se basant sur les résultats de l’écoulement de
Poiseuille).
Conclusion
Comme mentionné en introduction, ce cas test est assez sévère pour la formulation
SPH-VF complète à cause de la dissipation numérique à l’interface. Mais en restreignant
118
III.1.3. RÉSULTATS
la configuration avec seulement une surface libre, nous montrons que les majeures parties
composant la formulation SPH-VF permettent d’avoir des résultats satisfaisants : ainsi,
la version « SPH-VF (pi + pj ) » démontre les bonnes capacités du formalisme à fraction
de volume choisi. Et l’introduction d’un solveur de Riemann pour résoudre correctement
la partie hyperbolique (la version « SPH-VF sans limiteur ») permet de conserver des
résultats proches de son modèle équivalent sans fraction de volume (à savoir « SPH hybride
monofluide »). Mais toutes ces bonnes propriétés sont dégradées lorsque la correction par
limiteur est ajoutée, comme cela avait été montré par ailleurs sur les cas de tube à choc.
La partie « fautive » est donc parfaitement identifiée. Mais n’ayant pu trouver d’autre
moyen d’effectuer des simulations bifluides sans s’en passer, nous la conservons, comme
nous l’avons déjà mentionné.
119
CHAPITRE III.1. BALLOTTEMENT LINÉAIRE
120
Chapitre III.2
Instabilités de Rayleigh-Taylor
On étudie maintenant un cas de validation classique bi-dimensionnel, qui est celui des
instabilités de Rayleigh-Taylor. On notera que ce cas se place dans le cadre des équations
de Navier-Stokes. La tension superficielle joue en revanche un rôle faible dans cet écoule-
ment et on peut la négliger. Dans ce problème d’instabilités, l’interface entre deux fluides
visqueux différents nécessite d’être calculée précisément. En effet, initialement le fluide
lourd se situe dans la partie supérieure du domaine, au-dessus du fluide léger. Afin de se
retrouver dans une configuration stable, le fluide lourd va glisser sous le fluide léger en
créant des enroulements d’interface dus au cisaillement visqueux, dont la capture nécessite
d’avoir une description précise de l’évolution de l’interface entre les deux fluides.
III.2.1 Configuration
Le domaine du problème est rectangulaire (de hauteur 2m et de largeur 1m) et l’in-
terface est située à y = 1 − 0, 15 sin(2πx). Elle sépare le fluide lourd (de densité 1, 8kg.m−3 )
situé au dessus, et le fluide léger (ρ = 1kg.m−3 ) situé au dessous. Le nombre de Reynolds
basé sur la largeur du domaine est égal à 420.
La vitesse du son numérique est prise égale à 14m.s−1 pour le fluide léger et à 10m.s−1
pour le fluide lourd. L’équation d’état choisie est celle linéarisée pour les deux formulations.
Pour la formulation SPH-VF, il est appliqué une pression de fond correspondant à la
pression hydrostatique (soit p0 = 3, 6P a).
Une condition d’adhérence est appliquée sur les parois.
Dans cet écoulement, on peut supposer que les effets de tension superficielle sont
négligeables. Elle n’est pas modélisée par le schéma C.S.S.. Seule la formulation SPH-
Multiphase utilisera le paramètre I , pris égal à 0,02, pour éviter une inter-pénétration
des phases.
Pour cette dernière formulation, l’interface est repérée par les particules dont la valeur
propre de la matrice de renormalisation (calculée en ne considérant que les particules
voisines du même fluide) est supérieure à une certaine valeur seuil (ici égale à 0,6). Le
noyau d’interpolation est la gaussienne.
Pour la formulation SPH-VF, la notion d’appartenance à l’un ou l’autre fluide n’étant
présente qu’implicitement, l’interface devrait être repérée par l’isocontour de la fraction de
volume égal à 0,5. Cependant, les résultats donnés par la méthode SPH étant ponctuels
(définis sur les particules), la détermination des champs sur tout le domaine de calcul
(pour calculer les contours) ferait appel à une interpolation supplémentaire. Pour éviter
121
CHAPITRE III.2. INSTABILITÉS DE RAYLEIGH-TAYLOR
cela, nous avons tiré parti de la légère diffusion de l’interface au cours du temps et nous
n’affichons que les particules dont la fraction volumique est comprise entre 0,1 et 0,9
(choix totalement arbitraire mais suffisant pour avoir une approximation de la position
de l’interface). Le noyau choisi est la spline cubique (cf. paragraphe suivant pour le choix
de ce noyau).
1 1 1
y
y
0 0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x x
(a) SPH-Multifluide (b) Level-Set (c) SPH-VF
Fig. III.2.1 – Convergence spatiale pour les deux formulations SPH et comparaison à une
solution Level-Set. Instant t = 5s. Résolutions SPH : ● ∶ 75×150, ◆ ∶ 150×300, ◂ ∶ 300×600.
Résolutions Level-Set : 22 ∶ 78 × 156, 2 ⋅ 2 ∶ 156 × 312, 2 ∶ 312 × 624.
Les résultats sont regroupés sur la figure III.2.1. On peut noter que les ordres de
convergence semblent comparables (la convergence est peut-être un peu plus rapide en
SPH-Multifluide). Par contre, sur cette figure, on est de droite à gauche de plus en plus
précis, dans le sens où l’on est de plus en plus proche de la solution convergée alors qu’on
a des résolutions identiques. La méthode SPH-Multifluide est ainsi déjà convergée avec la
résolution moyenne, contrairement aux autres.
L’étude de la vorticité (cf. figure III.2.2), montre que celle-ci est générée à la fois à
l’interface et sur les parois, et est responsable de la complexité des enroulements créés. La
comparaison entre la formulation SPH-Multifluide et le Level-Set montre une similarité
122
III.2.2. RÉSULTATS NUMÉRIQUES
dans les formes de cette vorticité, et par conséquent de la forme de l’interface. Cependant,
la vorticité SPH est plus bruitée que la vorticité Level-Set, ce qui est lié à la répartition
désordonnée des particules.
L’analyse de l’intensité de la vorticité montre que, sur les deux formulations (SPH
et Différences Finies), la répartition est proche, même si, dans certaines régions (e.g.
au milieu du domaine), le Level-Set exhibe un champ plus fort (et qui est la cause de
l’enroulement plus intense qu’en SPH). Le Level-Set a aussi tendance à présenter une
vorticité plus diffusée dans l’espace : ceci peut être dû à la fonction Level-Set qui tend à
lisser les quantités entre les deux fluides, ou au fait que cette vorticité est probablement
dissipée plus rapidement en SPH.
(a) SPH-Multifluide
(b) Level-Set
123
CHAPITRE III.2. INSTABILITÉS DE RAYLEIGH-TAYLOR
1.5 1.5
1 1
y
0.5 0.5
0 0
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x x
(a) SPH-Multifluide (b) SPH-VF
Fig. III.2.3 – Importance du choix du noyau d’interpolation pour les deux formulations
SPH. Instant t = 5s. Résolution : 150 × 300. ● : spline cubique, ● : gaussienne
124
III.2.2. RÉSULTATS NUMÉRIQUES
125
CHAPITRE III.2. INSTABILITÉS DE RAYLEIGH-TAYLOR
α
1.5
0.95
0.9
0.85
0.8
0.75
0.7
0.65
1 0.6
y
0.55
0.5
0.45
0.4
0.35
0.3
0.25
0.5 0.2
0.15
0.1
0.05
0
0 0.5 1
x
(a) SPH-VF (b) SPH-VF triangulé (c) VOF
Fig. III.2.4 – Formulations SPH-VF et Volumes Finis VOF. Résolution 75 × 150. Instant
t = 5s.
Cela est d’autant plus notable que, pour la formulation VOF, la fraction de volume
étant convectée par une équation de transport propre, des techniques numériques sont
utilisées pour limiter la diffusion à l’interface à chaque pas de temps et reconstruire celle-
ci. Aucune de ces techniques n’est utilisée dans la formulation SPH-VF.
Conclusion
La simulation des instabilités de Rayleigh-Taylor avec les deux formulations SPH a
permis d’observer comparativement leur comportement sur un cas d’écoulement bidimen-
sionnel avec une interface pleinement développée. Il en ressort que la SPH-Multifluide
se montre beaucoup plus précise que la méthode SPH-VF, avec une convergence rapide.
La comparaison à des solutions de codes maillés classiques sur des résolutions identiques
permet en outre de conclure positivement quant à la validation de cette méthode sur ce
cas.
126
Chapitre III.3
Envahissement
Les cas abordés dans ce chapitre vont permettre de valider les formulations sur des
configurations où l’écoulement est plus violent. Les deux cas étudiés concernent des pro-
blèmes d’envahissement. On retrouve cette problématique dans les questions de surviva-
bilité des navires en cas d’avarie. Un des paramètres importants entrant alors en jeu est
le degré de confinement de la partie envahie par rapport à l’atmosphère environnante.
Afin de maı̂triser au mieux les conditions de réalisation de ces simulations, la géométrie
est simplifiée à l’extrême et aboutit à étudier l’envahissement d’un volume de type caisson.
Comme dans les études précédentes, nous nous limitons à l’approche bi-dimensionnelle de
ce problème.
III.3.1.1 Configuration
Comme présenté en introduction, la configuration choisie est relativement simple, cf.
schémas III.3.1(a) & III.3.1(b). D’un coté du domaine, il y a la partie que nous pouvons
appeler réservoir et qui contient le fluide « envahissant ». Séparée par une paroi d’une
certaine épaisseur se trouve, de l’autre coté du domaine, la zone prête à être envahie et
dont le confinement va pouvoir être modulé. En effet, tandis que le réservoir est toujours
laissé à pression atmosphérique, la partie de la zone envahie peut être plus ou moins fermée
avec des couvercles de taille variable (afin de faire varier la perte de charge représentant
le confinement plus ou moins imparfait d’un navire). Dans le cadre de cette étude nous
nous limiterons à deux particuliers : avec confinement total (zone envahie complètement
fermée) ou sans (zone envahie à la pression atmosphérique).
Les cotations du domaine sont mentionnées dans le tableau III.3.1.
Du point de vue numérique, les simulations étant conduites principalement en prenant
en compte les deux fluides, pour assurer la condition de pression atmosphérique on intro-
127
CHAPITRE III.3. ENVAHISSEMENT
X X X X
duit une surface libre air/vide au-dessus de l’air (d’où la position arbitraire de la partie
supérieure du domaine en milieu ouvert et dans la partie réservoir).
D’autre part, la configuration géométrique présente des angles aigus (du point de vue
du fluide modélisé) qui peuvent être problématiques dans l’implémentation algorithmique.
Il a été adopté la modélisation présentée par Le Touzé et alii (2006).
Les propriétés des fluides sont les suivantes : pour l’eau, la densité est de 1000kg.m−3
et la vitesse du son de 45m.s−1 . Pour l’air, la densité est prise égale à 1kg.m−3 et la vitesse
du son à 636m.s−1 .
Pour une raison de temps disponible, la formulation choisie a été SPH-Multifluide.
Pour que ce schéma reste stable numériquement, il a été choisi des viscosités dynamiques
issues de la relation suivante :
1
μ = αρch
8
où α est pris égal à 0,1. Les viscosités dynamiques ont donc été choisies telles que :
μX = 2, 25kg.m−1 .s−1 et μY = 2, 5.10−2 kg.m−1 .s−1 .
Sur ce cas les effets de tension superficielle vont être négligeables devant les effets iner-
tiels (W e = ρU 2 L/σ > 105 en se basant sur une vitesse maximale relevée numériquement
de 6m.s−1 , sur le coefficient de tension superficielle air-eau (σ = 0, 073N.m−1 ) et sur la
largeur de la zone envahie). Le coefficient I est donc pris égal à 0,08.
Exceptionnellement pour cette formulation, le noyau choisi a été la spline cubique pour
une question de coût de calcul. Pour simplifier l’implémentation, la résolution spatiale a
été choisie de telle façon que le rayon du support d’interaction soit strictement inférieur à
128
III.3.1. CAS DU LMF
l’épaisseur Eparoi de la paroi entre le réservoir et la zone envahie (de manière à ce que le
support d’interpolation ne soit pas à cheval sur les deux parties du domaine). La longueur
de lissage est donc prise égale à h = 4.10−3 m.
129
CHAPITRE III.3. ENVAHISSEMENT
Fig. III.3.2 – Zone envahie ouverte. Champ de densité. Instants dimensionnels. En haut :
formulation SPH-Multifluide avec un seul fluide modélisé. Au milieu : SPH-Multifluide
avec deux fluides modélisés. En bas : résultats expérimentaux.
130
III.3.1. CAS DU LMF
Fig. III.3.3 – Zone envahie fermée. Champ de densité. Instants dimensionnels. Colonnes
de gauche : formulation SPH-Multifluide avec deux fluides modélisés. Colonnes de droite :
résultats expérimentaux.
131
CHAPITRE III.3. ENVAHISSEMENT
Celle-ci pourrait être relevée en différents endroits sur les parois de la zone envahie pour
donner accès à la fois à la pression au sein de l’air et à la pression du liquide lorsque
celui-ci arrive sur ces parois.
Lcaisson
Hcaisson
Htirant eau
Hreservoir X
Dtrappe
Lreservoir
Fig. III.3.4 – Configuration du domaine.
Dans son article Ruponen utilise une méthode de correction de pression, basée sur la
résolution des équations de Bernoulli. Il se place dans une configuration tri-dimensionnelle
et ramène alors le problème à une résolution globale de la hauteur d’eau dans le caisson
en se basant sur les deux sections que sont les surfaces libres du réservoir et du caisson, en
voyant la trappe comme un simple orifice en paroi mince occasionnant une perte de charge
singulière de coefficient connu. Outre ce coefficient, les seules données qu’il se donne sont
le tirant d’eau du caisson (qui reste inchangé dans le temps dans son approximation),
132
III.3.2. CAS DE RUPONEN
De la même façon que dans le cas précédent, nous avons traité la paroi du caisson
séparant deux parties du domaine (réservoir et zone envahie) en imposant une certaine
épaisseur à celle-ci pour simplifier l’implémentation.
L’eau est caractérisée de la manière suivante : la densité est de 1000kg.m−3 et la vitesse
du son de 150m.s−1 . Sa viscosité est choisie telle que α = 0, 1 soit μX = 93, 75kg.m−1 .s−1 .
Le noyau utilisé est la spline cubique avec une longueur de lissage égale à h = 0, 05m.
Par ailleurs, on tire parti de la symétrie du domaine en ne simulant qu’un demi-domaine.
III.3.2.2 Résultats
Le caisson n’étant pas confiné dans cette configuration, on cherche à déterminer le
remplissage de celui-ci au cours du temps. Afin de comparer à la hauteur globale prédite
par Ruponen, on estime la hauteur d’eau moyenne dans le caisson (hauteur qui n’est pas
uniforme avec notre modélisation locale) à partir du volume d’eau entrant en comptant
simplement les particules entrantes. Les résultats sont présentés sur la figure III.3.5.
Nous pouvons observer que la tendance générale de la hauteur de remplissage est
correctement respectée. La méthode SPH-Multifluide prédit une évolution plus lente. En
fin de simulation, à t = 35s, l’équilibre n’est pas atteint entre la hauteur d’eau dans le
caisson et la hauteur de charge imposée correspondant au tirant d’eau (ici 4m). Cependant
ces premiers résultats sont à considérer avec précaution étant donné qu’aucune étude de
convergence spatiale n’a été menée et que la discrétisation est relativement grossière (6
particules dans la demi-largeur de la trappe). De plus, on rappelle que les hypothèses
étant très différentes, on ne s’attend pas à obtenir une évolution strictement identique du
remplissage.
De façon parallèle à la remarque de fin du paragraphe précédent, pour compléter la
comparaison entre les deux méthodes numériques, il serait intéressant d’étudier la pression
dans le cas d’un caisson confiné (cas abordé par Ruponen (2006)). Cela n’a pas encore été
effectué par manque de temps.
Conclusion
Ces deux cas tests sont intéressants dans la mesure où ils ont permis de confronter
la formulation SPH-Multifluide à des écoulements assez violents où des fragmentations et
133
CHAPITRE III.3. ENVAHISSEMENT
SPH-Multifluide
Fig. III.3.5 – Évolution temporelle de la hauteur d’eau dans le caisson. Comparaison aux
résultats issus de Ruponen (2006).
134
Chapitre III.4
Bulle isolée
Dans ce chapitre nous allons aborder des problèmes plus spécifiques et plus proches
du but de notre travail, l’étude de la séparation des fluides. Nous allons nous intéresser
à l’écoulement d’une bulle de gaz isolée et ascendante dans une colonne d’un liquide au
repos. Dans un premier temps, il s’agit d’une bulle avec un ratio de densité élevé (par
rapport au fluide environnant) mais avec des effets de tension superficielle faibles. La
comparaison est faite avec des résultats de la littérature. Dans un deuxième temps, on
étudie l’écoulement d’une bulle avec un ratio de densité faible et des effets de tension
superficielle importants. La convergence de la vitesse d’ascension est alors étudiée.
III.4.1 Sussman
III.4.1.1 Problème et configuration
Pour ce premier problème, nous nous plaçons dans une configuration qui a été étudiée
en SPH par Colagrossi et Landrini (2003), sur la base de travaux de Sussman et alii
(1994). Ces derniers ont été effectués par une méthode de type Volumes Finis Level-Set.
Il s’agit d’une bulle d’air remontant dans une colonne d’eau au repos. Le domaine est
représenté sur la figure III.4.1.
L’écoulement est caractérisé
√ par les nombres adimensionnels suivants : le nombre de
Reynolds (défini par Re = (2R) gρX /μX et valant 1000), et le nombre de Bond (défini par
3
Bo = 4ρX gR2 /σ et valant 200). Ce dernier caractérise les effets de tension superficielle (par
rapport aux effets d’inertie) et, ici, sa valeur élevée signifie que ces effets sont négligeables.
Le rayon de la bulle est pris égal à R = 0, 025m. L’accélération de la pesanteur est
prise égale à g = 9, 81m.s−2 . Une condition de glissement est imposée sur les parois, comme
dans Sussman et alii (1994).
Les deux fluides sont définis par les caractéristiques suivantes : pour l’eau, la densité
est de 1000kg.m−3 , la viscosité dynamique égale à 3, 5.10−2 kg.m−1 .s−1 et la vitesse du son
égale à 16m.s−1 . Pour le fluide gazeux on a : ρY = 1kg.m−3 , μY = 4, 5.10−3 kg.m−1 .s−1 et
cY = 505m.s−1 .
Pour la formulation SPH-Multifluide, l’équation d’état est celle de Tait (avec les coeffi-
cients polytropiques γX = 7 et γY = 1.4). Étant donné que les effets de tension superficielle
sont négligeables, ceux-ci ne sont pas modélisés par le schéma C.S.S. mais par la modéli-
sation alternative (cf. §II.2.2). Le coefficient I est pris égal à 0,08. Le noyau utilisé est la
gaussienne.
135
CHAPITRE III.4. BULLE ISOLÉE
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
10R xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
X
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Y
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2R
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
2R
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
6R
136
III.4.1. SUSSMAN
1
y/R
-1
1
y/R
-1
1
y/R
-1
137
CHAPITRE III.4. BULLE ISOLÉE
1.5
0.5
y/R
-0.5
-1
-1.5
-2
1.5
0.5
y/R
-0.5
-1
-1.5
-2
1.5
0.5
y/R
-0.5
-1
-1.5
-2
-2 -1 -2 -1 -2 -1
x/R x/R x/R
Fig. III.4.3 – Formulation SPH-VF en points bleus et carrés rouges pour la solution
Level-Set de Sussman et alii (1994).
138
III.4.1. SUSSMAN
Pour la formulation SPH-VF (cf. figure III.4.3), les résultats sont beaucoup moins en
accord avec ceux du Level-Set. Toutefois, de la même manière que sur le cas de validation
précédent (les instabilités de Rayleigh-Taylor), la solution ne semble pas être totalement
convergée à la plus forte résolution présentée ici. Les temps de calculs excessivement longs
de cette formulation nous ont amené à ne pas poursuivre l’étude de convergence en espace
au-delà de cette résolution.
Dans l’analyse des champs de vitesse et de pression (cf. figures III.4.4 & III.4.5) on
retrouve certaines tendances déjà mentionnées. Ainsi, le champ de vitesse dans le fluide
autour de la bulle en SPH-VF est bien plus faible qu’en SPH-Multifluide, conséquence
des déformations plus fortes avec cette dernière formulation. En revanche, le champ de
pression en SPH-VF apparaı̂t bien plus lisse spatialement qu’en SPH-Multifluide : ceci
est dû à l’emploi d’un solveur de Riemann. Il est à noter que des résultats semblables
pourraient être obtenus sur la formulation SPH-Multifluide en effectuant régulièrement
un lissage du champ de pression (par interpolation type MLS Belytschko et alii (1996)),
139
CHAPITRE III.4. BULLE ISOLÉE
cf. Colagrossi et Landrini (2003). Cela n’a pas été testé sur notre implémentation de la
formulation.
∑ mi u⃗i
i∈Y
u⃗G =
∑ mi
i∈Y
140
III.4.2. VITESSE TERMINALE ASCENDANTE
uc∞
0.2
h1
uyG
h2
0.1
h3
On note tout d’abord (cf. figure III.4.6) que, pour toutes les résolutions, la vitesse
converge en temps vers une certaine valeur, la vitesse terminale ascendante.
Nous avons cherché à déterminer la valeur de cette vitesse terminale à partir d’une
analyse théorique de l’écoulement. La simulation étant bidimensionnelle, la bulle plane
peut être assimilée à la section d’un cylindre infini dont on peut étudier l’effort de traı̂née
dans un écoulement à vitesse constante Lamb (1932) :
4πμX uc∞
Fc∞ =
1 Ruc∞ ρX
2 − γ − ln ( )
4μX
Pour cette configuration, la résolution numérique de cette équation donne pour solution :
uc∞ = 0, 22m.s−1 . Celle-ci est reportée sur la figure III.4.6.
Bien que la simulation numérique semble s’accorder avec cette vitesse terminale ana-
lytique, il faut noter plusieurs points essentiels sur la construction (et donc la validité) de
cette dernière. D’une part, le cylindre infini est supposé être un corps rigide donc l’écou-
lement se produisant à sa surface (condition d’adhérence sur un solide) est différent de la
condition existant à l’interface de la bulle simulée (condition d’adhérence mais entre deux
fluides). La circulation existant dans le gaz de la bulle va donc permettre un mouvement
de l’interface, ce qui va, à priori, avoir tendance à limiter les frottements et permettre à
la bulle de remonter plus vite que le cylindre.
Remarque :
141
CHAPITRE III.4. BULLE ISOLÉE
Pour être dans une configuration tout à fait équivalente au cylindre rigide, on aurait
pu simuler le mouvement libre d’un corps rigide.
D’autre part, le calcul effectué dans Lamb (1932) et présenté ici est valable pour de
faibles nombres de Reynolds (Re < 1), ce qui n’est pas tout à fait le cas dans la simulation
(Re = 5). Il serait nécessaire de prendre en compte l’asymétrie liée au décollement de
l’écoulement derrière le cylindre (cf. les corrélations données par Clift et alii (1978)).
Enfin, les effets de compressibilité du fluide gazeux sont importants (environ 20% de
variation du volume vair , cf. figure III.4.7(b)) : plus la bulle remonte, plus son volume
augmente (tout en conservant sa masse, acquise par construction du schéma). Alors la
différence entre la poussée d’Archimède et la force de pesanteur augmente (proportionnel-
lement à R2 ) tandis que les effets de traı̂née s’accentuent parallèlement. Néanmoins ces
efforts semblent s’équilibrer car la bulle n’accélère pas au cours du temps.
t = 1s 0.15
t = 3s
vair
0.1
n
y
t = 7s
0.05
0
-0.02 -0.01 0 0.01 0.02 0 2 4 6
x t
(a) Déformation de la surface libre à différents (b) Evolution du volume normalisé
0
n vair −vair
instants. Échelles identiques dans les deux direc- vair = 0
vair
de la bulle au cours du
tions. temps.
Nous pouvons observer sur la figure III.4.7(a) les déformations de la surface libre à
divers instants de remontée de la bulle. Celle-ci reste très circulaire et démontre l’efficacité
du schéma C.S.S. à reproduire les forts effets de tension superficielle (Bo = 0, 4).
Conclusion
Les écoulements qui ont été étudiés dans ce chapitre sont très importants dans la
mesure où leur validation permet de connaı̂tre l’aptitude des formulations à prédire des
situations qui se reproduiront dans un séparateur.
Toutefois, les cas présentés ne suffisent pas à balayer l’ensemble des configurations
d’écoulements de bulles pouvant réellement intervenir dans un séparateur : il serait néces-
saire de faire une étude plus exhaustive en faisant varier les caractéristiques de ces écoule-
ments de bulles isolées (ceci s’étudie classiquement dans le domaine We-Re). D’autre part,
il est à signaler que ce genre d’étude ne peut, à priori, être faite de manière satisfaisante
qu’en procédant à des simulations tridimensionnelles des écoulements. Cela permet alors
142
III.4.2. VITESSE TERMINALE ASCENDANTE
143
CHAPITRE III.4. BULLE ISOLÉE
144
Chapitre III.5
Séparateur eau-huile
Le but ultime de ces travaux est l’étude de la séparation de fluides, et plus particu-
lièrement de l’eau et de l’huile. Pour cela nous avons développé et confronté à des cas
de validation deux formulations dont l’une semble montrer plus de capacités pour pou-
voir simuler ce problème. Afin de tirer parti au mieux des caractéristiques favorables de
cette formulation SPH-Multifluide, nous avons construit un séparateur simplifié dont la
configuration est détaillée ci-dessous.
HY Y
HX
X
145
CHAPITRE III.5. SÉPARATEUR EAU-HUILE
d’intérêt d’un point de vue physique, nous avons, à la fois déplacé aléatoirement (selon
les 2 directions principales) le centre de chaque bulle autour de ces positions de référence,
et modulé le rayon de chacune d’entre elles de ±50% par rapport au rayon nominal.
Le domaine a été ensuite construit de manière à ce que des interactions puissent se
produire entre les bulles et les parois (avec L = HX = 26, 6R). Au dessus de la colonne d’eau
au repos contenant les bulles d’huile a été rajoutée une colonne d’huile (avec HX = 8, 4R)
aussi au repos, modélisant une part d’huile déjà séparée vers laquelle les bulles vont migrer.
Les fluides ont les caractéristiques suivantes : pour l’eau on a une densité égale à
ρX = 1000kg.m−3 , une vitesse du son de 7m.s−1 . Sa viscosité est choisie telle que α = 0, 08
soit μX = 7.10−3 kg.m−1 .s−1 . Pour l’huile on a choisi une densité égale à ρY = 800kg.m−3 , la
vitesse du son égale à 7, 8m.s−1 et une viscosité de μY = 6, 5.10−3 kg.m−1 .s−1 .
Le noyau choisi est la gaussienne avec une longueur de lissage égale à h = 10−4 m.
Le nombre de Reynolds de cet écoulement (cf. III.4.1.1) est de 40. Concernant les effets
de tension superficielle, nous ne les avons pas modélisés dans ce premier essai effectué
avant la validation du schéma C.S.S., seule la correction avec I = 0, 03 est appliquée.
Toutefois, en première approximation, on peut calculer le nombre de Bond moyen avec
un coefficient de tension superficielle eau-huile (σeau−huile = 2.10−2 N.m−1 ). On a alors
Boeau−huile = 4Δρ gR2 /σ ∼ 0, 4. En réalité, cet écoulement présente donc des effets de
tension superficielle importants. Les résultats obtenus ici seront donc à interpréter avec
précaution.
III.5.1.2 Résultats
Les résultats obtenus sont présentés sur la figure III.5.2. Qualitativement, on peut
mentionner le comportement cohérent de l’ensemble des bulles par rapport au compor-
tement observé avec une bulle isolée sans tension superficielle : d’importants « aplatisse-
ments » des bulles se produisent. Avec le nombre de Bond estimé pour l’écoulement réel,
les bulles devraient se déformer très peu tant qu’elles ne se touchent pas, surtout les plus
petites (cf. le cas de la bulle isolée dont on cherche à déterminer la vitesse terminale ascen-
dante). Ce n’est pas le cas ici, ce qui justifierait l’emploi du schéma C.S.S pour modéliser
des effets de tension superficielle.
Néanmoins, ces résultats sont encourageants dans la mesure où ils démontrent la
capacité de la formulation SPH-Multifluide à prendre en compte de manière relativement
précise un nombre élevé d’interfaces disjointes, pouvant se fragmenter et/ou se reconnecter,
comme ce serait le cas pour des nombre de Bond plus élevés. L’intérêt en est d’autant plus
grand dans la zone de fort mélange à l’interface des deux colonnes de fluide initialement
au repos, là où se déroule le processus même de séparation, où l’on peut observer à quel
point la non-miscibilité des fluides est bien respectée. Une telle simulation effectuée avec
une technique autorisant la diffusion de l’interface mènerait inévitablement à l’obtention
d’une zone de mélange épaisse.
Pour quelques instants choisis, les champs de vitesse et pression sont aussi présentés
sur la figure III.5.3. On remarque que le champ de pression est assez uniforme et dominé
par la pesanteur, l’écoulement étant assez lent. En revanche, on notera la complexité du
champ de vitesse dans le domaine.
Remarque :
146
III.5.1. CUVE FERMÉE
Fig. III.5.2 – Ensemble de bulles dans un séparateur simplifié. Champ de densité à divers
instants.
Le temps de calcul sur un processeur 2, 4GHz est de 18 jours pour 1, 32s de temps physique
simulé (171000 particules) avec un pas de temps de 1, 2 ∼ 10−5 s, soit 110000 itérations du
147
CHAPITRE III.5. SÉPARATEUR EAU-HUILE
schéma d’avance en temps Runge-Kutta 4e ordre. Des techniques spécifiques devront donc
être recherchées pour améliorer ces temps importants. On pourra notamment tenter de
s’affranchir de la condition CF Lvisq (cf. remarque au chapitre concernant la viscosité).
III.5.2.1 Configuration
Le système d’entrée/sortie développé par Oger et alii (2006) est adapté ici pour les
écoulements bifluides. Les modifications portent sur deux points : la zone d’entrée et
l’imposition des variables physiques.
Tout d’abord, l’objectif du système d’entrée est de permettre l’injection de deux fluides
mélangés sous forme de bulles. Comme précédemment, de façon à éviter des configura-
tions singulières, les bulles injectés vont être choisies avec un rayon et une position relative
dans la hauteur de la section d’entrée de manière aléatoire (en restant dans un certain
148
III.5.2. CUVE AVEC ENTRÉE/SORTIE
intervalle). Ainsi, pour faciliter la gestion de la position et du rayon des bulles par rap-
port au support d’interpolation des particules étant dans le domaine de calcul (zone entre
P osE − 3h et P osE sur la figure III.5.4), celles-ci sont créées dans une zone de « ger-
mination » (zone tampon en amont de la section d’entrée, entre P osT et P osE − 3h) et
advectées suivant la vitesse imposée jusqu’à l’entrée (il n’y a donc pas d’interpolation SPH
à proprement parler dans la zone tampon).
xxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
u⃗
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxx
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
149
CHAPITRE III.5. SÉPARATEUR EAU-HUILE
Le débit volumique imposé à l’entrée est égal à 2.10−4 m2 .s−1 (calcul 2D) pour l’eau, et
égal à 5 fois moins pour l’huile. Ainsi la vitesse du mélange à l’entrée est VsM = 24mm.s−1 .
Le noyau choisi est la gaussienne avec une longueur de lissage égale à h = 5.10−4 m (ce
qui correspond à environ 98 000 points de discrétisation).
Le nombre de Reynolds de cet écoulement (basé sur VsM et 2R) est de ∼ 1, 6. Le
nombre de Bond est égal à 0,2. Il est à noter que ce nombre adimensionnel est calculé
pour une bulle remontant dans un milieu infini au repos, ce qui n’est pas exactement le
cas ici : cependant,
√ les effets inertiels sont peu importants devant les effets de la gravité
(F r = VsM / g2R = 0, 17). Ainsi les effets visqueux et capillaires seront prépondérants par
rapports aux effets inertiels.
Y VsH
80r
X
VeM VsE
160r
III.5.2.2 Résultats
Les résultats sont présentés sur la figure III.5.6. Un détail du domaine de calcul est
aussi présenté sur cette même figure pour illustrer le fonctionnement du système d’entrée
sur ce cas-là. Qualitativement, on observe un écoulement qui est dominé par les effets de
gravité et qui se traduit par une remontée rapide des bulles après l’entrée. L’écoulement
des bulles est donc peu dominé par l’effet inertiel dû à la vitesse d’entrée.
Pareillement à la simulation dans la cuve fermé, le nombre de Bond de cet écoulement
est faible. Or ici le schéma C.S.S. a été utilisé. Il en résulte que les bulles conservent une
forme circulaire beaucoup plus longtemps au cours de leur remontée. Des irrégularités
à l’interface peuvent être cependant observées (cf. détail fig. III.5.6) : celle-ci peuvent
provenir d’une résolution spatiale insuffisante (ici on a R/h = 2 contre R/h = 10 dans
la précédente simulation). La discrétisation plus fine du domaine de calcul, qui apparaı̂t
ici nécessaire, se heurte aux dimensions même de ce grand domaine : une stratégie de
parallélisation algorithmique serait alors efficace pour résoudre ce genre de problème en
des temps de calcul raisonnables.
Néanmoins, de la même façon que précédemment, on peut souligner encore que la
formulation SPH-Multifluide prend en compte ces interfaces complexes sans diffusion nu-
mérique.
150
III.5.2. CUVE AVEC ENTRÉE/SORTIE
Fig. III.5.6 – Séparateur simplifié avec système d’entrée/sortie. Colonne de droite : zoom
sur l’entrée. Champ de densité aux instants t = 1s; 2s; 3s; 4s en partant du haut.
151
CHAPITRE III.5. SÉPARATEUR EAU-HUILE
Conclusion
Ce cas d’étude d’un séparateur simplifié conclut donc nos travaux en montrant les
capacités de la formulation SPH-Multifluide à simuler des problèmes de ce type. Comme
mentionné au chapitre précédent, il reste toutefois à approfondir les différents régimes
d’écoulements de bulles isolées et à s’intéresser aux phénomènes spécifiques lorsque plu-
sieurs bulles interagissent (coalescence, rupture) avant de pouvoir réellement valider le
comportement d’un ensemble plus important de bulles.
152
Conclusion
Les travaux présentés dans ce mémoire ont porté sur le développement et la validation
de formulations de la méthode SPH adaptées à la modélisation d’écoulements bifluides.
Ainsi, avec l’objectif final de modélisation du fonctionnement d’un séparateur eau-huile,
nous avons poursuivi le développement initié dans l’équipe par G. Oger dans une étude
préliminaire de la simulation d’écoulements bifluides avec la méthode SPH.
Après avoir identifié les faiblesses de cette formulation initiale, il a été proposé, en
collaboration avec A. Colagrossi, une nouvelle formulation plus à même de tenir compte
des problèmes numériques introduits par la modélisation de plusieurs fluides avec interface.
Parallèlement ont été aussi étendus les travaux effectués dans l’équipe sur les solveurs de
Riemann pour la méthode SPH, par J.-B. Deuff et P.-M. Guilcher, aux cas d’écoulements
bifluides avec cette technique. Ces deux formulations, respectivement dénommées dans le
manuscrit comme SPH-Multifluide pour la première et SPH-VF pour la seconde, ont alors
été développées de manière distincte.
Puis de nouvelles modélisations physiques ont été ajoutées aux formulations afin de
prendre en compte les effets de viscosité et les effets de tension superficielle. Ces schémas
additifs ont été implémentés pour chaque formulation (sauf la tension superficielle avec
SPH-VF) et testés sur des cas de validation spécifiques (écoulement de Poiseuille, oscilla-
tion d’une goutte en apesanteur). Si l’une des formulations (SPH-Multifluide) a montré de
bonnes capacités sur ces tests, l’autre a révélé un comportement beaucoup plus diffusif.
Enfin, d’autres cas de validation plus complets ont été simulés (ballottement linéaire,
instabilités de Rayleigh-Taylor, cas d’envahissement, écoulement de bulle isolée), s’atta-
chant à chaque fois à comparer les deux formulations aux solutions de référence. Globa-
lement, il a été retrouvé un meilleur comportement de la méthode SPH-Multifluide que
celle SPH-VF. Il faut néanmoins préciser que cette dernière formulation se base sur un
formalisme plus complexe. L’analyse des sources de différence aux solutions de référence
ainsi que les solution palliatives sont donc plus difficiles à mettre en œuvre. Néanmoins, on
a pu identifier le point dégradant excessivement sa précision : la correction ajoutée pour
permettre à la formulation de simuler de manière satisfaisante les écoulements bifluides
est clairement à l’origine de ce défaut. Aucune alternative efficace n’a été trouvée dans
cette étude pour s’en affranchir, mais la question reste ouverte.
La poursuite des cas-test avec la SPH-Multifluide a mené à démontrer finalement sa
capacité à simuler le problème cible de ce travail, la séparation eau-huile. On bénéfice là
de ce que cette formulation capture les formes complexes d’interface, incluant de multiples
fragmentations et reconnexions sans aucune diffusion numérique de l’interface.
Au moment de dresser le bilan de ces travaux, on peut donc dire que l’objectif initial
de développer une modélisation SPH bifluide apte à être utilisée pour réaliser la sépara-
153
CONCLUSION
tion eau-huile est globalement atteint. Toutefois le travail réalisé est loin d’être exhaustif
sur le sujet et bien des améliorations seront à y apporter afin d’augmenter le potentiel des
formulations proposées. Nous allons suggérer quelques pistes de travail.
Tout d’abord sur la modélisation des effets physiques introduits, ceux-ci se sont limités
à de faibles nombres de Reynolds pour des questions de stabilité. Il conviendra donc
d’étendre le champ d’application de la modélisation de la viscosité : ce travail est poursuivi
dans la thèse de M. De Leffe en cours. Pour les effets de tension superficielle, il sera
nécessaire de travailler sur les courants parasites et en adaptant les techniques existant
dans la littérature sur les méthodes maillées. Ayant observé qu’actuellement l’intensité
des courants ne convergeait pas en espace vers zéro, résoudre ce problème permettra de
conduire correctement des simulations à haute résolution.
La formulation SPH-Multifluide étant relativement robuste et ayant montré de bonnes
capacités, peu de corrections (hormis celles décrites ci-dessus) semblent à recommander
aujourd’hui dans le contexte de la séparation eau-huile. Bien sûr, elle pourra bénéficier des
améliorations qui seraient réalisées sur la formulation monofluide, en terme de précision
par exemple. Néanmoins, l’étude de cas où la dynamique est plus violente (impacts)
nécessitera d’introduire un lissage supplémentaire des champs (tel que pratiqué de son
coté par A. Colagrossi sur la divergence de la vitesse). L’étude fine du comportement
de la méthode en termes de pression dans l’eau et dans l’air sur des cas confinés type
envahissement devra aussi être menée.
Quant à la formulation SPH-VF, son potentiel réside dans la résolution précise de sa
partie hyperbolique (qui, au contraire, ne nécessite pas de lissage tel que celui proposé
ci-dessus). Néanmoins, la clé de son succès résidera dans la compréhension du problème
numérique justifiant l’introduction de la correction proposée. De plus, les temps de calculs
présentés dans les différents cas, qui peuvent paraı̂tre prohibitifs, sont à relativiser : des
gains sont à attendre en étudiant mieux l’influence du coefficient de la condition CFL (qui
pourrait être augmenté comme en formulation monofluide) et celle du nombre de voisins
(la précision en compensant alors leur plus faible nombre).
Enfin, concernant l’applicabilité de ce qui a été développé, l’inclusion relativement
aisée dans le code de calcul SPH-flow du laboratoire (tridimensionnel et parallélisé) per-
mettra d’étendre les possibilités d’étude et de validation. En ce qui concerne la séparation
eau-huile, on pourra, par exemple, valider le modèle sur des études d’écoulements de bulles
réalistes (tridimensionnels). Et, outre une exploration du domaine {Reynolds-Weber} et
son influence sur le comportement d’une bulle isolée, cela permettra de s’intéresser aussi
aux phénomènes de coalescence entre bulles ou de ruptures de celles-ci. L’étude d’un sé-
parateur sera alors beaucoup plus proche de la physique modélisée (avec possibilité de
pouvoir prendre en compte une géométrie réaliste de cet équipement).
154
Annexe A
On reprend ici la démonstration de Lamb (1932), qui détaille les modes propres d’os-
cillation d’une goutte bi-dimensionnelle en apesanteur et dans le vide, pour la compléter
du cas où un fluide est pris en compte dans l’environnement extérieur. Cette extension
est essentiellement basée sur un autre travail effectué par Lamb sur une goutte sphérique
dans un liquide environnant.
On suppose que la goutte plane, qui peut être assimilée à une section d’un cylindre
liquide infini (tel que l’a supposé Lamb), est de masse volumique ρi et qu’elle est entourée
d’un fluide de masse volumique ρe . La tension superficielle existant entre eux a pour
coefficient σ.
La goutte a un rayon moyen R. Les oscillations créées par l’excitation initiale (quel-
conque) vont déformer l’interface et on suppose que ces déformations ξ sont suffisamment
petites devant R. Le rayon local r peut donc se décomposer de la manière suivante :
r(θ) = R + ξ(θ).
Comme détaillé par Landau et Lifchitz (1989), on peut exprimer le rayon de courbure
κ en fonction de ces déformation :
1 1 ∂ 2ξ
κ = − 2 (ξ + 2 ) .
R R ∂θ
On introduit alors le potentiel des vitesse φ pour chaque fluide,
– pour celui intérieur :
A rn
φi = cos(nθ) cos(ωt + ϕ);
n Rn
– et pour celui extérieur :
A rn+1
φe = − cos(nθ) cos(ωt + ϕ)
n + 1 Rn+1
de telle sorte qu’ils puissent vérifier, sur l’interface r = R, la condition suivante :
∂ξ ∂φi ∂φe
= − = − .
∂t ∂r ∂r
155
ANNEXE A. PULSATION PROPRE D’UNE GOUTTE CYLINDRIQUE
pi − pe = σκ.
En introduisant les différents termes définis plus haut et en omettant les constantes (qui
ne contribuent pas à la détermination des modes propres), on a :
A A σ A
−ρi cos(nθ)ω sin(ωt+ϕ)−ρe cos(nθ)ω sin(ωt+ϕ) = − 2 (−1 + n2 ) cos(nθ) sin(ωt+ϕ),
n n R nR
relation valable à chaque instant t et pour tout angle θ. On en déduit donc la relation
suivante :
σ n(n + 1)(n2 − 1)
ω2 = 3 .
R (n + 1)ρi + nρe
Pour deux cas remarquables, cette relation se simplifie de la manière suivante :
– si le fluide extérieur est de densité négligeable (ρe ≪ ρi ), on retrouve la relation
démontrée par Lamb :
σ n(n2 − 1)
ω2 = 3 ;
R ρi
– si les deux fluides ont une densité identique (ρi = ρe ) :
σ n(n + 1)(n2 − 1)
ω2 = .
R3 (2n + 1)ρi
156
Bibliographie
Andrillon Yann. Simulation d’écoulements à surface libre par une méthode de capture
d’interface en formulation totalement couplée. Thèse de Doctorat, École Centrale de
Nantes, Mar 2004.
Basa Mihai, Quinlan Nathan J., et Lastiwka Martin. 2008. Robustness and accuracy of
sph formulations for viscous flow. Int. J. Numer. Meth. Fluids. Disponible à http://
dx.doi.org/10.1002/fld.1927.
Belytschko T., Krongauz Y., Dolbow J., et Gerlach C. 1998. On the comple-
teness of meshfree particle methods. International Journal of Numerical Me-
thods in Engineering, 43 :785–819. Disponible à http://dx.doi.org/10.1002/
(SICI)1097-0207(19981115)43:5%3C785::AID-NME420%3E3.0.CO;2-9.
Belytschko T., Krongauz Y., Organ D., Fleming M., et Krysl P. December 1996. Meshless
methods : An overview and recent developments. Computer Methods in Applied Me-
chanics and Engineering, 139(1-4) :3–47. Disponible à http://dx.doi.org/10.1016/
S0045-7825(96)01078-X.
Blondel Élise. Conception d’un protocole expérimental pour validation de simulation sph
d’envahissement. Master’s thesis, Université Pierre et Marie Curie, 2006.
Brackbill J. U., Kothe D. B., et Zemach C. June 1992. A continuum method for modeling
surface tension. Journal of Computational Physics, 100(2) :335–354. Disponible à
http://dx.doi.org/10.1016/0021-9991(92)90240-Y.
157
BIBLIOGRAPHIE
Clift R., Grace J. R., et E. Weber M. Bubbles, drops, and particles. 1978. Academic
Press, New-York. ISBN 0-12-176950-X.
Colagrossi A., Colicchio G., et Le Touzé D. Enforcing boundary conditions in sph appli-
cations involving bodies with right angles. In Proc. 2nd SPHERIC Workshop, 2007.
Colagrossi A., Le Touzé D., et Antuono M. 2008. Theoretical considerations on the free-
surface role in the smoothed-particle-hydrodynamics model. Physical Review E, 79,
Issue 5, 056701. Disponible à http://dx.doi.org/10.1103/PhysRevE.79.056701.
Colagrossi Andrea. A meshless lagrangian method for free-surface flows and interface
flows with fragmentation. PhD thesis, Università di Roma La Sapienza, Apr 2005.
Colicchio G. Violent disturbance and fragmentation of free surfaces. PhD thesis, University
of Southampton, 2004.
Cueille Pierre-Victor. Modélisation par SPH des phénomènes de diffusion présent dans
un écoulement fluide. Thèse de Doctorat, INSA Toulouse, Dec 2005.
Deuff Jean-Baptiste. Extrapolation au réel des mesures de pressions obtenus sur des cuves
modèle réduit. Thèse de Doctorat, École Centrale de Nantes, Oct 2007.
Di G. Sigalotti L., Klapp J., Sira E., Meleán Y., et Hasmy A. 1997. Sph simulations
of time-dependent poiseuille flow at low reynolds numbers. Journal of Computational
Physics, 191(2). Disponible à http://dx.doi.org/10.1016/S0021-9991(03)00343-7.
Doring Mathieu. Développement d’une méthode SPH pour les applications à surface libre
en hydrodynamique. Thèse de Doctorat, École Centrale de Nantes, Jun 2005.
Francois Marianne M., Cummins Sharen J., Dendy Edward D., Kothe Douglas B., Sicilian
James M., et Williams Matthew W. 2006. A balanced-force algorithm for continuous
and sharp interfacial surface tension models within a volume tracking framework. Jour-
nal of Computational Physics, 213(1) :141 – 173. ISSN 0021-9991. Disponible à http://
dx.doi.org/10.1016/j.jcp.2005.08.004.
Godunov S.K. 1958. A difference scheme for numerical solution of discontinuous solution
of hydrodynamic equations. Math. Sbornik, 47, 271-306. Translated by US Joint Publ.
Res. Service, JPRS 7226, 1969.
158
BIBLIOGRAPHIE
Grenier N., Antuono M., Colagrossi A., Le Touzé D., et Alessandrini B. 2009. An ha-
miltonian interface sph formulation for multi-fluid and free-surface flows. Journal of
Computational Physics. Disponible à http://dx.doi.org/10.1016/j.jcp.2009.08.
009.
Guilcher Pierre-Michel. Contribution au développement d’une méthode SPH pour la simu-
lation numérique des interactions houle-structure. Thèse de Doctorat, École Centrale
de Nantes, Oct 2008.
Harten Amiram. 1983. High resolution schemes for hyperbolic conservation laws. Journal
of Computational Physics, 49 :357–393. Disponible à http://dx.doi.org/10.1006/
jcph.1997.5713.
Harten Amiram et Hyman James M. 1983. Self adjusting grid methods for one-dimensional
hyperbolic conservation laws. Journal of Computational Physics, 50(2) :235–269. Dis-
ponible à http://dx.doi.org/10.1016/0021-9991(83)90066-9.
Harten Amiram, Lax Peter D., et Van Leer Bram. 1983. On upstream differencing and
godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25(1) :35–61.
Disponible à http://dx.doi.org/10.1137/1025002.
Hirt C. W. et Nichols B. D. 1981. Volume of fluid (vof) method for the dynamics of free
boundaries. Journal of Computational Physics, 39. Disponible à http://dx.doi.org/
10.1016/0021-9991(81)90145-5.
Hu X. Y. et Adams N. A. 2006. A multi-phase sph method for macroscopic and mesoscopic
flows. Journal of Computational Physics, 213 (2) :844–861. Disponible à http://dx.
doi.org/10.1016/j.jcp.2005.09.001.
Ivings M.J., Causon D.M., et Toro E.F. 1998. On riemann solvers for compres-
sible liquids. International Journal for Numerical Methods in Fluids, 28 :395–418.
Disponible à http://dx.doi.org/10.1002/(SICI)1097-0363(19980915)28:3<395::
AID-FLD718>3.0.CO;2-S.
Lafaurie B., Nardone C., Scardovelli R., Zaleski S., et Zanetti G. 1994. Modelling merging
and fragmentation in multiphase flows with surfer. Journal of Computational Physics,
113, Issue 1 :134–147. Disponible à http://dx.doi.org/10.1006/jcph.1994.1123.
Lamb Horace. Hydrodynamics. 1932. Dover, New York.
Landau L. et Lifchitz E. Physique théorique. Tome 6 : Mécanique des Fluides. 1989. Mir,
Moscou.
Lanson N. et Vila J. P. 2001. Meshless methods for conservation laws. Mathematics and
Computers in Simulation, 55 :493–501. Disponible à http://dx.doi.org/10.1016/
S0378-4754(00)00285-8.
Lanson Nathalie. Etude des méthodes particulaires renormalisées : applications aux pro-
blèmes de dynamique rapide. Thèse de Doctorat, INSA Toulouse, 2001.
Le Touzé D., Colagrossi A., et Colicchio G. Ghost technique for right angles applied to
the solution of benchmarks 1 and 2. In Proc. 1st SPHERIC Workshop, 2006.
159
BIBLIOGRAPHIE
Monaghan J. J. 1994. Simulating free surface flows with SPH. Journal of Computational
Physics, 110 :399–406. Disponible à http://dx.doi.org/10.1006/jcph.1994.1034.
Morris Joseph P. 2000. Simulating surface tension with smoothed particle hy-
drodynamics. International Journal of Numerical Methods in Fluids, 33 :333–
353. Disponible à http://dx.doi.org/10.1002/1097-0363(20000615)33:3%3C333::
AID-FLD11%3E3.0.CO;2-7.
Morris Joseph P., Fox Patrick J., et Zhu Yi. 1997. Modeling low reynolds number incom-
pressible flows using sph. Journal of Computational Physics, 136 :214–226. Disponible
à http://dx.doi.org/10.1006/jcph.1997.5776.
Nugent S. et Posch H. A. 2000. Liquid drops and surface tension with smoothed particle
applied mechanics. Physics Review E, 62. Disponible à http://dx.doi.org/10.1103/
PhysRevE.62.4968.
Oger G., Doring M., Alessandrini B., et Ferrant P. 2006. Two-dimensional sph simulations
of wedge water entries. Journal of Computational Physics, 213 (2) :803–822.
Prosperetti Andrea. 1980. Normal-mode analysis ofr the oscillations of a viscous liquid
drop in a immiscible liquid. Journal de Mécanique, 19(1) :149–182.
160
Randles P. W. et Libersky L. D. 1996. Smoothed particle hydrodynamics : some recent im-
provements and applications. Computer Methods in applied mechanics and engineering,
139 :375–408. Disponible à http://dx.doi.org/10.1016/S0045-7825(96)01090-0.
Raviart P.A. An analysis of particle methods. 1985, volume 1127 of Lecture Notes in Ma-
thematics. Springer, Berlin. Disponible à http://dx.doi.org/10.1007/BFb0074532.
Sussman M., Smereka P., et Osher S. 1994. A level set approach for computing solutions
to incompressible two-phase flow. Journal of Computational Physics, 114, 146-159.
Disponible à http://dx.doi.org/10.1006/jcph.1994.1155.
Takeda H., Miyama S. M., et Sekiya M. 1994. Numerical simulation of viscous flow
by smoothed particle hydrodynamics. Progress of Theoretical Physics, 92(5) :939–960.
Disponible à http://dx.doi.org/10.1143/PTP.92.939.
Toro E. F. Riemann solvers and numerical methods for fluid dynamics. 1997. Springer,
Berlin.
Van Leer B. 1979. Towards the ultimate conservative difference scheme. v. a second-
order sequel to godunovŠs method. Journal of Computational Physics, 32 :101Ű136.
Disponible à http://dx.doi.org/10.1016/0021-9991(79)90145-1.
Watkins S. J., Bhattal A. S., Francis N., Turner J. A., et Whitworth A. P. 1996. A
new prescription for viscosity in smoothed particle hydrodynamics. Astronomy and
Astrophysics Supplement Series, 119 :177–187. Disponible à http://dx.doi.org/10.
1051/aas:1996104.
Zhou Z.Q., De Kat J.O., et Buchner B. A nonlinear 3-d approach to simulate green water
dynamics on deck. In Proc. 7th Int. Conf. Num. Ship Hydrod., pages 5.1–1, 15, 1999.
BIBLIOGRAPHIE
Keywords : SPH, multiphase, interface, volume fraction , Riemann solver, bubble flows
162