Méthodes D'éléments Finis Pour Les Interactions Fluide-Structure
Méthodes D'éléments Finis Pour Les Interactions Fluide-Structure
Méthodes D'éléments Finis Pour Les Interactions Fluide-Structure
fluide-structure
Thèse
Aymen Jendoubi
Québec, Canada
Cette thèse concerne la modélisation des interactions fluide-structure et les méthodes numé-
riques qui s’y rattachent. De ce fait, la thèse est divisée en deux parties. La première partie
concerne l’étude des interactions fluide-structure par la méthode des domaines fictifs. Dans
cette contribution, le fluide est incompressible et laminaire et la structure est considérée rigide,
qu’elle soit immobile ou en mouvement. Les outils que nous avons développés comportent la
mise en oeuvre d’un algorithme fiable de résolution qui intégrera les deux domaines (fluide et
solide) dans une formulation mixte. L’algorithme est basé sur des techniques de raffinement
local adaptatif des maillages utilisés permettant de mieux séparer les éléments du milieu fluide
de ceux du solide que ce soit en 2D ou en 3D. La seconde partie est l’étude des interactions
mécaniques entre une structure flexible et un fluide incompressible. Dans cette contribution,
nous proposons et analysons des méthodes numériques partitionnées pour la simulation de
phénomènes d’interaction fluide-structure (IFS). Nous avons adopté à cet effet, la méthode
dite «arbitrary Lagrangian-Eulerian» (ALE). La résolution fluide est effectuée itérativement
à l’aide d’un schéma de type projection et la structure est modélisée par des modèles hyper
élastiques en grandes déformations. Nous avons développé de nouvelles méthodes de mouve-
ment de maillages pour aboutir à de grandes déformations de la structure. Enfin, une stratégie
de complexification du problème d’IFS a été définie. La modélisation de la turbulence et des
écoulements à surfaces libres ont été introduites et couplées à la résolution des équations de
Navier-Stokes. Différentes simulations numériques sont présentées pour illustrer l’efficacité et
la robustesse de l’algorithme. Les résultats numériques présentés attestent de la validité et
l’efficacité des méthodes numériques développées.
iii
Abstract
This thesis is concerned with the modeling of fluid-structure interactions (FSI) and the cor-
responding specific numerical methods. The thesis is divided into two principal parts. The
first part concerns the study of fluid-structure interactions using the fictitious domain method.
In this contribution, the fluid is incompressible and laminar and the structure is considered
rigid, whether stationary or moving. The tools we have developed include the implementa-
tion of a reliable resolution algorithm that incorporates both domains (fluid and solid) in a
common mixed formulation. The algorithm is based on adaptive local mesh refinement tech-
niques used to distinguish the elements in the fluid from those of the solid either in 2D or
3D. The second part is the study of the mechanical interactions between a flexible structure
and an incompressible fluid. In this context, we propose and analyze partitioned numerical
methods for simulating fluid-structure interaction phenomena (FSI). We adopt an "arbitrary
Lagrangian-Eulerian" (ALE) formulation for this purpose. The fluid resolution is performed
iteratively by means of a projection scheme and the structure is modeled by hyperelastic mod-
els in large deformations. We have introduced new mesh update methods to achieve large
deformation of the structure. Finally, a more complex strategy for FSI problem is proposed.
The turbulence and two-phase flows modelling are introduced and coupled to the resolution
of the Navier-Stokes equations for studying FSI problems. The numerical results presented
attest the validity and effiency of the proposed numerical methods developed.
v
Table des matières
Résumé iii
Abstract v
Remerciements xiii
Introduction 1
vii
5 Couplage multiphysique 123
5.1 Modélisation de la turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . 123
5.2 Modélisation des écoulements diphasiques . . . . . . . . . . . . . . . . . . . 133
5.3 Algorithme d’IFS turbulent . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.4 Algorithme d’IFS turbulent à surface libre . . . . . . . . . . . . . . . . . . . 144
5.5 Interaction entre un fluide et une structure en rotation avec surface libre . . 147
Conclusion 157
Bibliographie 161
viii
Liste des tableaux
3.1 Cas #1 : Qualité moyenne des éléments du maillage. Hext , Eext , Kext et Pext
présente respectivement les relèvements harmonique, élastique, krigeage et pa-
rabolique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.2 Cas #1 : Distribution de la qualité des éléments pour d = 0.18, 0.46 et 1.00. Le
nombre d’éléments retournés est indiqué entre parenthèses. . . . . . . . . . . . 80
3.3 Cas #2 : Qualité moyenne des éléments du maillage. . . . . . . . . . . . . . . . 82
3.4 Cas #2 : Distribution de la qualité d’éléments pour θ = 9°, 30°, 60°et 75°(d =
0.12, 0.40, 0.80 et 1.0 resp.). Le nombre d’éléments retournés est indiqué entre
parenthèses. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
3.5 Cas #3 : Qualité moyenne des éléments du maillage. . . . . . . . . . . . . . . . 86
3.6 Cas #3 : Distribution de la qualité des éléments pour d = 0.2, 0.4 et 1.0 . Le
nombre d’éléments retournés est indiqué entre parenthèses. . . . . . . . . . . . 86
3.7 Cas #4 : Qualité moyenne des éléments du maillage. . . . . . . . . . . . . . . . 88
3.8 Cas #4 : Distribution de la qualité des éléments pour d = 0.125, 0.50 et 1.0. Le
nombre d’éléments retournés est indiqué entre parenthèses. . . . . . . . . . . . 89
3.9 Cas #4 : Temps de calcul pour une itération. Tous les calculs ont été exécutés
sur un processeur Intel 64 bits de 3GHz typique . . . . . . . . . . . . . . . . . . 89
ix
Liste des figures
xi
3.5 Cas #1 : Maillages obtenus pour les différentes méthodes (d = 1.0) : à gauche
une vue rapprochée, à droite un agrandissement de l’extrémité de la barre . . . 81
3.6 Cas #2 : Géométrie du problème. . . . . . . . . . . . . . . . . . . . . . . . . . . 82
3.7 Cas #2 : Maillage initial et vue détaillée du profil d’aile. . . . . . . . . . . . . 82
3.8 Cas #2 : Maillages obtenus pour les différentes méthodes (θ = 60°). Au milieu
et à droite, une vue rapprochée des extrémités du profil. . . . . . . . . . . . . . 84
3.9 Mouvement rigide d’une sphère : géométrie considérée . . . . . . . . . . . . . . 85
3.10 Cas #4 : Géometrie et conditions aux limites. . . . . . . . . . . . . . . . . . . . 87
3.11 Cas #4 : Différents angles de vision du déplacement imposé à la structure pour
d = 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
3.12 Géométrie du problème de Turek et Hron (2006) . . . . . . . . . . . . . . . . . 94
3.13 Maillage du fluide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.14 Maillage de la structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
3.15 Champ de vitesses dans le cas test FSI2 . . . . . . . . . . . . . . . . . . . . . . 96
3.16 Champ de vitesses dans le cas test FSI3 . . . . . . . . . . . . . . . . . . . . . . 96
3.17 FSI 2 : déplacements en X et en Y du point A. . . . . . . . . . . . . . . . . . . 97
3.18 FSI 3 : déplacements en X et en Y du point A. . . . . . . . . . . . . . . . . . . 98
3.19 Représentation schématique du cas test des ondes de pression . . . . . . . . . . 98
3.20 Maillages fluide et structure non conformes . . . . . . . . . . . . . . . . . . . . 99
3.21 Propagation de l’onde de pression . . . . . . . . . . . . . . . . . . . . . . . . . . 100
3.22 Géométrie et conditions aux limites . . . . . . . . . . . . . . . . . . . . . . . . . 101
3.23 Maillages : domaine fluide (gauche) et domaine solide (droite) . . . . . . . . . . 102
3.24 Drapeau déformé avec lignes de courant du fluide pour 3 pas de temps différents 103
3.25 Déplacement horizontal du point A . . . . . . . . . . . . . . . . . . . . . . . . . 104
3.26 Problème du pain de gomme : géométrie et conditions aux limites . . . . . . . . 105
3.27 Géométrie et maillage de peau du pain gomme . . . . . . . . . . . . . . . . . . 105
3.28 Géométrie et maillage de peau du domaine fluide . . . . . . . . . . . . . . . . . 106
3.29 Déformée du pain de gomme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
3.30 Champ de vitesse obtenu dans un plan perpendiculaire au pain de gomme . . . 107
3.31 Vecteurs vitesse présentés dans une coupe transversale du domaine fluide . . . . 109
3.32 Répartition de CP sur la ligne centrale d’une roue en rotation, d’après Fackrell
et Harvey (1972) et Mears et al. (2002) . . . . . . . . . . . . . . . . . . . . . . . 110
3.33 Phénomènes caractéristiques de la roue en rotation, d’après Fackrell et Harvey
(1972) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
3.34 Géométrie du pneu considéré (fournie par Michelin) . . . . . . . . . . . . . . . 112
3.35 Vecteurs vitesse (eau et air) présentés autour du pneu . . . . . . . . . . . . . . 113
3.36 Vecteurs vitesse (eau seulement) présentés autour du pneu . . . . . . . . . . . . 113
3.37 Pression autour du pneu . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.1 Cascade d’énergie idéalisée des échelles de la turbulence (tirée de Lemay (2011)) 124
5.2 Couche limite en cisaillement : les solutions manufacturées . . . . . . . . . . . . 132
5.3 Géométrie et conditions aux limites pour l’écoulement dans un jet. . . . . . . . 137
5.4 Jet en 2D : évolultion de l’interface . . . . . . . . . . . . . . . . . . . . . . . . . 137
xii
5.5 Jet en 2D : Maillage adapté . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
5.6 Jet en 2D : solution finale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
5.7 Jet en 3D : conditions aux limites considérées . . . . . . . . . . . . . . . . . . . 139
5.8 Jet en 3D : évolution de l’interface . . . . . . . . . . . . . . . . . . . . . . . . . 140
5.9 Jet en 3D : Maillage adapté . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
5.10 Jet en 3D : solution finale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
5.11 Champs de vitesses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
5.12 Structure flexible déformée avec lignes de courant du fluide pour 4 pas de temps
différents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
5.13 Aire de contact entre un pneumatique et la chaussée en présence d’eau (tiré de
Michelin (2001)). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
5.14 Transfert des données entre les codes fluide et solide. . . . . . . . . . . . . . . . 150
5.15 Géométrie du domaine fluide. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
5.16 Coupe verticale dans le maillage fluide. . . . . . . . . . . . . . . . . . . . . . . . 152
5.17 Vecteurs vitesse du fluide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
5.18 Vecteurs vitesse du fluide selon une coupe longitudinale . . . . . . . . . . . . . 153
5.19 Courbe de vitesse suivant une ligne de coupe dans le fluide (eau). . . . . . . . . 154
5.20 Coupe de la pression sur la structure . . . . . . . . . . . . . . . . . . . . . . . . 155
5.21 Coupe longitudinale (y = 0) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
5.22 Vue de dessous et déplacement maximal . . . . . . . . . . . . . . . . . . . . . . 156
xiii
Remerciements
Tout d’abord, je tiens à remercier mon directeur de recherche, le professeur André Fortin pour
son soutien qu’il m’a offert tout au long de mes études graduées ainsi que sa disponibilité et
ses précieux conseils.
Je tiens également à remercier mon co-directeur, professeur Jean Deteix pour sa grande pa-
tience et son bon encadrement.
Je tiens aussi à remercier monsieur Patrice Hauret pour son immense encouragement à pour-
suivre mes études. J’en profite pour remercier la société MICHELIN et le CRSNG pour le
support financier dans le cadre de la chaire de recherche en calcul scientifique de haute per-
formance.
Mes remerciements vont aussi à toute ma famille : à mes frères et à ma soeur qui m’ont toujours
soutenu et accompagné tout au long de mes études. Une mention spéciale à ma femme Ons
et mon fils Adam pour leur amour, soutien et encouragements. Un merci spécial aussi à mon
père qui m’a donné un appui soutenu depuis le tout début.
Je ne voudrais surtout pas oublier, toute l’équipe du GIREF qui m’ont été d’un grand soutien.
Je tiens à remercier tous les gens qui de près ou de loin m’ont aidé à écrire cette thèse.
xv
Introduction
Les phénomènes d’interaction fluide-structure font partie de la vaste classe des problèmes
multi-physiques. Un grand nombre de situations font apparaître à ces phénomènes. Différents
régimes peuvent être distingués. Les applications en bio-mécanique font, en général, intervenir
un liquide et une structure déformable de masses volumiques proches : écoulements sanguins
dans les vaisseaux, déformation des globules rouges dans les capillaires, mais aussi la nage des
poissons. Dans le domaine du génie nucléaire, l’usure d’un faisceau tubulaire d’un échangeur
de chaleur, par instabilité sous écoulement, peut prendre à peine quelques secondes ; cet effet
de couplage est pris en compte de façon primordiale pour des raisons de sûreté des installations
de production d’énergie. La compréhension des effets de vibrations induites par écoulement
a initié de nombreuses campagnes expérimentales et justifie aujourd’hui le développement de
méthodes de calcul numérique en couplage fluide-structure. Dans le domaine du génie civil, on
cite fréquemment l’exemple de destruction du pont de Tacoma dont la compréhension a donné
lieu à une littérature scientifique abondante et qui illustre l’importance des effets d’interaction
fluide-structure. Dans le génie naval, l’exemple de l’impact hydrodynamique est l’un des effets
de couplage fluide-structure des plus parlants : il se pose pour des navires naviguant à grande
vitesse et dans des conditions de mer difficiles. Dans le domaine des pneumatiques, l’interaction
fluide-structure se manifeste dans l’étude du phénomène d’aquaplanage. L’aquaplanage est un
phénomène résultant de la perte de contact entre un pneu et la route lorsque le véhicule se
déplace à une certaine vitesse sur une route mouillée. C’est un problème important pour les
manufacturiers de pneus. Pour une vitesse donnée, l’interaction entre la flaque d’eau sur la
route et le pneu génère une pression importante à la surface du pneu dans la région de contact
entre le pneu et la réserve d’eau. Lorsque l’effort vertical généré par cette pression devient
1
supérieur au poids du véhicule, le contact entre ce véhicule et la route n’est plus maintenu et
l’aquaplanage se produit : l’adhérence entre le véhicule et la chaussée est perdue et la trajectoire
du véhicule n’est plus contrôlée. Ainsi, l’aquaplanage est responsable de la perte de contrôle
du véhicule. Ayant pour but d’augmenter la sécurité des usagers de la route, la simulation
numérique de l’interaction d’un pneu avec un liquide est un domaine de recherche actif dans
la dernière décennie. Il s’avère alors primordial d’établir les bases d’une modélisation fine de
l’interaction fluide-structure pour un pneu et de proposer de nombreux schémas numériques
à cet effet. Les principales difficultés de ce problème sont d’une part, le couplage entre la
structure et le fluide à l’interface ; d’autre part, la génération du maillage puisque la structure
est mobile, la gestion du temps de calcul et de la mémoire, etc. Toutefois, les comportements
non linéaires des fluides et des solides complexes soumis à de grands déplacements ainsi que
la déformation des interfaces induisent de nombreuses difficultés à surmonter, la résolution
analytique de tels problèmes étant impossible. En régime stationnaire, ces phénomènes d’IFS
sont étudiés depuis quelques années déjà mais restent d’une grande complexité et requièrent
d’importantes ressources informatiques. En régime instationnaire, les temps de calculs sont
démultipliés de sorte que ces phénomènes sont extrêmement complexes à modéliser et leur
compréhension reste encore délicate.
De nombreuses méthodes précises et robustes ont été développées au cours des dernières dé-
cennies afin de résoudre l’écoulement du fluide, d’une part, et la déformation de la structure,
d’autre part. Avec le développement de la puissance des calculateurs, la simulation numérique
du couplage de ces méthodes devient accessible. Cependant, il faut que ce couplage ne dégrade
pas la qualité des résultats engendrés par les méthodes dédiées aux fluides et aux solides. En
particulier, la simulation numérique se doit d’être stable pour obtenir un résultat de façon
assurée. Le fait que les méthodes dédiées aux fluides et aux solides prises séparément soient
stables n’assure en aucun cas la stabilité du système couplé. Par ailleurs, il faut assurer que
la méthode numérique obtenue est consistante avec les équations résolues afin que le résul-
tat numérique soit «proche», en un certain sens, de la solution physique du problème. On
est alors confronté aux problématiques suivantes : assurer la conservation de l’énergie lors
de la transmission des données entre les codes (fluide et solide) et gérer l’interface mobile
où s’échangent les données et où se crée le couplage. En théorie, on suppose que l’on a une
continuité des contraintes et des vitesses à l’interface entre le fluide et la structure. Ces condi-
tions sont très largement utilisées dans les différents modèles d’interaction. La difficulté est
de trouver une formulation adaptée pour la résolution des équations fluide-structure qui tient
compte des parois mobiles et permet de simuler des écoulements fortement convectifs. Deux
types d’approches peuvent être envisagées : une résolution explicite ou une résolution impli-
cite. La résolution explicite découplée consiste, en général, à résoudre le problème fluide sur
un domaine fixe, puis à intégrer une seule fois les efforts du fluide sur la structure pour en
déduire le déplacement. Cette approche ne donnera un résultat raisonnable que dans le cas où
le mouvement de la structure est extrêmement lent par rapport au temps caractéristique de
2
l’évolution du fluide. Cette approche est instable si le fluide est incompressible. Au contraire,
nous souhaitons étudier des problèmes où la dynamique de la structure est comparable à celle
du fluide. Dans ce cas, une approche totalement implicite permettant une résolution complète
du système est utilisée.
Une difficulté majeure tient au fait que le domaine de résolution du fluide varie en fonction
du temps suivant le déplacement du solide avec, dans les cas qui nous intéressent, de grandes
déformations de ce domaine. Trois principaux types d’approches sont envisageables pour traiter
la déformation du domaine fluide : une approche multiphasique, une approche avec déformation
de maillage ou une approche par domaine fictif. L’approche multiphasique (lagrangienne ou
eulérienne) serait limitée au cas où le solide et le fluide seraient décrits par les mêmes équations
avec des paramètres physiques variables affectés à chaque phase et advectés au cours du
mouvement de l’interface. En général, les méthodes eulériennes (le domaine de calcul est fixe
et ne suit pas le matériau) sont préférées pour le calcul de l’écoulement fluide, afin d’éviter
de trop grandes distorsions des mailles de calcul. Au contraire, le solide est plutôt simulé avec
des méthodes lagrangiennes (le domaine se déforme en suivant le déplacement matériel), qui
permettent de suivre les discontinuités matérielles avec précision. Les approches de déformation
de maillage permettent d’assurer une transition lisse entre ces deux types de modélisation.
L’exemple typique d’une telle approche est la méthode «Arbitrary Lagrangian – Eulerian»
(ALE) : le domaine fluide est déformé de façon à suivre l’interface fluide-solide. Le couplage
de ces différentes formulations rend l’étude de ces systèmes complexes, tant d’un point de
vue numérique que mathématique. La déformation du maillage induit des termes additionnels
dans la résolution des équations du fluide. Cependant, des difficultés apparaissent lorsque le
déplacement de la structure est trop important : le domaine fluide se déforme fortement ce qui
nécessite l’introduction de nouvelles méthodes de mise à jour du maillage fluide (voir chapitre
3). La troisième approche est la méthode des domaines fictifs. Cette méthode a été développée
pour traiter les cas des interactions entre un fluide incompressible et une structure rigide (voir
chapitre 2).
Dans le cadre de cette thèse, pour traiter le cas des structures flexibles, on a fait le choix d’une
formulation ALE pour le fluide. Cette formulation permet de gérer des maillages mobiles avec
des déplacements d’amplitude modérée et/ou importante, c.à-d. qu’elle permet de résoudre
les équations du problème fluide dans un domaine mobile s’adaptant aux déformations de la
structure. Nous avons utilisé une formulation lagrangienne pour la structure, où la simulation
numérique est basée sur la méthode des éléments finis. Le couplage entre les codes fluide et
structure se fait de la façon suivante : pour un pas de temps fixé et dans une grande boucle de
point fixe, on calcule dans un premier temps, les chargements fluides exercés sur la structure.
Ensuite, ce chargement est utilisé pour le calcul mécanique et on en déduit le déplacement de
la structure utilisée comme condition limite du calcul fluide à l’itération suivante du même
point fixe. On itère entre les deux solveurs fluide et solide jusqu’à convergence du point fixe
3
et ensuite on passera au prochain pas de temps. Nous notons la présence d’un mécanisme de
mouvement de maillage fluide à chaque itération du point fixe.
— le développement d’une méthode de domaine fictif avec des conditions aux limites pré-
cises pour le cas des structures rigides.
Dans le premier chapitre, nous présenterons les modèles fondamentaux et les méthodes numé-
riques développées pour la simulation des écoulements fluides, dans un premier temps, et des
équations d’élasticité linéaire, dans un second temps.
Dans le troisième chapitre, le schéma de couplage est étendu au cas d’un solide déformable.
Ce chapitre contient des discussions sur les formulations numériques des interactions fluide-
structure (monolitique vs partitionnée). Après une introduction des différentes méthodes de
couplage fluide-structure, nous présenterons en détail la méthode ALE. Un algorithme de
couplage est alors mis au point pour la mise en oeuvre de cette interaction fluide-structure. Cet
algorithme, basé sur un schéma implicite, permet le transfert de champs de façon interactive.
Nous avons détaillé les techniques de couplage utilisées dans les présentes études, l’algorithme
de couplage et aussi des exemples tests. Ces exemples tests concernent les écoulements de
fluide dans des conduites en présence d’obstacles. Les obstacles sont considérés comme solides
déformables. Notre approche a été testée et se compare avantageusement aux résultats de la
littérature sur plusieurs problèmes tests.
Au chapitre 4, nous présenterons une comparaison de la méthode des domaines fictifs avec
l’approche ALE dans le cadre des obstacles rigides. Les tests de vérification permettant d’ap-
précier la performance et la fiabilité des deux méthodes y seront également présentés.
4
Au chapitre 5, une stratégie de complexification du problème d’IFS a été définie. La modélisa-
tion de la turbulence et des écoulements à surfaces libres ont été introduites. Nous avons fixé
notre choix de modèle en considérant les équations de Navier-Stokes moyennées basées sur le
modèle URANS (Unsteady Reynolds Averaged Navier Stokes). En présence d’une surface libre,
nous avons introduit la méthode des surfaces de niveau très largement utilisée pour ce type de
problème. Ainsi, nous présenterons plusieurs algorithmes de couplage multi-physique et nous
montrons la robustesse des méthodes de couplage développées sur des résultats numériques
présentant un intérêt pour les applications industrielles futures visées.
Une synthèse des différents résultats et discussions obtenus dans cette thèse sera donnée dans
la conclusion.
5
Chapitre 1
Nous présenterons, dans ce chapitre, les équations écrites du point de vue de la mécanique des
milieux continus et utilisées dans les cas d’un fluide et d’un solide. Les formulations faibles
associées à ces équations seront construites. Nous allons utiliser la formulation eulérienne au
niveau du fluide. Cette approche est l’une des plus couramment utilisées en mécanique des
fluides que ce soit en éléments finis, en différences finies ou en volumes finis. Elle permet d’ob-
tenir des solutions numériques pour des maillages fixes dans l’espace. Néanmoins, le domaine
solide est décrit généralement dans un cadre lagrangien où l’on suit une particule dans son
mouvement et on regarde sa position à un instant quelconque. La particule est identifiée par
sa position initiale et non pas par sa vitesse.
Ce chapitre est articulé autour de deux grandes parties. La première partie est une présentation
des équations fluides et des méthodes de résolution par la méthode des éléments finis, adoptées
dans cette thèse. La deuxième partie sera consacrée à la présentation des équations pour le
solide et des méthodes de résolution qui s’y rattachent. Notons que le couplage entre ces
équations sera effectué au chapitre 3.
Un fluide est un objet très déformable qui n’a pas de forme propre. La modélisation des fluides
7
a été introduite par Euler au 17ième siècle puis complétée par Navier et Stokes au 18ième siècle,
notamment par la prise en compte de la viscosité. Les équations régissant ce type d’écoulement
portent le nom d’équations de Navier–Stokes. Ces équations sont à la base de la mécanique
des fluides incompressibles et proviennent des ingrédients suivants :
— les axiomes de la mécanique : conservation de la masse et de la quantité de mouvement,
— la contrainte d’incompressibilité,
— la loi de comportement rhéologique des fluides newtoniens.
Les axiomes de la mécanique stipulent que la masse de tout volume de particules transportées
par le mouvement est constante et que la dérivée temporelle de sa quantité de mouvement est
égale à la résultante des forces appliquées.
La contrainte d’incompressibilité dit que l’on peut considérer le volume de tout domaine
transporté par l’écoulement comme étant constant. C’est une approximation qui est valide
lorsque la vitesse de l’écoulement est petite devant la vitesse du son dans le fluide.
La loi de comportement rhéologique des fluides newtoniens demande plus d’attention et sera
décrite à la section suivante.
σ = 2µγ̇(u) − pI
8
Remarque :
Le mouvement d’un fluide est gouverné par l’équation de Navier–Stokes qui s’écrivent par :
Cette équation est valable dans l’ensemble du domaine Ω, occupé par le fluide, et elle doit être
complétée par une condition aux limites sur la frontière ∂Ω :
u = uD sur ∂Ω (1.2)
Notons qu’il existe d’autres conditions aux limites qui peuvent être appliquées ; mais pour
alléger l’explication sans perte de généralités, nous présentons seulement la condition 1.2.
Remarque :
De façon générale, l’équation (1.1) doit être aussi complétée par l’équation locale de conser-
vation de la masse. Dans le cas compressible,
∂t ρ + ∇ · (ρ u) = 0. (1.3)
Dans ce cas, il faut également, pour fermer le problème, prescrire la relation thermodynamique
qui relie la pression p à la masse volumique ρ et possiblement à d’autres paramètres tels
que la température, le champ magnétique, etc. Dans cette thèse, notre étude sera consacrée
uniquement aux fluides incompressibles. La masse volumique ρ est constante et uniforme.
Quant à l’équation locale de conservation de la masse (1.3), elle se simplifie sous la forme :
∇ · u = 0. (1.4)
Dans ce cas, la pression n’est plus connue à travers sa relation p(ρ) (nous sommes dans la limite
où cette relation est très raide de sorte que p évolue de façon significative, tandis que ρ peut être
considérée comme étant essentiellement constante). On peut voir p comme un multiplicateur de
Lagrange associé à la propriété d’incompressibilité (1.4). De ce fait, en appliquant l’opérateur
divergence à l’équation (1.1) et en utilisant la condition d’incompressibilité (1.4), on obtient :
∆p = −ρ∇ · (u · ∇ u) + ∇ · f (1.5)
9
∆u, on obtient
∂ ∂ ∂
∇ · ∆u = ∆u1 + ∆u2 + ∆u3
∂x ∂y ∂z
∂ ∂ 2 u1 ∂ 2 u1 ∂ 2 u1 ∂ ∂ 2 u2 ∂ 2 u2 ∂ 2 u2
= + 2 + 2 + + 2 + 2 (1.6)
∂x ∂ 2 x ∂ y ∂ z ∂y ∂ 2 x ∂ y ∂ z
2 2 2
∂ ∂ u3 ∂ u3 ∂ u3
+ + 2 + 2
∂z ∂ 2 x ∂ y ∂ z
2 2 2
∂ ∂ ∂
= 2
+ 2 + 2 ∇·u
∂ x ∂ y ∂ z
= ∆ (∇ · u) = 0.
On voit donc apparaître un problème central dans la mise en œuvre des méthodes numériques
pour la résolution des équations de Navier-Stokes incompressibles. En effet, l’équation de Pois-
son (1.5) sur p est une équation elliptique qui doit être munie de conditions de bord pour que le
problème soit bien posé, c.-à-d., pour que la solution soit unique. Cette condition peut être par
exemple soit de type Dirichlet ou de type Neumann. Or, on ne sait pas écrire proprement ces
conditions, dans ce cas. La pression p n’ayant pas de condition de bord naturelle (physique),
les modèles introduisent alors, des conditions aux limites arbitraires pour effectuer ce calcul.
Si cela permet de calculer un multiplicateur de Lagrange et de forcer l’incompressibilité dans
une certaine mesure, la valeur de la pression reste inévitablement entachée par ces conditions
arbitraires.
En principe la condition de bord sur p est définie par l’équation de Navier-Stokes (1.1) elle-
même, qui prescrit l’ensemble du gradient de pression ∇p au bord, et il est facile de se
convaincre qu’imposer ∇p au bord, sous réserve de régularité suffisante, c’est à la fois im-
poser :
• une condition de dérivée tangentielle sur ∂Ω
10
des champs de vitesse qui sont la solution des équations de NS. On comprend, dès lors, que
les algorithmes garantissent mal que le champ de vitesse reste dans l’espace de ces champs
compatibles et qu’il n’est plus possible de contrôler certaines erreurs au bord. Ce phénomène
est particulièrement visible lorsqu’on utilise les méthodes de types pas fractionnaires, comme
par exemple l’algorithme de projection de Chorin-Temam qui sera décrit plus tard.
∂u
ρ ∂t + u · ∇u − 2µ div(γ̇(u)) + ∇p = f dans Ω
(1.9)
= 0 dans Ω
∇·u
A ces équations, il faut ajouter les conditions aux limites appropriées suivantes :
∂u 2
∂t + u · ∇u + ∇p − Re ∇ · γ̇(u) = f
(1.12)
∇·u = 0
où Re est le nombre adimentionnel de Reynolds qui est défini par :
U Lρ
Re =
µ
où U et L sont respectivement des vitesses et longueurs caractéristiques.
Puisque, dans le cadre de notre étude, on s’intéresse à des cas instationnaires, on ne peut pas
∂u
négliger le terme et on cherchera alors à obtenir une approximation de la fonction u(t)
∂t
11
pour t0 ≤ t ≤ tf . En utilisant l’approche par semi discrétisation proposée par Lions (1969),
Quarteroni et Valli (1999), on considère une partition de l’intervalle de temps [t0 ,tf ] et on
considère, en premier lieu, un schéma d’Euler implicite (voir Fortin et Garon (2013)) :
∂ui n un − un−1 ∆t ∂ 2 ui n
(t ) = i i
+ (ξ ) avec ξ n ∈ [tn−1 , tn ] (1.13)
∂t ∆t 2 ∂t2 i
Le système 1.15 est non linéaire par la présence du terme un · ∇un . L’utilisation des méthodes
de type points fixes, telle que la méthode de Newton, est envisageable (voir par exemple
Quarteroni et Valli (2008), Boffi et al. (2013)). Notons, toutefois, qu’il est plus efficace d’utiliser
une extrapolation linéaire de Lagrange de (tn−2 , uin−2 ) à (tn−1 , un−1
i ) pour linéariser le terme
un · ∇un . En effet, puisqu’on a
∂ 2 ui n
uni = 2un−1 − un−2 + (µ )(∆t)2 avec µni ∈ [tn−2 , tn ] (1.16)
i i
∂t2 i
alors, on remplace un · ∇un par (2un−1 − un−2 ) · ∇un . Mentionnons que l’approche favorisée
par Turek (1999) est basée sur ce type d’extrapolation. Le problème à résoudre à chaque pas
de temps est maintenant linéaire et s’écrit :
3 n 2 1
u + 2un−1 − un−2 · ∇un − ∇ · γ̇(un ) + ∇pn = f + 4un−1 − un−2
2∆t Re 2∆t
∇ · un = 0
(1.17)
Pour obtenir la formulation variationnelle, on multiplie le système (1.17) respectivement par
des fonctions test v f ∈ V et q ∈ Q où V = (HΓ1 (Ω))d (d est la dimension de l’espace) et
12
Q = L2 (Ω), on obtient par conséquent, la formulation variationnelle du problème 1.17
Z
3 n f 2 n f f n−1 n−2 n
f
u ·v + γ̇(u ) : γ̇(v ) − p∇·v + (2u −u ) · ∇u · v dv =
Ω 2∆t Re
Z Z
1
Z
f f
f · v dv + t · v dS + (4un−1 − un−2 ) · v f dv
2∆t
Ω ΓN Ω
Z
q∇ · un dv = 0
Ω
(1.18)
La discrétisation par éléments finis du problème de Navier-Stokes suit les étapes habituelles
de discrétisation. Soit τh (Ω) une discrétisation de Ω de paramètre h. Nous considérons, dans
ce travail, des maillages constitués d’éléments notés K. En dimension 2, nous considérons les
triangles et les quadrilatères, tandis qu’en dimension 3, nous nous limiterons aux tétraèdres et
aux hexaèdres. On choisit deux espaces Vh ⊂ V et Qh ⊂ Q de dimension finie et on approxime
(u,p) par des fonctions (uh ,ph ) ∈ Vh × Qh . Le problème discrétisé du problème 1.18 peut
s’écrire :
f f f f
a(uh ,v h ) + b(v h ,ph ) = (f,v h ) ∀v h ∈ Vh
(1.19)
b(uh ,qh ) = 0 ∀qh ∈ Qh
Sous les hypothèses d’existence et d’unicité des solutions (uh ,ph ) du problème discret de
Navier–Stokes (voir Boffi et al. (2013) pour plus de détails), on part des solutions initiales
(un−1 ,pn−1 ) et (un ,pn ) et on cherche une correction (δuh ,δph ) de sorte que (un +δuh ,pn +δph )
soit une solution du problème de Navier–Stokes. On obtient :
Z Z
f
an (δuh ,v h ) − f
δph ∇ · v h dv = − Rn (v fh ) dv
Ω Ω
(1.21)
Z Z
qh ∇ · δuh dv = − qh ∇ · un dv
Ω Ω
13
Z
3 δu f 2
an δu,v f γ̇(δu) : γ̇(v f ) + ((2un−1 − un−2 ) · ∇)δu · v f
= v + dv
Ω 2 ∆t Re
Z
3 n f 2
Rn (v fh ) n f f
= u ·v + γ̇(u ) : γ̇(v ) − p∇·v dv
Ω 2∆t Re
Z Z Z
n−1 n−2 n f f
t · v f dS
+ (2u −u ) · ∇u · v dv − f · v dv −
Ω Ω ΓN
Z
1
(4un−1 − un−2 ) · v f dv
−
2∆t Ω
(1.22)
Posons :
nu
XD
K
αjK ΨK
u(x)|K ' uh (x)|K = u (x) = u,j
j=1
npD
X
K
pK K
p(x)|K ' ph (x)|K = p (x) = j ψp,j
j=1
p
u,j sont des fonctions de base vectorielles et que les nD fonctions de base en
Notons que les ΨK
pression et les nuD fonctions de base en vitesse sont généralement de natures différentes (Fortin
et Garon (2013)). Les approximations en vitesse et en pression ne peuvent pas être choisies
d’une manière arbitraire puisqu’elles doivent vérifier la condition inf-sup. Le lecteur peut se
référer à Boffi et al. (2013) pour plus de détails. Il existe plusieurs choix possibles. Cependant,
le choix doit être en accord avec le phénomène à l’étude. Il s’agit, ici, d’interactions fluide-
structure exigeant de la précision au passage de l’information à l’interface fluide solide. Ainsi,
un élément compatible tel que l’élément Mini décrit dans Arnold et al. (1984), pour lequel la
discrétisation de la pression est linéaire tandis que la vitesse est linéaire mais enrichie d’un
nœud de calcul au centre de l’élément, ne sera pas considéré. En effet, cet élément converge à
l’ordre 1 en espace et donc il est peu précis. De plus, il ne permet pas un transfert adéquat
des contraintes fluides vers la structure dans les cas de problèmes d’IFS. En se basant sur
14
Figure 1.1 – Élément de Taylor-Hood P2 − P1 (O(h2 )) en 2D et 3D.
le contexte IFS, nous avons choisi l’élément de Taylor-Hood souvent appelé P2 − P1 où les
u,j sont quadratiques et les ψp,j linéaires sur un élément (voir Brezzi et Fortin (1991) et
ΨK K
Boffi et al. (2013)). Pour cet élément, les nœuds de calcul sont situés aux sommets (vitesse et
pression) et aux milieux d’arêtes (vitesses seulement) (voir la figure 1.1) et on obtient avec cet
élément une convergence à l’ordre 2 en espace (voir Boffi et al. (2013) pour plus de détails).
L’élément P2 − P1 en dimension 3 est représenté aussi à la figure 1.1. Ces approximations
en vitesse et en pression vérifient bien la condition inf-sup et sont dites compatibles. Bien
entendu, d’autres choix sont possibles. Cependant, l’approche que nous présenterons exige, en
plus de la précision au passage des contraintes fluides et structures à l’interface, une continuité
de la discrétisation en pression.
Dans cette section, nous nous consacrons à la description des solveurs utilisés pour résoudre
numériquement les équations de Navier–Stokes. Nous avons considéré trois types de résolution
pour ces équations. La première est la résolution par un solveur direct. La deuxième est
itérative de type point-selle. Quant à la troisième, c’est la résolution par une méthode de
projection.
15
Résolution par une méthode directe
Dans cette partie, on s’intéresse à une brève description d’un solveur itératif pour résoudre
numériquement les équations de Navier-Stokes.
Afin de réduire les coûts inhérents aux éléments quadratiques, nous avons utilisé une méthode
de résolution multi-niveaux basée sur la hiérarchie naturelle entre les éléments finis linéaires et
quadratiques, d’où le nom de méthode hiérarchique. Elle possède plusieurs points en commun
avec les méthodes multi-grilles mais a l’avantage de s’appliquer aux géométries complexes
et aux maillages non structurés. L’utilisation de cette méthode comme préconditionneur à
une méthode de Krylov (ici un préconditionneur variable) permet d’obtenir une méthode très
efficace.
Le système 1.23 peut s’écrire matriciellement sous forme d’un système linéaire non symétrique
de type point-selle à deux inconnues u et p. La résolution efficace de ce système joue un rôle
majeur dans le traitement numérique des équations de Navier-Stokes. Pour cela, nous avons
utilisé un préconditionneur à droite de format triangulaire par bloc. Pour rendre ce précon-
ditionneur efficace, nous avons fait appel à trois ingrédients : l’ajout du terme r∇(div(u))
aux équations continues de Navier-Stokes où sa forme faible s’écrit par rdiv(u)div(v) (voir
Raviart (1981) et Franca et Frey (1992)), une résolution efficace en vitesse par la méthode
hiérarchique et un préconditionneur additif pour le complément de Schur (voir Cahouet et
Chabard (1988) et Turek (1999)). Pour plus de détails sur cette méthode de résolution, le
lecteur peut se référer à El maliki (2007), El maliki et Guénette (2010).
Confrontés aux défis numériques lors de la résolution des équations de NS, Chorin et Temam
ont produit, dans les années soixante séparément, une série d’articles qui introduisaient et
analysaient la méthode de projection (voir Chorin (1968) et Temam (1968)). Cette méthode
a été introduite dans le but de résoudre efficacement les équations de NS. Elle est basée
sur des schémas en temps fractionnaire qui découple le problème de convection-diffusion de
16
la contrainte d’incompressibilité, avec correction de la vitesse et de la pression. Ainsi, elle
devint très populaire dans les applications d’écoulements visqueux et incompressibles, voire
pour les écoulements à faible nombre de Mach. En effet, d’un point de vue mathématique, la
dynamique de ces écoulements (à faible nombre de Mach) est similaire à celle des écoulements
incompressibles. D’une part, dans l’équation de conservation de la quantité de mouvement,
la masse volumique est indépendante de la pression dynamique. D’autre part, le champ de
vitesse satisfait une contrainte équivalente à la condition de divergence nulle des équations
incompressibles.
Ces méthodes procèdent dans une première étape par le calcul d’un champ de vitesse intermé-
diaire ũ en utilisant l’équation de quantité de mouvement. Ce champ ne vérifie pas la condition
d’incompressibilité ∇ · ũ = 0, et il est donc nécessaire de le projeter dans un espace où les
champs de vecteur sont à divergence nulle pour obtenir un+1 . Le schéma le plus simple, qui
est le schéma explicite d’ordre 1 en temps, se présente sous sa forme semi-discrète :
Étape 1 :
ũ un
ρ + ρ(un · ∇) un − µ∆un = f + ρ (1.24)
∆t ∆t
Étape 2 :
∆t
ũ = un+1 + ∇pn+1
ρ (1.25)
n+1
∇·u = 0.
Un schéma implicite d’ordre 1 en temps peut aussi être écrit sous la forme :
Étape 1 :
un
ρ ũ + ρ(un · ∇) un − µ∆ũ = f +ρ
∆t ∆t (1.26)
ũ = uD sur ∂Ω.
Étape 2 :
∆t
ũ = un+1 + ∇pn+1
ρ (1.27)
n+1
∇·u = 0.
La présence de la condition aux limites sur ũ pour résoudre (1.26) constitue une très importante
différence entre la résolution des problèmes (1.24)-(1.25) et (1.26)-(1.27). Dans ces schémas,
on calcule la pression en résolvant une équation de Poisson provenant de l’incompressibilité
de un+1 . Ce terme de pression intervient dans l’équation de quantité de mouvement comme
un multiplicateur de Lagrange assurant que la vitesse vérifie la contrainte d’incompressibilité.
Ces problèmes découlent des équations (1.25) ( ou bien (1.27)), en leur appliquant l’opérateur
divergence pour obtenir :
ρ
∆pn+1 = ∇ · ũ (1.28)
∆t
17
Cette dernière équation est similaire au problème (1.5) : La résolution numérique de ce pro-
blème de Poisson passe, donc, par l’introduction de conditions aux limites qui sont plus ou
moins arbitraires, de type dérivée tangentielle (1.7) ou de type Neumann (1.8). Afin d’expli-
quer simplement ce qui se passe sans perte de généralité, prenons une condition de Dirichlet
homogène sur la vitesse, c.-à-d. u = 0 sur le bord ∂Ω. Combinant cela avec les conditions
(1.7) et (1.8), on obtient :
∇p · τ = µ ∆u · τ + f · τ , (1.29)
∇p · n = µ ∆u · n + f · n. (1.30)
∆t
ũ · n = un+1 · n + ∇pn+1 · n, (1.31)
ρ
Or, ceci est incohérent avec la condition (1.30). Cette inconsistance (due aux conditions aux
limites) provoque naturellement des erreurs, soit sur le calcul de la pression, soit sur la vitesse.
Ces erreurs sont localisées au bord du domaine dans une couche limite. On consultera les
articles suivants : Gresho et Sani (1987), Weinan et Jian-Guo (1995), E et Liu (1996), Dagan
(2003) pour plus de détails.
Dans E et Liu (1996), la localisation de cette erreur au niveau de la couche limite est expliquée
grâce à une étude de ses modes propres. En effet, en combinant les équations (1.26), (1.27)
et (1.32) sans le terme convectif, afin que les explications soient plus claires, on obtient le
problème suivant : (
∆pn+1 − ∆t∆2 pn+1 = 0
(1.33)
∇pn+1 · n = 0 sur ∂Ω.
∆pn+1 = 0. (1.34)
Dans Dagan (2003), Weinan et Jian-Guo (1995), les auteurs montrent que ce problème est
singulier et que la couche limite est d’ordre O(∆t1/2 ) dans l’approximation de la pression.
Dans Patankar et Sharma (2005), une méthode de projection a été présentée afin de remédier
à ces inconsistances numériques. Cette méthode consiste à décomposer un pas de temps en
une étape de prédiction, sans prise en compte de la contrainte, puis à construire le champ à
l’étape suivante par projection restreinte aux particules sur l’espace des mouvements rigides.
18
Les champs construits sont discontinus aux interfaces, sans briser la convergence en norme L2
en espace.
Les améliorations des schémas proposées par Chorin et Temam demeurent un enjeu essentiel
dans la résolution des équations de Navier Stokes. Beaucoup de travaux sur les méthodes à
pas fractionnaires proposent des solutions aux problèmes cités précédemment selon différents
principes, en fixant de manière exacte les conditions aux limites :
• Soit en vitesse. Le schéma le plus connu est celui de Leveque et Oliger (1983), Kim et
Moin (1985)
• Soit en pression, voir les travaux de Orszag et al. (1986), Karniadakis et al. (1991),
• Soit en proposant des schémas corrigés. La version la plus précise connue actuellement
est la méthode de correction de pression proposée par Goda (1979), van Kan (1986) et
connue sous le nom de schéma incrémental de correction de pression .
19
— Étape 4 : mise à jour de la vitesse
2δt
uk+1 = ũk+1 − ∇ϕk+1 , (1.38)
3ρ
Remarques :
— Dans ce cas, on a des conditions aux limites artificielles
— Shen (1996), Guermond (1997, 1999) et Guermond et Quartapelle (1998) ont montré
que le schéma incrémental de correction de pression est d’ordre maximal égal à 1 en
temps. Notons que nous avons pris, dans ce cas, un schéma de différence arrière d’ordre
2 pour la discrétisation temporelle du terme instationnaire de Navier-Stokes. Ce schéma
incrémental de correction de pression nous fait perdre un ordre.
||u − uexacte ||H 1 (Ω)d + ||p − pexacte ||L2 (Ω) ≤ Cδt (1.39)
L’idée proposée par Timmermans et al. (1996) consiste à utiliser l’identité suivante :
∆u = ∇(∇ · u) − ∇ × (∇ × u).
Remarques :
— Cette méthode implique une équation consistante pour la pression telle que :
∆pk+1 = ∇ · f k+1 ,
k+1
k+1 k+1
∇p ·n = f −∇×∇×u · n,
20
Dans ce travail, nous allons utiliser deux schémas principaux pour la résolution des équations
de Navier Stokes. Dans le chapitre 2, nous allons utiliser la méthode de type point-selle alors
que, dans le chapitre 3, nous allons considérer le schéma incrémental de correction de pression.
Les résultats numériques des deux chapitres montrent bien que ces méthodes fonctionnent bien
dans le cas de la résolution des problèmes d’interaction fluide-structure.
Dans un premier temps, nous considérons le cas de la structure rigide considérée par plusieurs
auteurs tels que Glowinski et al. (2000), Su et al. (2007) et Valette et al. (2007). Plus récem-
ment, mentionnons que Jendoubi (2010) et aussi Ilinca et Hétu (2011) ont utilisé ce type de
structures pour étudier l’interaction d’un corps rigide avec un fluide.
La modélisation des structures rigides peut aussi tenir compte de la dynamique propre du
solide en lui associant une masse (et donc une inertie). Le profil subit alors, par conséquent,
des mouvements de corps rigide soumis à la pression du fluide. On parle, alors, d’aéroélasticité
(Causin et al. (2005)).
Dans un deuxième temps, nous considérons le modèle des structures flexibles qui tient compte
des déformations structurelles et introduit des équations des milieux continus régissant ces
déformations. Il existe de nombreuses équations dans la mécanique des milieux continus per-
mettant de modéliser les déplacements et les déformations d’une structure. On peut, toutefois,
les classer en 3 catégories (sans parler de la plasticité bien au-delà du champ d’investigation
de notre étude) :
1. Élasticité linéaire,
2. Hyper-élasticité,
3. Hyper-élasticité finie.
Ainsi, un modèle d’élasticité linéaire suffit pour modéliser les structures en petits déplacements.
Dans l’hypothèse des petites perturbations, les modèles d’élasticité linéaire peuvent être facile-
ment mis en œuvre. Nous citons plusieurs auteurs qui ont considéré cette catégorie (Pederzani
et Haj-Hariri (2006), Borazjani et Sotiropoulos (2008), Huang et Sung (2009), Van Loon et al.
(2006)). Toutefois, dès qu’on a de grands déplacements, un modèle hyperélastique est néces-
saire. De nombreuses lois de comportement non linéaires existent mais c’est celle de Saint-
Venant Kirchhoff qui est la plus répandue dans le cadre de la mécanique des solides (Turek
et Hron (2006), Gerbeau et Vidrascu (2003), Von Scheven (2009), Étienne et Pelletier (2004),
Huber et al. (2004)). Les autres modèles hyperélastiques tels que le modèle Mooney-Rivlin
21
(Lian et Shyy (2005)) et le modèle Néo-Hookéen (Wood et al. (2008)) sont encore peu utili-
sés dans le domaine des IFS. Dans ce cadre, nous pouvons citer par exemple Huang et Sung
(2009) qui ont considéré des structures flexibles modélisées par la dynamique des solides dans
des problèmes en biologie.
Notons la déformation du solide par T s : Ω̂s × [0,T ] → Ωs (t) alors, on introduit le gradient de
s déf déf s
déformation correspondant par F̂ (x̂,t) = ∇x̂ T s (x̂,t) et son déterminant Jˆs = dét (F̂ (x̂,t)).
s déf
Le déplacement du domaine est donné par d̂ (x̂,t) = T s (x̂,t) − x̂ et sa vitesse est notée par
s
ûs (x̂,t) = ∂t d̂ (x̂,t).
Ayant adopté une formulation lagrangienne, le comportement de la structure est ainsi régi par
les équations d’équilibre exprimées sur la configuration initiale non déformée Ω̂s ce qui permet
de définir le problème de la structure comme suit :
s s
Trouver le déplacement d̂ = d̂ (x̂,t) : Ω̂s × [0,T ] → Ωs (t) telle que
s
s ∂ û
s s
0 ∂t − ∇ · Π̂ (d̂ ) = ρ̂0 r̂ dans Ω̂ ,
ρ̂ s s s
(1.44)
s
= ûs dans Ω̂s ,
∂t d̂
déf
où ρ̂s0 = Jˆs ρ̂s , ρ̂s est la densité du solide et r̂ s est une force de volume donnée. Le tenseur
s s s
Π̂ = Π̂ (d̂ ) est appelé le premier tenseur de Piola-Kirchoff (voir Bonet et Wood (2008)) qui
est relié au tenseur des contraintes de Cauchy par :
s déf s
Π̂ = Jˆs σ̂ s (F̂ )−T , (1.45)
22
et les conditions aux limites suivantes
s s s
Π̂ n̂s = Jˆs ||(F̂ )−T n̂s ||ĥ sur Γ̂sN (1.47)
s s
d̂ = d̂D sur Γ̂sD (1.48)
s
où d̂D et ĥ sont des fonctions données, n̂s est la normale unitaire à ∂ Ω̂s et Γ̂sD ∪ Γ̂sN = ∂ Ω̂s .
L’équation 1.48 impose les déplacements sur une partie du bord Γ̂sD ⊂ ∂ Ω̂s alors que l’équation
1.47 impose la contrainte superficielle sur Γ̂sN ⊂ ∂ Ω̂s , écrite sous forme lagrangienne qui est
équivalente à la condition :
σ s ns = hs sur ΓsN
écrite en description eulérienne.
Beaucoup de lois de comportement peuvent être utilisées pour un solide. Dans ce qui suit,
nous fournissons la formulation générale d’un matériau hyperélastique. Pour introduire cette
formulation, on définit le deuxième tenseur de Piola-Kirchhoff S qui est un tenseur symétrique
donné par :
déf s s s
S = (F̂ )−1 Π̂ = Jˆs (F̂ )−1 σ̂ s (F̂ )−T , (1.49)
Modèle de Saint-Venant-Kirchoff
Nous ferons un rapide survol sur les caractéristiques de ce modèle. Pour plus de détails, le
lecteur se refèrera à Bonet et Wood (2008). Le modèle de Saint-Venant-Kirchoff est la généra-
lisation aux grandes déformations du modèle utilisé en élasticité linéaire (petites déformations).
Le potentiel d’énergie prend la forme (voir par exemple Bonet et Wood (2008)) :
1
Ψ = λ(trE)2 + µ(E : E) (1.51)
2
Les paramètres constants λ et µ sont les coefficients de Lamé définis par :
E Eν
µ= et λ = (1.52)
2(1 + ν) (1 + ν)(1 − 2ν)
où E est le module d’élasticité (ou module d’Young) et ν est le coefficient de Poisson du
matériau. E est le tenseur de déformation de Green-Lagrange défini par :
1
E = (C − I) (1.53)
2
23
À partir des relations 1.50, 1.51, 1.52 et 1.53, on vérifie facilement que :
∂Ψ
S = λ(trE)I + 2µE = (1.54)
∂E
On définit le tenseur d’élasticité d’ordre 4 par :
∂S
C= (1.55)
∂E
qui a pour composantes :
I ik I jl + I il I jk
Cijkl = λI ij I kl + 2µ (1.56)
2
et qui permet d’obtenir une «généralisation» de la loi de Hooke : la forme linéarisée de la
conservation de la quantité de mouvement aura la forme d’une loi de Hooke (voir Fortin et
Garon (2013) pour plus de détails).
p = −k(J − 1) (1.59)
S = S 0 − pJC −1 (1.60)
Pour plus d’informations sur les matériaux en grandes déformations, on réfère le lecteur à For-
tin et Garon (2013) et Bonet et Wood (2008).
24
1.2.3 Formulation éléments finis
Pour la structure, on définit l’espace suivant :
s déf s s d s 1 s d s
V̂ = v̂ : Ω̂ → R ,v̂ ∈ [H (Ω̂ )] , v̂ |Γ̂s = 0 , (1.62)
D
En multipliant l’équation du solide 1.44 par v̂ s ∈ V̂ s et en intégrant par parties tout en tenant
compte des conditions aux limites 1.47 et 1.48 on obtient :
∂ ûs s
Z Z Z Z
s
ρ̂s · v̂ dx̂ + Π̂ : ∇x̂ v̂ s dx̂ = ρ̂s r̂ s v̂ s dx̂ + Jˆs ||(F̂ )−T n̂s ||ĥ · v̂ s dγ̂ (1.63)
Ω̂s ∂t Ω̂s Ω̂s s
Γ̂N
N
X
dˆh (t,x̂) = d̂i (t)ϕi (x̂).
i=1
on obtient le système non linéaire d’équations différentielles ordinaires (EDO) donné par :
¨
Ms d̂s + Ks (d̂s ) = F s , t>0 (1.64)
avec :
— d̂s = [dˆ1 (t),.....,dˆN (t)]T est le vecteur solution
Z
— (Ms )ij = ρ̂s ϕi ϕj dx̂ est la matrice masse
Ω̂s
Z
— (Ks (d̂s ))i = Π̂ : ∇x̂ ϕi dx̂ est le terme non-linéaire de raideur
Ω̂s
1.2.4 Discrétisation
On peut utiliser un solveur d’équations différentielles ordinaires pour résoudre le système 1.64.
Le schéma populaire pour résoudre ce type de problème est celui de Newmark (voir Bathe
(1996) et Zienkiewicz et Taylor (1989)) décrit dans ce qui suit.
n n
On définit dns ≈ ds (tn ), ḋs ≈ ḋs (tn ), d̈s ≈ d̈s (tn ) et deux paramètres γ ∈ [0,1] et β ∈ [0, 12 ].
Le système :
n
Ms d̈s + Ks (dns ) = F ns , à tn
n n−1 n n−1
ḋs = ḋs + ∆t(γ d̈s + (1 − γ)d̈s ), (1.65)
Remarque : Pour des systèmes linéaires de deuxième ordre, le système 1.65 est incondition-
nellement stable pour β ≥ 1
4 et précis au second ordre pour γ = 1
2
25
n
Dans un schéma de Newmark, nous pouvons utiliser les extensions pour dns et ḋs pour exprimer
n
l’accélération d̈s en fonction du déplacement dns telle que :
n 1
d̈s = dn − η n
β∆t2 s
où
1 n−1 1 − 2β n−1
ηn = (dn−1 + ∆tḋs ) + d̈s
β∆t2 s 2β
De cette manière, le système algébrique 1.64 peut être écrit en fonction du déplacement uni-
quement :
1
Ms dns + Ks (dns ) = F ns + Ms η n (1.66)
β∆t2
Ainsi, à chaque pas de temps, nous résolvons un système nonlinéaire du vecteur déplacement
dns . Pour traiter cette non-linéarité, la méthode de Newton semble la plus appropriée pour la
résolution et elle est fréquemment utilisée.
Remarques :
26
Chapitre 2
Ce chapitre a fait l’objet d’une publication sous le titre "An immersed boundary method for fluid
flows around rigid objects." International Journal Numerical Methods in Fluids, 75(1) :63–80
(2014).
La méthode des domaines fictifs est souvent utilisée dans la résolution numérique d’équations
aux dérivées partielles dans un domaine avec des obstacles mobiles. Elle a été étudiée notam-
ment par : Peskin (2002), Glowinski et al. (2000), Ilinca et Hétu (2011), Su et al. (2007) et
beaucoup d’autres, au cours des dernières années.
En particulier, cette méthode est applicable pour l’interaction fluide-structure. En effet, le
problème à résoudre est étendu du domaine fluide au domaine total fluide-solide, en prolon-
geant le champ de vitesse au solide. Un problème fluide est ainsi résolu sur tout le domaine
fictif fluide-solide. De cette manière, la condition de rigidité d’obstacles n’est plus explicite
et il convient de la rajouter dans la formulation faible. Le calcul du champ de vitesse se fait
sur le domaine total, qui a une géométrie simple et fixe. Le maillage est, donc, fixe et ne doit
être calculé qu’une seule fois (Valette et al. (2007)), ce qui n’est pas le cas pour la résolution
des problèmes d’écoulement de fluide avec un obstacle (corps rigide) mobile. Ces derniers pro-
blèmes peuvent induire, d’une part, un déplacement de maillages à chaque pas de temps de
calcul ; d’autre part, des difficultés d’utilisation de maillages réguliers puisque le domaine à
mailler peut être de géométrie complexe (Valette et al. (2007)).
27
La méthode des domaines fictifs a été introduite par Peskin (2002) qui l’a appliquée avec
succès pour simuler les valves cardiaques ainsi que dans beaucoup de problèmes du domaine
biomédical. Peskin a considéré que la géométrie complexe de la valve dans le maillage consi-
déré cartésien peut être simulée par la génération d’un champ de force externe afin d’imiter
la frontière immergée.
Dans le même cadre, Van Loon et al. (2006) ont proposé une méthode de domaine fictif
étendue avec un algorithme local d’adaptation de maillage pour fournir la flexibilité nécessaire
à l’égard du mouvement et des déformations de la valve et assurer la capacité d’imposer la
pression exercée par le solide sur le fluide.
Su et al. (2007) ont employé dans leur formulation un mélange de variables eulériennes et
lagrangiennes, où la frontière solide est représentée en description lagrangienne et exerce des
forces dans le domaine fluide eulérien. Les interactions entre les variables lagrangiennes du
solide et les variables eulériennes du fluide définies sur le maillage eulérien fixe sont données
par une simple fonction delta discrétisée. Les forces de frontière sont distribuées au maillage
eulérien en utilisant la fonction delta discrétisée. Si on modélise la position de la frontière
solide en description lagrangienne par X(s) pour 0 ≤ s ≤ Lf ront , où s est le paramètre
de la configuration de référence de la frontière, alors la vitesse du fluide sur la frontière est
U (X(s),t). On obtient le champ de force eulérien dans tout le domaine solide-fluide par la
relation :
Z Lf ront
f (x,t) = F (X(s),t)δ(x − X(s))ds
0
Ainsi, la force agissant sur le fluide f est due à la force F de la frontière immergée (domaine
solide) qui sera calculée, d’abord, à partir de la relation entre les vitesses lagrangienne du
solide et eulérienne du fluide.
Noor et al. (2009) ont développé aussi la méthode des domaines fictifs mais en utilisant la
méthode de forçage direct (direct-forcing method). Une force virtuelle est ajoutée à l’équation
de Navier-Stokes. La frontière immergée (immersed boundary) est modélisée par une force
singulière qui est incorporée dans un terme de forçage, f ∗ , dans les équations de Navier-
Stokes. L’interaction entre le fluide et la frontière immergée est liée, grâce à la diffusion de la
force singulière du maillage lagrangien au cartésien et l’interpolation de la vitesse du maillage
cartésien au lagrangien en utilisant une fonction discrétisée delta. Cette méthode utilise un
terme de forçage déterminé par la différence entre les vitesses interpolées sur les points de la
frontière et les vitesses limites souhaitées. Elle est également connue sous le nom de méthode
de forçage discret, car le forçage peut être soit explicitement ou implicitement appliqué aux
équations de Navier-Stokes.
Dans leur formulation, Noor, Chern et Horng ajoutent une force virtuelle pour les équations
28
de Navier-Stokes pour manifester l’interaction entre le fluide et le solide :
∂u 1
∂t + u · ∇u − Re ∇· 2γ̇(u) − pI = f∗
∇·u = 0
où f ∗ est le terme de la force virtuelle adopté pour cette méthode. Cette force sert à immobiliser
ou déplacer le solide. Pour satisfaire la condition aux limites du non-glissement appliquée aux
solides, la force agissant sur le solide devrait s’assurer que la vitesse du fluide u est égale à la
vitesse du solide us au (n + 1)-ième pas de temps. De là, la force virtuelle est proportionnelle
à la différence entre la vitesse du solide au (n + 1)-ième pas de temps et la vitesse du fluide
au n-ième pas de temps. Cette force existe sur le corps solide et elle est nulle ailleurs, de telle
façon que :
un+1 − un
f ∗(n+1) (x) = η(x)
s
∆t
Le solide est identifié par une fonction η nommée fonction volume du solide (VOS). η désigne
une fraction de solides au sein d’une cellule. Elle est égale à 1 et 0 respectivement pour les
cellules des solides et des fluides, et égale à une fraction sur les cellules de la frontière du
solide-fluide. Le fluide occupe virtuellement le domaine du solide avec un champ de vitesse
bien déterminé. Le champ de pression dans le solide obéit à la même équation régissant la
pression extérieure.
Huang et Sung (2009) ont aussi adopté la méthode des domaines fictifs pour la simulation
de l’interaction d’un fluide avec un solide déformable. Dans leur formulation, ils prennent
en compte à la fois la dynamique du fluide et du solide et déterminent avec précision le
déplacement libre de la frontière du solide flexible. Un maillage lagrangien est utilisé pour
discrétiser le domaine solide et l’équation du mouvement du solide est déterminée par la
méthode de la dérivation de l’énergie et est résolue par une méthode itérative. L’interaction
fluide-structure est formulée par la méthode de forçage direct comme dans Noor et al. (2009),
c’est-à-dire par l’ajout d’une force dans les équations de Navier-Stokes. Cette force est obtenue
sur le maillage lagrangien du solide par l’utilisation de l’équation de mouvement du solide et
elle se diffuse dans tout le domaine du fluide par l’intermédiaire d’une fonction delta de Dirac.
Avec cette méthode, Huang et Sung ont simulé beaucoup de problèmes en biologie et ont
obtenu de bons résultats.
Plusieurs autres auteurs ont adopté la méthode des domaines fictifs parmi lesquels Van Loon
et al. (2004), Baaijens (2001) et Van Loon et al. (2007). Dans sa formulation, Baaijens (2001)
a étudié l’interaction fluide-structure par la méthode des domaines fictifs en modélisant le
fluide par les équations de Navier-Stokes en description eulérienne et le solide en description
Lagrangienne par un modèle néo-hookéen. Un multiplicateur de Lagrange est associé à la
formulation faible des domaines fluide et solide. Le couplage entre le fluide et le solide est
assuré par l’intermédiaire d’un multiplicateur de Lagrange pour assurer la contrainte us −u = 0
29
sur la frontière immergée. Par conséquent, il y aura un terme de plus dans les équations de
Navier-Stokes comme suit :
Z n
−
Z
u u 2
ρ w+ γ̇(u) : γ̇(w) − p∇·w + (u · ∇u) · w dv = f w dv
∆t Re
Ω Ω
Z
γ · w dS
+
Γ1
Z
q∇ · u dv
= 0
Ω
Z
λ · (us − u) ds = 0
Γ1
où γ est le multiplicateur de Lagrange et λ est la fonction poids associée. γ peut être interprétée
comme étant la force de surface exercée sur le fluide et le solide à la frontière Γ1 pour maintenir
le couplage entre le fluide et le solide. La résolution de ce problème est réalisée d’une manière
itérative par la méthode de Newton en tenant compte du couplage assuré par la dernière
équation du système précédent à chaque itération de la boucle de Newton.
Dans cette méthode, la position de la frontière Γ n’est pas précise, et donc, la convergence
quadratique du schéma du Newton n’est pas obtenue (Baaijens (2001)).
Tout comme Glowinski et al. (1994), Baaijens (2001) a proposé une méthode des domaines
fictifs basée sur un multiplicateur de Lagrange. Cette méthode a été étendue afin d’obtenir
des solutions par éléments finis pour le déplacement des corps rigides par l’utilisation de
maillages quelconques (le lecteur pourra se référer à Glowinski et al. (1999)). La méthode
consiste à imposer faiblement des contraintes cinématiques au fluide dans les régions occupées
par les corps rigides et ceci, par l’intermédiaire des multiplicateurs de Lagrange. Un problème
d’écoulement sur l’ensemble du domaine est, ainsi, résolu. D’autres améliorations de cette
méthode sont développées dans Glowinski et al. (2000).
On présente, maintenant, la méthode des domaines fictifs que nous avons adoptée pour le
calcul numérique d’interaction fluide-structure. L’approche proposée repose sur une méthode
de calcul des champs de vitesse et de pression du fluide en présence d’un objet rigide placé
comme obstacle mais qui est à priori fictif (Ilinca et Hétu (2011)). L’interaction de la frontière
du solide avec le fluide est assurée par des forces locales appliquées au fluide sur tous les points
occupés par le solide. Ces forces locales sont dues à l’application des contraintes cinématiques.
Ces contraintes sont imposées par l’application des conditions de Dirichlet. De ce fait, deux cas
de conditions de Dirichlet se présentent selon le cas d’étude. Si l’étude est faite sur un solide
immobile, on applique dans ce cas une condition de Dirichlet présentant une vitesse nulle sur
tous les points occupés par le solide. Par contre, si le solide est mobile, on applique comme
conditions, celles de Dirichlet sur tous les nœuds du solide illustrées par la vitesse définie par
le mouvement de ce dernier. Ainsi, on prolonge le champ de vitesse au solide et on résout
30
le problème sur tout le domaine fluide-solide, tout en appliquant, selon le cas, les conditions
aux limites nécessaires (condition d’adhérence) sur le solide. Pour assurer l’application de
cette condition et pour bien délimiter le volume ainsi que la frontière du solide, nous avons
développé une méthode de découpage du maillage présentée dans la section suivante.
Les maillages structurés permettent une discrétisation simple et précise des équations, mais
ils ne sont pas, en général, adaptés au traitement des cas industriels mettant en jeu des
obstacles ou des interfaces de formes complexes. Afin d’améliorer la prise en compte d’interfaces
complexes et assurer la précision des résultats donnés par la méthode des domaines fictifs, nous
avons développé une méthode de découpage de maillage.
Considérons le cas d’une fonction f dont l’isovaleur zéro définie un profil NACA. Un agran-
dissement du maillage autour de l’aile est présenté dans la figure 2.1. En parcourant toutes
les arêtes du maillage, on cherche à l’aide d’une méthode de bissection une intersection des
arêtes du maillage considéré avec l’isovaleur zéro de la fonction f ainsi donnée. Une fois un
point d’intersection trouvé, on ajoute un nœud au maillage en ce point c’est-à-dire sur l’arête
contenant le point d’intersection et on construit une coquille autour de l’arête du maillage mis
en jeu. On raffine l’arête en deux et on coupe les éléments voisins. On répète la dernière étape
sur toutes les arêtes contenant les points d’intersection du maillage avec la fonction f . L’ajout
de nœuds sur l’interface ainsi que le raffinage des arêtes sont représentés dans la figure 2.2 et
le maillage ainsi construit, dans ce cas, est représenté dans la figure 2.3.
Notons que la technique est la même pour le découpage d’autres formes, il faut juste donner
la fonction f définissant la géométrie de l’obstacle qu’on veut découper.
Remarque
Lors du découpage du maillage, il faut prendre certaines précautions pour assurer un bon fonc-
tionnement. En effet, dans la méthode de bissection utilisée, nous avons ajouté une condition
sur la recherche des points d’intersection. Cette dernière consiste à négliger tous les points
d’intersection dont la distance par rapport aux extrémités des arêtes est de l’ordre de 10−6
31
∂Ωr
∂Ωr
∂Ωr
32
multiplié par la longueur de l’arête en question. Cela évite de créer des éléments de taille très
petite. En effet, si l’intersection de la structure rigide avec une arête du maillage est située à
proximité d’un nœud existant, un élément avec au moins une très petite arête peut être créé.
De plus, nous avons envisagé une autre condition qui consiste à ne pas avoir des éléments
d’aires très faibles après raffinement d’arêtes. Ceci est réalisé en calculant l’aire de tous les
éléments ainsi construits après découpage du maillage (c’est-à-dire après la modification de la
coquille construite au préalable) et en éliminant les éléments dont l’aire est inférieure à 10−6
multiplié par l’aire totale de la coquille. En effet, ce choix a été effectué suite à plusieurs tests
effectués sur différents maillages ; nous avons testé le cas où on élimine les éléments dont l’aire
est de l’ordre de 10−7 ou 10−8 ou 10−9 fois l’aire totale de la coquille et nous avons constaté
que le maillage comporte des éléments d’aires très faibles. Par contre, en éliminant les éléments
dont l’aire est inférieure à 10−6 fois l’aire totale de la coquille, on aboutit à de bons maillages.
Comme troisième condition, nous avons vérifié l’existence d’éléments superposés ou retournés
(d’aire négative). Ces éléments peuvent survenir suite au raffinement d’arêtes après l’ajout de
nœuds. En gros, on s’assure que le maillage est valide.
Ces trois conditions ont été ajoutées après beaucoup de tests sur différentes fonctions définies
et sur différents maillages.
Malgré l’introduction des critères présentés à la remarque précédente, des éléments de mauvaise
qualité peuvent appraître dans le maillage.
Rappelons que pour un élément tétraédrique K avec des longueurs d’arêtes li , i = 1,..,6 et un
volume VK , la qualité des éléments peut être définie comme :
12 VK
Q3D
K = √ √ .
2 l1 l2 l3 l4 l5 l6
Il peut être facilement vérifié que cette quantité est égale à 1 pour un tétraèdre régulier dont
toutes les arêtes sont de même longueur. Une définition similaire existe pour un triangle dans
le cas 2D :
27(l1 + l2 − l3 )(l1 + l3 − l2 )(l2 + l3 − l1 )
Q2D
K = .
(l1 + l2 + l3 )3
ce qui donne évidemment 1 pour un triangle équilatéral.
Améliorer la qualité des éléments signifie donc, essayer de se rapprocher le plus possible des
triangles équilatéraux ou des tétraèdres réguliers. La qualité d’un maillage (ou d’un sous-
maillage) est définie comme étant la qualité minimale de ses éléments.
33
Pour améliorer la qualité des éléments, nous introduisons la notion de retournement d’arêtes.
En effet, notre méthode d’adaptation de maillage (voir Bois et al. (2012a,b)) est basée sur
le calcul précis d’un estimateur d’erreur conduisant à des opérations locales sur le maillage :
raffinement d’arêtes, retournement d’arêtes, déraffinement et déplacement de sommets. Dans
le contexte des méthodes des domaines fictifs, l’estimation de l’erreur et le remaillage de
l’ensemble du domaine à chaque pas de temps serait probablement prohibitif.
Un aspect fondamental de notre approche est le retournement d’arêtes, qui est illustré dans le
cas 2D à la figure 2.4. Pour une arête donnée, nous construisons un sous-maillage contenant
ses deux triangles adjacents. Des deux sous-maillages illustrés à la figure 2.4, nous choisissons
celui dont les éléments sont de plus haute qualité. Notons qu’il existe des situations où le
retournement d’arêtes ne peut pas être fait. Ce sera le cas, si l’aire des sous-maillages change
par exemple.
Le cas tridimensionnel est beaucoup plus complexe et mérite quelques commentaires. Pour
une arête donnée e avec deux noeuds a et b (voir figure 2.5), on construit premièrement le
sous-maillage Ωe contenant tous les tétraèdres partageant l’arête en question. Nous examinons
ensuite, tous les nœuds dans ce sous-maillage sauf a et b :
Ne = {xi ∈ Ωe |xi 6= a et b}
Ceci est illustré à la figure 2.5 où l’arête e est supposée perpendiculaire au plan de la figure
(avec des nœuds a et b de chaque côté, superposés sur la figure). Afin de simplifier l’expli-
cation, nous considérons dans la figure le cas particulier où tous les nœuds Ne se trouvent
dans le même plan. Ce n’est pas le cas en général, mais la procédure est la même. À partir
de l’un des nœuds dans Ne noté X dans la figure, un maillage triangulaire en forme d’étoile
est créé par l’ajout de nouvelles arêtes entre ce nœud et tous les autres (lignes en pointillés
sur la figure 2.5). Dans le cas 2D, cette opération est évidente, car il y a seulement deux
nœuds à gauche et une seule arête peut être ajouté. En 3D, à partir de chacun des nouveaux
triangles, deux tétraèdres sont, ensuite, construits en les reliant respectivement à a (au-dessus
du plan) et b (en dessous du plan). Dans le cas illustré, six tétraèdres seraient créés. Différents
maillages peuvent être obtenus selon le choix du premier noeud X. Toutes ces différentes
configurations sont explorées et la qualité moyenne de chaque sous-maillage est calculée. La
configuration ayant la plus haute qualité d’éléments est conservée. Toutes les précautions sont
prises pour éviter la création d’éléments non valides (avec volume négatif par exemple). Un
critère important est que les volumes initial et final des sous-maillages devraient être les mêmes.
Le retournement d’arêtes a l’avantage d’améliorer la qualité des éléments sans créer de nou-
veaux nœuds. Les nouvelles arêtes ainsi créées sont à leur tour découpées si nécessaire créant
ainsi un peu plus de nœuds.
Toute la procédure est illustrée à la figure 2.6 pour un objet rigide circulaire (cylindrique).
34
b
a b a b
X X
Les éléments avec des arêtes situées à l’interface de l’objet rigide sont illustrés en rouge,
avant et après le retournement d’arêtes, montrant clairement une amélioration de leurs formes
géométriques.
Ainsi, l’algorithme complet nécessite un maillage de base de départ (constant) noté MBase et
peut être décrit comme suit :
Algorithme
1. Pour un pas de temps donné tn :
1.1 En utilisant la position connue de l’objet rigide au temps tn , effectuer le découpage
du maillage de base Mbase . On notera le maillage, ainsi, obtenu par Mtn .
1.2 Appliquer le retournement d’arêtes et un nouveau découpage pour améliorer la
qualité des éléments du maillage Mtn .
35
Y
VU Z X
$1
$1
Figure 2.6 – Éléments autour d’un cercle avant et après retournement d’arêtes
1.3 Réinterpoler les champs de vitesses et de pression du pas de temps précédent (un−1 ,
un−2 et pn−1 ) sur le nouveau maillage Mtn (notons que cette étape peut se faire
seulement sur les nouveaux noeuds),
1.4 Appliquer les bonnes conditions aux limites au temps tn . En particulier, les condi-
tions aux limites de Dirichlet sont imposées à l’intérieur de l’objet rigide et à tous
les nœuds du bord, y compris les nouveaux nœuds introduits par la méthode de
découpage.
1.5 Résoudre les équations de Navier-Stokes et calculer la solution (un ,pn ),
1.6 Retour au maillage initial MBase .
2. Passer au pas de temps suivant.
On notera que, à chaque pas de temps, la technique de découpage (voir étape 1.1) est effectuée
à partir d’un maillage de base Mbase . Dans le cas contraire, le maillage obtenu serait de plus
en plus fin à chaque pas de temps (et éventuellement complètement déformé) dans toute la
région où l’objet se déplace. Évidemment, si l’objet rigide est immobile, le découpage et le
36
retournement d’arêtes sont effectués une seule fois. Une fois le maillage Mtn créé (étape 1.3),
les solutions de vitesse et de pression des pas de temps précédents (un−1 , un−2 et pn−1 ) doivent
être réinterpolées sur le nouveau maillage. Notons que la plupart des nœuds sont communs
aux deux maillages. Par conséquent, l’interpolation à ces nœuds est triviale et sans erreur
numérique additionnelle.
Pour les nouveaux nœuds, l’interpolation est effectuée d’une manière classique. En utilisant
les coordonnées d’un nouveau nœud dans le maillage Mtn , on détermine dans quel élément
de l’ancien maillage il est situé, puis il suffit d’utiliser les fonctions de base de Lagrange pour
interpoler correctement.
L’application de cet algorithme et les résultats numériques obtenus par la méthode des do-
maines fictifs (qu’on notera par Méthode IB « Immersed Boundary » par la suite) seront
détaillés dans la prochaine section.
37
immédiate. Les coefficients de traînée et de portance sont définis par :
F1 F2
CD = 1 2 2
et CL = 1 2 2
2 ρU πR 2 ρU πR
24
1 + 0.15 Re0.687 (2.2)
CD =
Re
La géométrie utilisée est illustrée dans la figure 3.9 qui n’est pas à l’échelle. On impose à
l’entrée et aux frontières latérales un écoulement uniforme u = (1,0,0) et à la sortie on a :
σ · n = 0. Une condition de non-glissement a été appliquée sur la sphère.
Sortie σ · n =
(0,0,0)
u = (1,0,0)
u = (0,0,0)
35
u = (1,0,0)
(-10,10,-10)
20
x2 x
1
20
(-10,-10,-10) (-10,-10,10) x3
Entrée
38
Les simulations ont été faites en commençant par un maillage grossier mais sans retournement
d’arêtes. La valeur calculée du coefficient de traînée CD = 56 était très loin de la valeur
empirique CD = 247 donnée par la relation (2.2). Ceci peut être expliqué par la mauvaise
qualité des éléments dans la région autour de la sphère. Le retournement d’arête donne une
meilleure qualité d’éléments autour de la sphère et un calcul plus précis du coefficient de
traînée, soit CD = 232. Ce simple exemple montre que la méthodologie globale fonctionne
bien mais montre aussi la nécessité du retournement d’arêtes. En poursuivant le calcul sur le
maillage plus fin (avec retournement d’arêtes), nous avons obtenu une valeur plus précise du
coefficient de traînée, soit CD = 250. Des résultats similaires ont été obtenus pour Re = 0.5
tel qu’illustré dans le tableau 2.1.
Le même problème a été aussi résolu en utilisant la méthode avec un véritable obstacle (mé-
thode BF). C’est-à-dire que l’obstacle est effectivement représenté dans le maillage et n’est
donc pas fictif. C’est cette approche qui nous servira de solution de référence. On a considéré
2 maillages constitués respectivement de 25 470 et 128 161 éléments (un nombre d’éléments
similaire à celui qu’on a utilisé pour la méthode IB). Les résultats sont aussi présentés dans le
tableau 2.1 et montrent que les valeurs calculées des coefficients de traînée sont très compa-
rables à celles données par la formule 2.2. Le raffinement de maillage minimise l’erreur pour
les 2 méthodes considérées. La méthode IB produit, ainsi, des valeurs de traînée très simi-
laires en précision à celles obtenues par la méthode classique BF. Plus de détails sont élaborés
dans Jendoubi et al. (2014).
Cas 2D
Nous étudions maintenant le cas d’un écoulement autour d’un cylindre. On considère un cy-
lindre de diamètre 1 et centré à l’origine dans le domaine : −10 ≤ x ≤ 25 et −10 ≤ y ≤ 10 pour
éviter l’influence indue des conditions aux limites en entrée et sortie. À l’entrée, on impose
u = (1,0) de même que sur les frontières inférieures et supérieures, et on impose une condition
naturelle nulle en sortie, soit σ · n = (0,0). De plus, on applique une condition de Dirichlet
sur le cylindre, soit u = (0,0), et ceci comme l’indique la figure 2.8.
Le maillage considéré (voir la figure 2.9) est constitué de 12 704 éléments de Taylor-Hood,
et présente un découpage du cylindre centré à l’origine pour bien délimiter le domaine de
l’obstacle et pouvoir y appliquer la condition d’adhérence. On résout le problème de Navier-
39
Stokes instationnaire pour un pas de temps égal à 0.1 et pour différentes valeurs du nombre
de Reynolds : Re = 1, Re = 45 et Re = 100. On s’intéresse plus précisément au dernier cas où
Re = 100. Une fois les effets transitoires disparus, on présente à la figure 2.10 les vecteurs de la
composante Ux de la vitesse autour du cylindre. On remarque le développement de plusieurs
zones de recirculation (tourbillons).
Pour ce problème, la condition d’adhérence est bien respectée (vitesse nulle à l’intérieur du
cylindre (voir figure 2.11)), c’est-à-dire qu’aucun vecteur vitesse ne se développe à l’intérieur
du cylindre, et ceci est dû à notre ajout de nœuds sur l’interface (méthode de découpage).
On conclut qu’il existe une instabilité développée à l’arrière du cylindre qui est, en fait, l’allée
de Von Karman classique. Ceci est conforme aux résultats de Fortin et Garon (2013) et En-
gelman et Jamnia (1990). En effet, Re = 100 est supérieur au nombre critique Re ≈ 45 pour
que l’écoulement devienne instationnaire, ce qui justifie l’instabilité du problème et le déve-
loppement de l’allée de von Karman comme l’indique la figure 2.11.
On trace, à la figure 2.12, l’évolution de la composante Uy du vecteur vitesse en un point
A(3,0) du domaine situé à l’arrière du cylindre (voir la figure 2.8) dans les 400 derniers pas de
temps de notre simulation. Par ailleurs, on montre que l’écoulement est périodique de période
6.3 en accord avec Fortin et Garon (2013). Puisque le pas de temps de notre calcul est de 0.1
alors, 63 pas de temps sont utilisés pour couvrir une période d’oscillations, ce qui est largement
suffisant.
u = (1,0)
(25,10)
u = (0,0)
u = (1,0) σ·n=0
A(3,0)
(-10,-10)
u = (1,0)
Figure 2.8 – Conditions aux limites pour l’étude de l’écoulement autour d’un cylindre.
Cas 3D
Un écoulement autour d’un cylindre de section transversale circulaire est maintenant considéré.
La géométrie et les conditions aux limites de ce problème sont illustrées dans la figure 2.13.
40
Y
VU Z X
Tout comme dans Schäfer et Turek (1996), on applique une condition de Dirichlet sur le cy-
lindre et sur les frontières inférieures et supérieures du domaine, soit u = (0,0,0). La condition
aux limites de Dirichlet d’entrée est donnée par :
72
u(0, x2 , x3 ) = x2 x3 (H − x2 )(H − x3 ), 0, 0
H4
ce qui donne une vitesse moyenne U = 2. La hauteur et la largeur du canal sont de H = 4,1
unités, et le diamètre du cylindre est fixé à l’unité (D = 1). La fonction distance est, donc,
donnée par ϕ(x1 , x2 , x3 , t) = x21 + x22 − 1/4. Le nombre de Reynolds est basé sur le diamètre
du cylindre et la vitesse moyenne à l’entrée. Nous considérons les écoulements à Re = 20 à
des fins de comparaison.
41
Figure 2.10 – Vitesse Ux autour du cylindre (Re = 100).
42
Figure 2.12 – Comportement de la composante Uy en un point derrière le cylindre.
Sortie σ · n =
(0,0,0)
u = (0,0,0)
u = (0,0,0)
19.5
u = (0,0,0)
1.6
xe
(0,H,0)
1.0
1.5 xa 4.5
4.1
x2 x
4.1 1
(0,0,0) (0,0,H) x3
Entrée
43
sont situés de part et d’autre du cylindre avec les coordonnées xa = (4.5,2.0,2.05) et xe =
(5.5,2.0,2.05)) tel qu’illustré à la figure 2.13.
Dans cette section, on compare aussi les résultats obtenus avec la méthode IB, la méthode BF
et les résultats numériques de Schäfer et Turek (1996) qui étaient basés, aussi, sur un maillage
BF. Le tableau 2.2 résume les différents résultats. Comme on peut le constater, les méthodes
considérées conduisent à des résultats très similaires.
Nous avons montré dans ces deux premières sections les résultats obtenus au cas où les ob-
jets seraient statiques et auraient permis de vérifier notre approche. Dans ce qui suit, nous
considérons les objets mobiles.
Dans un premier exemple des objets en mouvement, l’objet rigide rotatif est un quadrifolium
(également connu sous le nom de trèfle à quatre feuilles). L’objet peut être exprimé avec des
coordonnées polaires (au temps t = 0 d’où l’indice supérieur 0 ) par :
0
x3 (θ) = (α + β cos(4θ)) cos(θ),
θ ∈ (0,2π)
0
x2 (θ) = (α + β cos(4θ)) sin(θ)
où α = 1
5 et β = 10 .
1
En coordonnées cartésiennes, ceci correspond à :
5
q
0 2 2
(x02 )2 (x03 )2 0 2 0 0 (2.3)
− (x2 ) + (x3 ) 2 2
0.3 − (x2 ) + (x3 ) = 0
4
qui est présenté à la figure 2.14. La largeur (dans la direction x1 ) du quadrifolium est fixée à
0,04 unités.
Tel qu’illustré à la figure 2.15, l’objet est en rotation dans le sens des aiguilles d’une montre
autour de l’axe x1 avec une vitesse angulaire ω. Pour un temps donné t et un point (x2 ,x3 )
donné, on doit déterminer si (x2 ,x3 ) est à l’intérieur ou non de l’obstacle. La frontière du
quadrifolium est donnée par une fonction distance définie par :
5
q
0 2 2
(x02 )2 (x03 )2 0 2 0 0
ϕ(x1 , x2 , x3 , t) = − (x2 ) + (x3 ) 2 2
0.3 − (x2 ) + (x3 ) = 0
4
44
Figure 2.14 – Quadrifolium ou trèfle à quatre feuilles
où
x03 = x3 cos(ωt) − x2 sin(ωt)
(2.4)
x02 = x3 sin(ωt) + x2 cos(ωt)
Cette dernière expression est une rotation dans le sens contraire des aiguilles d’une montre du
point (x3 , x2 ) pour le remettre à sa position d’origine. L’équation (2.3) peut être utilisée pour
déterminer la position de la frontière du quadrifolium. À partir de (2.4), on obtient facilement :
Ce problème peut être considéré comme une approximation d’une hélice en rotation et est
présenté pour illustrer que la méthode fonctionne bien dans le cas des objets en mouvement.
Cette hélice est en mouvement de rotation dans un canal de diamètre D = 1. Une analyse
dimentionnelle du problème montre qu’il existe deux nombres adimentionels qui sont : le
ρU D ωD
nombre de Reynolds Re = et le nombre de Strouhal St = , tous les deux définis
µ U
en utilisant la vitesse moyenne à l’entrée U , le diamètre D du cylindre et la vitesse angulaire
ω de rotation de l’objet. Le nombre de Strouhal représente le rapport du temps d’advection
et du temps caractéristique de l’instationnarité, c’est à dire qu’il représente le rapport entre
45
les forces inertielles dues à l’instationnarité de l’écoulement et les forces d’inertie dues aux
variations de la vitesse, point par point, dans l’écoulement.
Bien que le solide soit en rotation, ce problème peut encore être résolu en utilisant à la fois la
méthode BF et la méthode IB si les systèmes de coordonnées appropriés sont bien choisis. La
géométrie et les conditions aux limites sont illustrées pour les deux méthodes à la figure 2.15.
Méthode IB
Dans ce cas, on utilise un système de coordonnées statique. La vitesse à l’entrée est donnée
par
2 2 2 2
(2.5)
uentrée = R − x2 − x 3 , 0, 0
R2
où R est le rayon du cylindre (R = 0.5), conduisant à une vitesse moyenne U = 1 et St = ω.
Une condition de non-glissement est imposée sur le cylindre (ucyl = (0,0,0)) et l’objet rigide
est en rotation dans le sens des aiguilles d’une montre avec une vitesse angulaire St :
ur = St(0, −x3 , x2 )
Méthode BF
Dans ce cas, on considère un système de coordonnées en rotation avec l’hélice. L’objet rigide
n’est pas en mouvement (ur = (0,0,0)), mais c’est le cylindre qui est en rotation dans le sens
contraire des aiguilles d’une montre
2 2 2 2
uentrée = R − x2 − x3 , 0, 0 − St(0, −x3 , x2 )
R2
Notons que puisqu’on est dans un système de coordonnées en rotation, les termes de force
de Coriolis doivent être ajoutés aux équations de mouvement (Olshanskii et al. (2010), Petcu
et al. (2004)) comme suit
3 un + 2un−1 − un−2 · ∇un − 2 ∇ · γ̇(un ) + ∇pn + 2 (w × un ) =
2∆t Re
1 n−1 n−2
f+ 4u −u − w × (w × r)
2∆t
où w = St(1,0,0) est l’axe du vecteur rotation et r = (x1 , x2 , x3 ) est le vecteur position.
Les calculs ont été effectués pour Re = 50 et St = 10 à la fois pour les méthodes IB et
46
x2
BF : ucyl = −St(0, −x3 ,x2 )
x1 Sortie
x3 IB : ucyl = (0,0,0)
a2
BF : ur = (0,0,0) σ · n = (0,0,0)
IB : ur = St(0, −x3 ,x2 )
xc
Entrée
a1
1
2 7
BF : ucyl = −St(0, −x3 ,x2 )
IB : ucyl = (0,0,0)
BF. La méthode IB a été appliquée sur un maillage tétraédrique non structuré ayant 115 372
éléments (20 038 sommets et 136 629 arêtes). La méthode BF a également été appliquée sur
un maillage tétraédrique non structuré ayant 104 001 éléments (18 607 sommets et 125 033
arêtes). Les sections transversales des deux maillages, sur un plan passant par l’hélice, sont
présentées dans la figure 2.16. Notez que ce ne sont pas des maillages triangulaires réguliers
mais seulement les intersections de maillages tridimensionnels avec le plan. Le maillage de la
méthode IB est illustré après retournement d’arêtes et est un peu plus raffiné dans la région
où se trouve l’hélice. Rappelons également que tous les nœuds à l’intérieur de l’hélice seront
éliminés des systèmes linéaires, réduisant ainsi la charge de calcul.
Le pas de temps a été fixé à π/180. Par conséquent, une rotation complète est effectuée tous
les 36 pas de temps. Pour la méthode IB, cela signifie que 36 maillages doivent être construits
à partir du maillage de base et, ensuite, réutilisés pour les pas de temps ultérieures. Une fois
les effets transitoires disparus, on obtient la vitesse longitudinale u1 tracée le long de l’axe
x1 dans la figure 2.17. Nous avons également présenté à la figure 2.18, la comparaison de la
pression le long de l’axe du cylindre pour les deux méthodes considérées.
Il est également possible de comparer les vitesses transversales (u2 ) en utilisant la relation
(2.6), mais pas sur l’axe des x1 , car elles sont nulles. La figure 2.19 montre la comparaison
entre les vitesses transversales tracées entre a1 = (−2,0.2,0) et a2 = (5,0.2,0). Dans tous les
cas, l’accord est très satisfaisant. Les solutions numériques obtenues avec les deux méthodes
sont similaires. C’est, aussi, le cas pour les coefficients de traînée et de portance calculés, qui
sont donnés dans le tableau 2.3.
47
Figure 2.16 – Maillages de l’hélice pour la méthode BF (à gauche) et la méthode IB (à droite)
u1 Méthode IB
Méthode BF
x1
u1 Méthode IB
Méthode BF
1.5 2.5 x1
Vue agrandie dans l’intervalle [1.5, 2.5] proche de l’hélice
Figure 2.17 – Vitesses longitudinales au long de l’axe x1 : comparaison entre les méthodes
BF et IB
48
Méthode IB
P Méthode BF
x1
P Méthode IB
Méthode BF
1.5 2.5 x1
Vue agrandie dans l’intervalle [1.5, 2.5] proche de l’hélice
u3
Méthode IB
Méthode BF
x1
u3
Méthode IB
Méthode BF
1.5 2.5 x1
Vue agrandie dans l’intervalle [1.5, 2.5] proche de l’hélice
49
profil est décrit à la figure 2.21. La position du point xc (t) = (xc1 (t), xc2 (t)) dans la figure 3.6
et l’angle θ(t) du profil sont donnés par (voir Cori (2011)) :
θ(t) = θ0 sin(ωt)
h(t) = h0 sin(ωt)
xc1 (t) = 0
c
x2 (t) = h(t)
uy = 0
30
10
u = (1,0) xc 20 σ · n = (0,0,0)
1
3 1
y
x
uy = 0
Figure 2.20 – Domaine de définition pour la modélisation de l’écoulement autour d’un profil
oscillant
Le mouvement est composé d’une rotation avec un angle θ(t), suivi par une translation (bat-
tement) d’amplitude (xc1 (t), xc2 (t)). Ainsi, à un temps donné t, sa frontière est donnée par la
fonction distance suivante :
" r
0.15c 2
0
x01 − 1/3
x1 − 1/3
ϕ(x1 , x2 , t) = (x02 )2
− 0.2969 − 0.1260
0.2 c c
0 2 0 3 0 4 #2
x1 − 1/3 x1 − 1/3 x1 − 1/3
−0.3537 + 0.2843 − 0.1015 =0
c c c
50
Axe x
de rotation
p
h0
θ0
U∞
Battement
où, pour revenir à la position initiale, nous enlevons la translation et tournons avec un angle
−θ.
x01 = (x1 − xc1 (t)) cos(θ(t)) + (x2 − xc2 (t)) sin(θ(t))
(2.8)
x02 = −(x1 − xc1 (t)) sin(θ(t)) + (x2 − xc2 (t)) cos(θ(t))
En inversant cette dernière expression, nous obtenons
ce qui donne la position d’un point sur le profil à tout moment. En prenant la dérivée, on
trouve l’expression de la vitesse qui doit être imposée sur ce point
x01 (t) = θ0 (t) −x01 sin(θ(t)) − x02 cos(θ(t))) + (xc1 )0 (t) = −θ0 (t)(x2 − xc2 ) + (xc1 )0 (t)
x02 (t) = θ0 (t) x01 cos(θ(t)) − x02 sin(θ(t))) + (xc2 )0 (t) = θ0 (t)(x1 − xc1 ) + (xc2 )0 (t)
Comme dans Cori (2011), nous considérons le cas où la longueur de la corde c = 1, l’amplitude
de battement h0 = 1, l’amplitude de la rotation est définie par θ0 = π/3 et la fréquence
d’oscillation est 0.18. Par conséquent, nous avons ω = 2π × 0.18 = 1.13 et une période
ρU∞ c
correspondant à T = 5.55. Nous considérons l’écoulement à Re = 1100 où Re = µ .
Il s’agit d’un problème difficile, car le domaine de calcul est très grand par rapport à la dimen-
sion réelle du profil oscillant. Le maillage initial (avant découpage et retournement d’arêtes)
51
est présenté à la figure 2.22 et contient 86 580 nœuds, un nombre d’éléments comparable à
celui des maillages utilisés dans Cori (2011) (120 000 nœuds ) et Kinsey et Dumas (2008) (plus
de 200 000 nœuds).
Le pas de temps a été fixé à 0,05 et il faut donc exactement 111 pas de temps pour couvrir
une oscillation complète du profil d’aile. Par conséquent, 111 maillages ont été créés pour un
seul cycle et utilisés de façon répétée.
Notons qu’à partir d’un champ de vitesse nul, les effets transitoires ont disparu après deux
périodes complètes et un écoulement bien établi a été obtenu. On définit le coefficient de
traînée comme suit :
F1
CD = 1 2
2 ρU∞ c
qui est une quantité critique dans ce type d’applications d’ingénierie. Le coefficient de portance
calculé sur une période complète est présenté à la figure 2.23 et comparé avec les résultats
de Cori (2011). Comme on peut le constater, l’accord est très satisfaisant.
Enfin, nous présentons, à la figure 2.24, les isolignes de la composante u1 de la vitesse ; derrière
le profil d’aile aux différents temps, à travers une période T . Une allée de von Kármán peut
être facilement visualisée.
52
CL Travail courant
Cori et al. 2011
t
T
Dans notre dernier test de validation, nous présentons un exemple illustrant la flexibilité de
notre méthode IB. Ce problème ne peut pas être résolu en utilisant la méthode BF. Nous
considérons une géométrie similaire à celle de la section 2.3.3 mais nous ajoutons une seconde
hélice contrarotative avec une vitesse angulaire différente. Le domaine de calcul et les condi-
tions aux limites pour ce problème sont présentés à la figure 2.25. Les calculs ont été effectués
pour Re = 50 et pour deux vitesses angulaires différentes : ω = 1 Hz pour la première hélice
et ω = 10 Hz pour la seconde. Le pas de temps a été fixé à dt = π/180.
La figure 2.27 présente la vitesse transversale u2 entre les mêmes points a1 et a2 utilisés dans
la section 2.3.3. Les vitesses transversales induites par les objets en rotation sont clairement
visibles pour un facteur d’environ 10 entre les deux hélices.
53
t = T /2
t = 3T /4
t=T
54
Sortie
a2
σ · n = (0,0,0)
Entrée
a1 1
1 x2
2 7
x1
x3
u2 Méthode IB
x1
55
2.4 Conclusion
Dans ce chapitre, nous avons présenté une méthode des domaines fictifs où un maillage de base
est modifié localement (à travers le découpage et le retournement d’arêtes) afin de s’adapter
à la frontière de la structure rigide. Notons que plus le maillage de base est fin, plus la
reconstruction de l’objet rigide est précise . Le retournement d’arêtes est utilisé pour améliorer
la qualité des éléments et améliorer, par conséquent, la précision des quantités telles que les
coefficients de traînée et de portance. Les résultats numériques montrent que la méthode
fonctionne remarquablement bien dans diverses situations, que ce soit dans le cas des structures
rigides immobiles ou/et mobiles. Plusieurs cas tests numériques présentés dans ce chapitre ont
été publiés dans Jendoubi et al. (2014).
Dans le chapitre suivant, nous présenterons une autre approche de modélisation de l’IFS où
la structure sera dans ce cas flexible. Une comparaison des deux méthodes dans le cas des
structures rigides sera faite dans le chapitre 4.
56
Chapitre 3
La principale difficulté liée à l’interaction d’un fluide et d’une structure réside dans le cou-
plage de modèles décrits de manière différentes. En effet, le fluide est classiquement décrit en
formulation eulérienne et la structure élastique en formulation lagrangienne. Le couplage de
ces différentes formulations rend l’étude de ces systèmes très complexe, tant d’un point de vue
numérique que mathématique. Les méthodes de couplage fluide-structure se divisent en deux
catégories.
La première consiste à un couplage où les maillages sont mobiles. Dans ce contexte, Donea
et al. (2004) présentent la méthode ALE (Arbitrary Lagrangian Eulerian) basée sur un com-
promis entre les descriptions lagrangienne et eulérienne. Cette méthode consiste à déformer le
domaine fluide à partir d’une configuration de référence, ce qui introduit des termes complé-
mentaires liés au mouvement du maillage de calcul, dans les équations aux dérivées partielles
du modèle fluide. Une autre avancée significative a été réalisée par Dunne (2007) qui pré-
sente une méthode eulérienne-lagrangienne. Un maillage solide lagrangien est superposé à un
maillage fluide eulérien fixe. Le principe est, alors, de pénaliser la différence de vitesse entre
l’interface solide et le fluide. Cette approche a été aussi utilisée dans Aquelet et al. (2006) et
consiste à considérer les noeuds à l’interface solide comme «esclaves» de noeuds «maîtres»
constitués de particules fluides virtuelles (voir Aquelet et al. (2006) pour plus de détails).
La deuxième catégorie de couplage est le couplage à maillages fixes. Dans cette catégorie,
une première avancée significative a été réalisée par Peskin (2002) qui permet de considérer
l’ensemble fluide/solide comme un unique milieu continu. Des multiplicateurs de lagrange sont
introduits pour calculer les forces élastiques appliquées au solide. Dans le cas particulier d’une
membrane élastique immergée dans un fluide, une autre méthode a été proposée par Cottet
et Maitre (2004, 2006). Elle consiste à représenter l’interface entre les deux milieux par une
fonction courbe de niveau «level set» qui est advectée par la vitesse du fluide. La nouveauté
57
de la méthode réside dans la prise en compte d’une partie de l’élasticité de la structure à l’aide
de cette fonction level set. Le couplage fluide-structure est ainsi réduit à une formulation
complètement eulérienne.
Dans ce travail, on a fait le choix d’une formulation ALE dans le fluide et d’une formulation
lagrangienne dans la structure. En effet, puisqu’il s’agit de modéliser par éléments finis une
structure hyperélastique en grandes déformations, il semble naturel de considérer une descrip-
tion lagrangienne (voir chapitre 1). De plus, une formulation purement eulérienne ne permet
pas de suivre efficacement les phénomènes à l’interface fluide-solide et soulève la question de
la loi de comportement à imposer à l’interface. L’approche ALE est construite pour éviter ces
inconvénients : l’écoulement fluide est calculé sur un domaine qui est déformé de façon à suivre
le mouvement de l’interface (lagrangien près du solide) et dont la vitesse de déformation ne
suit pas nécessairement celle du fluide à l’intérieur du domaine.
Dans ce chapitre, nous décrivons les méthodes numériques les plus couramment utilisées pour
effectuer des simulations d’interaction fluide-structure. Nous nous intéressons, ensuite, à la
description de la formulation ALE et à la formulation variationnelle du problème couplé
fluide/solide. Enfin, nous décrivons différentes méthodes de relèvements permettant la mise à
jour du maillage à chaque pas de temps et nous introduisons le relèvement parabolique. Pour
vérifier numériquement notre approche, nous considérons plusieurs tests numériques d’IFS.
Pour clore ce chapitre, nous présentons une étude de deux applications industrielles.
La première méthode est le couplage dit monolithique ou encore direct. Cette méthode consiste
à résoudre simultanément le système couplé du fluide et du solide en discrétisant les équations
du mouvement dans chaque milieu. Le système matriciel est, ensuite, fermé par l’application
des conditions aux limites à l’interface fluide-structure. Ainsi, cette méthode consiste à écrire
une formulation complète couplant toutes les inconnues du problème et permet de résoudre
les deux problèmes (fluide et solide) dans un cadre unifié qui exige implicitement que les
conditions aux limites d’interaction fluide-structure soient satisfaites d’une façon naturelle.
Huber et al. (2004) ont adopté cette méthode. Dans leurs approches, les formes discrètes
58
de la conservation de la masse fluide et la conservation de la quantité de mouvement pour
l’ensemble du système fluide-structure sont formulées dans un seul système d’équations et
résolues avec le mouvement des mailles du fluide dans une boucle d’itérations unique. Notons
qu’une discrétisation éléments finis stabilisée est appliquée pour discrétiser toutes les équations
modèles, conduisant à une discrétisation uniforme de l’ensemble du système.
Plus récemment, Étienne et Pelletier (2004) ont développé une méthode monolithique de
couplage d’IFS. La méthode numérique proposée repose sur une formulation directe couplée
à des intégrateurs en temps d’ordre élevé (Runge-Kutta/BDF). Les équations incompressibles
de Navier-Stokes pour le fluide, hyperélastiques de Saint-Venant Kirchhoff pour la structure,
de Newton pour la masse ponctuelle (équations de déplacement de corps rigide) et d’équilibre
pour les termes de couplage forment un large système monolithique à résoudre. L’approche
IFS monolithique utilise des noeuds coïncidants sur les interfaces fluide-structure afin que les
efforts, les déplacements et les vitesses soient évalués au même endroit en un temps identique.
Le problème global est résolu de manière implicite grâce à une approche éléments finis de
Newton-Raphson utilisant un pseudo-solide. Plus de détails sont donnés dans Cori (2011).
Dans ce même contexte, on trouve dans la littérature une autre façon pour assurer un tel
couplage. Greenshields et Weller (2005) ont procédé par une méthode monolithique qui consiste
à résoudre les équations de conservation dans tout le domaine par une formulation unifiée.
Cette formulation s’obtient d’une loi constitutive basée sur un paramètre de phase avec lequel
on distingue les deux milieux : fluide et solide. Notons qu’un seul système d’équations de
quantité de mouvement et de continuité est dérivé, gouvernant à la fois le fluide et le solide et
résolu avec un maillage unique en utilisant les schémas de discrétisation de volume finis. La
méthode est validée en simulant les oscillations dynamiques d’une poutre élastique encastrée.
Dunne (2007) propose une formulation monolithique totalement eulérienne pour la résolu-
tion des problèmes d’interaction fluide-structure. Au lieu de changer le cadre de référence du
problème fluide exactement comme dans le méthode ALE présentée ultérieurement, le cadre
de référence de la structure a été changé pour correspondre à la configuration eulérienne du
fluide. Ainsi, toutes les variables de la structure ont été formulées dans un cadre eulérien. Il a
été nécessaire d’introduire une variable qui contient soit la position initiale ou le déplacement
des points matériels. L’avantage majeur de cette formulation complètement eulérienne est que
le maillage global du domaine est fixe et le couplage est monolithique : le fluide et le solide
sont résolus dans un seul bloc et la méthode des positions initiales a permis de déformer le
maillage à chaque pas de temps. L’inconvénient de cette méthode est le coût de calcul puisqu’il
y a 4 variables indépendantes à traiter partout dans le domaine global du fluide et du solide.
Notons que la méthode monolithique a été adoptée, aussi, par Cottet et Maitre (2004, 2006).
L’idée consiste à capturer l’interface par une surface de niveau (level set) qui contient des in-
formations sur l’étirement des membranes solides dans un fluide. La fonction level set remplace
la méthode des positions initiales présentée dans Dunne (2007).
59
Dans ce même contexte, Legay et Belytschko (2006) ont décrit une méthode d’interaction
fluide-structure où l’interface et les surfaces libres sont aussi définies par une surface de ni-
veau. Le fluide est classiquement décrit en formulation eulérienne et la structure élastique en
formulation lagrangienne. La méthode level set fournit un cadre compact pour écrire la forme
faible et pour le développement et le calcul des équations discrétisées.
La deuxième famille des méthodes de couplage d’IFS présentée dans la littérature est le cou-
plage dit partitionné. Cette méthode permet la résolution du problème fluide et du problème
solide séparément. Le couplage est assuré par un transfert d’informations entre les solveurs
(du fluide et du solide) pour assurer la continuité du problème global à l’interface. Générale-
ment, cela se fait en échangeant, en alternance, des informations à chaque pas de temps. Cette
méthode est, donc, caractérisée par la présence de deux solveurs à l’inverse de la méthode
monolithique.
Une manière simple de réaliser un couplage partitionné est d’effectuer une seule fois l’échange
d’informations entre les solveurs du fluide et du solide dans un même pas de temps. Ce cou-
plage partitionné est nommé méthode explicite et non itérative. Matthies et Steindorf (2003)
affirment que ce type de méthode n’est pas d’une grande utilité, surtout, quand il s’agit d’un
couplage fort entre le fluide et le solide. En effet, les invariants globaux tels que la conser-
vation de l’énergie seront presque certainement violés par ce type de méthode. Autrement
dit, le schéma partitionné explicite devient généralement instable quand les forces qui causent
la déformation solide proviennent essentiellement du champ de pression et du cisaillement du
fluide à l’interface (voir par exemple Matthies et Steindorf (2003) et Rugonyi et Bathe (2001)).
Burman et Fernández (2007) ont essayé de stabiliser la méthode explicite avec la méthode de
Nitsche mais les applications restent limitées.
L’alternative à ce couplage explicite est de considérer un couplage partitionné implicite (Mat-
thies et Steindorf (2001)). Le couplage implicite consiste à itérer plusieurs fois dans un seul
pas de temps pour effectuer l’échange d’information à l’interface solide-fluide tel que l’indique
la figure 3.1.
Satisfaire les conditions à l’interface implicitement est souvent appelé un couplage fort. En
60
outre, ce schéma partitionné implicite permet une meilleure stabilité grâce à la conservation
de l’énergie comme l’ont confirmé Mouro et LeTallec (2001). Dans ce contexte, le couplage
partitionné implicite est souvent approprié pour traiter des problèmes faisant intervenir des
interactions fortes. Notons que le schéma de couplage explicite est, généralement, approprié
pour les couplages faibles entre le fluide et la structure tandis que le schéma implicite est
toujours valide. Des comparaisons de la méthode partitionnée implicite avec la méthode mo-
nolithique ont été proposées par Degroote et al. (2009) et Degroote et al. (2010) qui ont obtenu
les mêmes résultats en appliquant les deux approches pour des problèmes similaires.
Dans cette thèse, puisqu’il s’agit de modéliser, par éléments finis, un couplage fort tridimen-
sionnel entre une structure hyperélastique en grandes déformations et un fluide newtonien
incompressible, nous avons décidé d’adopter la méthode partitionnée, étant donnée la taille
remarquable des problèmes tridimensionnels à coupler. De plus, cette méthode permet de
résoudre à l’aide des codes de calculs dédiés à la physique du fluide d’une part, et du so-
lide, d’autre part. Pour chaque milieu, un code de calcul spécifique est utilisé, ce qui permet
d’utiliser toutes les fonctionnalités du solveur solide déjà développé au GIREF.
En général, les formulations lagrangienne et eulérienne présentent chacune des forces et des
faiblesses. La méthode ALE a été développée afin de combiner les avantages de ces deux
formulations. Elle permet d’exprimer les deux physiques dans leur cadre naturel. Pour cette
raison, au lieu de placer le solide dans un cadre eulérien comme dans Dunne (2007), une
alternative, qui est la modélisation par une approche ALE, est utilisée par la majorité des
auteurs tel que dans Dettmer et Peric (2006).
Les méthodes ALE ont été proposées, pour la première fois, dans les méthodes de différences
finies et de volumes finis. Les développements originaux ont été faits, entre autres, par Noh
(1964), Franck et Lazarus (1964) et Hirt et al. (1997). La méthode a, ensuite, été adoptée
dans le contexte des éléments finis et les premières applications se trouvent dans le travail
de Donea et al. (1977), Belytschko et Kennedy (1978) et Hughes et al. (1981) pour résoudre
les problèmes d’écoulements à surface libre et d’interaction fluide-structure. La méthode ALE
est une méthode hybride combinant les descriptions lagrangienne et eulérienne permettant
de déformer le maillage de manière adéquate, et éviter, ainsi, les nombreuses interpolations
intervenant dans d’autres méthodes. Au niveau de la mise en oeuvre, elle permet d’exploiter
des codes pré-existants.
Cette approche consiste à mettre en place un cadre (fixe) de référence qui est transformé à
61
chaque fois au domaine physique désiré. Les équations de mouvement sont, ensuite, transfor-
mées à chaque instant t, d’une configuration de référence Ω̂f vers une configuration courante
Ωf (t) (Ω̂f = Ωf (0)) par la transformation ALE suivante :
A : Ω̂f × R+ → Ωf (t),
(x̂,t) → x = A(x̂,t),
At : Ω̂f → Ωf ,
x̂ → x
Remarque :
Le choix de la transformation est en général arbitraire, à condition que la loi donnée pour le
mouvement de la frontière du domaine soit respectée. En conséquence, une approche classique
consiste à construire la transformation A à partir de l’évolution de la frontière ∂Ωf (t) du
domaine fluide Ωf (t) (voir section 3.5).
déf
Remarquant que la transformation ALE At (x̂) = A(x̂,t) représente la déformation du do-
maine à n’importe quel temps t > 0, on peut définir la vitesse de maille correspondante :
déf
ŵ(x̂,t) = ∂t A(x̂,t)
déf f
Jˆf (x̂,t) = dét F̂ (x̂,t)
Pour une fonction donnée fˆ : Ω̂f × R+ → R définie dans le domaine de référence ALE, sa
description eulérienne est donnée par :
f (x,t) = fˆ(A−1 f
t (x),t)), ∀x ∈ Ω (t), t > 0;
Inversement
fˆ(x̂,t) = f (At (x̂),t), ∀x̂ ∈ Ω̂f .
62
Ainsi, pour un point donné x du domaine Ωf (t), on a w(x,t) = ŵ(x̂,t) où w(x,t) est la vitesse
de maille dans la configuration courante.
Remarque :
Dans le cas général, la vitesse de maille du domaine est différente de la vitesse du fluide, c’est
à dire
w(x,t) 6= u(x,t)
Nous avons maintenant à calculer les termes intervenant dans l’équation de conservation de
la quantité de mouvement sachant que le maillage se déforme. Les termes comportant des
dérivées spatiales intervenant dans les équations ne posent pas de problème car le maillage est
fixé à chaque pas de temps. La difficulté réside dans le calcul de l’accélération qui fait intervenir
une dérivée temporelle de la vitesse. Pour cela, on a besoin d’introduire la relation entre la
dérivée temporelle eulérienne pour un champ eulérien g, notée ∂t g et la dérivée correspondante
dans le sens ALE notée ∂t g|A . Cette relation est donnée dans Nobile (2001) et Fernández et
Gerbeau (2009) telle que
∂g ∂g
= w · ∇g +
∂t A ∂t
Par conséquent, il a été choisi de développer une formulation ALE. Cette approche permet de
modéliser de grandes déformations tout en conservant un maillage de bonne qualité tout au
long du calcul, si elle est couplée à des techniques de mise à jour de maillage adéquates. Ces
techniques seront expliquées plus tard.
Notons l’interface entre le fluide (indexé par f ) et le solide (indexé par s) par Σ(t) = Ωf (t) ∩
Ωs (t) au temps t et par Σ̂ = Ω̂f ∩ Ω̂s à la configuration initiale. La continuité géométrique
63
(non-glissement) impose que :
uf = us sur Σ(t) (3.2)
ce qui signifie que la frontière du fluide doit suivre les déplacements de la structure à la
configuration actuelle (pas de glissement).
La conservation du moment à l’interface est imposée par la condition :
ce qui demande une égalité des tractions, soit, en d’autres termes, un équilibre de contraintes.
Ainsi, les tractions du fluide sont projetées dans la configuration de référence avec un signe
opposé puisque les vecteurs normaux du fluide et de la structure ont des directions opposées.
Dans la section suivante, nous reformulons ces conditions de couplage pour avoir le problème
couplé de l’IFS.
Pour le fluide, puisque Ωf (t) est un domaine mobile, nous considérons la formulation ALE
pour les équations de Navier-Stokes 3.1. Pour la structure, on considère les équations élasto-
dynamiques définies au chapitre 1 dans les équations 1.44. Les deux problèmes sont complétés
des conditions initiales et aux limites données dans 1.10 et 1.11.
D’un point de vue géométrique, pour assurer la compatibilité entre les deux formulations
(fluide et structure) on s’assure que le volume de contrôle de fluide suit le mouvement de
l’interface c’est à dire :
A(x̂,t) = T s (x̂,t), sur Σ̂ (3.4)
64
où T s : Ω̂s × [0,T ] → Ωs (t) est la déformation du solide définie au chapitre 1.
s
Puisque le mouvement de la structure est décrit par le déplacement d̂ , alors il est logique
f
de définir la transformation ALE (A) à partir du déplacement du volume de contrôle d̂ :
Ω̂f × R+ → Rd défini par :
f def
d̂ (x̂,t) = A(x̂,t) − x̂, ∀x̂ ∈ Ω̂f
On note que cette équation définit un couplage géométrique entre les deux problèmes. Finale-
ment, en dérivant l’équation 3.5 par rapport à t on a :
Remarque :
f
Le déplacement du volume de contrôle d̂ (et donc A) peut être choisi arbitrairement. En
s
réalité, il peut être n’importe quel relèvement «extension» de d̂ |Σ̂ sur Ω̂f . On notera par la
suite cette opération par :
f s
d̂ = Ext(d̂ |Σ̂ ) (3.7)
À partir de l’équation 3.7, la configuration courante du domaine fluide Ωf (t) est paramétrée
s
par la transformation ALE A(x̂,t) = x̂ + Ext(d̂ |Σ̂ )
f
Ωf (t) = At (Ω̂f ) = (I Ω̂f + d̂ )(Ω̂f )
D’un point de vue mécanique, les conditions de couplage définies dans la section 3.2.2 sont
applicables. L’équation 3.3 s’écrit dans une description lagrangienne comme suit :
f
Π̂ · n̂s + Jˆf σ̂ f · (F̂ )−T · nf = 0, sur Σ̂ (3.8)
Nous avons toutes les équations maintenant pour définir le problème d’interaction fluide-
structure (voir par exemple dans Fernández et Gerbeau (2009) et Fernández et al. (2009)) :
Trouver la vitesse du fluide uf : Ωf × R+ → R, la pression p : Ωf × R+ → R, le déplacement
s
de la structure d̂ : Ω̂s × R+ → R tel que
65
— Problème du fluide :
f
∂u
ρf
+ (u − w) · ∇u − 2µdiv(γ̇(uf )) + ∇p = f dans Ωf (t)
f f
∂t A
∇ · uf = 0 dans Ωf (t)
(3.9)
uf = ufD sur ΓfD
= g sur ΓfN
f
σ · nf
— Problème de la structure :
s
∂ 2 d̂ s s
s
− ∇ · Π̂ (d̂ ) = ρ̂s0 r̂ s dans Ω̂s ,
ρ̂0
∂t 2
s
= ûs dans Ω̂s ,
∂t d̂
(3.10)
s s
sur Γ̂sD
d̂ = d̂D
s
Π̂ · n̂s = Jˆs ||(F̂ )−T n̂s ||ĥ sur Γ̂sN
— Conditions de couplage :
f s f
d̂ = Ext(d̂ |Σ̂ ), ŵ = ∂t dˆf dans Ω̂f , Ωf (t) = (I Ω̂f + d̂ )(Ω̂f )
us = w sur Σ(t), (3.11)
Π̂ · n̂s + Jˆf σ̂ f · (F̂ f )−T · nf = 0, sur Σ̂
Remarque :
Puisqu’il s’agit dans notre cas d’utiliser la méthode partitionnée définie dans 3.1, résoudre le
problème couplé (3.9 - 3.11), revient à résoudre un problème de type Dirichlet-Neumann. En
effet, une condition aux limites de Dirichlet est imposée à l’interface du fluide (3.112 ) alors
que, pour le problème de la structure, une condition de Neumann est appliquée à la frontière
du solide (3.113 ). Pour plus de détails, le lecteur peut se référer à Quarteroni et Valli (1999)
et Causin et al. (2005).
D’un point de vue énergétique, le problème couplé 3.9 - 3.11 satisfait le bilan énergétique
suivant
ρf f 2 ρs ˆ s 2 ρf f 2
Z Z Z Z Z
d s f 2
|u | + |ḋ | + W (E(d̂ )) + 2µ|γ̇(u )| = g − |u | n · uf
dt Ωf (t) 2 Ω̂s 2 Ω̂s f
Ω (t) Γ f 2
| {z } | {z } | N {z }
Énergie totale Dissipation visqueuse Puissance externe
(3.12)
66
Pour la preuve, le lecteur peut se référer à Fernández et Gerbeau (2009). Ici, on observe
que la dissipation dans le système est produite par la viscosité du fluide et que la puissance
échangée par le fluide et la structure est équilibrée exactement à l’interface. Cet équilibre est
une conséquence directe de la réalisation des conditions de couplage (3.112 ) et (3.113 )
Notons que contrairement aux fonctions tests v̂ f ∈ V̂ f et q̂ ∈ Q̂ qui sont définies dans le
domaine de référence fixe Ω̂f , les fonctions v f ∈ V f et q ∈ Q dépendent du temps puisqu’elles
sont définies dans le domaine courant Ωf . Notons aussi que puisque v̂ f et q̂ sont indépendantes
du temps t, alors on peut écrire
∂v f
∂q
=0 ; =0 (3.18)
∂t A ∂t A
En multipliant les équations fluide 3.9 par (v f ,q) ∈ V0f × Q et en intégrant par parties compte
tenu des conditions aux limites on obtient, ∀(v f ,q)
f
Z Z
f ∂u f f f f f
f
ρ
∂t · v dx + f ρ ((u − w) · ∇u ) · v dx
Ω (t) A Ω (t)
(3.19)
Z Z Z
+ (2µγ̇(uf ) − pI) : ∇v f dx − g · v f dγ + q∇ · uf dx = 0
Ωf (t) ΓfN Ωf (t)
67
Ainsi, la formulation variationnelle du problème du fluide s’écrit :
∀(v f ,q) ∈ V0f × Q et t ∈ R+ , trouver la vitesse uf (x,t) ∈ H 1 (Ωf (t))d et la pression p(x,t) ∈
L2 (Ωf (t)) satisfaisant uf = w sur Σ(t) et uf = ufD sur ΓfD et telle que
f
Z Z
f ∂u f
ρ · v dx + ρf ((uf − w) · ∇uf ) · v f dx
f
Ω (t) ∂t
A f
Ω (t)
(3.20)
Z Z Z
+ (2µγ̇(uf ) − pI) : ∇v f dx + q∇ · uf dx = g · v f dγ
Ωf (t) Ωf (t) ΓfN
∂uf
Z
Remarque : En utilisant un changement de variables dans ρ · v dx et en
f f
f
Ω (t) ∂t A
combinant avec 3.18, il est possible de montrer que
f
Z Z Z
f ∂u d
ρ · v f
dx = ρ f f
u · v f
− ρf (∇ · w)uf · v f dx (3.21)
f
Ω (t) ∂t
A dt f
Ω (t) f
Ω (t)
En remplaçant l’équation 3.21 dans 3.20, on obtient une autre forme de la formulation varia-
tionnelle (voir par exemple Fernández (2001) et Astorino (2010))
s déf
V = s s
v (t) : Ω (t) → R ,v = v̂ ◦ d s s
Tt−1 , v̂ s ∈ V̂ s
, (3.23)
En multipliant l’équation du solide 3.10 par v̂ s ∈ V̂ s et en intégrant par parties tout en tenant
compte des conditions aux limites 3.103 et 3.104 , et de la condition à l’interface 3.113 , on
obtient Z 2 d̂s Z
s∂ s
ρ̂ · v̂ dx̂ + Π̂ : ∇x̂ v̂ s dx̂
Ω̂s ∂t2 Ω̂s
(3.24)
Z Z
s −T s −T
= Jˆs ||(F̂ ) n̂s ||ĥ · v̂ s dγ̂ + Jˆs σ̂ f · (F̂ ) · n̂s · v̂ s dγ̂
Γ̂sN Σ̂
En remplaçant la deuxième équation de 3.10 dans 3.24 et en utilisant la transformation entre
la configuration de référence et la configuration courante sachant que nf = −ns , on obtient
la formulation variationnelle du problème de la structure qui s’écrit
s
Trouver d̂ ∈ V̂ s tel que ∀v̂ s ∈ V̂ s on a
Z 2 s Z
s ∂ d̂ s
ρ̂ 2
· v̂ dx̂ + Π̂ : ∇x̂ v̂ s dx̂ =
Ω̂s ∂t Ω̂s
(3.25)
Z Z
s
Jˆs ||(F̂ )−T n̂s ||ĥ · v̂ s dγ̂ − (σ f · nf ) · v s dγ
Γ̂sN Σ(t)
68
Formulation variationnelle du problème couplé
Nous avons maintenant tous les outils pour définir la formulation variationnelle du problème
couplé de fluide-structure. Le problème 3.9, 3.10 et 3.11 s’écrit
∀(v f ,q,v̂ s ) ∈ V0f × Q × V̂ s ,
s
trouver (uf ,p,d̂ ) ∈ H 1 (Ωf (t))d × L2 (Ωf (t)) × V̂ s satisfaisant uf = w sur Σ(t) et uf =
ufD sur ΓfD telle que
— Problème du fluide
∂uf
Z Z
ρf · v f
dx + ρf ((uf − w) · ∇uf ) · v f dx
Ωf (t) ∂t A f
Ω (t)
(3.26)
Z Z Z
+ (2µγ̇(uf ) − pI) : ∇v f dx + q∇ · uf dx = g · v f dγ
Ωf (t) Ωf (t) ΓfN
— Problème de la structure
Z 2 d̂s Z
s∂ s
ρ̂ · v̂ dx̂ + Π̂ : ∇x̂ v̂ s dx̂
Ω̂s ∂t2 Ω̂s
(3.27)
Z Z
s −T
= Jˆs ||(F̂ ) n̂s ||ĥ · v̂ s dγ̂ − (σ f · nf ) · v s dγ
Γ̂sN Σ(t)
— Conditions de couplage
f s f
d̂ = Ext(d̂ |Σ̂ ), ŵ = ∂t dˆf dans Ω̂f , Ωf (t) = (I Ω̂f + d̂ )(Ω̂f ) (3.28)
69
3.5.1 Revue de la littérature
Parmi les nombreuses méthodes introduites au cours des dernières années, beaucoup sont ba-
sées sur la détermination du mouvement de maillage par une équation différentielle partielle
en utilisant des déplacements au bord, comme conditions aux limites. Une méthode équipo-
tentielle proposée par Winslow (1963) considère les lignes de maillage comme deux ensembles
entrecroisés d’équipotentielles. En utilisant la méthode de Winslow (1963), Godunov et Pro-
kopov (1972) ont conçu un algorithme pour générer des maillages pour des problèmes où les
changements dans les données des conditions aux limites sont reflétées dans les changements
du maillage. Astorino (2010), Fernández (2001) et Nobile (2001) ont utilisé les relèvements
harmoniques qui consistent à résoudre l’équation de Laplace pour mettre à jour le maillage
fluide à chaque pas de temps, lors de la résolution des problèmes de fluide-structure. Brackbill
et Saltzman (1982) ont aussi résolu l’équation de Laplace mais avec des termes non homo-
gènes pour optimiser la régularité, l’orthogonalité et la variation du volume des cellules. Wang
et McLay (1986) proposent la résolution d’équations différentielles du quatrième ordre (opé-
rateur de Bilaplacien) pour le traitement des mouvements du maillage. Chen et al. (1988)
ont proposé la résolution des équations d’élasticité modifiées dans lesquelles le jacobien des
transformations est exclu descalculs, introduisant, ainsi, un effet variable de raideur dans le
domaine de calcul (voir aussi Johnson et Tezduyar (1994) et Tezduyar et al. (1996)). Dans le
même contexte, Stein et al. (2004) ont conçu une méthode où ils considèrent que les éléments
les plus petits, généralement placés à proximité des surfaces solides, sont plus raidis que les
grands. Chiandussi et al. (2000) montrent que les critères géométriques ne sont pas efficaces
si on les compare avec les critères basés sur des paramètres structuraux tels que le champ de
déformation ou la densité d’énergie de déformation. Bottasso et al. (2005) ont étudié le pro-
blème de déformation de maillages non structurés en utilisant des systèmes pseudo-structurels
à paramètres localisés. Dans leur méthode, une analogie par un modèle à ressort est utilisé. En
effet, un réseau de base d’arêtes du maillage est complété par un ensemble supplémentaire de
ressorts linéaires qui s’opposent au retournement des éléments. La même méthode est explorée
en 3D par Degand et Farhat (2002). Ces derniers utilisent l’analogie à un système de ressorts
pour mettre à jour la position des maillages dynamiques non structurés.
Rendall et Allen (2008) ont présenté un schéma d’interpolation multivariée, utilisant des
fonctions de base radiale (RBF), qui se traduit par une formulation totalement unifiée de
l’interpolation fluide-structure et des problèmes de mouvement de maillages (voir aussi Rendall
et Allen (2009a,b)). Les applications des méthodes RBF sont présentées dans Ahrem et al.
(2006) et dans Beckert et Wendland (2001). Une comparaison utile des méthodes RBF, d’une
part, et d’autres méthodes d’échange d’informations entre les solveurs fluides et structures,
d’autres part, est établie par de Boer et al. (2005) pour un écoulement quasi 1D. La théorie
générale des méthodes RBF est présentée par Buhmann (2003) et Wendland (2005).
Dans ce qui suit, on présentera quatre méthodes de relèvements pour déplacer les maillages et
70
on en comparera l’efficacité. On s’intéresse à la capacité de ces méthodes à soutenir de grands
déplacements.
Une des formes, les plus classiques et plus utilisées dans la construction d’une extension des
déplacements de bord, consiste à résoudre une équation de Laplace (ou Poisson) pour chaque
composante de la position des nœuds. Cette méthode, initialement proposée par Winslow
(1963) et nommée «relèvement harmonique», est parfois appelée «génération elliptique de
maillage». Elle a été utilisée par plusieurs auteurs (voir Brackbill et Saltzman (1982), Fernández
(2001), Astorino (2010)). Malgré la simplicité de cette technique, elle présente un inconvénient
important car, dans certains cas, la méthode «pousse» les nœuds intérieurs à l’extérieur du
domaine, produisant, ainsi, des éléments non valides. Pour garantir un maillage résultant
valide, un coefficient de diffusion κ est introduit, conduisant à la formulation :
f
−∇ · (κ(x)∇d̂ ) = 0 dans Ω̂f ,
f
d̂ · nf = 0 sur ∂ Ω̂f \Σ̂, (3.29)
d̂f s
sur Σ̂,
= d̂
où les bords du domaine fluide appartenant à ∂ Ω̂f \Σ̂ sont autorisés à glisser le long des
parois. Il convient de noter que l’imposition d’une condition de Dirichlet sur le bord permet
f
un découplage des composantes de d̂ et une réduction de la taille du système discret (voir par
exemple Fernández et Gerbeau (2009)). Dans ce cas, le problème 3.29 est résolu séparément
pour chaque composante du déplacement. Cependant, même si dans certains cas cela pourrait
être possible, puisque nous bloquons tout mouvement le long de ∂ Ω̂f \Σ̂, cela va produire
des maillages contenant des éléments allongés et éventuellement invalides. La condition de
glissement proposée, ici, a tendance à réduire la détérioration des éléments au prix d’un système
plus important à résoudre. Bien entendu, si le domaine a une forme rectangulaire, nous pouvons
f
fixer l’axe de telle sorte qu’un découplage total des composantes de d̂ est possible, même pour
la condition de glissement.
En ce qui concerne le coefficient de diffusion κ, dans le cadre des éléments finis, deux approches
sont fréquemment utilisées. L’une en prenant en compte la distance par rapport à la frontière
en déplacement. Quand à la seconde, elle prend en compte la taille des éléments (permettant
une plus grande déformation pour un plus grand élément). Suivant les travaux de Wick (2011),
nous présentons, ici, les résultats obtenus en utilisant un coefficient de diffusion en fonction
du volume de l’élément en 3D (ou aire de l’élément en 2D).
71
En désignant par T le maillage des éléments finis, la diffusion est définie par :
M 2
M h
h = max he , κ(x) = , x ∈ τe ,
τe ∈T he
où he est le volume de l’élément τe ∈ T . Avec une telle diffusion, les éléments de petites tailles,
généralement près de la frontière mobile auront un coefficient de diffusion relativement élevé,
conduisant à de petites déformations, alors que les grands éléments, où une grande déformation
est moins susceptible de produire des maillages invalides, auront un coefficient de diffusion plus
petit.
Encore une fois, une condition aux limites de glissement est préférée à un blocage total des
déplacements, car elle empêche la déformation des éléments proches de la frontière de Ω̂f . Tel
est le cas du relèvement harmonique : pour retarder le plus possible la production de maillages
non valides, les deux paramètres du matériau virtuel dépendent de la taille élémentaire. Nous
proposons la relation suivante
M 2 2
8 h he
E(x) = 10 , ν(x) = 0.45 − 0.15 , x ∈ τe .
he hM
Notons que si ν est au plus égal à 0.45, il assure la validité de la formulation des déplacements.
En utilisant une valeur plus grande pour ce coefficient (près de 0.5), on se verrait obligé d’uti-
liser une formulation mixte (déplacement-pression) pour tenir compte de l’incompressibilité.
Certes, il est tentant de compter sur l’incompressibilité afin de préserver l’aire des éléments
(volume en 3D), mais le coût de calcul d’une telle approche serait prohibitif, comme indiqué
dans Wick (2011). Notez que dans cette approche, nous ne pouvons pas résoudre le problème
72
composante par composante. Cette méthode est nettement plus coûteuse (d’un point de vue
informatique) que le relèvement harmonique.
Notons que X = (X1 ,X2 ,X3 ) et X j = (X1j ,X2j ,X3j ). Cela donne un système linéaire de la
forme :
73
où Aij = g(kX i − X j k) et dépend, donc, de la distance entre les nœuds des frontières grâce à
la fonction radiale g. La matrice B est une matrice 4 × n composée pour chaque colonne par
1, suivis par les coordonnées de chaque nœud. Comme dans le cas de l’extension harmonique,
les déplacements peuvent être calculés pour chaque composante séparément (en utilisant la
même matrice exacte). Il est facile de montrer que le krigeage peut récupérer exactement les
mouvements des corps rigides. C’est en effet, l’une des caractéristiques intéressantes dans le
contexte de déplacement des nœuds à l’interface fluide-structure. Dans ce cas, les coefficients
αi disparaissent tous et le mouvement exact est donné par la partie linéaire.
Dans sa forme originale, le krigeage nécessite la solution d’un système linéaire dense, limi-
tant, ainsi, son application aux problèmes bidimensionnels ou à de petits problèmes en trois
dimensions (voir par exemple Olivier (2014) pour plus de détails).
Pour pallier ce problème, un certain nombre de stratégies a été introduit. La première est de
négliger la contribution de la ime ligne de la matrice A des nœuds situés trop loin de X i et de
récupérer ainsi une matrice creuse. Le choix de la distance maximale est, cependant, délicat.
Comme dans Rendall et Allen (2009a), il est également possible d’envisager un sous-ensemble
approprié des nœuds de la frontière, réduisant ainsi la taille de la matrice globale. Le choix
des nœuds est fait de façon à minimiser l’erreur d’interpolation.
telle que proposée dans Beckert et Wendland (2001) qui rend la matrice A définie positive.
Aucune différence significative n’a été observée par rapport à notre fonction g utilisée dans ce
travail.
Dans ce travail, nous voulons comparer les résultats finaux de la méthode RBF, à savoir les
maillages obtenus suite à de grands déplacements imposés. Compte tenu de l’écart impor-
tant de cette méthode par rapport aux autres (toutes basées sur la solution d’un problème
d’éléments finis), l’évaluation de la performance de calcul est en dehors du cadre du présent
document. Par conséquent, nous avons utilisé la méthode RBF dans sa forme la plus simple,
sans chercher à optimiser son efficacité numérique. De part notre connaissance de la méthode,
74
comme mise en œuvre pour nos besoins, elle ne devrait pas présenter des différences drama-
tiques par rapport aux résultats obtenus grâce à des implémentations plus numériquement
efficaces, telles que celles proposées par Rendall et Allen (2009a).
Nous proposons une méthode de relèvement basée sur l’équation de la chaleur en utilisant un
temps artificiel et arbitraire.
f
∂ d̂ f
dans Ω̂f ,t ≥ 0
− ∇ · (κ(x)∇d̂ ) = 0
∂t
f (3.34)
d̂ · nf = 0 sur ∂ Ω̂f \Σ̂,
f
s
d̂ = d̂ sur Σ̂,
L’idée de départ de cette méthode était de garder la simplicité du relèvement harmonique tout
en améliorant son comportement et évitant parallèlement la génération d’éléments invalides
près de la frontière mobile. La solution de l’équation de la chaleur a de nombreuses propriétés
analytiques : régularité, principe du maximum, comportement asymptotique, ceci contribue
à l’intérêt d’une telle solution comme relèvement. En exploitant ces propriétés, le relèvement
produit par cette équation converge vers le relèvement harmonique quand t tend vers l’infini.
Le coefficient de diffusion κ(x) peut jouer un rôle similaire à celui du coefficient de diffusion
dans le relèvement harmonique (ou le coefficient d’élasticité pour le relèvement élastique). Il
peut être construit sous la forme d’une fonction du volume des éléments (aire en 2D), ou une
fonction de la distance par rapport à la frontière mobile. Il est intéressant de remarquer que
dans les cas présentés ici, il n’y avait aucun gain observable à l’utilisation d’une définition plus
élaborée de κ(x) et il a été tout simplement mis à 1. La simplicité relative de cette approche
et sa capacité à traiter les mouvements de grandes amplitudes rendent cette méthode très
attrayante et numériquement efficace.
En ce qui concerne les conditions aux limites sur ∂ Ω̂f \Σ̂, les mêmes remarques que pour
le relèvement harmonique peuvent être faites : l’imposition d’une condition de Dirichlet sur
f
∂ Ω̂f \Σ̂ permet un découplage des composants de d̂ et une réduction de la taille du système
discret. Notons que résoudre, composante par composante, permet de réduire, d’un facteur de
3, la taille du système discret associé. Toutefois, cela se fait au détriment de la qualité des
maillages près de ∂ Ω̂f \Σ̂. Nous préférons, donc, imposer une condition de glissement présentée
dans (3.34). Comme pour les relèvements harmoniques, si le domaine de calcul est de forme
rectangulaire ou prismatique rectangulaire, la condition de glissement peut être découplée et
la taille du système sera réduite, en conséquence.
75
Résolution du relèvement parabolique
Différents schémas d’intégration en temps peuvent être utilisés pour résoudre le système (3.34).
Quelques remarques sont nécessaires concernant notre choix. La précision de la solution ap-
prochée du système (3.34) n’est pas indispensable et un schéma en temps de premier ordre
(Euler) est largement suffisant et mènera à un faible coût de calcul.
Deux méthodes ont été explorées, un schéma d’Euler implicite et un schéma d’Euler expli-
cite avec une matrice masse condensée. Bien que cette dernière méthode n’ait pas été retenue
pour le présent travail, dans un contexte plus pratique (par exemple dans le cas des vraies
interactions fluide-structure), cette méthode pourrait être extrêmement efficace : l’utilisation
d’une matrice masse condensée conduit à un algorithme de résolution explicite (pas de système
linéaire à résoudre). Ce schéma est difficile à utiliser suite à sa stabilité conditionnelle. De part
notre expérience, seul un petit nombre de pas de temps peut être utilisé. Dans le contexte
actuel où de grandes amplitudes de déplacements sont imposées, cette méthode conduit inévi-
tablement à des maillages invalides. Dans un contexte plus réaliste, où le temps physique est
faible (ce qui conduit à de petits déplacements des frontières), cette approche pourrait être
utilisée. Dans ce cas, ce serait la méthode la plus efficace car elle fait impliquer seulement des
opérations simples sur les vecteurs. Comme pour le schéma d’Euler explicite, cela impliquerait
à résoudre un système algébrique (impliquant la matrice masse), conduisant à des coûts infor-
matiques comparables à ceux du schéma implicite. Pour cette raison et étant donné qu’il est
inconditionnellement stable, le schéma préféré est le schéma implicite, permettant l’utilisation
d’un pas de temps de longueur arbitraire.
La résolution du système (3.34) sur une longue période de temps (éventuellement avec de
nombreux pas de temps) conduirait à des maillages très proches de ceux obtenus avec le
relèvement harmonique. Il serait, en effet, un moyen très inefficace pour obtenir un relèvement
harmonique. Dans les exemples numériques présentés dans ce travail, la situation est quelque
peu artificielle : nous imposons instantanément un très grand déplacement. Pour limiter le
nombre de pas de temps (pour des raisons de calcul), nous devons utiliser une longueur de pas
suffisante pour que l’effet de diffusion ait lieu, sans état stationnaire. D’après la définition du
coefficient de diffusion, nous avons
D2
κ=1=
τ
où τ est le temps nécessaire à la diffusion pour qu’elle se propage sur une distance D. Dans
la pratique, D peut être choisie comme la valeur maximale du déplacement imposé, donnant
ainsi τ . Même pour les grands déplacements D, un seul pas de temps, de longueur ∆t = D2 ,
était suffisant pour obtenir des résultats acceptables (dans tous les cas, un deuxième pas de
temps a donné un gain modeste dans la qualité du maillage). Notons qu’en utilisant plusieurs
76
petits pas de temps, nous obtenons des résultats similaires, mais ceci augmenterait le coût de
calcul.
Nous sommes à la recherche d’une méthode, ayant la capacité de produire des maillages de
qualité acceptable pour de nombreux pas de temps (pour des phénomènes instationnaires ou
quasi-statiques). De manière équivalente, nous cherchons à montrer, à partir de ces cas d’étude,
la capacité de chaque relèvement à maintenir le plus possible un déplacement de la frontière.
On s’intéressera dans ces tests à deux aspects : l’amplitude maximale de la déformation avant
l’apparition d’éléments invalides et la qualité générale des maillages obtenus. En effet, la qualité
des éléments doit être maintenue pour garantir que le solveur éléments finis puisse converger,
spécialement quand on utilise des méthodes itératives. La validité du maillage est cependant
la mesure principale, puisqu’un maillage valide peut toujours être amélioré, ce qui n’est pas le
cas d’un maillage invalide.
Il existe un certain nombre de définitions de la qualité d’un élément et cette notion est souvent
utilisée dans les méthodes d’adaptation de maillage. La qualité d’un élément est généralement
un nombre compris entre 0 pour un élément dégénéré et 1 pour un triangle équilatéral ou de
tétraèdre. Dans le cas bidimensionnel, on utilise la définition proposée par Bank et Xu (1996) :
√
4 3|T |
Q(T ) = 2
l1 + l22 + l32
où |T | est la valeur absolue de l’aire du triangle T et li , 1 ≤ i ≤ 3 sont ses longueurs d’arêtes.
En dimension trois, on définit
12 |T |
Q(T ) = √ √
2 l1 l2 l3 l4 l5 l6
où |T | est maintenant la valeur absolue du volume du tétraèdre T et li , 1 ≤ i ≤ 6 sont ses
longueurs d’arêtes.
Les quatre cas illustrent les avantages et les limites des méthodes décrites dans la section 3.5.
Ils impliquent, tous, un paramètre d (0 ≤ d ≤ 1) qui contrôle l’amplitude du mouvement. Tous
les tests (sauf cas # 1) peuvent être trouvés dans la littérature d’interaction fluide-structure.
Le profil aérodynamique mobile (cas # 2), est un problème bidimensionnel classique (voir par
exemple Johnson et Tezduyar (1994); Lin et al. (2014)). La sphère en mouvement (cas # 3),
est également traitée dans Lin et al. (2014). Enfin, le drapeau flottant (cas # 4) a été inspiré
des travaux de Schäfer et Turek (1996) et est également considéré dans Wick (2011).
77
Comme nous le verrons, le relèvement harmonique est la méthode la moins exigeante en termes
de coût de calcul, suivie par la méthode parabolique. Le relèvement élastique est la méthode la
plus coûteuse en mémoire et en temps de calcul. Puisque cette conclusion était prévisible, nous
présentons les valeurs relatives au coût computationnel seulement pour le dernier cas (voir le
tableau 3.9), qui est le plus exigeant. Rappelons que nous ne recueillons pas des indicateurs de
performance numériques pour la méthode de krigeage, puisque notre mise en œuvre de cette
méthode est clairement sous-performante d’un point de vue programmation.
Nous présentons un cas bidimensionnel comprenant un déplacement imposé d’une tige rigide
placée dans un canal. Le canal est de hauteur h = 1 et de longueur l = 6. La tige rigide est de
hauteur H = 0.45 et d’épaisseur e = 0.1. Notons que cette tige est fixée à sa partie inférieure
par une condition de Dirichlet nulle, soit U s = 0, voir figure 3.3.
h
H
Figure 3.3 – Translation horizontale imposée à une tige élastique, déplacement maximal de
la tige (à droite).
78
d Hext Eext Kext Pext
0.09 0.94 0.94 0.94 0.94
0.18 0.94 0.94 0.93 0.94
0.28 0.92 0.93 0.92 0.93
0.46 — 0.90 0.89 0.91
0.64 — 0.88 0.85 0.90
0.82 — 0.85 0.80 0.88
1.00 — — 0.76 0.81
Table 3.1 – Cas #1 : Qualité moyenne des éléments du maillage. Hext , Eext , Kext et Pext
présente respectivement les relèvements harmonique, élastique, krigeage et parabolique
s 0.55y 2
d̂ = d
H2
où d est une constante variant entre 0.0 (position initiale) et 1.0. Pour une valeur donnée de d,
la déformation totale est imposée en une seule étape, et donc le problème devient plus difficile
lorsqu’on augmente d.
Les calculs sont effectués à partir d’un maillage triangulaire non structuré, ayant 4 475 éléments
(2 373 nœuds et 6 847 arêtes). Le maillage initial a été obtenu en utilisant GMSH (voir
Remacle et Geuzaine (2009)) et est illustré à la figure 3.4. La qualité moyenne des éléments
est de 0.94. Le déplacement maximal est obtenu lorsque d = 1 et y = H et donc D = 0.55
(voir section 3.5.5). Le pas de temps artificiel pour le relèvement parabolique a, donc, été fixé
à ∆t = (0.55)2 ' 0.3. Dans le tableau 3.1, nous présentons les résultats de chaque méthode
de relèvement en utilisant différentes valeurs de d.
Le symbole «—» indique un maillage invalide (une aire d’un triangle négative). De ce tableau,
nous pouvons facilement conclure que le relèvement harmonique impose une sérieuse limitation
de l’amplitude du déplacement. Des maillages invalides sont obtenus à partir de d = 0.25. Pour
une application en fonction du temps, ceci correspond à une limitation de la longueur du pas
de temps ou à un remaillage plus fréquent par rapport aux autres méthodes. Notons que le
krigeage et le relèvement parabolique fonctionnent très bien et maintiennent une bonne qualité
moyenne d’éléments.
Dans le tableau 3.2, nous présentons pour chaque méthode la distribution de la qualité
des éléments pour d = 0.18, 0.46 et 1.0 (le nombre d’éléments invalides est indiqué entre
parenthèses). La distribution initiale de la qualité des éléments est également indiquée dans
chaque cas puisque nous commençons toujours par le même maillage. Il est facile de voir que le
maillage initial est essentiellement constitué de triangles presque équilatéraux. Pour d = 0.18,
le relèvement harmonique ne produit aucun élément invalide, mais génère 25 éléments invalides
79
d Qualité Initiale (%) Hext (%) Eext (%) Kext (%) Pext (%)
0.18 invalide (0) (0) (0) (0) (0)
0.00 - 0.25 0 0 0 0 0
0.25 - 0.50 0 0 0 0 0
0.50 - 0.75 0.11 0.54 0.53 1.02 0.26
0.75 - 1.00 99.89 99.46 99.47 98.98 99.74
0.46 invalide (0) (1) (0) (0) (0)
0.00 - 0.25 0 — 0 0 0
0.25 - 0.50 0 — 0.10 0.47 0.04
0.50 - 0.75 0.11 — 6.40 10.80 3.60
0.75 - 1.00 99.89 — 93.50 88.73 96.36
1.00 invalide (0) (25) (8) (0) (0)
0.00 - 0.25 0 — — 2.9 1.6
0.25 - 0.50 0 — — 14.4 6.9
0.50 - 0.75 0.11 — — 19.3 13.0
0.75 - 1.00 99.89 — — 63.4 78.4
Table 3.2 – Cas #1 : Distribution de la qualité des éléments pour d = 0.18, 0.46 et 1.00. Le
nombre d’éléments retournés est indiqué entre parenthèses.
pour d = 1.0. Le relèvement parabolique ne produit aucun élément invalide pour toutes les
valeurs de d et présente le pourcentage le plus élevé de la qualité des éléments.
À la figure 3.5, nous présentons les maillages obtenus en utilisant les différentes méthodes
de relèvements pour d = 1.0. Pour les relèvements harmoniques et élastiques, les éléments
retournés peuvent être vus sur la pointe de la barre et les maillages sont donc invalides. Ces
éléments retournés commencent à apparaître respectivement pour d > 0.44 et d > 0.78. Le
relèvement parabolique et le relèvement par krigeage fournissent des maillages valides pour
tous d ≤ 1.0.
Nous illustrons, au prochain cas, l’efficacité et les limites de chacune de ces méthodes dans un
deuxième problème en 2D.
Nous étudions, maintenant, le mouvement de rotation d’un profil d’aile rigide NACA0012. La
géométrie du problème et le mouvement du profil sont décrits dans la figure 3.6. La déformation
du maillage est produite par une rotation de l’aile d’un angle d’attaque initial de 0° jusqu’à
une position d’angle θc = 75°. En notant les coordonnées dans la configuration initiale Γ0 par
(x01 ,x02 ), la nouvelle position du profil avec une rotation de θ = d θc sera
80
Relèvement harmonique
Relèvement élastique
Relèvement parabolique
Figure 3.5 – Cas #1 : Maillages obtenus pour les différentes méthodes (d = 1.0) : à gauche
une vue rapprochée, à droite un agrandissement de l’extrémité de la barre
81
30
10
θ
20
1
3 c=1
y
Les calculs ont été effectués pour un maillage triangulaire non structuré, ayant 5 678 éléments
(2 921 nœuds et 8 599 arêtes). Ce maillage initial et une vue à proximité du profil d’aile sont
présentés dans la figure 3.7. Sa qualité moyenne d’éléments est de 0.95. Pour ce problème, le
déplacement maximal est observé au niveau du bord de fuite du profil et est approximativement
égal à D = 0.8. Le pas de temps artificiel du relèvement parabolique a donc été fixé à ∆t =
(0.8)2 .
82
d Qualité Initiale (%) Hext (%) Eext (%) Kext (%) Pext (%)
0.12 invalide (0) (0) (0) (0) (0)
(9°) 0.00 - 0.25 0 0 0 0 0
0.25 - 0.50 0 0 0 0 0
0.50 - 0.75 0.18 0.19 0.3 0.16 0.08
0.75 - 1.00 99.82 99.1 99.7 99.84 99.92
0.40 invalide (0) (0) (0) (0) (0)
(30°) 0.00 - 0.25 0 0 0 0 0
0.25 - 0.50 0 0.02 0.04 0 0
0.50 - 0.75 0.18 2.94 2.50 0.48 0.74
0.75 - 1.00 99.82 97.04 97.46 99.52 99.26
0.80 invalide (0) (40) (0) (0) (0)
(60°) 0.00 - 0.25 0 — 0.04 0 0
0.25 - 0.50 0 — 1.70 0 8.58
0.50 - 0.75 0.18 — 10.30 4.33 13.95
0.75 - 1.00 99.82 — 87.96 91.67 77.47
1.0 invalide (0) (168) (1) (0) (0)
(75°) 0.00 - 0.25 0 — — 0 5.5
0.25 - 0.50 0 — — 0.04 15.5
0.50 - 0.75 0.18 — — 9.73 4.0
0.75 - 1.00 99.82 — — 90.3 75.0
Table 3.4 – Cas #2 : Distribution de la qualité d’éléments pour θ = 9°, 30°, 60°et 75°(d = 0.12,
0.40, 0.80 et 1.0 resp.). Le nombre d’éléments retournés est indiqué entre parenthèses.
Dans le tableau 3.3, nous présentons la qualité moyenne des éléments du maillage initial et
des maillages obtenus pour les différentes méthodes de relèvements à différentes valeurs de
d. La distribution de la qualité des éléments du maillage pour chaque méthode est, à son
tour, présentée dans le tableau 3.4 pour d = 0.12 (9°), d = 0.4 (30°), d = 0.8 (60°) et
d = 1.0 (75°). Comme on peut le voir, des éléments retournés sont présents pour le relèvement
harmonique à d = 0.8 et pour le relèvement élastique à d = 1.0. En fait, les relèvements
harmoniques et élastiques permettent seulement un maximum de rotation (en une seule étape)
d’environ θ = 46°(d = 0.61) et θ = 72°(d = 0.96), respectivement. Une perte importante de
la qualité d’éléments est observable (voir le tableau 3.4) pour d = 0.8 et d = 1.0 pour les
relèvements : élastique, par krigeage et parabolique (ce dernier souffre le plus). Même si le
relèvement parabolique a produit des maillages valides jusqu’à d = 1.0, il s’agit d’un des
rares cas où ce relèvement semble incapable d’être aussi efficace que la méthode de krigeage.
L’utilisation d’un coefficient variable (coefficient κ dans (3.34)) et/ou de multiples pas de temps
n’ont pas augmenté la qualité du maillage. Bien que le maillage obtenu par le relèvement
parabolique puisse être considéré acceptable, nous concluons que pour ce cas test, seul le
relèvement par krigeage fonctionne adéquatement. Le krigeage produit une bonne distribution
de la qualité des éléments pour toutes les valeurs de d avec plus de 90% de bons éléments pour
83
Relèvement harmonique
Relèvement élastique
Relèvement parabolique
Figure 3.8 – Cas #2 : Maillages obtenus pour les différentes méthodes (θ = 60°). Au milieu
et à droite, une vue rapprochée des extrémités du profil.
Nous présentons à la figure 3.8 les maillages obtenus après une rotation de θ = 60°. Il s’agit
d’une situation typique qui a déjà été prise en compte dans Kinsey et Dumas (2008) et Cori
(2011). Le maillage obtenu avec le relèvement harmonique est invalide. Le relèvement élastique
est capable de produire des maillages valides pour θ ≤ 72°, comme le montre l’agrandissement
du bord d’attaque à la figure 3.8, tandis que le relèvement parabolique est allé jusqu’à θ ≤ 86°.
84
Pour ce problème spécifique, la méthode de krigeage semble être la meilleure méthode et est
capable de traiter une rotation allant jusqu’à 120°(toutefois la qualité du maillage se dégrade).
Le cas suivant traité est un cas en 3D et consiste à appliquer un déplacement vertical à une
sphère.
Nous étudions maintenant un mouvement rigide d’une sphère. Les configurations du problème
sont illustrées à la figure 3.9.
20
35
x2
x1
x3 20
La sphère de diamètre unitaire est placé à la position (0,0,0). Le but de cette section est de
déplacer verticalement la sphère dans le domaine. La condition de déplacement appliquée sur
la sphère est donnée par
s
d̂ = −d (0, 5, 0)
Les calculs ont été faits pour un maillage du domaine constitué de 6 386 éléments P2 tétra-
édriques, soit 1 529 sommets, 8 866 arêtes et 13 725 faces. La qualité moyenne des éléments est
de 0.58. Même si la qualité initiale peut être améliorée, ce test nous permettra aussi d’illustrer
la capacité des différents relèvements face à un maillage dont la qualité initiale est pauvre.
Pour ce problème, le déplacement maximal est obtenu quand d = 1 et donc D = 5 (voir
section 3.5.5). Le pas de temps artificiel a, donc, été fixé à ∆t = (5)2 .
Le tableau 3.5 résume, encore une fois, les résultats obtenus avec les différentes approches
pour les valeurs de d entre 0.2 et 1. Les éléments retournés apparaissent rapidement pour
85
d Hext Eext Kext Pext
0.2 0.58 0.58 0.58 0.58
0.4 0.58 — 0.57 0.58
0.6 — — 0.56 0.57
0.8 — — — 0.57
1.0 — — — 0.56
d Qualité Initiale (%) Hext (%) Eext (%) Kext (%) Pext (%)
0.2 invalid (0) (0) (0) (0) (0)
0.00 - 0.25 3.5 3.5 3.5 3.5 3.2
0.25 - 0.50 20.7 20.8 21.0 21.4 21.1
0.50 - 0.75 72.4 72.5 72.3 71.8 72.3
0.75 - 1.00 3.4 3.2 3.2 3.3 3.4
0.4 invalide (0) (0) (1) (0) (0)
0.00 - 0.25 3.5 3.5 — 3.5 3.1
0.25 - 0.50 20.7 21.1 — 22.3 21.2
0.50 - 0.75 72.2 72.5 — 71.3 72.3
0.75 - 1.00 3.4 3.2 — 2.9 3.4
1.0 invalide (0) (8) (10) (4) (0)
0.00 - 0.25 3.5 — — — 4.1
0.25 - 0.50 20.7 — — — 24.9
0.50 - 0.75 72.4 — — — 68.4
0.75 - 1.00 3.4 — — — 2.6
Table 3.6 – Cas #3 : Distribution de la qualité des éléments pour d = 0.2, 0.4 et 1.0 . Le
nombre d’éléments retournés est indiqué entre parenthèses.
toutes les méthodes sauf le relèvement parabolique qui permet le déplacement maximal. La
distribution de la qualité des éléments est détaillée dans le tableau 3.6.
De manière assez surprenante, ce test semble extrêmement difficile pour toutes les méthodes et
révèle une sensibilité relativement élevée par rapport au maillage initial pour les relèvements
harmoniques, élastiques et par krigeage. La mauvaise qualité du maillage initial est clairement
en cause et un meilleur maillage initial pourrait probablement changer ces résultats. Il est
intéressant d’observer que la distribution (voir tableau 3.6) est relativement stable. Il n’y a
pas une détérioration nette de la qualité d’éléments jusqu’à l’apparition soudaine d’éléments
retournés. Le relèvement parabolique est la seule méthode capable de soutenir le déplacement
maximal imposé.
86
Figure 3.10 – Cas #4 : Géometrie et conditions aux limites.
Dans cette section, nous présentons, en trois dimensions, un autre cas d’une membrane souple
en mouvement dans un fluide. Le problème peut être considéré comme un modèle simplifié pour
l’interaction d’un drapeau dans le vent. La géométrie du problème et toutes les dimensions
pertinentes sont données à la figure 3.10.
où d est une constante qui varie entre 0.125 et 1.0. Cette déformation importante, imposée à
la structure, est illustrée (pour d = 1.0) à la la figure 3.11.
Le déplacement maximal pour ce cas, obtenu quand d = 1, est obtenu pour D = 0.036 (voir
section 3.5.5). Le pas de temps artificiel a, donc, été fixé à ∆t = (0.036)2 .
Les calculs ont été effectués à partir d’un maillage tétraédrique non structuré, ayant 30 571
éléments (6 692 noeuds et 40 335 arêtes). La qualité moyenne des éléments du maillage initial
est de 0.65. Encore une fois, il s’agit d’un maillage de faible qualité initiale.
Le tableau 3.7 résume les différents résultats. Comme on peut le voir, les méthodes de relève-
ments harmonique et élastique sont plutôt limités à de petites déformations (d < 0.375). Pour
le relèvement par krigeage, des éléments retournés (invalides) apparaissent pour d > 0.50.
87
Figure 3.11 – Cas #4 : Différents angles de vision du déplacement imposé à la structure pour
d = 1.0.
88
d Qualité Initiale (%) Hext (%) Eext (%) Kext (%) Pext (%)
0.125 invalide (0) (0) (0) (0) (0)
0.00 - 0.25 0.41 0.42 0.41 0.42 0.40
0.25 - 0.50 3.43 3.43 3.48 3.33 3.11
0.50 - 0.75 89.10 89.11 89.00 89.19 89.34
0.75 - 1.00 7.06 7.04 7.05 7.06 7.15
0.50 invalide (0) (2) (0) (0) (0)
0.00 - 0.25 0.41 — 0.45 0.44 0.42
0.25 - 0.50 3.43 — 3.82 3.79 3.61
0.50 - 0.75 89.10 — 88.92 88.80 88.91
0.75 - 1.00 7.06 — 6.81 6.97 7.06
1.0 invalide (0) (29) (11) (24) (0)
0.00 - 0.25 0.41 — — — 0.54
0.25 - 0.50 3.43 — — — 4.68
0.50 - 0.75 89.10 — — — 88.62
0.75 - 1.00 7.06 — — — 6.16
Table 3.8 – Cas #4 : Distribution de la qualité des éléments pour d = 0.125, 0.50 et 1.0. Le
nombre d’éléments retournés est indiqué entre parenthèses.
Table 3.9 – Cas #4 : Temps de calcul pour une itération. Tous les calculs ont été exécutés
sur un processeur Intel 64 bits de 3GHz typique
La méthode la plus efficace pour ce problème, est encore une fois, le relèvement parabolique
car il conduit à des maillages valides pour d = 1.0. Il est à noter que le relèvement parabolique
semble peu sensible à la piètre qualité initiale du maillage.
Le tableau 3.8 montre la distribution de la qualité d’éléments pour d = 0.125, 0.50 et 1.0.
Les résultats présentés montrent, encore une fois, que le relèvement parabolique est la seule
méthode capable de traiter les grands déplacements (pour d = 1.0).
Une dernière série de tests avec d = 1.25 a été effectuée. Des maillages invalides ont été
produits par toutes les méthodes. Toutefois, le nombre d’éléments retournés de chacune des
méthodes était en accord avec les résultats présentés dans le tableau 3.8. Dans ce cas, le nombre
d’éléments invalides était repectivement 40, 16, 28 et 2 pour les relèvements harmonique,
élastique, par krigeage et parabolique. Ceci sous-entend que le relèvement parabolique peut
produire des maillages valides pour les valeurs de d plus grandes que 1.
89
mémoire utilisée pour chaque méthode. Ceci est fait au tableau 3.9 pour le cas test #4 ; des
résultats semblables ont été obtenus pour les trois premiers cas (mais non présentés). Nous
n’avons pas indiqué le temps de calcul pour la méthode de krigeage, puisque notre implémen-
tation de la méthode n’est pas du tout optimale et la comparaison serait donc inappropriée.
Le relèvement parabolique produit des maillages valides à un coût très faible (comparable au
relèvement harmonique) et apparaît comme la meilleure méthode pour la majorité des cas.
À partir des études comparatives présentées, nous constatons que la méthode de relèvement
parabolique est la plus efficace, pour mettre à jour les maillages non structurés, dans la mesure
où elle semble peu sensible à la qualité du maillage. La méthode de relèvement parabolique
sera utilisée pour éliminer la nécessité de remaillage dans la solution des problèmes couplés de
fluide-structure. Les résultats numériques présentés plus tard illustrent cette idée.
Un système partitionné est dit faiblement ou fortement couplé. D’une part, un système parti-
tionné est faiblement (ou peu) couplé lorsque les conditions de couplage ne sont pas exactement
satisfaites à chaque pas de temps. D’autre part, un système partitionné est fortement couplé
lorsque les conditions de couplage sont exactement vérifiées c.à-d. elle sont appliquées avec
une grande précision à l’aide de sous-itérations, effectuées à chaque pas de temps entre les
deux solveurs. Le couplage faible cause des problèmes du point de vue numérique. En effet, le
développement d’algorithmes stables faiblement couplés est très difficile.
Dans Causin et al. (2005), un modèle simplifié d’IFS, représentant l’interaction entre un fluide
potentiel et un solide élastique linéaire, est étudié. Le non respect de la conservation à l’inter-
face fait apparaître un phénomène nommé «effet de masse ajoutée». Cet effet prend en compte
l’inertie virtuelle ajoutée à la structure, en raison de l’interaction avec le fluide. Dans le cas ex-
plicite, le couplage fluide-structure peut être condensé en une action supplémentaire de masse
sur la structure (d’où la terminologie «effet de la masse ajoutée»). Des difficultés numériques
s’en suivent et peuvent être expérimentées pour des valeurs croissantes de ce terme inertiel,
qui dépend du rapport de densité ρs /ρf et de la géométrie du domaine. Il en ressort une
instabilité des approches explicites. En résumé, l’effet de masse ajoutée est essentiel dans le
cadre d’un couplage explicite, ce qui n’est pas toujours le cas du couplage totalement implicite
(voir Fernández (2001)).
90
Naturellement, les schémas partitionnés ne sont pas à l’abri des difficultés liées au rapport de
densité. Pour plus de détails, le lecteur est référé à Causin et al. (2005) mais, cela ne nous
empêche pas de citer, brièvement, quatre détails essentiels liés à ce sujet :
— Pour une géométrie donnée, les instabilités numériques peuvent être expérimentées dans
les algorithmes de couplage partitionné dès que le rapport de densité ρs /ρf est inférieur
à un certain seuil.
— Pour un rapport de densité donné, les instabilités numériques peuvent être expérimen-
tées dans les algorithmes de couplage partitionné dès que la longueur du domaine est
supérieure à un certain seuil.
— Pour une géométrie donnée, une augmentation du nombre de sous-itérations peut être
observée dans les algorithmes partitionnés fortement couplés (schéma de couplage im-
plicite) dès que le rapport de densité ρs /ρf diminue.
— Pour un rapport de densité donné, une augmentation du nombre de sous-itérations peut
être observée dans les algorithmes à couplage implicite dès que la longueur des domaines
augmente.
Ainsi, les algorithmes de couplage partitionné implicite ne sont pas à l’abri des instabilités
numériques surtout si ρs ρf . Dans Causin et al. (2005), une étude approfondie de ce
phénomène a été étudiée. En effet, l’algorithme de couplage implicite ne convergera plus (ou
très mal) si ρs ρf parce qu’il y a une condition à respecter sur un facteur de relaxation.
Cette condition se traduit par :
2(ρs + aδt2 )
0<w< (3.35)
ρs + ρf µmax + aδt2
où w est le facteur de relaxation défini par
s ˆ s
d̂k+1 = wd̃sk+1 + (1 − w)d̂k
µmax est la plus grande valeur propre de la masse ajoutée et est une longueur caractéristique
du domaine.
Remarque :
— La condition 3.35 limite le facteur de relaxation avec une borne supérieure dépendant
des densités et du pas de temps. Dans le cas où w = 1, il suffit de développer un peu
plus la formule 3.35 et démontrer que ρf µmax − ρs < aδt2 . Donc, le schéma est contraint
à un certain seuil de δt et à des densités du solide et du fluide bien déterminées.
— Lorsque le rapport des densités solide-fluide est grand (c.à-d. ρs ρf ), la condition 3.35
n’est plus active, on peut utiliser la relaxation pour faire une extrapolation (accélération
de Aitken) ou la désactiver (w = 1).
— Si le rapport est petit (c.à-d. ρs ρf ), alors la borne devient «active» et le facteur de
relaxation peut être plus petit que 1 (dépend du pas de temps) : si le pas de temps tend
vers zéro, alors, le facteur w devient proportionnel au rapport des densités.
91
Pour notre cas d’étude, nous considérons un schéma implicite puisqu’il présente plus d’avan-
tages (convergence et stabilités numériques) que les deux schémas explicite et semi-implicite.
Étant données les applications considérées plus tard, la densité de la structure ρs est, en gé-
néral, supérieure ou égale à celle du fluide ρf , soit (ρs ≥ ρf ), par conséquent, nous notons
l’absence de l’effet de masse ajoutée décrit dans Causin et al. (2005). Les instabilités numé-
riques liées à cet effet sont aussi négligeables. Le lecteur peut se référer à Fernández et Gerbeau
(2009) pour plus de détails.
s s
ufn,0 = ufn−1 , pfn,0 = pfn−1 , d̂n,0 = d̂n−1 , p̂sn,0 = p̂sn−1 , ŵn,0 = ŵn−1 , Ωfn,0 = Ωfn−1
• Point fixe pour j ≥ 1
1.1 Résoudre les équations de Navier-Stokes sous forme ALE dans Ωfn,j−1 (équations
3.9) avec la condition à l’interface du fluide (3.112 ) par la méthode de projection ba-
sée sur un schéma incrémental de correction de pression 1.1.4 et calculer la solution
(ufn , pfn ) et donc σ fn,j ,
1.2 Résoudre les équations structurelles 3.27 avec la condition de Neumann appliquée
s
à la frontière du solide (3.113 ) et calculer la solution d̂n,j , p̂sn,j ,
1.3 Appliquer la méthode de relèvement parabolique (3.5.5) et résoudre :
f s
d̂n,j = Ext(d̂n,j |Σ̂ ) dans Ω̂f ,
1.6 Transformer les champs de vitesses et de pression du pas de temps précédent (un−1 ,
un−2 et pn−1 ) et la vitesse de maille ŵn,j sur le nouveau maillage Ωfn,j ,
s s
• Arrêt de point fixe si kd̂n,j − d̂n,j−1 k < ε , kufn,j −ufn,j−1 k < ε et kpfn,j −pfn,j−1 k < ε
92
2. Passer au pas de temps suivant.
Remarque :
— Dans nos calculs, nous avons choisi un critère d’arrêt du point fixe ε = 10−8 suite à
plusieurs tests numériques.
— Le passage des efforts s’effectue par une projection de toutes les contraintes sur tous
les points d’intégration des éléments ou encore les points de Gauss des éléments des
maillages fluide et solide.
— Bien que nous n’ayons pas eu besoin de le faire, nous notons que de nombreux travaux
ont été axés sur le développement de techniques pour accélérer la convergence dans les
schémas de couplage implicite. Un certain nombre de stratégies de point fixe ont été
proposées : une formule d’accélération Aitken est utilisée dans Wall et al. (1999) et Mok
et al. (2001) et des conditions aux limites spéciales de «transpiration» sont utilisées
dans Deparis et al. (2003) pour éviter le calcul des matrices de fluide dans chaque sous-
itération. Heil (2004) et Michler et al. (2005) ont proposé des techniques d’accélération
en utilisant des espaces de Krylov. Dans Raback et al. (2001), les auteurs proposent
une méthode de relaxation de la contrainte d’incompressibilité fluide afin d’accélérer la
convergence des itérations de point fixe. Plus récemment, des procédures de couplage fort,
basées sur un couplage Robin-Robin, ont également été étudiées, comme par exemple
dans Badia et al. (2008) pour les stratégies itératives du point fixe et dans Badia et al.
(2009) pour les méthodes de Krylov.
Une fois la formulation variationnelle du problème couplé et l’algorithme de résolution pré-
sentés, nous allons appliquer nos outils pour résoudre des problèmes classiques ainsi que de
nouveaux problèmes d’IFS.
93
Ce problème implique un écoulement incompressible autour d’un cylindre sur lequel une plaque
flexible est ajoutée sur le côté aval. Le cylindre induit des tourbillons de von Karman qui,
dans ce cas précis, font obligatoirement osciller la plaque flexible. Ce problème de référence
est présenté, ici, dans un but de validation. La géométrie du problème est décrite à la figure
3.12.
H h
A
`
Nous avons utilisé les mêmes dimensions que celles de Turek et Hron (2006) : H = 0.41 m,
L = 2.5 m, h = 0.02 m, et l = 0.35 m. Le cylindre de rayon de 0.05 m est centré à (0.2, 0.2).
À des fins de comparaison, nous considérons la position initiale d’un point A(0.6,0.2) tel
qu’illustré à la figure 3.12. Un profil parabolique est imposé à l’entrée du domaine fluide
y(H − y) 4
uf (0,y) = 1.5Ū H 2
= 1.5Ū y(0.41 − y)
(2) 0.1681
tel que la vitesse moyenne à l’entrée est Ū et la vitesse maximale de l’écoulement est 1.5Ū .
On impose une condition naturelle nulle en sortie, soit σ · n = (0,0) et des conditions de
non-glissement sur l’ensemble des bords du domaine fluide, y compris sur le cylindre.
Turek et Hron (2006) présentent différentes analyses de ce cas. Le premier cas est un test CFD
qui consiste à valider les valeurs des coefficients de traînée et de portance. Le deuxième cas est
un test CSM «Computational Structural Mechanics» qui consiste en une étude mécanique de
la structure. On a appliqué la force gravitationnelle sur la structure, dans un écoulement en
arrêt. Par ailleurs, on a étudié la déformation de la barre, en réponse à la force appliquée. Ces
deux premiers tests ont été effectués par application de notre algorithme 3.6 et on a validé les
différentes valeurs obtenues dans Turek et Hron (2006). On ne présente pas les résultats, car
ils sont moins pertinents dans le contexte actuel.
En outre, Turek et Hron (2006) proposent 3 cas d’étude d’IFS, mais nous nous intéressons aux
deux derniers cas instationnaires qui sont FSI2 et FSI3.
94
Le fluide est Newtonien incompressible et la structure est modélisée par un modèle hyperélas-
tique de Saint-Venant-Kirchhoff, tel que défini à la section 1.2.2. Dans ces deux cas tests, les
paramètres du problème sont :
Pour les deux cas d’étude, le maillage du fluide considéré (voir la figure 3.13) est constitué de
3 196 éléments triangulaires, soit 1 760 sommets et 4 956 arêtes. Nous adoptons des éléments
P2 −P1 . Notons que, dans ce cas, le nombre de degrés de liberté est comparable à celui proposé
par Turek et Hron (2006). Le maillage de la structure (voir la figure 3.14) est constitué de
742 éléments P2 triangulaires, soit 448 sommets et 1 189 arêtes. Nous avons considéré une
approximation en espace, par éléments finis, permettant ainsi des maillages (fluide et structure)
non-conformes.
Les simulations ont été faites avec les paramètres donnés dans Turek et Hron (2006) pour un
pas de temps fixé à dt = 0.001 et nous présentons, aux figures 3.15 et 3.16, les champs de
vitesses obtenus pour différents pas de temps dans les deux cas (FSI2 et FSI3).
95
Figure 3.15 – Champ de vitesses dans le cas test FSI2
96
Figure 3.17 – FSI 2 : déplacements en X et en Y du point A.
Pour caractériser davantage les résultats obtenus, nous comparons nos résultats avec ceux
obtenus par Turek et Hron (2006). À cet effet, nous présentons respectivement, aux figures
3.17 et 3.18, l’évolution du déplacement du point A dans le cas de FSI2 et FSI3. Comme on
peut le constater les résultats sont très satisfaisants.
Ce premier cas test nous a servi comme un test de validation. Les résultats numériques obtenus
montrent l’efficacité de l’algorithme 3.6.
97
Figure 3.18 – FSI 3 : déplacements en X et en Y du point A.
98
Pentrée
H
considéré pour cet exemple, est un rectangle (H = 1cm, L = 6cm) présenté à la figure 3.19. Le
fluide est incompressible laminaire et la structure est modélisée par un modèle hyperélastique
de Saint-Venant-Kirchhoff, tel que défini dans la section 1.2.2. Le fluide est initialement au
repos ; une pression a été imposée à l’entrée, pendant 0.005 secondes. Les paramètres physiques
du problème sont :
Notons que les paramètres, choisis pour ce cas test, sont des valeurs réalistes dans un contexte
d’écoulement sanguin. En fait, ce problème peut être interprété comme une interaction entre
le sang et une artère. On parle, aussi, de tension artérielle car cette pression est aussi la force
exercée par le sang sur la paroi des artères.
Les simulations ont été faites pour un pas de temps fixe dt = 0.0001. Pour les faire, on a pris en
considération, d’une part, un maillage de fluide constitué de 7 328 éléments P2 triangulaires :
soit 3 849 sommets et 11 176 arêtes présenté en noir à la figure 3.20 ; d’autre part, un maillage
de solide constitué de 1 780 éléments P2 triangulaires : soit 1 148 sommets et 2 926 arêtes
présenté en rouge à la figure 3.20. Notons que les deux maillages ne sont pas conformes (voir
figure 3.20).
Bien que cette étude soit qualitative, elle nous a permis de vérifier notre algorithme de couplage
d’IFS dans le cas en 2D. On présentera à la prochaine section des résultats pour un cas en 3D.
99
Agrandissement à l’interface fluide structure
Afin de faciliter les premières étapes des calculs d’interaction fluide-structure, nous ne com-
mençons pas à t = 0, mais plutôt nous effectuons un calcul de fluide seul de t = −2 s à t = 0.
Nous avons considéré une vitesse à l’entrée qui augmente suivant l’équation
1 t+1
Ventrée = sin π +1 .
2 2
L’interaction fluide-structure commence à t = 0.
L’écoulement du vent est modélisé par un fluide Newtonien alors que le «drapeau» d’épaisseur
0.0006 m est modélisé par un modèle de Saint-Venant Kirchhoff. Les paramètres physiques de
ce problème sont donnés dans le tableau suivant :
Les simulations ont été faites pour un pas de temps fixe dt = 0.01 par la prise en considération
de deux maillages pour le fluide et de deux autres pour la structure. Pour le fluide, nous avons
considéré, en premier lieu, un maillage grossier constitué de de 2 777 éléments P2 tétraédriques :
soit 559 sommets et 3 515 arêtes. En second lieu, nous avons considéré un maillage plus fin
100
Figure 3.21 – Propagation de l’onde de pression
101
Figure 3.22 – Géométrie et conditions aux limites
constitué de de 30 741 éléments P2 , soit 6 756 sommets et 40 633 arêtes. Pour la structure,
nous avons aussi considéré deux maillages : le premier est grossier et constitué de 345 éléments
P2 tétraédriques : soit 138 sommets, 617 arêtes alors que le deuxième maillage était plus fin
et constitué de 1428 éléments P2 : soit 556 sommets et 2 537 arêtes. Les différents maillages
utilisés pour les calculs sont présentés à la figure 3.23.
Il est assez difficile de choisir le résultat le plus pertinent pour les problèmes d’écoulements tri-
dimensionnels, et encore plus pour les problèmes d’interaction fluide-structure. Par conséquent,
afin d’obtenir une image qualitative des résultats calculés, nous présentons à la figure 3.24 la
déformée du drapeau pour 3 pas de temps différents soit t = 4.12s, t = 7.40s et t = 11.52s.
Nous présentons aussi, à la figure 3.25, le déplacement horizontal du point A de coordonnées
(0.1, 0.055, 0.07), montrant le mouvement du drapeau (voir figure 3.22). En raison de la
charge symétrique dans l’écoulement laminaire, la simulation de la structure ne présente pas
de déformation dans les premières 2.5 s. Après l’initiation de la formation de tourbillons, le
drapeau est déformé uniformément dans la direction y positive. Ce résultat est en accord avec
les résultats Von Scheven (2009).
Dans cette section, nous allons considérer deux exemples d’application de notre algorithme
d’IFS. Le premier exemple consiste à simuler numériquement l’interaction d’un pain de gomme
d’un pneu avec un écoulement laminaire. Le deuxième exemple consiste à simuler le phéno-
102
Figure 3.23 – Maillages : domaine fluide (gauche) et domaine solide (droite)
mène d’aquaplanage, c’est à dire l’interaction d’un pneu avec un écoulement laminaire, ce qui
correspondrait à une voiture se déplaçant très lentement.
Pain stationnaire
Dans un premier temps, nous avons étudié l’interaction d’un pain de gomme avec un fluide
incompressible. La géométrie et les conditions aux limites sont présentées à la figure 3.26.
Le maillage du pain de gomme a été fourni par notre partenaire industriel, soit la Manufac-
ture Française des Pneumatiques Michelin. Il s’agit d’un maillage constitué de 2 992 éléments
P2 hexaédriques, soit 3 876 sommets, 10 681 arêtes et 9 798 faces. La figure 3.27 présente la
géométrie et le maillage de peau du pain. Le maillage du fluide présenté à la figure 3.28 a,
aussi, été fourni et présente 114 414 éléments P2 tétraédriques.
103
t = 4.12s
t = 7.40s
t = 11.52s
Figure 3.24 – Drapeau déformé avec lignes de courant du fluide pour 3 pas de temps différents
104
Figure 3.25 – Déplacement horizontal du point A
Sortie σ · n =
(0,0,0)
pain
uf = (1,0,0)
150
uf = (1,0,0)
40
y
60
x
uf = (1,0,0) z
Entrée
105
Figure 3.27 – Géométrie et maillage de peau du pain gomme
Remarque : Notons que les deux maillages, fluide et structure, ne sont pas conformes puisque
celui du fluide est tétraédrique alors que celui du pain est hexaédrique. Les paramètres du
fluide et de la structure sont donnés par
Notons que les paramètres utilisés, dans cette partie, ne sont pas réalistes. En effet, le but de
ce test est de valider l’algorithme et le processus de couplage fluide-structure implémenté dans
le code MEF++ du GIREF (1996). Les simulations ont été effectuées pour un pas de temps
fixe dt = 0.1.
Nous observons, à la figure 3.29, la déformation du pain de gomme. Notons que ce problème
est stationnaire et que la norme euclidienne du déplacement maximal du pain de gomme
106
Figure 3.29 – Déformée du pain de gomme
est de l’ordre de 2.85 × 10−5 m. Nous présentons, à la figure 3.30, une coupe dans le plan
perpendiculaire au pain pour illustrer les vecteurs de vitesses dans ce plan.
Pour valider numériquement notre approche de mise à jour de la géométrie fluide, une deuxième
simulation a été effectuée pour l’écoulement autour du pain de gomme. Dans ce cas, nous
considérons, en premier lieu, un mouvement de rotation du pain dans le domaine fluide ; en
second lieu, un mouvement de translation verticale.
Le but de ce test est d’étudier la capacité de l’approche adoptée dans le cas d’une structure
rigide en mouvement. En effet, l’objectif principal de ce test consiste à déplacer le pain de
gomme dans le domaine du fluide avec le minimum de remaillage possible. Les paramètres
du fluide sont les mêmes utilisés dans la section précédente alors que pour la structure nous
avons augmenté les valeurs λs = 1012 kg/ms2 et µs = 0.4 × 1012 kg/ms2 pour rigidifier le pain.
Nous illustrons, à la figure 3.31, deux coupes du maillage du fluide présentant, deux positions
différentes du pain de gomme dans le domaine fluide. Comme on peut le constater dans cette
figure, une zone de recirculation du fluide se forme à l’intérieur du pain de gomme.
Ce cas test illustre bel et bien que notre approche de relèvement parabolique, définie à la
section 3.5.5, nous a permis d’effectuer les deux mouvements du pain de gomme. Notons que
107
Figure 3.30 – Champ de vitesse obtenu dans un plan perpendiculaire au pain de gomme
dans ce cas, seule la méthode de relèvement parabolique permet d’effectuer une mise à jour
de la géométrie fluide. Toutes les autres méthodes retournent des éléments et ne permettent
pas au pain de gomme de se déplacer verticalement, ni de tourner.
L’un des phénomènes les plus étudiés par les manufacturiers de pneus est l’aquaplanage. C’est
un phénomène qui se produit quand au moins une roue d’un véhicule perd de l’adhérence par
rapport au sol, en glissant complètement sur une surface aqueuse (flaques, chaussée mouillée,
etc.). Pour simuler, donc, ce phénomène, l’étude de l’interaction entre un fluide et un pneu
s’avère essentielle. Cette étude réalisée dans le cadre d’une collaboration avec Michelin est,
en effet, consacrée essentiellement à l’étude de l’influence des contraintes d’un écoulement
laminaire sur les pneumatiques.
Les premières études aérodynamiques, menées sur l’écoulement qui se développe autour des
roues de voitures, ont été réalisées au début des années 1970. L’expérience de référence est
celle de Fackrell et Harvey (1972) qui, à l’époque, fut la seule à considérer une roue en rota-
tion en contact complet avec un sol en translation. Fackrell et Harvey (1972) et Mears et al.
(2002) ont effectué des mesures de pression statique et des visualisations par fumée sur des
108
Translation verticale (du haut vers le bas) du pain de gomme
Figure 3.31 – Vecteurs vitesse présentés dans une coupe transversale du domaine fluide
109
roues présentant différents profils de pneumatiques et différentes largeurs. Leur objectif est de
fournir des données de référence pour une configuration relativement simple. Les roues étu-
diées présentaient toutes des surfaces lisses et non-déformées. Des essais menés sur différentes
configurations de roues fixes ou en rotation ont notamment permis d’obtenir la répartition des
coefficients de pression (portance) CP autour de la roue sur la ligne centrale (voir la figure
3.32).
Figure 3.32 – Répartition de CP sur la ligne centrale d’une roue en rotation, d’après Fackrell
et Harvey (1972) et Mears et al. (2002)
Les résultats de cette expérience, souvent pris comme référence, donnent des renseignements
importants sur la physique de l’écoulement se développant autour de la roue. Ils permettent,
en particulier, de mettre en évidence les conséquences principales de la rotation. Le point
d’arrêt en amont (CP = 1 autour de 0°) se trouve légèrement abaissé par rapport à la roue
fixe. Entre 0°et 90°, on observe une baisse de pression due à l’accélération de l’écoulement ;
puis, à l’approche de la surface de contact (second point d’arrêt à 90°), l’écoulement ralentit de
110
nouveau. Cependant, dans le cas de la roue en rotation, les CP obtenus dépassent amplement
la valeur de 1, limite théorique pour un écoulement sans apport d’énergie. Pour des angles
supérieurs à 90°, c’est-à-dire en aval de la zone de contact, la répartition de CP suit, ensuite, un
«palier» correspondant à la pression statique dans la zone massivement décollée à l’arrière de la
roue, seule zone, fortement instationnaire autour de la ligne médiane. La nouvelle augmentation
de CP marque l’apparition du décollement sur le haut de la roue.
À partir de ces données fournies par Fackrell et Harvey (1972) et Mears et al. (2002), nous
considérons un fluide incompressible laminaire (à bas Reynolds) et un pneu écrasé au sol
et comprenant 2 sillons longitudinaux sur la bande de roulement et plusieurs rainures tels
qu’illustrés à la figure 3.34.
111
Figure 3.34 – Géométrie du pneu considéré (fournie par Michelin)
Rivlin tel que présenté au chapitre 1. Les simulations ont été faites pour un pas de temps
fixe dt = 0.01. La vitesse du fluide, à l’entrée, est de 1.0 m s−1 c’est à dire 3.6 km h−1 . Le
sol est en translation avec la même vitesse de rotation de pneu, soit Ωr = 3.33 rad s−1 et la
hauteur de l’eau est de 1 cm. Les caractéristiques du fluide sont respectivement 1000 kg m−3
et 1,1768 kg m−3 pour les densités de l’eau et de l’air ; et de 10−3 P a s et 1.85 10−5 P a s,
pour les viscosités dynamiques dans l’eau et l’air.
À la figure 3.35, nous présentons les vecteurs vitesse de l’écoulement autour du pneu. Nous
constatons, d’une part, la présence du phénomène de «jetting», que nous rappelons sa définition
par l’apparition d’une forte surpression en amont de la zone de contact, décrit dans Fackrell
et Harvey (1972) et Mears et al. (2002) ; d’autre part, l’apparition d’une forte surpression en
amont du pneu. Ceci engendre une forte accélération en aval du pneu. Une représentation des
isolignes de vitesse de l’eau autour du pneu est présentée à la figure 3.36.
Le phénomène d’aquaplanage est caractérisé par la présence d’une pression dite «pression
112
u Magnitude
0 0.0816
113
d’aquaplanage» à la surface du pneumatique et la présence d’une zone de contact entre le pneu
et la réserve d’eau formée par l’accélération en aval du pneu. Lorsque l’effort vertical engendré
par cette pression devient supérieur à la masse du véhicule, le contact entre ce véhicule et
la route n’est plus maintenue et l’aquaplanage se produit : l’adhérence entre le véhicule et la
route est perdue et la trajectoire du véhicule n’est plus contrôlée. De la littérature, il est bien
connu que le lien entre la vitesse du véhicule, notée Vh et la pression d’aquaplanage, notée Ph
suit, au premier ordre, la relation (voir Tunnicliffe et al. (2006))
Ph = KVh2 (3.38)
La figure 3.37 illustre une coupe longitudinale au centre du pneu. À partir de cette figure,
nous constatons que les valeurs simulées de la pression maximale sur le pneu sont de l’ordre
de 512 P a et ceci est conforme à la loi présentée à la relation 3.38. En effet
Remarque :
On note que ce premier test a été effectué à des buts de validation du couplage fluide structure
et que les résultats obtenus confirment que notre algorithme fonctionne bien même dans les cas
tridimensionnels de grandes tailles. D’autres applications industrielles plus complexes seront
présentées plus tard.
Notre approche nous a permis de résoudre différents problèmes. Nous nous sommes assurés
du bon fonctionnement de notre algorithme à partir des comparaisons avec des cas tests tels
que celui de Turek et al. (2010). Nous avons aussi présenté d’autres résultats numériques en
2D et 3D.
Dans le prochain chapitre, nous allons effectuer une comparaison de cette approche ALE avec
la méthode des domaines fictifs présentée, au chapitre 2, dans le cas des structures rigides
immobiles et en mouvement.
114
Figure 3.37 – Pression autour du pneu
115
Chapitre 4
On a présenté dans les deux chapitres précédents, deux modèles pour la modélisation de
l’interaction fluide-structure. Le premier est celui par la méthode des domaines fictifs et le
deuxième est l’approche ALE. Une étude comparative des deux approches s’avère intéressante
pour valider numériquement et comparer l’efficacité des deux méthodes. Nous ferons l’étude
des écoulements des fluides présentant des objets rigides en mouvement. Bien sûr, nombreux
sont les cas où on peut appliquer ces méthodes et nous allons en étudier deux. Le premier
consiste en l’étude d’un écoulement d’un fluide autour d’un cylindre oscillant horizontalement.
Le deuxième cas d’étude de cette section est celui de l’écoulement autour du profil oscillant
présenté à la section 2.3.4.
x = 0.14D sin(2πf t)
c
2πfc
y=0
Ainsi, les conditions aux limites appliquées au fluide sur le cylindre sont : Ux = 0.14 D cos(2πfc t)
et Uy = 0, où l’amplitude de l’oscillation est égale à 0.14 et D est le diamètre du cylindre. Le
117
nombre de Reynolds a été fixé à 100 (Re = 100).
Nous avons appliqué l’approche des domaines fictifs présentée au chapitre 2. On rappelle
l’algorithme utilisé pour la résolution de ce problème.
Algorithme
1.3 Réinterpoler les champs de vitesse et de pression du pas de temps précédent (un−1 ,
un−2 et pn−1 ) sur le nouveau maillage Mtn (notons que cette étape peut se faire
seulement sur les nouveaux nœuds),
1.4 Appliquer les bonnes conditions aux limites au temps tn . En particulier, les condi-
tions aux limites de Dirichlet sont imposées à l’intérieur de l’objet rigide et à tous
les nœuds du bord, y compris les nouveaux nœuds introduits par la méthode de
découpage.
Notons que le découpage à chaque pas de temps se fait sur un maillage de baseMBase , c’est-à-
dire qu’après chaque résolution à un pas de temps tn−1 , on revient au maillage de base pour
le découper à l’instant tn . Dans le cas contraire, le maillage obtenu serait de plus en plus fin
à chaque pas de temps (et éventuellement complètement déformé) dans toute la région où
l’objet se déplace. Évidemment, si l’objet rigide est immobile, le découpage et le retournement
d’arêtes sont effectués une seule fois. Une fois le maillage Mtn créé (étape 1.3), les solutions
de vitesse et de pression des pas de temps précédents (un−1 , un−2 et pn−1 ) doivent être ré-
interpolées sur le nouveau maillage. Notez que la plupart des nœuds sont communs aux deux
maillages. Par conséquent, l’interpolation à ces nœuds est triviale.
118
4.1.2 Approche ALE
Le même problème a été étudié par l’approche ALE présentée au chapitre 3. Notons dans ce
cas qu’il n’y aura pas d’équations structurelles à résoudre puisque la structure est rigide. On
peut, donc, remplacer les équations 3.10 du problème couplé 3.9, 3.10 et 3.11 par l’application
de la condition aux limites Ux = 0.14 D cos(2πfc t). On présente, alors, l’algorithme utilisé
pour la résolution de ce problème
Algorithme
4.1.3 Comparaison
Nous pouvons voir, à la figure 4.1, les positions occupées par le cylindre, par rapport à la
coordonnée (x = 10, y = 10, z = 0) qui est la position du centre du cylindre à l’instant
t = 0. Dans les deux cas présentés, on conclut qu’il existe une instabilité développée à l’arrière
du cylindre qui est, en fait, l’allée de von Karman. Suite à de nombreux pas de temps pour
éliminer les effets transitoires, nous obtenons un régime périodique pour les deux approches
(voir figure 4.2).
Pour mieux valider les deux méthodes, on s’intéresse à une étude CFD en comparant certains
coefficients adimensionnels obtenus de cet écoulement. Nous définissons alors les coefficients
suivants :
— CDp est le coefficient de traînée dû à la seule contribution de la pression défini par
Z
p nx dA
CDp =
1 2
ρU D
2 moy
119
Approche IB Approche ALE
120
Figure 4.2 – Composante Uy en un point quelconque en arrière du cylindre.
Table 4.1 – Comparaison : coefficients de traînée et de portance pour les maillages grossiers
et (fins)
Notons que nous avons effectué le calcul sur des maillages plus raffinés (' 70 000 éléments
pour les deux cas) avec un pas de temps dt = 0.01. Les résultats sont présentés aussi dans
le tableau 4.1 et montrent que les valeurs des coefficients de traînée et de portance sont très
comparables à celles données par Rajani et al. (2009). Ainsi, Le raffinement des maillages
minimise l’erreur pour les 2 méthodes considérées.
121
Dans ce premier cas d’étude, les résultats obtenus dans les deux approches (IB et ALE)
concordent assez bien avec les données de la littérature. Le deuxième cas d’étude dans ce
chapitre sera le cas d’un écoulement d’un fluide autour d’un profil oscillant.
122
Remarque : Pour l’approche IB, notons que, dans la région du mouvement du profil, le
maillage doit être très fin pour pouvoir appliquer la méthode de découpage du maillage et
reproduire la frontière du profil.
Ici, encore, une attention particulière est portée au calcul du coefficient de portance qui est
donné par :
F2
CL = 1 2
2 ρU∞ c
où F2 est la force transversale appliquée par l’écoulement sur le profil. Le coefficient de portance
calculé sur une période complète est présenté à la figure 4.4 et comparé avec les résultats de Cori
(2011). Rappelons ici que la méthode utilisée par Cori (2011) est une méthode monolithique
de couplage d’IFS qui repose sur une formulation directe couplée à des intégrateurs en temps
d’ordre élevé (voir la section 3.1 pour plus de détails sur cette méthode). Comme on peut le
constater, l’accord est très satisfaisant.
CL Approche IB
Cori et al. 2011
Approche ALE
t
T
4.3 Conclusion
Les approches ALE et IB présentent des résultats numériques satisfaisants pour les deux cas
tests de solide rigide présentés dans ce chapitre. On note que pour les deux cas d’étude, les
123
résultats sont très comparables. Néanmoins, compte tenu du cas industriel traité dans cette
thèse et surtout la complexité de la structure géométrique du pneu, la méthode des domaines
fictifs ne sera pas utilisée pour le traitement de l’IFS dans le cas des structures flexibles. En
effet, il est très difficile d’appliquer la méthode de découpage et reproduire la frontière exacte
d’un pneu surtout avec la présence des sillons et des sculptures. On réfèrera au cas du profil
NACA, présenté à la figure 4.3, où il a été primordial d’avoir un maillage très fin pour pouvoir
appliquer la méthode de découpage du maillage et reproduire correctement la frontière du
profil. Ceci ne sera pas possible dans le cas du pneu vue la taille remarquable des problèmes
tridimensionnels et du coût du calcul. On note que l’approche ALE produit des résultats
similaires à ceux de l’approche IB et que les résultats présentés dans le cas des structures
flexibles concordent bien avec la littérature. Dans le prochain chapitre, nous étudions une
physique plus complexe du problème d’IFS où nous considérons un couplage multiphysique de
l’IFS avec les écoulements turbulents et à surfaces libres.
124
Chapitre 5
Couplage multiphysique
Jusqu’à présent, nous avons présenté des résultats d’interaction fluide-structure en ne considé-
rant que les écoulements laminaires. Nous allons entamer, maintenant, une étude des phéno-
mènes plus complexes tels que la modélisation de la turbulence et le traitement des écoulements
à surfaces libres. Ces deux phénomènes seront modélisés et combinés avec l’IFS pour présenter
un algorithme d’IFS plus complexe et plus représentatif des interactions réelles.
Ainsi, dans la plupart des applications industrielles de la mécanique des fluides, l’état de
l’écoulement est de nature turbulente. Ce régime turbulent apparaît lorsque le nombre de
Reynolds est suffisamment grand. Il se caractérise par la perte de tout élément d’organisation et
le développement de toute une hiérarchie de mouvements qui induisent un état nécessairement
tridimensionnel instationnaire.
125
généralement associées à des éléments structurels de vorticité). Les plus grands tourbillons sont
d’une dimension comparable à l’étendue de la région turbulente, dénotée l et appelée échelle
intégrale tandis que les plus petits sont de dimension telle que la viscosité devient un agent
effectif de transport. Elle diffuse, en conséquence, les gradients de vitesse. L’échelle de longueur
associée à la diffusion visqueuse, également associée à la dissipation d’énergie turbulente par
effet visqueux, est dénotée ici η et est généralement appelée échelle de Kolmogorov. Entre ces
deux extrêmes, l’écoulement turbulent contient une multitude de tourbillons de dimensions
caractéristiques intermédiaires. Le schéma de la figure 5.1 illustre ce concept.
Figure 5.1 – Cascade d’énergie idéalisée des échelles de la turbulence (tirée de Lemay (2011))
Le concept de cascade d’énergie stipule que pour les écoulements turbulents, à partir des
grandes échelles (échelles intégrales), les non-linéarités des équations de Navier-Stokes génèrent
des fluctuations d’échelles de plus en plus petites jusqu’à ce que la viscosité se dissipe et
transforme, en chaleur interne, les mouvements des plus petites échelles (échelles dissipatives,
voir 5.1). Cette cascade est représentée par une zone dite inertielle où le spectre d’énergie
présente une décroissance en -5/3 en fonction du nombre d’ondes. On rappelle bien sûr que
toute échelle spatiale l est reliée à un nombre d’ondes κ par l ∼ κ−1 .
126
Direct Numerical Simulation (DNS)
La méthode de simulation des grandes échelles dite LES consiste à ne résoudre qu’une partie
du spectre turbulent, les échelles les plus énergétiques, en introduisant un filtre dit de sous-
maille. Seules les petites échelles, dont le comportement est plus universel, sont modélisées si
bien qu’elle permettent de réduire les coûts de calcul en augmentant notamment la taille de la
maille nécessaire. Notons que la taille du filtre est liée à la finesse de la discrétisation spatiale.
Donc, plus le maillage est fin, plus la simulation est fidèle. En réduisant la gamme des échelles
résolues, l’intérêt majeur de la méthode LES est de pouvoir s’appuyer sur un maillage plus
grossier que la DNS, tout en garantissant un haut niveau de fidélité. Elle reste, cependant,
relativement coûteuse pour les problèmes industriels.
Les méthodes de simulations instationnaires URANS et stationnaires RANS sont basées sur
un principe proche de celui de la LES : réduire le coût des calculs en ne résolvant qu’une partie
de l’écoulement. Le plus souvent, on divise l’écoulement en mouvement moyen et mouvement
fluctuant ; chaque grandeur physique est représentée, suivant la règle de décomposition dite
de Reynolds, par la somme d’une quantité moyenne et d’une valeur fluctuante :
U = Ū + u et p = p̄ + p0
127
∂ Ū i ∂ Ū i 1 ∂ p̄ ∂ ∂ Ū i
+ Ū j + − ν + τij = 0 (5.2)
∂t ∂xj ρ ∂xi ∂xj ∂xj
où τij = −ui uj est le tenseur de Reynolds qui permet de modéliser l’effet de toutes les échelles
tourbillonnaires, constituant les parties fluctuantes, sur le mouvement moyen résolu. La demi-
trace du tenseur ui uj définit l’énergie cinétique turbulente :
1 1
k = ui ui = (u21 + u22 + u23 ) (5.3)
2 2
Les tensions de Reynolds apparaissent comme des inconnues supplémentaires qui font que le
système formé des équations 5.1, 5.2 soit un système ouvert. On a, alors, plus d’inconnues
qu’il n’y a d’équations pour résoudre les équations gouvernant l’écoulement moyen. Précisons
que le problème de fermeture se pose pour tous les écoulements turbulents. Les techniques qui
en permettent la résolution font appel à des hypothèses supplémentaires destinées à obtenir
un système comportant autant d’équations qu’il y a d’inconnues. Ces hypothèses, nommées
hypothèses de fermeture, ont pour rôle de fournir un moyen de déterminer les tensions de
Reynolds le plus adéquatement possible. La plupart des méthodes de fermeture des équations
du mouvement moyen (5.2) sont basées sur le concept de viscosité tourbillonnaire. Il s’agit
d’une démarche qui établit une relation directe entre les tensions turbulentes ui uj et les
vitesses moyennes Ū i . Cette hypothèse, introduite par Boussinesq en 1877, se traduit en
URANS par la relation :
2 ∂ Ū i ∂ Ū j
τij + kδij = νt + (5.4)
3 ∂xj ∂xi
Méthodes hybrides
Les méthodes hybrides RANS/LES sont généralement utilisées pour accélérer la convergence
et diminuer les coûts des calculs LES. Il s’agit, ici, d’utiliser la LES dans la partie du domaine
de calcul où le coût reste raisonnable et de la coupler à une approche de type RANS pour la
résolution des couches limites dont le coût est particulièrement élevé. La difficulté principale
réside, alors, dans le critère de transition entre les résolutions RANS et LES. Soulignons
de manière particulière la méthode DES (Detached Eddy Simulation) proposée par Spalart
(2009) qui traite les régions proches des couches limites d’une manière RANS et le reste de
l’écoulement, d’une manière LES.
Le choix du type de simulations utilisées dans le cadre de ce projet s’est porté sur l’approche
instationnaire URANS présentée en section 5.1.3. Il se justifie par des raisons à la fois phy-
siques, numériques et économiques :
128
— D’un point de vue physique, le sillage des corps émoussés tels que les roues de voitures
est le siège d’une instationnarité naturelle qu’il semble important de caractériser dans un
travail dont l’objectif principal est la compréhension de l’écoulement et des phénomènes
physiques mis en jeu. La résolution instationnaire apporte des informations temporelles
non seulement nécessaires à la description de la dynamique des sillages mais également
très utiles pour la compréhension des phénomènes d’interaction entre les tourbillons et
avec le pneumatique.
— D’un point de vue numérique, l’instationnarité naturelle de l’écoulement requiert l’utili-
sation de simulations aptes à reproduire les phénomènes instationnaires principaux tels
que celui des lâchers tourbillonnaires. L’utilisation de simulations stationnaires RANS
se heurte à des problèmes numériques dans le sens où il est impossible de converger vers
une solution moyenne sans l’ajout artificiel de viscosité numérique via un relâchement
du maillage ou par l’utilisation des schémas d’intégration dissipatifs. En outre, il est éta-
bli que, pour les écoulements naturellement instationnaires, les méthodes instationnaires
fournissent également des résultats moyennés d’une meilleure précision, notamment pour
la prévision des efforts aérodynamiques exercés sur l’objet.
— Enfin, dans un cadre industriel, le coût de simulations de type Direct Navier-Stokes
(DNS), Large-Eddy Simulations (LES) ou même hybrides RANS/LES reste encore pro-
hibitif. Les études CFD ont donc été restreintes à des cas de simulations URANS basées
sur des modèles de turbulence à 1 ou 2 équations.
Bien qu’il existe plusieurs modèles à 2 équations pour la modélisation de la turbulence, nous
nous limitons à l’application du modèle k − . Les équations moyennées de Reynolds, qui
modélisent les écoulements instationnaires de fluides incompressibles en régime turbulent avec
l’hypothèse de Boussinesq, s’écrivent :
∂uf
T
ρ
+ ρuf · ∇uf = −∇p + ∇ · [(µ + µT )(∇uf + ∇uf )] + ρf,
∂t
(5.5)
∇ · uf
= 0,
où ρ est la densité, u est le vecteur vitesse, p est la pression, f représente une force volumique
et µ définit la viscosité dynamique du fluide. La viscosité turbulente µT définie par :
k2
µT = ρ Cµ . (5.6)
est calculée à l’aide du modèle k − proposé par Launder et Spalding (1972). Le système est
complété par les équations de transport de l’énergie cinétique turbulente k, et par celles de
son taux de dissipation suivantes (voir Launder et Spalding (1972) pour plus de détails)
129
∂k
µT
ρ + ρuf · ∇ k = ∇ · µ + σk ∇ k + µT P (uf ) − ρ,
∂t
(5.7)
∂
µT 2
ρ + ρuf · ∇ = ∇· µ+ ∇ + C1 k µT P (uf ) − C2 ρ k ,
∂t σ
Les valeurs des paramètres utilisés C1 , C2 , Cµ , σ k , σ sont celles suggérées par Launder et
Spalding (1972) :
Cµ C1 C2 σk σ
0.09 1.44 1.92 1.0 1.3
Les modèles de turbulence à deux équations ont été et sont toujours largement employés
dans le calcul des écoulements turbulents. Nous savons que les deux équations de transport
permettent d’évaluer des variables de turbulence pour calculer, ensuite, le niveau de la viscosité
turbulente.
Par définition, les variables de turbulence sont toujours positives. Bien que les équations
de transport qui modélisent leur comportement admettent des solutions strictement positives
(voir Mohammadi et Pironneau (1993)), rien ne garantit que les solutions numériques le seront
aussi. En effet, numériquement les solutions en k et en peuvent prendre des valeurs négatives
ou nulles. Pour contourner ce problème, plusieurs approches ont été utilisées (voir par exemple
Kuzmin et Mierka (2006)). Dans ce travail, nous allons utiliser une formulation logarithmique
en k et .
Formulation logarithmique en k et
Pour aboutir à une formulation logarithmique des problèmes 5.7 en k et en , il faut d’abord
réécrire certains termes dans les équations de transport en utilisant la formule de la viscosité
tourbillonnaire 5.6 (voir par exemple Ilinca et Pelletier (1998) et Ilinca et al. (1998))
k2
− ρ = −ρ2 Cµ et C1 µT P (u) = ρC1 Cµ kP (u) (5.8)
µT k
130
Par combinaison des équations 5.7, 5.8 et 5.9, on obtient la formulation logarithmique suivante :
∂K
µT µT
ρ + ρuf · ∇ K = ∇ · µ + σk ∇K + µ + σk |∇K|2 + FK1 P (uf ) − FK2 ,
∂t
∂E
µT µT
+ ρuf · ∇ E = ∇· µ+ ∇E + µ + |∇E|2 + FE1 P (uf ) − FE2 ,
ρ
σ σ
∂t
(5.10)
où
eK eK
FK1 = µT e−K , FK2 = −ρ2 Cµ , FE1 = C1 µT e−K et FE2 = −ρ2 C2 .
µT µT
Plus de détails sur l’obtention de la forme logarithmique sont donnés par Wane (2012). La
résolution de K et E garantit que les variables k, et
uf
ψ̂K = ψK + τ f · ∇ψK ,
||u ||
uf
ψ̂E = ψE + τ f · ∇ψE
||u ||
où le paramètres τ dépend essentiellement de la taille de l’élément et de la diffusion. En
multipliant les équations 5.10 par ces fonctions tests, on obtient la formulation variationnelle
des équations 5.10 donnée aussi dans Wane (2012)
Z Z Z
∂K µT
ρ ψ̂K dv + ρuf · ∇ K ψ̂K dv = µ+ (∇K · ∇K)ψ̂K dv+
Ω ∂t Ω Ω σk
uf
XZ Z
µT µT
∇· µ+ ∇K τ f · ∇ψK dv − µ+ ∇K · ∇ψK dv (5.12)
K σk ||u || Ω σk
K∈Ω
Z
+ FK1 P (uf ) + FK2 ψ̂K dv
Ω
131
et
Z Z Z
∂E µ
ρ ψ̂E dv + ρuf · ∇ E ψ̂E dv = µ + T (∇E · ∇E)ψ̂E dv+
Ω ∂t Ω Ω σ
uf
XZ Z
µT µT
∇· µ+ ∇E τ f · ∇ψE dv − µ+ ∇E · ∇ψE dv (5.13)
K σ ||u || Ω σ
K∈Ω
Z
+ FE1 P (uf ) + FE2 ψ̂E dv
Ω
132
a. Résoudre les équations de Navier-Stokes 5.5 par la méthode de projection et calculer
uf et p
b. Résoudre l’équation 5.12 pour K et mise à jour de K
c. Résoudre l’équation 5.13 pour E et mise à jour de E
d. Mise à jour de µT
e. Retour à l’étape c si le critère n’est pas satisfait ou dans le cas contraire passer à
l’étape f
f. Retour à l’étape a si le critère de convergence globale n’est pas satisfait ou dans le
cas contraire, arrêter les calculs.
1+r 1−r σy
ux (x,y) = U1 + erf − ,
2 2 x
1−r 1 σy 2
uy (x,y) = U1 √ exp − ,
2 σ π x
p(x,y) = 0, (5.14)
σy 2
k(x,y) = k0 ck + exp − ,
x
0 σy 2
(x,y) = ck + exp − ,
x x
avec les valeurs de constantes :
10−4
U1 = 1.0 , r = 0.3 , σ = 13.5 , ck = , ρ = 1.0 , µ = 10−4 ,
k0
343 σ 343 σ2
k0 = U12 (1 − r) √ , 0 = Cµ U13 (1 − r)2
75 000 π 22 500 π
Les conditions aux limites de type Dirichlet sont appliquées pour toutes les variables sur toute
la frontière du domaine. La solution exacte, ci-dessus, est substituée dans les équations de
133
Navier-Stokes et de transport de la turbulence. Concernant les termes sources appropriés, ils
sont déterminés par l’utilisation du le logiciel Maple.
Ce problème a été étudié aussi dans la thèse de Wane (2012) pour la validation de l’adaptation
de maillage. Dans notre cas, pour vérifier la convergence de notre algorithme, nous avons
considéré deux maillages : le premier est un maillage grossier constitué de 4 032 éléments ; le
deuxième est un maillage plus fin, constitué de 15 130 éléments. Une convergence quadratique a
été observée pour toutes les variables. Nous présentons, à la figure 5.2, les différentes solutions
de ce problème.
Notons que les résultats fournis à la figure 5.2 concordent bien avec les résultats de Wane (2012)
et Turgeon et al. (2004). Cet exemple nous a permis de valider numériquement notre algorithme
de résolution des écoulements turbulents. À la prochaine section, nous allons présenter un autre
type de couplage des équations de Navier-Stokes, issu d’autres phénomènes physiques.
134
5.2 Modélisation des écoulements diphasiques
5.2.1 Écoulements à surface libre
On retrouve des problèmes à surfaces libres dans de nombreuses applications. On peut citer
par exemple le remplissage de moules dans l’industrie des polymères où on doit prédire la
position d’un front de matière qui doit remplir entièrement le moule. Lors de la coextrusion
de plusieurs polymères, il y a présence de plusieurs couches de polymères séparées par des
surfaces libres. Les surfaces libres sont présentes aussi dans la nature sous forme de gouttes
de pluie, de vagues sur l’océan, etc.
Nous avons considéré, dans ce travail, la méthode des surfaces de niveau («level-set») initiale-
ment introduite par Osher et Sethian (1988). Cette méthode est très populaire pour tous les
types de problèmes à surface libre stationnaire ou instationnaire. La résolution d’une équation
de convection permet alors de prédire les mouvements de l’interface dans un champ de vitesse
donné. Les travaux de Fortin et Benmoussa (2006) et Benmoussa et Fortin (2006) sur la dé-
formation de gouttelettes en cisaillement ou en élongation, que ce soit dans le cas 2D ou 3D
sont des exemples d’applications de ces méthodes.
Le signe permet de distinguer adéquatement les deux milieux fluides. La valeur de cette fonc-
tion ϕ est, ensuite, convectée avec le fluide. L’équation d’évolution s’écrit de la façon suivante :
∂ϕ
+ uf · ∇ϕ = 0 (5.16)
∂t
avec la condition initiale ϕ(x,y,0) = ϕ0 (x,y).
Les propriétés du fluide sont évaluées de la façon suivante : la fonction distance ϕ est négative
dans le fluide 1, positive dans le fluide 2 et nulle à l’interface. La densité et la viscosité du
135
fluide se traduisent par :
ρ(ϕ) = ρ1 + H(ϕ)(ρ2 − ρ1 ) (5.17)
où H est la fonction de Heaviside qui est nulle si ϕ < 0 et qui vaut 1 ailleurs. L’application
de cette méthode à la simulation d’écoulements à surfaces libres n’est pas immédiate. En
outre, les fonctions de Heaviside doivent être régularisées pour permettre la convergence de la
méthode numérique. En général, on remplace la fonction de Heaviside, sur un petit intervalle
au voisinage de zéro, par une fonction sinusoïdale donnée dans 5.19. Plus de détails sont donnés
dans les travaux de Benmoussa (2006).
0 si ϕ ≤ −,
1 ϕ 1 πϕ
H(ϕ) ' F (ϕ) = 1 + + sin( ) si − ≤ ϕ ≤ , (5.19)
2 π
si ≤ ϕ
1
L’écoulement dans ce cas se modélise par les équations de Navier-Stokes, qui s’écrivent en
suivant la surface libre :
f
∂u
ρ(ϕ) + u · ∇u − 2µ(ϕ)div(γ̇(uf )) + ∇p = ρ(ϕ)f dans Ωf (t)
f f
∂t
(5.20)
∇ · uf = 0 dans Ωf (t)
Il est bien connu que la solution d’une équation hyperbolique telle que l’équation 5.16 est
difficile et requiert une technique de stabilisation. Dans ce travail, on utilise la même technique
de stabilisation pour les équations de transport de k et de de la section 5.1.3. On considère
alors la méthode SUPG pour éviter les oscillations numériques et on définit les fonctions tests
SUPG notées par Φ̂
uf
Φ̂ = Φ + τ f · ∇Φ ,
||u ||
La formulation variationnelle de l’équation 5.16 est
Z
∂ϕ f
+ u · ∇ϕ Φ̂dv = 0 (5.21)
Ω ∂t
Des polynômes quadratiques P2 ont été utilisés pour la discrétisation en espace de cette équa-
tion hyperbolique 5.21 et une différence arrière du second ordre totalement implicite (O(∆t)2 )
a été utilisée pour la dérivée en temps. Si la solution ϕ(x,y,t(n) ) de l’équation 5.16 au temps
t = t(n) est notée par ϕ(n) , la descrétisation temporelle s’écrit
136
la condition initiale ϕ0 (x,y) est souvent choisie comme une distance signée. C’est la base de
la méthode des surfaces de niveau introduite par Sethian (1996). Cela signifie que l’interface
coïncide avec la surface de niveau zéro de la fonction de niveau ϕ0 et que si ϕ0 (x,y) = 1 le
point (x,y) est, alors, à une distance unitaire à partir de l’interface. Dans le cas contraire,
si ϕ0 (x,y) = −1 le point (x,y) est, donc, à une distance unitaire dans la direction opposée.
La construction de la fonction distance ϕ0 (x,y) est relativement facile dans la majorité des
applications.
L’équation 5.16 transporte l’interface à la bonne vitesse. Néanmoins, malgré le fait que la
fonction ϕ0 (x,y) est initialement une fonction distance, la solution ϕ(x,y,t) pour t > 0 ne
conservera pas nécessairement cette propriété. Pour de grands pas de temps, de forts gradients
apparaissent près de l’interface, provoquant des oscillations et des solutions irrégulières. Par
conséquent, la fonction de distance doit être réinitialisée périodiquement (voir Sussman et al.
(1994) pour plus de détails) .
Il est facile de vérifier que la fonction de distance doit satisfaire l’équation suivante :
Comme nous l’avons mentionné, la fonction ϕ(n) n’est pas nécessairement une fonction distance
et si elle n’est pas rénitialisée périodiquement, des oscillations apparaissent à l’interface. Pour
remédier à ce problème, on remplace ϕ(n) par une fonction distance notée ϕ̂, obtenue par la
résolution de l’équation différentielle suivante :
∂ ϕ̂
= sign(ϕ(n) ) 1 − |∇ϕ̂| (5.23)
∂ t̂
avec ϕ̂(x,y,0) = ϕ(n) (x,y). Dans cette équation, t̂ est un temps artificiel et sign(ϕ(n) ) est une
fonction signée définie par :
−1 si ϕ < 0,
(n)
sign(ϕ(n) ) = 0 si ϕ(n) = 0, (5.24)
si ϕ(n) > 0.
1
Remarques :
— Notons bien qu’à la convergence de l’équation 5.23 dans le cas stationnaire, |∇ϕ̂| = 1 ;
par conséquent, ϕ̂ est une fonction distance.
— Puisque la fonction sign(ϕ(n) ) n’est pas différentiable, on utilise en pratique une fonction
lisse définie par :
−1 si ϕ(n) ≤ −,
ϕ 1 πϕ
sign(ϕ(n) ) = + sin( ) si − ≤ ϕ(n) ≤ , (5.25)
π
si ≤ ϕ(n) .
1
137
Cette régularisation produit une fonction différentiable et une interface d’épaisseur 2.
En pratique, on choisit la valeur de faible, de l’ordre de 10−2 . Cette même régularisation
a été utilisée dans Sussman et Fatemi (1999). On référera le lecteur à Benmoussa (2006)
pour la résolution de l’équation 5.23.
Le système global non linéaire de résolution d’écoulement à surface libre est résolu par l’algo-
rithme suivant :
1. Initialiser les variables uf , p et ϕ
2. Jusqu’à la convergence de toutes les variables :
Remarque :
Dans l’algorithme précédent, les étapes 2a et 2b peuvent être résolues d’une manière différente.
En effet, nous pouvons transporter ϕ par le champ ũ qui est un champ à divergence non nulle
(voir 1.1.4) et effectuer, ensuite, une étape de correction sur le champ de vitesse u. Plus de
détails sur cette méthode sont donnés dans Deteix et al. (2014).
On s’intéresse, dans cette section, à la présentation des résultats numériques pour les écoule-
ments à surface libre. L’exemple présenté consiste en l’étude d’un jet d’un fluide incompressible
et a pour but de valider numériquement l’algorithme présenté à la section précédente. On trai-
tera ce problème dans le cas 2D et aussi dans le cas 3D.
Cas 2D
138
(2,4) (8,4)
uy = 0
σ·n=0
u = (0,0) y=1
(0,1)
u = (1 − y 2 ,0)
(0,0)
Figure 5.3 – Géométrie et conditions aux limites pour l’écoulement dans un jet.
fixée à = 0.25. La position initiale de l’interface entre ces deux fluides est donnée par
ϕ(y) = y − 1 = 0 (5.26)
La forme de l’interface à différents instants (par exemple, t = 4s, 8s, 12s, 30s) est tracée sur
la figure 5.4. On peut y voir la perturbation de l’interface qui se propage jusqu’à la sortie.
Notons que le diamètre d’un jet de fluide Newtonien à la sortie d’un écoulement de Poiseuille
est peu différent du diamètre du tube. À très faible Reynolds, on observe un gonflement
de l’ordre de 19%. Nous avons utilisé l’adaptation de maillage hiérarchique pour améliorer
la solution (voir Fortin (2001) et Bois (2012) pour plus de détails). Le maillage obtenu après
adaptation est donné à la figure 5.5. La figure 5.6 présente la solution obtenue après adaptation
de maillage. Nous avons obtenu un gonflement de 19%. Ce résultat est en accord avec ceux
obtenus par Mitsoulis (1999). Le gonflement s’explique ici, par le réarrangement du profil de
vitesse de l’écoulement purement Newtonien.
Cas 3D
Nous considérons maintenant le cas du jet en 3D. La géométrie et les conditions aux limites
sont données à la figure 5.7. Dans ce cas aussi, nous considérons les paramètres de l’eau et de
l’air, présentés à la section précédente. La position initiale de l’interface entre ces deux fluides
est donnée par
ϕ(y,z) = y 2 + z 2 − 1 = 0 (5.27)
Nous montrons à la figure 5.8 la forme de l’interface à différents instants (par exemple, t =
4s, 10s, 20s, 60s) sur une coupe suivant l’axe des z du domaine. Notons qu’à l’instant t=60 s,
l’écoulement est stabilisé et nous obtenons un gonflement d’environ 12%.
139
t= 4s t= 8s
t= 12s t= 30s
140
Figure 5.6 – Jet en 2D : solution finale
u = (1 − y 2 − z 2 ,0,0)
141
t= 4s t= 10s
t= 20s t= 60s
142
Encore une fois, nous avons utilisé l’adaptation de maillage hiérarchique pour améliorer la
solution. Le maillage obtenu après adaptation est donné à la figure 5.9. La figure 5.10 présente
la solution obtenue après adaptation de maillage.
Ce résultat est conforme à celui obtenu par Middleman et Gavis (1961) et Nickell et al. (1974).
En effet, Middleman et Gavis (1961) ont montré expérimentalement qu’un jet formé à partir
d’un fluide Newtonien pouvait connaître une expansion de son diamètre en fonction de la
valeur du Reynolds Re associée à l’écoulement au sein de l’injecteur. L’injecteur utilisé dans
leurs expériences est un tube de diamètre constant (de l’ordre de la centaine de microns)
et de quelques centimètres de long. Le nombre de Reynolds Re est défini ici à partir du
diamètre de l’injecteur et de la vitesse moyenne au sein de l’injecteur. Ils démontrent qu’à
faible Reynolds, un gonflement de l’ordre de 12% est observé en sortie. Ce gonflement s’explique
par la réorganisation du profil des vitesses qui impose généralement ∂uz
∂z 6= 0. L’équation de
continuité en coordonnées cylindriques implique alors
∂u u ∂uz
+ =− 6= 0
∂r r ∂z
D’où, l’existence d’une composante radiale de vitesse u, positive ou négative, qui modifie les
dimensions de l’écoulement et implique un gonflement en sortie.
143
Figure 5.10 – Jet en 3D : solution finale
Les deux tests présentés ici, mettent en évidence notre approche pour la résolution des écou-
lements à surfaces libres, que ce soit en 2D ou en 3D. Nous pouvons affirmer que l’algorithme
de résolution fonctionne bien et donne des résultats qui concordent avec ceux de la littérature.
Une fois les algorithmes 5.1.5 et 5.2.2 validés numériquement, nous en considérons le couplage
avec celui d’IFS.
5.3.1 Algorithme
144
s s
ufn,0 = ufn−1 , pfn,0 = pfn−1 , Kn,0 = Kn−1 , En,0 = En−1 , d̂n,0 = d̂n−1 , p̂sn,0 = p̂sn−1 , ŵn,0 =
ŵn−1 , Ωfn,0 = Ωfn−1
• Point fixe pour j ≥ 1
1.2 Résoudre les équations structurelles 3.27 avec la condition de Neumann appliquée
s
à la frontière du solide (3.113 ). et calculer la solution d̂n,j , p̂sn,j
1.3 Appliquer la méthode de relèvement parabolique (3.5.5) et résoudre :
f s
d̂n,j = Ext(d̂n,j |Σ̂ ) dans Ω̂f ,
1.6 Transformer tous les champs du pas de temps précédent (un−1 , un−2 , pn−1 , Kn−1
et En−1 ) et la vitesse de maille ŵn,j sur le nouveau maillage Ωfn,j (attribuer tous
les champs aux noeuds déplacés),
s s
• Arrêt de point fixe si kd̂n,j − d̂n,j−1 k < ε , kufn,j − ufn,j−1 k < ε , kpfn,j − pfn,j−1 k < ε
, kKn,j − Kn,j−1 k < ε et kEn,j − En,j−1 k < ε
Remarque :
Pour notre cas d’étude, nous considérons un schéma implicite puisqu’il présente plus d’avan-
tages (convergence et stabilités numériques). Cet algorithme de couplage partitionné implicite
145
n’est pas à l’abri des instabilités numériques surtout si ρs ρf (voir Causin et al. (2005) pour
plus de détails). Étant données les applications considérées plus tard, la densité de la structure
ρs est, en général, supérieure ou égale à celle du fluide ρf , soit (ρs ≥ ρf ) ; par conséquent,
nous notons que l’effet de masse ajoutée, décrit dans Causin et al. (2005), est négligeable.
Dans notre cas, le rapport de masse volumique permet d’éviter l’utilisation d’une relaxation
telle que décrit au chapitre 3.
La vitesse à l’entrée du canal a été fixée à U0 = 0.1 m/s, l’énergie cinétique turbulente initiale
est donnée par
k0 = 0.05 · U02
Les valeurs initiales de k et de ont été choisies telles que présentées dans Mahbubar et al.
(2007). Pour ce test, nous avons repris le maillage du fluide constitué de 3 196 éléments trian-
gulaires, soit 1 760 sommets et 4 956 arêtes. Nous adoptons des éléments P2 − P1 . Le maillage
de la structure est constitué de 742 éléments P2 triangulaires.
Les simulations ont été faites pour un pas de temps fixé à dt = 0.001 et nous présentons à la
figure 5.11 les champs de vitesses obtenus pour différents pas de temps.
146
Figure 5.11 – Champs de vitesses
Nous allons maintenant chercher à réaliser des simulations instationnaires turbulentes pour des
écoulements à surface libre en interaction avec une structure flexible en grandes déformations.
Pour ce faire, l’approche consiste à combiner l’algorithme de résolution des problèmes d’IFS
3.6 avec l’algorithme de résolution des écoulements turbulents 5.1.5 et aussi l’algorithme de
résolution des écoulements à surface libre 5.2.2. Par rapport à la section précédente, la diffé-
rence se situe au niveau de la viscosité et de la densité qui varient pour prendre en compte
les propriétés physiques des différents fluides (l’eau et l’air pour les problèmes qui nous inté-
ressent). On doit, de plus, résoudre une équation de transport supplémentaire pour ϕ pour
déterminer l’espace occupé par chaque milieu.
147
5.4.1 Algorithme
1. Pour un pas de temps donné t = tn :
• Initialisation des inconnues du problème d’IFS turbulent à surface libre :
s s
ufn,0 = ufn−1 , pfn,0 = pfn−1 , Kn,0 = Kn−1 , En,0 = En−1 , ϕn,0 = ϕn−1 , d̂n,0 = d̂n−1 , p̂sn,0 =
p̂sn−1 , ŵn,0 = ŵn−1 , Ωfn,0 = Ωfn−1
• Point fixe pour j ≥ 1
1.1
a. Calculer la distribution initiale de µT à l’aide de 5.11,
b. Résoudre les équations de Navier-Stokes 5.5 par la méthode de projection et calculer
uf et p
c. Résoudre l’équation 5.21 pour ϕ et mise à jour de ϕ
d. Résoudre l’équation 5.12 pour K et mise à jour de K
e. Résoudre l’équation 5.13 pour E et mise à jour de E
f. Mise à jour de µT
f
g. Calculer la solution (ufn ,pfn ) et donc σn,j
1.2 Résoudre les équations structurelles 3.27 avec la condition de Neumann appliquée
s
à la frontière du solide (3.113 ). et calculer la solution d̂n,j , p̂sn,j
1.3 Appliquer la méthode de relèvement parabolique (3.5.5) et résoudre :
f s
d̂n,j = Ext(d̂n,j |Σ̂ ) dans Ω̂f ,
148
5.4.2 Résultats numériques
Les simulations présentées ici illustrent la validité du couplage entre les différents opérateurs
mis en jeu dans cette interaction fluide-structure. Il est à noter que nous n’avons pas réussi à
trouver des études mettant en jeu de telles simulations dans la littérature. Pour cette raison,
nous considérons le problème turbulent tridimensionel d’une structure flexible en interaction
avec l’eau et l’air. La géométrie considérée est décrite à la section 3.7.3. La vitesse à l’entrée
du canal a été fixée à U0 = 0.5 m/s, l’énergie cinétique turbulente initiale est k0 = 0.05U02 et
la dissipation d’énergie intiale est 0 = k01.5 /0.01. La position initiale de l’interface entre l’eau
et l’air est donnée par
ϕ(z) = z − 0.055 = 0
Les simulations ont été faites pour un pas de temps fixe dt = 0.01 en considérant un maillage
fluide constitué de 30 741 éléments P2 et un maillage de la structure, constitué de 1 428
éléments P2 aussi. Ces maillages sont présentés à la figure 3.23. Nous présentons à la figure
5.12 la déformée du drapeau pour 4 pas de temps différents soit t = 1.12s, t = 3.40s, t = 5.21s
et t = 7.03s.
5.5.1 Aquaplanage
Par son effet lubrifiant, la présence éventuelle d’eau sur la chaussée possède une forte influence
sur les performances d’adhérence du pneumatique. La route y contribue par sa pente, son
dévers, sa texture et sa structure (enrobés drainants, par exemple). Le pneumatique y contribue
149
t = 1.12s
t = 3.40s
t = 5.21s
t = 7.03s
Figure 5.12 – Structure flexible déformée avec lignes de courant du fluide pour 4 pas de temps
différents
150
notamment par la forme de l’aire de contact, la forme de ses sculptures, et sa sensibilité aux
efforts de pression dans l’aire de contact.
Dans l’aire de contact d’un pneumatique sur une chaussée mouillée, on peut distinguer trois
mécanismes qui apparaissent progressivement de l’avant à l’arrière (voir figure 5.13) :
— à l’entrée du contact, l’eau est progressivement évacuée par les sculptures du pneuma-
tique et les aspérités de la chaussée. La hauteur d’eau est supérieure à 0,5 mm et des
phénomènes hydrodynamiques se produisent. Dans cette zone appelée zone de pénétra-
tion, l’adhérence est quasi-nulle. C’est généralement dans cette zone que prend nais-
sance l’hydroplanage. L’hydroplanage, communément appelé aquaplaning, est une perte
d’adhérence due à une quantité d’eau importante sur la chaussée ainsi qu’une vitesse
élevée du véhicule circulant sur la chaussée. Sous l’effet de la vitesse de déplacement,
le pneumatique génère une mise en pression de l’eau. Lorsque cette pression devient
supérieure à la pression moyenne d’appui du pneumatique sur la chaussée (de l’ordre de
2 bars pour une voiture), le pneu ne peut plus repousser l’eau et il se retrouve soulevé
(voir Michelin (2001) pour plus de détails).
— à l’arrière de la zone de pénétration se trouve une zone de transition. Une partie de l’eau
ayant été évacuée à l’avant du contact, le pneumatique entre progressivement en contact
avec certaines aspérités du revêtement de la chaussée. La hauteur d’eau est comprise
entre quelques micromètres et 0,5 mm.
— à l’arrière de l’aire de contact se trouve la zone de contact humide ou sèche où peut
persister la présence d’un film d’eau résiduel discontinu. C’est souvent dans cette zone
que se manifeste le viscoplanage. Il s’agit une perte d’adhérence se déroulant dans la
zone de transition où la hauteur d’eau, présente entre le pneumatique et la surface de la
chaussée, est très faible. Ce film d’eau interrompt les liaisons moléculaires entre la gomme
et le revêtement de la chaussée empêchant le phénomène d’adhésion de se produire. La
présence de ce film d’eau résulte de l’existence d’efforts de cisaillement visqueux dans le
film d’eau.
L’augmentation de la vitesse génère, tout d’abord, une projection d’eau qui se propage dans
les sculptures et sous le pneu. De ce fait, la force de levage sur le pneu augmente, ce qui
réduit le contact avec le sol. Ensuite, la friction de l’écoulement dans et autour du pneu ainsi
que la turbulence jouent un rôle important. La turbulence influence le pneu en le déformant.
Tout écoulement devient turbulent au-dessus d’une certaine vitesse. Finalement, les aspects
d’écoulement multiphases eau-air ont un impact non négligeable sur l’aquaplanage.
151
Figure 5.13 – Aire de contact entre un pneumatique et la chaussée en présence d’eau (tiré
de Michelin (2001)).
Les étapes de transfert des données entre le code fluide et le code structure se déroulent de la
manière suivante (voir la figure 5.14) :
La procédure de couplage est basée sur un algorithme implicite. Cet algorithme permet le
passage des efforts par une projection de toutes les contraintes sur tous les points d’intégration
des éléments ou encore les points de Gauss des éléments des maillages fluide et solide. Ceci
constitue une difficulté majeure, étroitement liée aux maillages fluide et solide.
Les techniques développées dans ce travail prennent en compte le calcul de l’interaction tel que
présenté, le caractère turbulent du fluide avec une surface libre et le comportement dynamique
de la structure en grandes déformations, avec les contributions exposées dans les chapitres
précédents. Le cas de l’aquaplanage nécessite, en outre, le traitement simultané du contact
152
Figure 5.14 – Transfert des données entre les codes fluide et solide.
qui sort du cadre de ce travail. Dans ce contexte, nous présentons un cas dont le contact est
absent. Il s’agit schématiquement d’une pompe entraînant un film d’eau.
La configuration du problème choisi est illustrée par la figure 5.15. L’écoulement d’eau possède
une épaisseur initiale de 8 mm ; il est turbulent et possède une surface libre avec l’air situé au
dessus. Au centre du domaine, un cylindre déformable en grandes déformations est entraîné
en rotation par son axe central.
L’une des difficultés principales est le maillage du domaine fluide. En effet, la peau du solide
doit coîncider géométriquement avec le domaine fluide c.à-d. qu’aucun espace n’est toléré entre
les domaines fluide et solide. Avec les outils disponibles dans GMSH (voir Remacle et Geuzaine
(2009)), nous avons élaboré un maillage très fin dans la zone occupée par l’eau et un autre
plus grossier dans la zone de l’air. La figure 5.16 illustre cette idée.
Notons que le maillage, ainsi obtenu, est constitué de 1 242 822 éléments P2 tétraédriques, soit
322 771 sommets et 1 779 014 arêtes.
Le fluide est Newtonien incompressible modélisé par un modèle URANS (voir 5.1.3) et la
structure est modélisée par un modèle hyperélastique de Mooney-Rivlin incompressible (voir
1.2.2). La vitesse à l’entrée U0 est égale à 22.88 m/s c.à-d. 80 km/h (on impose à l’entrée une
condition de Dirichlet sur le vecteur vitesse U = (22.8, 0, 0)), le nombre de Reynolds basé sur
la hauteur d’eau initiale et la vitesse à l’entrée est de 1.83 × 106 .
153
Figure 5.15 – Géométrie du domaine fluide.
154
L’énergie cinétique turbulente initiale est donnée par
k0 = 0.05 · U02
Les simulations ont été effectuées pour un pas de temps fixé à dt = 0.0001. Sur le sol, une
vitesse longitudinale est imposée, cohérente avec la vitesse de rotation de la structure, soit Ωr =
76 rad s−1 ; la hauteur d’eau est de 8 mm. Les caractéristiques du fluide sont respectivement :
1000 kg m−3 et 1,1768 kg m−3 pour les densités de l’eau et de l’air. Quant aux viscosités
dynamiques dans l’eau et l’air, elles sont respectivement 10−3 Pa s et 1.85 10−5 Pa s.
155
À la figure 5.17, nous présentons les vecteurs vitesses de l’écoulement autour de la structure.
Nous constatons l’apparition d’une forte surpression en amont, qui engendre une forte accé-
lération en aval (voir figure 5.18). Cette forte accélération est présentée à la figure 5.19 où
nous présentons les valeurs de la vitesse du fluide sur une ligne de coupe dans le domaine (voir
figure 5.19).
On observe aussi un bourrelet d’eau à l’aval, avec un effet de pompe accélératrice de l’écoule-
ment (voir figure 5.21).
Pour valider davantage notre simulation, nous avons calculé la portance exercée par le fluide
sur la structure. La valeur trouvée est de 470 N. Ce résultat rejoint les résultats obtenus par
Wane (2012).
156
u Norme de u
Figure 5.19 – Courbe de vitesse suivant une ligne de coupe dans le fluide (eau).
157
P
P
158
Figure 5.21 – Coupe longitudinale (y = 0)
159
Conclusion
Les modèles d’interaction fluide-structure développés dans cette thèse font appel à des mé-
thodes développées à partir d’éléments finis. Nous nous sommes concentrés sur l’analyse nu-
mérique et le développement d’algorithmes efficaces pour la résolution de ces problèmes.
Bien qu’il existe une multitude de méthodes et de modèles pour résoudre l’IFS, nous avons
considéré deux méthodes de résolution selon le type de structure considérée : structure rigide
ou structure flexible. Dans le cas des structures rigides, une nouvelle approche d’étude de
l’interaction fluide-structure est présentée. Nous avons adopté une méthode de domaines fictifs
bien qu’il existe plusieurs stratégies pour aborder ce problème. L’approche est basée sur une
méthode d’éléments finis qui impose avec précision les conditions aux limites à l’interface
fluide-solide. L’approche domaine fictif permet d’éviter le remaillage complet du domaine de
calcul, en présence de variations géométriques ou de changement de positions des obstacles.
Son principal avantage réside dans le fait qu’un un seul maillage est généré et par la suite
modifié d’une manière simple.
Nous avons aussi développé une statégie adaptative permettant de rajouter sur un maillage
quelconque des noeuds à l’interface fluide-obstacle pour mieux appliquer les conditions aux
limites. Nous avons fait plusieurs tests de la méthode sur différents problèmes dans le cas
où les obstacles rigides sont fixes et/ou mobiles. Notre approche a été testée et comparée
avantageusement avec les résultats de la littérature sur plusieurs problèmes tests. Les données
de la littérature et les courbes simulées concordent assez bien et ont été publiées dans Jendoubi
et al. (2014).
Dans le cas des structures flexibles, la méthode arbitraire lagrangienne-eulérienne (ALE), les
équations de Navier-Stokes et de la dynamique d’un solide élastique nous ont aidé à établir le
modèle avec les dernières approches issues des méthodes d’éléments finis. Les géométries du
fluide et du solide sont considérées comme deux domaines complémentaires qui se déforment
au fur et à mesure que le fluide et le solide agissent l’un sur l’autre. Sans modifier la topologie
intrinsèque du maillage de chacun de ces deux domaines, les noeuds peuvent subir un mou-
vement propre sans passer par un remaillage systématique. Dans ce cadre, on a perfectionné
la gestion de la déformation du maillage due aux grands déplacements de la structure via le
relèvement parabolique (voir Jendoubi et al. (2016) pour plus de détails). Plusieurs investiga-
161
tions ont été réalisées et des résultats intéressants ont été discutés. Après avoir présenté les
solveurs numériques utilisés pour résoudre les domaines fluides et solides ainsi que la méthode
de mise à jour du maillage, un algorithme totalement implicite et partitionné de couplage
d’IFS est présenté. Des résultats de diverses simulations numériques impliquant des problèmes
d’IFS sont présentés et il est démontré que la méthode est fiable, même pour des problèmes
d’interaction en grandes déformations en 3D.
Nous avons présenté d’une part, une comparaison de la méthode des domaines fictifs avec l’ap-
proche ALE dans le cadre d’obstacles rigides ; d’autre part, des tests de validation permettant
d’apprécier la performance et la fiabilité des deux méthodes sont présentés.
Nous nous sommes aussi intéressés à la simulation numérique de phénomènes réels d’inter-
action entre un fluide et un pneumatique. Pour cela, une stratégie de complexification du
problème d’IFS a été définie. La modélisation de la turbulence et des écoulements à surfaces
libres a été introduite. Nous avons fixé notre choix de modèle en considérant les équations
de Navier-Stokes moyennées basées sur le modèle URANS. En présence d’une surface libre,
nous avons introduit la méthode des surfaces de niveau. Ainsi, nous avons présenté plusieurs
algorithmes de couplage multi-physique et nous avons développé des algorithmes en langage
C++ dans le code générique MEF++ du GIREF. Nous avons montré la robustesse des mé-
thodes de couplage développées sur des résultats numériques, présentant un intérêt pour des
applications industrielles. Les applications traitées ont montré la possibilité de construire un
modèle d’interaction fluide-structure industriel réaliste. Les résultats obtenus sont comparés
avantageusement avec ceux de la littérature.
Pour conclure cette thèse, nous présentons certains aspects de ce travail offrant des perspectives
de développement :
— La modélisation du fluide serait plus réaliste en prenant en compte des modèles dits
de seconde génération comme par exemple le modèle R.S.M. ou des modèles de type
multi-échelles comme par exemple la simulation des gros tourbillons L.E.S. pour la mo-
délisation de la turbulence. En effet, l’avantage des modèles basés sur l’hypothèse de
type Boussinesq, comme par exemple le modèle URANS utilisé dans cette thèse, est la
facilité de mise en oeuvre numérique. De façon générale, si les résultats obtenus ne sont
pas toujours parfaits, ils sont rarement très fortement erronés. Si ce concept est bien
adapté au calcul de couches limites, il existe cependant des écoulements pour lesquels
l’hypothèse de viscosité tourbillonnaire n’est pas applicable.
— Nous n’avons pas utilisé d’adaptation de maillage dans nos calculs d’IFS, mais nous
sommes persuadés que cela augmenterait beaucoup la précision des résultats. Passer
à des maillages adaptés permet d’obtenir de nouveaux maillages permettant de mieux
coller aux phénomènes physiques modélisés. Ceci, en effet, permettera de réduire les
temps de calcul et d’améliorer la précision des résultats numériques.
162
— Nous avons considéré, dans ce travail, une approche lagrangienne totale pour la résolution
des équations structurelles. Puisque dans certains cas, le solide peut subir de grands
déplacements, il sera intéressant d’adopter une méthode de Lagrangien actualisé, pour
permettre plus de flexibilité dans la gestion des mouvements complexes de la structure.
Cette méthode utilise le maillage déformé au temps précédent pour générer le nouveau
maillage sur la configuration déformée. Il s’agit d’une approche incrémentale qui permet
d’exprimer les variables cinématiques et statiques par rapport à la configuration déformée
antérieure. L’avantage de cette description par rapport à l’approche lagrangienne totale
est qu’on peut remailler des configurations intermédiaires pour avoir plus de précision
sur les solutions calculées. Plus de détails à ce sujet sont élaborés dans Léger et al.
(2014).
— Nous avons considéré dans ce travail, un couplage implicite d’IFS. Il s’agit d’un schéma
fortement couplé et il est, en général, beaucoup plus coûteux en temps de calcul que la
somme des temps de calcul du fluide et du solide, car la résolution itérative du problème
de point fixe demande de résoudre les problèmes fluide et solide à chaque itération.
La stabilité est, en général, assurée sous des conditions moins contraignantes que pour
les schémas explicites ; cependant, ces propriétés rendent la simulation de problèmes
d’interaction fluide-structure souvent trop coûteux pour des problèmes réalistes. C’est
pour rendre de telles simulations possibles qu’il serait également intéressant de passer à
des schémas semi-implicites.
— L’étude de nouvelles géométries de pneumatiques, incluant par exemple différents types
de sculptures transverses, nécessitera l’utilisation d’autres moyens numériques en termes
de maillage et de conditions aux limites afin de prendre en compte la rotation et la
déformation des éléments.
— La prise en compte de nouvelles caractéristiques des roues (pneumatiques et jantes). La
complexification de la géométrie des jantes avec la prise en compte de leur ouverture
côté intérieur est susceptible d’agir sur l’écoulement interne au passage de la roue. Cela
ajouterait un degré de réalisme aux simulations et augmenterait leur précision et leur
coût de calcul.
— Un axe de recherche reste encore à approfondir pour de futures simulations numériques :
la prise en compte des problèmes d’élasticité en grandes déformations avec le phénomène
de contact. En effet, nous avons exclu dans cette étude les phénomènes mécaniques
entrant en jeu dans un contact roue-sol en milieu naturel. Lorsqu’un véhicule évolue
sur un terrain qui peut se déformer, sa dynamique est influencée par la nature et le
profil du sol. Afin d’avoir un diagnostic plus fiable des situations à risques (phénomène
d’aquaplannage), le modèle doit prendre en compte les déformations du sol et de la
roue. L’enfoncement de la roue dans le sol engendre des forces résistives susceptibles de
modifier le comportement du pneumatique.
163
Bibliographie
Amestoy, P., Guermouche, A., l’Excellent, J.-Y. et Pralet, S. (2006). Hybrid sche-
duling for the parallel solution of linear systems. Parallel Computing, 32(2):136–156.
Amestoy, P. R., Duff, I. S., Koster, J. et L’Excellent, J.-Y. (2001). A fully asynchronous
multifrontal solver using distributed dynamic scheduling. SIAM Journal on Matrix Analysis
and Applications, 23(1):15–41.
Arnold, D. N., Brezzi, F. et Fortin, M. (1984). A stable finite element for the Stokes
equations. Calcolo, 21:337–344.
Bach, P. et Hassager, O. (1985). An algorithm for the use of the Lagrangian specification in
Newtonian fluid mechanics and applications to free-surface flow. Journal of Fluid Mechanics,
152:173–190.
165
Bank, R. E. et Xu, J. (1996). An algorithm for coarsening unstructured meshes. Numerische
Mathematik, 73(1):1–36.
Batina, J. T. (1990). Unsteady Euler airfoil solutions using unstructured dynamic meshes.
AIAA Journal, 28(8):1381–1388.
Benmoussa, K. (2006). Une méthode d’éléments finis adaptative pour les problèmes à surfaces
libres instationnaires. Thèse de doctorat, Université Laval, Québec, Canada.
Blom, F. J. (2000). Considerations on the spring analogy. International Journal for Numerical
Methods in Fluids, 32(6):647–688.
Boffi, D., Brezzi, F. et Fortin, M. (2013). Mixed Finite Element Methods and Applications,
volume 44 de Springer Series in Computational Mathematics. Springer, Berlin, Heidelberg.
Bois, R., Fortin, M. et Fortin, A. (2012a). A fully optimal anisotropic mesh adaptation
method based on a hierarchical error estimator. Computer Methods in Applied Mechanics
and Engineering, 209-212:12–27.
Bois, R., Fortin, M., Fortin, A. et Couët, A. (2012b). High order optimal anisotropic
mesh adaptation using hierarchical elements. European Journal of Computational Mecha-
nics/Revue Européenne de Mécanique Numérique, 21(1-2):72–91.
Bonnet, D. M. et Frangi, A. (2007). Analyse des solides déformables par la méthode des
éléments finis. European Journal of Computational Mechanics, 16(5):667–668.
166
Borazjani, I. et Sotiropoulos, F. (2008). Numerical investigation of the hydrodynamics of
carangiform swimming in the transitional and inertial flow regimes. J. Exp. Biol., 211(10):
1541–1558.
Bottasso, C., Detomi, D. et Serra, R. (2005). The ball-vertex method : a new simple
spring analogy method for unstructured dynamic meshes. Comput. Meth. Appl. Mech.
Eng., 194(39-41):4244–4264.
Brackbill, J. et Saltzman, J. (1982). Adaptive zoning for singular problems in two dimen-
sions. J. Comput. Phys., 46(3):342–368.
Brezzi, F. et Fortin, M. (1991). Mixed and Hybrid Finite Element Methods, volume 15 de
Springer Series in Computational Mathematics. Springer-Verlag, New York.
Burman, E. et Fernández, M. (2007). Stabilized explicit coupling for fluid structure inter-
action using Nitsche’s method. C. R. Acad. Sci. Paris, Ser. I, 345:467–472.
Cahouet, J. et Chabard, J. (1988). Some fast 3D finite element solvers for thegeneralized
Stokes problem. International Journal for Numerical Methods in Fluids, 8.
Causin, P., Gerbeau, J.-F. et Nobile, F. (2005). Added-mass effect in the design of partitio-
ned algorithms for fluide-structure problems. Comput. Mth. App. Mech. Engrg., 194:4506–
4527.
Chiandussi, G., Bugeda, G. et Onate, E. (2000). A simple method for automatic update
of finite element meshes. Commun. Numer. Meth. Engng., 16:1–19.
Chorin, A. (1968). Numerical solution of the Navier–Stokes equations. Math. Comp., 22:745–
762.
167
Cottet, G.-H. et Maitre, E. (2004). A level-set formulation of immersed boundary methods
for fluid-structure interaction problems. C. R. Acad. Sci. Paris, 338:581–586.
Cottet, G.-H. et Maitre, E. (2006). A level-set method for fluid-structure interactions with
immersed surfaces. Mathematical Models and Methods in Applied Sciences, 16:415–438.
Dagan, A. (2003). Numerical consistency and spurious boundary layer in the projection
method. Computers & Fluids, 32(9):1213 – 1232.
de Boer, A., Bijl, H. et Van Zuijlen, A. (2005). Comparing different methods for the
coupling of non-matching meshes in fluid–structure interaction computations. AIAA Paper,
(AIAA-2005-4620).
De Boer, A., Van der Schoot, M. et Bijl, H. (2007). Mesh deformation based on radial
basis function interpolation. Comp. & Struct., 85:784–795.
Deteix, J., Jendoubi, A. et Yakoubi, D. (2014). A coupled prediction scheme for solving the
Navier–Stokes and heat equations. SIAM Journal of Numerical Analysis, 52(5):2415–2439.
Donea, J., Fasoli-Stella, P. et Giuliani, S. (1977). Lagrangian and Eulerian finite element
techniques for transient fluid-structure interaction problems. In Trans. 4th Int. Conf. on
Structural Mechanics in Reactor Technology, San Francisco.
168
Dunne, T. (2007). Adaptive finite element approximation of fluid-structure interaction based
on Eulerian and arbitrary Lagrangian-Eulerian variational formulations. Thèse de doctorat,
University of Ruprecht Karls, Heidelberg.
Fackrell, J. et Harvey, J. (1972). The Flow Field and Pressure Distribution of an Isolated
Road Wheel. IC aero report.
Farhat, C., Degand, C., Koobus, B. et Lesoinne, M. (1998). Torsional springs for two-
dimensional dynamic unstructured fluid meshes. Computer Methods in Applied Mechanics
and Engineering, 163:231–245.
Fernández, M., Formaggia, L., Gerbeau, J.-F. et Quarteroni, A. (2009). The deri-
vation of the equations for fluids and structure. In Formaggia, L., Quarteroni, A. et
Veneziani, A., éditeurs : Cardiovascular Mathematics, volume 1 de MS&A, pages 77–121.
Fortin, A. et Garon, A. (2013). Les éléments finis : de la théorie à la pratique. 400 pages,
Notes de cours.
Fortin, M. (2001). Anisotropic mesh adaptation through hierarchical error estimators, pages
53–65. Nova Science Publishers, Inc., Commack, NY, USA.
169
Franca, L. et Frey, S. (1992). Stabilized finite element methods : II. the incompressible
Navier–Stokes equations. Computer Methods in Applied Mechanics and Engineering, 99(2):
209 – 233.
Glowinski, R., Pan, T.-W., Hesla, T.and Joseph, D. et Periaux, J. (1999). A distribu-
ted Lagrange multiplier/fictitious domain method for flows around moving rigid bodies :
application to particulate flows. International Journal for Numerical Methods in Fluids,
30:1043–1066.
Glowinski, R., Pan, T.-W., Hesla, T., Joseph, D. et Periaux, J. (2000). A fictitious
domain approach to the direct numerical simulation of incompressible viscous flow past
moving rigid bodies : Application to particulate flow. J.Comput.Phys., pages 169–363.
Glowinski, R., Pan, T.-W. et Periaux, J. (1994). A fictitious domain method for external
incompressible viscous flow modeled by Navier-Stokes equations,. Comput. Methods Appl.
Mech. Engrg, 112:133–148.
Goda, K. (1979). A multistep technique with implicit difference schemes for calculating two-
or three-dimensional cavity flows. J. Comput. Phys, 30:76–95.
Gong, D., Chang, J. et Wei, C. (2011). An adaptive method for choosing center sets of
RBF interpolation. Journal of Computers, 6(10):2112–2120.
Guermond, J.-L. (1997). Un résultat de convergence d’ordre deux pour l’approximation des
équations de Navier-Stokes par projection incrémentale. C. R. Acad. Sci. Paris Sér. I Math.,
325(12):1329–1332.
170
Guermond, J.-L. (1999). Un résultat de convergence d’ordre deux en temps pour l’approxi-
mation des équations de Navier-Stokes par une technique de projection incrémentale. M2AN
Math. Model. Numer. Anal., 33(1):169–189.
Heil, M. (2004). An efficient solver for the fully coupled solution of large-displacement fluid-
structure interaction problems. Computer Methods in Applied Mechanics and Engineering,
193(1-2):1 – 23.
Ilinca, F. et Hétu, J.-F. (2011). A finite element immersed boundary method for fluid flow
around rigid objects. International Journal for Numerical Methods in Fluids, 65(7):856–875.
Ilinca, F., Hétu, J.-F. et Pelletier, D. (1998). A unified finite element algorithm for
two-equation models of turbulence. Computers & Fluids, 27(3):291–310.
Ilinca, F. et Pelletier, D. (1998). Positivity preservation and adaptive solution for the k-
turbulence model. AIAA J., 36(1):44–50.
Jendoubi, A. (2010). Interaction fluide-structure par la méthode des domaines fictifs. Mé-
moire de D.E.A., Département de génie mécanique, Université Laval.
171
Jendoubi, A., Deteix, J. et Fortin, A. (2016). A simple mesh-update procedure for moving
boundary flows and fluid–structure interaction. Accepté à Computers & Structures.
Karniadakis, G., Israeli, M. et Orszag, S. (1991). High-order splitting methods for the
incompressible Navier–Stokes equations. Journal of Computational Physics, 97(2):414 – 443.
Krige, D. (1951). A statistical method for mine variation problems in the Witwatersrand.
Journal of Chemistry and Metallurgy of the Mining Society of South Africa, 52:119–139.
Lin, T. J., Guan, Z. Q., Chang, J. H. et Lo, S. H. (2014). Vertex-ball spring smoothing :
An efficient method for unstructured dynamic hybrid meshes. Computers & Structures,
136:24–33.
172
Lions, J.-L. (1969). Quelques méthodes de résolution des problèmes aux limites non linéaires.
Dunod.
Léger, S., Fortin, A., Tibirna, C. et Fortin, M. (2014). An actualized Lagrangian for-
mulation with automatic remeshing for large deformation problems. International Journal
for Numerical Methods in Engineering, 100(13):1006–1030.
Matheron, G. (1973). The intrinsec ramdom functions and their applications. Advances in
Applied Probability, 5:439–468.
Mears, A., Dominy, R., Sims-Williams, D. W. et Wohlmuth, B. (2002). The air flow
about an exposed racing wheel. Rapport technique, Durham University.
Mok, D., Wall, W. et Ramm, E. (2001). Accelerated iterative substructuring schemes for
instationary fluid-structure interaction. Computational Fluid and Solid Mechanics, pages
1325–1328.
Mouro, J. et LeTallec, P. (2001). Fluid structure interaction with large structural displa-
cement. Comput. Meth. Appl. Mech. Eng, 190:3039–3067.
173
Nickell, R. E., Tanner, R. I. et Caswell, B. (1974). The solution of viscous incompressible
jet and free surface flows using finite element methods,. J. Fluid Mech., 65:189–206.
Noor, D. Z., Chern, M.-J. et Horng, T.-L. (2009). An immersed boundary method to solve
fluid-solid interaction problems. Comput Mech, 44:447–453.
Patankar, N. A. et Sharma, N. (2005). A fast projection scheme for the direct numerical
simulation of rigid particulate flows. Communications in Numerical Methods in Engineering,
21(8):419–432.
Pederzani, J. et Haj-Hariri, H. (2006). A numerical method for the analysis of flexible bo-
dies in unsteady viscous flows. International Journal for Numerical Methods in Engineering,
68(10):1096–1112.
Petcu, M., Temam, R. et Wirosoetisno, D. (2004). Existence and regularity results for
the primitive equations in two space dimensions. Commun. Pure Appl. Anal., 3(1):115–131.
174
Raback, P., Ruokolainen, J., Lyly, M. et Järvinen, E. (2001). Fluid-structure interac-
tion boundary conditions by artificial compressibility. In ECCOMAS Computational Fluid
Dynamics Conference, Swansea, UK.
Raviart, P. (1981). Les méthodes d’éléments finis en mécanique des fluides. Eyrolles.
Rendall, T. C. S. et Allen, C. B. (2009a). Efficient mesh motion using radial basis functions
with data reduction algorithms. Journal of Computational Physics, 228:6231–6249.
Rugonyi, S. et Bathe, K. (2001). On finite element analysis of fluid flows fully coupled with
structural interactions. Comput. Model.Simulat.Eng, 2:195–212.
Sethian, J. (1996). Level Sets Methods and Fast Marching Methods. Numéro 3 de Cam-
bridge Monograph on Applied and Computational Mathematics. Cambridge University
Press, Cambridge.
Shen, J. (1996). On error estimates of the projection methods for the Navier-Stokes equations :
second-order schemes. Math. Comp., 65(215):1039–1065.
175
Stein, K., Tezduyar, T. et Benney, R. (2004). Automatic mesh update with the solid-
extension mesh moving technique. Comput. Meth. Appl. Mech. Eng., 193:2019–2032.
Su, S.-W., Lai, M.-C. et Lin, C.-A. (2007). An immersed boundary technique for simulating
complex flows with rigid boundary. Computers and Fluids, 36:313–324.
Sussman, M., Smereka, P. et Osher, S. (1994). A Level Set Approach for computing
solutions to incompressible two-phase flows. J. Comp. Phys., 114:146–159.
Tezduyar, T., Aliabadi, S.and Behr, M., Johnson, A., Kalro, V. et Litke, M. (1996).
Flow simulation and high performance computing. Comput. Mech., 18:397–412.
Turek, S. (1999). Efficient Solvers for Incompressible Flow Problems, volume 6 de Lecture
Notes in Computational Science and Engineering. Springer-Verlag, Berlin.
Turek, S., Hron, J., Razzaq, M., Wobker, H. et Schäfer, M. (2010). Numerical bench-
marking of fluid-structure interaction : A comparison of different discretization and solution
approaches. In Bungartz, H.-J., Mehl, M. et Schäfer, M., éditeurs : Fluid Structure
Interaction II, volume 73 de Lecture Notes in Computational Science and Engineering, pages
413–424. Springer Berlin Heidelberg.
176
Turgeon, É., Pelletier, D. et Borggard, J. (2004). A general continuous sensibility
equation formulation for the κ − formulation. Int. J. Comput. Fluid dynamics, 18(1):29–
46.
Valette, R., Bruchon, J., Digonnet, H., Laure, P., Leboeuf, M., Silva, L., Vergnes,
B. et Coupez, T. (2007). Méthodes d’interaction fluide-structure pour la simulation multi-
échelles des procédés de mélange. Mécanique et Industries, 8:251–258.
van Kan, J. (1986). A second-order accurate pressure-correction scheme for viscous incom-
pressible flow. SIAM J. Sci. Statist. Comput., 7(3):870–891.
Van Loon, R., Anderson, P., de Hart, J. et Baaijens, F. P. (2004). A combined fic-
titious domain/adaptative meshing method for fluid-structure interaction in heart valves.
International Journal for Numerical Methods in Fluids, 46:533–544.
Van Loon, R., Anderson, P., van de Vosse, F. et Sherwin, S. (2007). Comparaison of va-
rious fluid-structure interaction methods for deformable bodies. Computers and Structures,
85:833–843.
Wall, W. A., Mok, D. et Ramm, E. (1999). Partitioned analysis approach of the transient
coupled response of viscous fluids and flexible structures. Solids, structures and coupled
problems in engineering, proceedings of the European conference on computational mechanics
ECCM, 99.
Wane, B. (2012). Adaptation de maillages et méthodes itératives avec application aux écou-
lements à surfaces libres turbulents. Thèse de doctorat, Département de mathématiques et
de statistiques, Université Laval, Québec, Canada.
Wang, H. et McLay, R. (1986). Automatic remeshing scheme for modeling hot forming
process. J. Fluid. Eng., 108:465–9.
177
Wendland, H. (2005). Scattered Data Approximation (1st edn). Cambridge Monographs on
Applied and Computational Mathematics. Cambridge University Press, Cambridge.
Wick, T. (2011). Fluid-structure interactions using different mesh motion techniques. Com-
put. Struct., 89(13-14):1456–1467.
Wood, C., Gil, A., Hassan, O. et Bonet, J. (2008). A partitioned coupling approach for
dynamic fluid-structure interaction with applications to biological membranes. International
Journal for Numerical Methods in Fluids, 57(5):555–581.
Zienkiewicz, O. et Taylor, R. (1989). The finite element method : basic formulations and
linear problems, volume 1. McGraw-Hill, London, 4 édition.
178