Couplage Fluide Structure pour la simulation numérique

des écoulements fluides dans une conduite à parois

rigides ou élastiques, en présence d’obstacles ou non.
Lahcen Ait Moudid

To cite this version:

Lahcen Ait Moudid. Couplage Fluide Structure pour la simulation numérique des écoulements
fluides dans une conduite à parois rigides ou élastiques, en présence d’obstacles ou non.. Sciences
de la Terre. Université d’Artois, 2007. Français. <tel-00341094>

HAL Id: tel-00341094

[Link]
Submitted on 24 Nov 2008

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Université d’Artois
Faculté des Sciences Jean Perrin

Sujet de thèse :

Couplage Fluide-Structure pour la simulation numérique des

écoulements fluides dans une conduite à parois rigides ou élastiques, en
présence d’obstacles ou non.

Thèse présentée par :

Lahcen AIT MOUDID

Évaluée par un jury composé des pérsonnes suivantes :

M. Aziz Hamdouni Professeur, Univ. La Rochelle Rapporteur

Mme. Elisabeth Longatte HDR, EDF Chatou Rapporteur
M. Brahim Khelifa Professeur, Univ. Artois Président
M. Mhamed Souli Professeur, USTL Lille1 Examinateur
M. Christian Mathieu Professeur, Univ. Artois Examinateur
M. Hakim Naceur HDR, UTC Compiègne Examinateur
M. Andrew Parry Dr. Ingénieur, Schlumberger Clamart Examinateur
M. Pasquale Mammone Professeur, Univ. Artois Codirecteur
M. Abdellatif Ouahsine Professeur, UTC Compiègne Directeur de thèse

Thèse acceptée le : 24 octobre 2007

 REMERCIEMENTS

Cette thèse a été effectuée à la fois au sein de l’Université d’Artois et aussi en

collaboration avec l’Université de Technologie Compiègne. Je tiens tout d’abord à
remercier Abdellatif Ouahsine, Profeseur à l’université de Compiègne qui a accepté
de diriger mes recherches. Ses grandes qualités humaines et scientifiques ainsi que
sa disponibilité ont été essentielles à la réalisation de cette thèse.
Je remercie Madame elisabeth Longatte, HDR à EDF Chatou, et Aziz Hamdouni,
Professeur à l’université de La Rochelle, pour avoir accepté de rapporter sur ma
thèse. Je remercie aussi Monsieur Christian Mathieu, Professeur à l’université
d’Artois, pour sa rigueur et pour sa présence aux moments difficiles. Je suis très
honoré que Monsieur Brahim Khelifa, Professeur à l’université d’artois, ait accepté
de présider ce jury.
Mes vives remerciements vont également à Monsieurs les Professeurs Mhamed
Souli, Andrew Parry et Pasquale Mammone pour leur aide, les multiples dicus-
sions très fructueuses et pour avoir acceptés de juger mon travail. Je remercie
également Monsieur Hakim Naceur, Profeseur à l’université de Compiègne, pour
son aide et ses conseils. Je l’en remercie sincèrement pour l’intérêt qu’il a manifesté
pour mon travail et pour avoir accepté de le juger.
Mes remerciement vont également à toute ma famille, et en particulier à mes par-
ents, mes frères, mon épouse et ma fille Salma qui m’ont toujours soutenu dans les
moments difficiles.
Je remercie également Monsieur Adlanne Sayede, Maı̂tre de conférences à l’université
d’Artois, pour son aide et ses conseils. Que mes amis, Sami Kaidi, Bouras Raneme...
soient ici vivement remerciés pour leur fructueuse collaboration et leur aide matérielle
pour la réalisation de ces travaux.
 RÉSUMÉ

La simulation numérique de l’interaction fluide-structure par la méthode des éléments

finis a été étudiée dans le cadre des équations de Navier-Stokes pour un fluide
visqueux newtonien incompressible en interaction avec un solide élastique.
La formulation Euler-Lagrange Arbitraire (ALE) a été utilisée, en considérant un
maillage dynamique, où le solide est décrit par une formulation lagrangienne et
le fluide par une formulation eulérienne. Le modèle fluide est validé en con-
sidérant des cas tests académiques et concernent les cas : de la marche, de la
cavité, de l’écoulement autour d’un cylindre, etc... Le modèle solide est validé en
considérant le cas d’une poutre encastrée-libre et encastrée-encastrée, le cas d’un
cylindre soumis à son poids propre et le cas d’une arche elliptique.
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 explicite, permet
le transfert de champs de façon interactive. L’efficacité de la méthode ALE et du
couplage fluide-structure a été évaluée en considérant plusieurs cas tests: solide
immergé dans un canal, où s’écoule un fluide en écoulement transitoire ou station-
naire, écoulement d’un fluide dans une canalisation à parois élastiques, etc...
Les résultats montrent que ce couplage explicite-interactif permet d’utiliser un
maillage et un schéma différent pour le fluide et la structure, et jouit de l’avantage
de ne pas utiliser de grosses matrices de stockage des données.
Mots-clés :
Fluide, Equation de Navier-Stokes, Structure, Eléments finis de coque, formulation
ALE, interaction fluide-structure.
 ABSTRACT

The numerical simulation of fluid-structure interaction by the finite element method

has been studied in the context of The Navier-Stokes equations for an incompress-
ible Newtonian viscous fluid in interaction with an elastic solid.
The formulation Arbitrary Lagrangian-Eulerian (ALE) was used, by considering a
dynamic grid, where the solid is described by a Lagrangian formulation and the fluid
by a Eulerian formulation. The fluid model validation is carried out by considering
academic tests and concerns: flow around a step, the Cavity, the flow around a
cylinder... The solid model is validated by considering the case of a cantilever
beam and fixed beam, the case of a cylinder subjected to its actual weight and the
case of an elliptical arch.
An algorithm of coupling is then developed for the implementation of this inter-
action fluid-structure. This algorithm, based on an explicit schema, allows the
transfer of fields in an interactive way. The effectiveness of the method ALE and
the coupling fluid-structure was evaluated by considering several cases tests: solid
immersed in a channel, where runs out a transitory or stationary fluid flow, flow of
a fluid in a drain with elastic walls...
The results show that this coupling explicit-interactive makes it possible to use a
different mesh and schemas for the fluid and the structure, and enjoys the advantage
of not using large matrices of data storage.
Keywords :
Fluid, Navier-Stokes equation, Structure, shell finite element model, ALE formu-
lation, fluid-structure interaction.
 vi

TABLE DES MATIÈRES

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

CHAPITRE 1: ETAT DE L’ART . . . . . . . . . . . . . . . . . . . . . 4

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Fluide-structure et problèmes associés . . . . . . . . . . . . . . . . . 6
1.3 Méthodes de couplage Fluide-structure . . . . . . . . . . . . . . . . 7
1.3.1 Méthodes numériques . . . . . . . . . . . . . . . . . . . . . . 7
1.3.2 Méthodes analytiques . . . . . . . . . . . . . . . . . . . . . . 10
1.4 Formulation ALE des équations de Navier-Stokes . . . . . . . . . . 11

CHAPITRE 2: RÉSOLUTION DES ÉQUATIONS D’ÉCOULEMENT

DE FLUIDES PAR ELÉMENTS FINIS . . . . . . . 13
2.1 Introduction : Méthode des éléments finis . . . . . . . . . . . . . . 13
2.2 Eléments d’analyse fonctionelle . . . . . . . . . . . . . . . . . . . . 14
2.3 Mise en oeuve de la méthode des éléments finis . . . . . . . . . . . . 16
2.3.1 Formulation variationnelle forte . . . . . . . . . . . . . . . . 16
2.3.2 Formulation variationelle faible . . . . . . . . . . . . . . . . 16
2.3.3 Conditions aux limites . . . . . . . . . . . . . . . . . . . . . 17
2.4 Les équation du mouvement des fluides incompressible . . . . . . . 18
2.4.1 Équation de continuité . . . . . . . . . . . . . . . . . . . . . 18
2.4.2 Conservation de la quantité de mouvement . . . . . . . . . . 18
2.5 Ecriture variationnelle faible dans le cas de fluides incompressibles . 19
2.5.1 Discrétisation par éléments finis . . . . . . . . . . . . . . . . 20
2.6 Méthodes de stabilité numérique . . . . . . . . . . . . . . . . . . . . 22
 vii

CHAPITRE 3: MODÈLE FLUIDE : TRAITEMENT NUMÉRIQUE

ET VALIDATION . . . . . . . . . . . . . . . . . . . . . 24
3.1 Résolution numérique . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.1.1 Algorithmes de résolution : Discrétisation Temporelle . . . . 24
3.1.2 Linéarisation : Méthode itérative de Newton-Raphson . . . 26
3.1.3 Stabilité de la solution : Condition de Ladysenskaya-Brezzi-
Babuska (LBB) . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.4 Méthode SUPG (Streamline Upwind Petrov-Galerkin) . . . . 29
3.2 Validation du module fluide . . . . . . . . . . . . . . . . . . . . . . 30
3.2.1 Cas test de Poiseuille . . . . . . . . . . . . . . . . . . . . . . 31
3.2.2 Cas test d’écoulement autour d’un cylindre . . . . . . . . . . 33
3.2.3 Cas test d’une Cavité carrée . . . . . . . . . . . . . . . . . . 40
3.2.4 Cas test d’écoulement autour d’un obstacle . . . . . . . . . . 43

CHAPITRE 4: MODÈLE SOLIDE : TRAITEMENT NUMÉRIQUE

ET VALIDATION . . . . . . . . . . . . . . . . . . . . . 49
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
4.2 Synthèse bibliographique . . . . . . . . . . . . . . . . . . . . . . . . 50
4.3 Description des éléments finis utilisés dans le modèle structure . . . 51
4.3.1 Eléments de coques de type Kirchhoff discret . . . . . . . . . 51
4.3.2 Eléments de coques de type Reissner-Mindlin . . . . . . . . 56
4.3.3 Elément à deux dimention en coordonnées cartisiennes, Elément
triangulaire à trois noeuds (T3) . . . . . . . . . . . . . . . . 58
4.4 Résultats numériques : tests de validation . . . . . . . . . . . . . . 62
4.4.1 Modèles éléments finis . . . . . . . . . . . . . . . . . . . . . 62
4.4.2 Panneau cylindrique soumis à son poids propre. . . . . . . . 63
4.4.3 Poutre 2D encastré-libre, soumise à un effort concentré : . . 69
 4.4.4 Poutre 2D en flexion encastré-encastré. . . . . . . . . . . . . 72
4.4.5 Déformation d’une arche elliptique . . . . . . . . . . . . . . 75

CHAPITRE 5: COUPLAGE FLUIDE-STRUCTURE . . . . . . . . 78

5.1 Formulations eulériennes, lagrangiennes et mixtes . . . . . . . . . . 78
5.1.1 Formulation eulérienne . . . . . . . . . . . . . . . . . . . . . 78
5.1.2 Formulation lagrangienne . . . . . . . . . . . . . . . . . . . 79
5.1.3 Formulation ALE . . . . . . . . . . . . . . . . . . . . . . . . 80
5.2 Procédure du couplage . . . . . . . . . . . . . . . . . . . . . . . . . 86
5.2.1 Conditions de couplage . . . . . . . . . . . . . . . . . . . . . 86
5.2.2 Algorithme de couplage . . . . . . . . . . . . . . . . . . . . 87
5.3 Résultats de la procédure de couplage . . . . . . . . . . . . . . . . . 88
5.3.1 Cas test d’un canal plan en présence d’un obstacle vertical . 89
5.3.2 Écoulement dans un canal plan en présence de deux obstacles 92
5.3.3 Conduite en forme de T, avec une paroi flexible . . . . . . . 96

CONCLUSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

BIBLIOGRAPHIE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
 ix

LISTE DES FIGURES

3.1 Domaine de calcul pour le problème de Poiseuille et Maillage (2000

éléments). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.2 Isovaleurs de vitesse u dans le canal. . . . . . . . . . . . . . . . . . 32
3.3 Isovaleurs de pression correspondantes. . . . . . . . . . . . . . . . . 32
3.4 Pression en fonction de x au milieu du canal, comparer avec la pres-
sion analytique et stabiliser par SUPG. . . . . . . . . . . . . . . . . 33
3.5 Domaine de calcul et conditions aux limites. . . . . . . . . . . . . . 34
3.6 Maillage du domaine, 11000 noeuds et 23000 éléments. . . . . . . . 34
3.7 Lignes de courant pour Re = 100. . . . . . . . . . . . . . . . . . . . 35
3.8 Isovaleurs de la vitesse u. . . . . . . . . . . . . . . . . . . . . . . . . 35
3.9 Isovaleurs de la vitesse v. . . . . . . . . . . . . . . . . . . . . . . . . 36
3.10 Isovaleurs de la pression. . . . . . . . . . . . . . . . . . . . . . . . . 36
3.11 Lignes de courant pour Re = 100. Résultats Filippone [48]. . . . . . 36
3.12 Allure de vitesse et de la pression en amont de l’obstacle (couleur
jaune et rouge) et en aval de l’obstacle (couleur verte). . . . . . . . 37
3.13 Maillage homogène, 2689 noeuds. . . . . . . . . . . . . . . . . . . . 37
3.14 Champ de vitesses, Re = 100. . . . . . . . . . . . . . . . . . . . . . 38
3.15 Champ de vitesses, Re = 400. . . . . . . . . . . . . . . . . . . . . . 38
3.16 Champ de vitesses, Re = 1000. . . . . . . . . . . . . . . . . . . . . . 38
3.17 Champ de vitesses, Re = 2500. . . . . . . . . . . . . . . . . . . . . . 39
3.18 Isovaleurs vitesses, step = 60, step = 65. . . . . . . . . . . . . . . . 39
3.19 Isovaleurs vitesses, step = 70, step = 75. . . . . . . . . . . . . . . . 39
3.20 Géométrie de la cavité et le maillage (30 × 30). . . . . . . . . . . . 40
3.21 A gauche : vitesse u. A droite : vitesse v pour Re = 1. . . . . . . . 41
3.22 A gauche : vitesse u. A droite : vitesse v pour Re = 100. . . . . . . 41
 x

3.23 A gauche : vitesse u(y). A droite : vitesse v(x) pour Re = 100. . . . 42

3.24 Vitesse en fonction de nombre de Reynolds. . . . . . . . . . . . . . 42
3.25 Maillage de domaine (4000 noeuds et 8000 éléments) et conditions
aux limites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.26 Pression obtenue par notre modèle, Re = 1. . . . . . . . . . . . . . 44
3.27 Pression obtenue par TDYN, Re = 1. . . . . . . . . . . . . . . . . . 44
3.28 Vitesse u obtenue par notre modèle, Re = 10. . . . . . . . . . . . . 45
3.29 Vitesse u obtenue par TDYN, Re = 10. . . . . . . . . . . . . . . . . 45
3.30 Les lignes de courant obtenues par notre modèle Re = 10. . . . . . . 45
3.31 Les lignes de courant obtenues par TDYN Re = 10. . . . . . . . . . 45
3.32 Domaine et conditions aux limites. . . . . . . . . . . . . . . . . . . 46
3.33 Les lignes de courant, Re = 100. . . . . . . . . . . . . . . . . . . . . 46
3.34 Les lignes de courant, Re = 250. . . . . . . . . . . . . . . . . . . . . 47
3.35 Les lignes de courant, Re = 500. . . . . . . . . . . . . . . . . . . . . 47
3.36 Coupe de la composante horizontale de la vitesse, nos résultats. . . 48
3.37 Coupe de la composante horizontale de la vitesse, résultats Kouakou
[96]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.1 Géométrie et variables nodales de l’élément DKQ2. . . . . . . . . . 52

4.2 Géométrie et variables nodales de l’élément DKT 15. . . . . . . . . 53
4.3 Elément de coque courbe Q4γ20. . . . . . . . . . . . . . . . . . . . 57
4.4 Elément de référence à gauche et l’élément réel à droite . . . . . . . 58
4.5 Panneau cylindrique soumis à son poids propre. . . . . . . . . . . . 64
4.6 Type de maillage utilisé. . . . . . . . . . . . . . . . . . . . . . . . . 65
4.7 Convergence du déplacement W au bord libre B. Eléments à 3 noeuds. 65
4.8 Convergence du déplacement W au centre C. Eléments à 3 noeuds. 66
4.9 Convergence du déplacement W au bord libre B. Eléments à 4 noeuds. 67
4.10 Convergence du déplacement W au centre C. Eléments à 4 noeuds. 67
 xi

4.11 Déplacement V le long de DA. . . . . . . . . . . . . . . . . . . . . . 68

4.12 Déplacement W le long de CB. . . . . . . . . . . . . . . . . . . . . 68
4.13 Domaine et maillage (707 noeuds, 600 éléments). . . . . . . . . . . 69
4.14 Domaine déformé, f est la flèche). . . . . . . . . . . . . . . . . . . . 69
4.15 Vecteurs déplacements vy en mm. . . . . . . . . . . . . . . . . . . . 70
4.16 L’évolution de la solution y = f (x). . . . . . . . . . . . . . . . . . . 71
4.17 Domaine et maillage (707 noeuds et 600 éléments). . . . . . . . . . 72
4.18 Vecteurs déplacements et domaine déformé (×200). . . . . . . . . . 74
4.19 L’évolution de la solution calculée par notre modèle et la solution
analytique. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.20 Géométrie et conditions aux limites d’après G. Dhatt [41]. . . . . . 75
4.21 Maillage : 181 noeus et 150 éléments. . . . . . . . . . . . . . . . . . 76
4.22 Vecteurs déplacements (mm) suivant l’axe y et maillage déformé
(facteur= 200). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.23 Les déplacements suivant l’axe y pour N1 , N2 et N3 . . . . . . . . . . 77

5.1 Evolution du maillage suivant la formulation utilisée. . . . . . . . . 81

5.2 Transfert des données entre les deux codes. . . . . . . . . . . . . . . 87
5.3 Algorithme d’interaction fluide-structure. . . . . . . . . . . . . . . . 88
5.4 Domaine et conditions aux limites. . . . . . . . . . . . . . . . . . . 89
5.5 Maillage : 793 noeuds et 1400 éléments. . . . . . . . . . . . . . . . . 89
5.6 Zoom sur l’interface fluide-structure. . . . . . . . . . . . . . . . . . 89
5.7 Champ de pression le long de la conduite. . . . . . . . . . . . . . . 90
5.8 Déformation de la structure, E = 2.104 M P a. . . . . . . . . . . . . . 90
5.9 Les lignes de courant, Re = 1. . . . . . . . . . . . . . . . . . . . . . 91
5.10 Module de Young= 104 M P a. . . . . . . . . . . . . . . . . . . . . . 91
5.11 Domaine de calcul et conditions aux limites. . . . . . . . . . . . . . 92
5.12 Pression le long du canal : rouge à travers le premier obstacle et
jaune à travers le deuxième obstacle. . . . . . . . . . . . . . . . . . 93
5.13 Champs des pressions, Modude de Young E = 4.104 M P a. . . . . . 94
5.14 Déformation des structures avec un modude de Young E = 4.104 M P a. 94
5.15 Champs des pressions, Modude de Young E = 104 M P a. . . . . . . 94
5.16 Déformation des structures avec un modude de Young E = 104 M P a. 95
5.17 Lignes de courant, Modude de Young E = 4.104 M P a et Re  1. . . 95
5.18 Domaine et conditions aux limites. . . . . . . . . . . . . . . . . . . 96
5.19 Champ de pression (Pa), ( E = 104 M P a). . . . . . . . . . . . . . . 96
5.20 Vitesse u(m/s) de l’écoulement (E = 104 M P a). . . . . . . . . . . . 97
5.21 Vitesse v(m/s) de l’écoulement ( E = 104 M P a). . . . . . . . . . . . 97
5.22 Lignes de courants ( E = 104 M P a). . . . . . . . . . . . . . . . . . . 97
5.23 Déformation de la paroi ( E = 104 M P a). . . . . . . . . . . . . . . . 98
5.24 Déformation de la paroi ( E = 75.103 M P a). . . . . . . . . . . . . . 98
5.25 Déformation du paroi ( E = 50.103 M P a). . . . . . . . . . . . . . . 98
5.26 Vecteurs déplacements (E = 104 M P a). . . . . . . . . . . . . . . . . 99
5.27 Vecteurs déplacements (E = 75.103 M P a). . . . . . . . . . . . . . . 99
5.28 Vecteurs déplacements (E = 50.103 M P a). . . . . . . . . . . . . . . 99
5.29 Positions des coupes longitudinales. . . . . . . . . . . . . . . . . . . 100
5.30 Evolution de la vitesse u(x), en fonction de la position de la coupe. 100
5.31 Evolution de la vitesse v(x), en fonction de la position de la coupe. 101
 xiii

LISTE DES NOTATIONS ET DES SYMBOLES

I tenseur identité
J matrice de préconditionnement de Jacobi
M matrise masse
η viscosité cinématique pour un fluide newtonien
Ω domaine physique
Γ frontière du domaine physique
ρ densité
(ξ, η, ζ) coordonnées cartésiennes pour le passage à l’élément de réference
∇ opérateur gradient
k.k norme vectorielle
t quantité dépendant du temps
D condition de Dirichlet en vitesse
N condition de Neumann en vitesse
d dimension d’un espace
e quantité évaluée sur un élément
max quantité maximale
t variable temps
x, y coordonnées cartésiennes
u, v vitesse du fluide
p pression du fluide
µ viscosité dynamique du fluide
ndle nombre de degrés de libérté par élément
ndlt nombre de degrés de libérté total
ν coefficient de Poisson
 INTRODUCTION

Le couplage fluide-structure apparaı̂t systématiquement dès qu’une structure

vibrante est immergée dans un fluide au repos ou vise versa. Il est dû à l’accéleration
des particules du fluide, au voisinage de la structure. La variation de pression qui en
résulte induit sur la structure, une force d’inertie qui agit comme une masse ajoutée
fixée au système mécanique, et se déplaçant solidairement avec lui. En présence
d’un fluide visqueux, un amortissement est induit par la viscosité et s’ajoute à
l’amortissement propre de la structure.
Avec le développement des codes de calcul en mécanique des fluides et l’augmen-
tation des ressources informatiques, la simulation numérique devient une alterna-
tive intéressante et complémentaire pour l’étude des couplages fluide-structure.
En effet on peut envisager de reproduire numériquement les couplages fluide-
structure étudiés expérimentalement et analytiquement. On est alors confronté
aux problématiques suivantes : assurer la conservation de l’énergie lors de la trans-
mission des données entre les codes et gérer l’interface mobile où s’échangent les
données et où se crée le couplage. Les couplages fluide-structure s’opèrent au niveau
de l’interface fluide-structure. En théorie on suppose que l’on a continuité des con-
traintes et des vitesses à l’interface entre le fluide et la structure. 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 forte-
ment convectifs. La formulation Lagrangian-Eulérian Arbitraire (ALE) est une so-
lution qui permet de gérer des maillages mobiles avec des déplacements d’amplitude
modérée. Le couplage entre codes fluide et structure se fait de la façon suivante
: on calcule dans un premier temps, les chargements fluides exercés sur la struc-
ture. 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 à
 2

l’itération suivante. Toutefois la gestion des conditions limites en particulier pour

la prise en compte de parois mobiles avec des écoulements turbulents reste un
problème complexe.
Dans le cadre de cette thèse, le fluide est considéré incompressible et le solide
est un corps solide élastique. Pour le solide, la cinématique de son mouvement peut
être décomposée en une translation (deux ou trois degrés de liberté suivant que le
problème est plan ou tri-dimensionnel) et/ou une rotation (un ou trois degrés de
liberté suivant la dimension du problème) [60, 84, 95].
Une multitude de méthodes et de modèles est utilisée pour résoudre l’interaction
fluide-structure. Ceci révèle l’extrême complexité de ce type de problème, malgré
l’ajout de certaines simplifications. L’interaction entre le fluide et le solide fait in-
tervenir la continuité de la vitesse et la condition d’équilibre mécanique à l’interface.
Ces conditions sont très largement utilisées dans les différents modèles d’interaction.
On a fait le choix d’une formulation ALE pour le fluide et d’une formulation la-
grangienne pour la structure, où la simulation numérique est basée sur la méthode
des éléments finis.
Cette interaction est basée sur une procédure de couplage instantané, où à
chaque itération la pression obtenue de modèle numérique fluide est utilisée comme
sollicitation externe pour le modèle numérique solide. A cet effet, un modèle de cou-
plage a été développé pour le pilotage instantané des deux modèles fluide-structure.
Un des avantages de cette procédure de couplage instantané est l’optimisation
du temps de calcul et de mémoire, contrairement au couplage fort, qui nécessite une
grande capacité de stockage pour les matrices masse et rigidité. L’algorithme pro-
posé permet aussi facilement de modifier séparément les codes fluide ou structure
afin d’envisager des lois de comportement différentes.
Enfin l’utilisation de la méthode ALE (Arbitrary Lagrangian-Eulerian method)
permet de résoudre les équations du problème fluide dans un domaine mobile
s’adaptant aux déformations de la structure.
La structure de cette de cette thèse est la suivante. Dans le premier chapitre,
nous présentons l’état de l’art sur l’interaction fluide-structure ainsi qu’une revue
bibliographique sur le sujet. Le deuxième chapitre est consacré à la formulation
mathématique de la méthode des éléments finis, son principe et quelques méthodes
de stabilité utilisés dans la programmation du modèle fluide-structure. Dans le
troisième chapitre, des tests de validation numérique de modèle fluide proposé sont
effectués et discutés. Ils concernent des ”Benchmark” habituellement utilisés en
CFD. Dans le chapitre 4 une description détaillée du modèle structure est présentée.
Des tests de validation de ce modèle, qui concernent plus particulièrement des co-
ques minces sont donnés. Au chapitre 5 la couplage fluide-structure est discuté.
Il contient des discutions sur les formulations numériques des interactions fluide-
structure. Après une introduction des différentes formulations eulérienne, lagrang-
ienne et Lagrangian-Eulérian Arbitraire (ALE), nous avons détaillés les techniques
de couplage utilisées dans les présentes investigations, l’algorithme de couplage et
enfin des exemples tests. Ces exemples tests concernent les écoulements de flu-
ide dans des conduite en présence ou pas d’obstacles. Les obstacles sont considérés
comme solides déformables. Les canalisations sont également considérées aussi bien
comme des structures rigides ou flexibles.
Une synthèse des différents résultats et discutions obtenus dans cette thèse sont
données dans la conclusion.
 4

CHAPITRE 1

ETAT DE L’ART

1.1 Introduction

L’interaction fluide-structure s’intéresse au comportement d’un système con-

stitué par deux entités mécaniques considérées comme distinctes : une structure
mobile (rigide ou déformable) et un fluide (en écoulement ou au repos) autour ou
à l’intérieur de la structure. L’évolution de chacune des deux entités dépendant de
celle de l’autre, un phénomène de couplage apparaı̂t. Plus précisément, le mouve-
ment de la structure est influencé par l’écoulement du fluide à travers les efforts
transmis à l’interface, et réciproquement, le mouvement de la structure influence
l’écoulement du fluide par les déplacements de l’interface qui entraı̂ne le fluide dans
son mouvement.
Dans la nature, nombre de situations relèvent de ce type de problème, [36,99,110].
Nous nous intéressons au cas où la structure subit en présence de l’écoulement flu-
ide des pressions qui modifie sa géométrie. Ce type de situation intervient par
exemple dans les cas suivants :
hautes cheminées, gratte-ciels, ponts suspendus, câbles de téléphérique vibrant
sous l’effet du vent, ailes d’avion, pâles d’hélice, aubages de turbine, tuyauteries
de réfrigération, réservoirs subissant l’effet de ballottement de la surface libre du
liquide, etc.
En pratique, on suppose que la force fluide exercée sur une structure peut se
décomposer en trois contributions :
- une force, dite aléatoire ou force indépendante du mouvement : elle désigne
la force exercée par le fluide sur la structure fixe. Elle est due à la perte de charge
et aux frottements visqueux au niveau de la paroi fixe. Elle peut être stationnaire
 5

ou présenter un caractère aléatoire lors de lâcher de tourbillons ou en présence

d’écoulements turbulents ou diphasiques.
- une force fluide-structure : elle apparaı̂t lorsqu’on a une structure vibrante
dans un fluide au repos sans écoulement propre.
- une force fluide-élastique : cette force est due au mouvement de la structure
dans un écoulement fluide.
Le couplage fluide-élastique est un couplage dynamique non-conservatif. Les
flux d’énergie échangés entre le fluide et la structure dépendent de la vitesse de
convection de l’écoulement moyen. Or, dans certaines configurations au delà d’une
certaine vitesse caractéristique de l’écoulement moyen, le fluide transfère à la struc-
ture plus d’énergie qu’elle ne peut en dissiper. Le système devient alors ”instable”
: la structure se met à osciller fortement et peut éventuellement impacter les struc-
tures avoisinantes. Ce type de comportement, qualifié de départ en ”instabilité
fluide-élastique”, est en général suivi d’un endomagement [16].
De telles vibrations peuvent endommager la structure concernée. En effet, il
se peut, qu’à partir d’un seuil de vitesse du fluide, des mécanismes de couplage
instables apparaissent, où une partie de l’énergie cinétique du fluide se transforme
en énergie vibratoire de la structure. On observe alors une instabilité dite de
flottement particulièrement dangereuse pour la structure, pouvant mener ainsi à la
ruine de cette dernière.
Parmi les systèmes propices à de telles instabilités, les lignes de tuyauteries
de refrigération dans les réacteurs nucléaires ont été les premières étudiées, [36].
Paı̈doussis, [ 110,97,13] traite du problème du flottement d’un tube encastré en
écoulement interne. Dans le domaine de l’aéronautique, on peut citer l’accident
du chasseur F117 lors d’une rencontre aérienne en 1997 aux États Unis. L’origine
de cet accident a été attribué au flottement d’un aileron, problème générique de
cet avion qui semblait pourtant avoir été résolu par le renforcement des struc-
 6

tures. La présence d’un écoulement peut aussi conduire à des instabilités dites de
flambage, d’intérêt particulier en biomécanique pour les problèmes de contraction
de vaisseaux sanguins ou d’asthme[80,130]. Dans ces phénomènes, l’instabilité est
d’origine géométrique et est induite par la rotation de l’interface et de la charge
appliquée.

1.2 Fluide-structure et problèmes associés

Il existe une grande variété de phénomènes d’interaction fluide-structure. On

peut distinguer formellement les problèmes instationnaires et stationnaires. Ces
derniers peuvent apparaı̂tre paradoxaux, comme la forme stationnaire d’une voile
dans un vent permanent, ou celle d’un ballon de baudruche après gonflage. On
peut surtout distinguer les problèmes par les équations régissant chacun des sous
systèmes fluide ou structure. La distinction est moins naturelle quand on regarde les
équations régissant la structure (modèles linéaires, en petits déplacements ou non,
en petites déformations ou non, etc...) Par contre, une classification en fonction de
modèle retenu pour le fluide prend plus de sens.
- Génie Civil : un des premiers grands secteurs où les simulations numériques
et expérimentales d’interactions fluide-structure sont apparues est le Génie Civil,
notamment pour les structures très élancées (ponts suspendus) et faisceaux de
câbles.
- Génie Maritime: dans ce domaine, de nombreuses études numériques concer-
nent différents aspects des interactions fluide-structure. La plupart des études font
jouer un rôle particulier à la surface libre de la mer. Une des sollicitations les plus
importantes pour des constructions portuaires ou des structures pour l’extraction
off-shore du pétrole est simplement la houle (ondes sans déplacement de matière à
la surface libre).
- Hydrodynamique navale : est un autre champ d’application essentiel. Il s’agit
 7

de concevoir et d’optimiser des formes de coques et d’autres éléments sur différents

véhicules.
- Génie Biomédical : concerne l’interaction entre des tissus biologiques. Dans
de nombreux cas, les nombreux fluides du corps humains (sang, liquides oculaires,
etc...) peuvent interagir avec les environnants. Ceux-ci étant déformables, de
véritables interactions fluide-structure peuvent apparaı̂tre naturellement (veines
collabables, rupture d’anévrisme, etc...). Pour une meilleure compréhension de
l’apparition et de la rupture d’anévrisme, certaines équipes ont déjà proposé des
chaı̂nes plus ou moins automatiques allant de l’imagerie médicale sur le patient
lui-même, à l’extraction des parois des vaisseaux [49].
- Génie Industriel : de nombreuses études sont menées dans ce domaine. Elles
concernent : les structures très fines et très élancées comme les câbles ou les pan-
neaux, ou les réservoirs. Des calculs ont par exemple été menés sur des lignes haute
tension (interaction avec le vent, chargement de neige ou de glace, etc...) ou plus
simplement sur des câbles.

1.3 Méthodes de couplage Fluide-structure

1.3.1 Méthodes numériques

La structure est, en général, résolue par une formulation classique en éléments

finis. Le problème du couplage fluide-structure réside dans la recherche d’une so-
lution numérique pour le problème fluide. Une fois, cette solution trouvée, il faut
la coupler à la formulation en éléments finis pour la structure afin de résoudre
le problème d’interaction. Les méthodes numériques doivent être conçues pour
suivre les déformations du fluide. Dans un premier temps, la formulation Lagrang-
ienne ou ALE en éléments finis. L’interaction fluide-structure est alors gérée par
des méthodes de contact. Dans ce cas, cette approche est limitée par la distor-
 8

sion des mailles. Des méthodes de remaillage peuvent être envisagées mais leur
coût CPU est très important. Ensuite, les solutions apportées par la méthode des
éléments frontières (BEM) ou par la méthode SPH seront décrites: elles permet-
tent de s’affranchir du maillage et donnent des résultats très proches des solutions
physiques.

[Link] La méthode du ’Lagrangien augmenté’ (pénalisation)

Le problème majeur des équations de Navier-Stokes est le traitement du cou-

plage entre la pression et la vitesse. La difficulté réside dans le fait de calculer
un champ de pression résultant d’un écoulement incompressible avec V qui vérifie
l’équation de continuité, le tout sans disposer d’équation traduisant de manière
explicite les évolutions de la pression. En incompressible, plusieurs méthodes de
résolution existent, dont la relativement récente méthode dite du ”Lagrangien aug-
menté”. Cet algorithme a été élaboré dans un premier temps dans le but de résoudre
le couplage vitesse-pression induit par les équations de Navier-Stokes écrites pour
un écoulement incompressible. Il repose sur une méthode de minimisation sous la
contrainte de l’équation de continuité, où la pression, découplée par rapport à la
vitesse, apparaı̂t comme un multiplicateur de Lagrange [51]. La contrainte est en
fait directement introduite dans l’équation du mouvement sous la forme du terme
de pénalisation qui couple les différentes composantes de la vitesse. Le calcul du
couple vitesse-pression est effectué par l’algorithme itératif d’Uzawa [51, 58]. Les
itérations de lagrangien augmenté (notées k) sont répétées jusqu’à ce que la valeur
moyenne de la divergence de la vitesse, dans l’ensemble du domaine, soit suff-
isamment petite. En notant K le nombre maximal d’itérations nécessaires pour
satisfaire le précédent critère.
 9

[Link] Méthode des éléments frontières (BEM method)

Longuet-Higgins et Cokelet [78] furent les premiers à introduire cette technique

connue aussi sous le nom de méthode Euler-Lagrange mixte. Des variantes de
cette méthode ont été appliquées à une large gamme de problèmes 2D et 3D par
de nombreux auteurs tels que Faltinsen [47], Vinje et Brevig [127], Cointe [38],...
La présentation de cette approche est détaillée dans les travaux de Zhao et Faltin-
sen. Ils présentèrent, dans [138], une méthode numérique permettant de calculer la
répartition de pression sur le dièdre à partir d’une méthode non-linéaire d’éléments
frontières. Cette méthode numérique s’applique à des problèmes linéaires ellip-
tiques dont les équations partielles peuvent être reformulées sous formes intégrales.

[Link] Méthode SPH

En 1977, les travaux de Lucy et Monaghan [79] donnent lieu à une nouvelle
méthode numérique: la méthode SPH (Smoothed Particle Hydrodynamics). Cette
méthode particulaire s’affranchit d’une connectivité géométrique, contrairement
aux formulations classiques. Les domaines physiques sont discrétisés par des partic-
ules qui suivent le mouvement du matériau. Pour le tossage, la formation du jet ne
pose ainsi plus de problèmes de modélisation. Cette méthode fut d’abord appliquée
en astrophysique, [79] et [57], avant d’être utilisée par la suite en Mécanique, [85].
Les particules représentent les points d’interpolation où les inconnues mécaniques
du problèmes sont déterminées. La méthode SPH se caractérise par la masse des
particules et par deux grandeurs de discrétisation, la distance interparticulaire et
la longueur de lissage. Le principe est d’approximer un champ Φ sur un domaine
Ω par un ensemble de points de discrétisation appelés particules.
 10

[Link] Méthode VOF

Cette méthode a été introduite la première fois par Hirt et Nichols [89] dans
un code appelé SOLA-VOF pour traquer les interfaces matérielles dans une grille
Eulérienne. La méthode est liée à la résolution des équations d’Euler ou de Navier-
Stokes par un schéma de différences finies ou de volumes finis. Le principe du
suivi d’interface est de déterminer la quantité de fluide dans chaque volume de
contrôle (cellule) en évaluant une quantité appelée fraction volumique. Si cette
dernière vaut 1 alors la cellule est remplie de fluide ; si elle vaut 0, la cellule
est soit vide, soit remplie d’un autre fluide. Les fractions volumiques comprises
entre 0 et 1 déterminent la position de l’interface matérielle ou de la surface libre.
En connaissant la distribution de la fraction volumique d’un fluide, il est possible
de déterminer les limites physiques du domaine fluide dans un maillage. Pour
cela, Young [136] propose une technique pour calculer précisément la position de
l’interface dans une cellule en utilisant les fractions volumiques des cellules voisines.
La méthode VOF a les trois caractéristiques nécessaires à une technique de suivi des
interfaces matérielles : un schéma localisant l’interface, un algorithme déterminant
précisément la position de l’interface dans le maillage et un algorithme appliquant
les conditions aux limites sur cette frontière.

1.3.2 Méthodes analytiques

Pour présenter les deux approches analytiques dominantes permettant la résolution

du problème d’impact, on peut effectuer un rapide parallèle entre elles en définissant
leurs limites [137]. Chacune de ces méthodes a ses propres caractéristiques :
1) La méthode des similitudes respecte les conditions de surface libre et elle
n’est applicable que pour des dièdres entrant à vitesse constante dans le fluide.
2) La méthode asymptotique est limitée à des angles d’incidence petits, inférieurs
à 30◦ d’après [138]. Cette dernière méthode est le fruit d’une succession d’articles
 11

dont le premier est [128]. Von karman fut, en effet, le pionnier dans l’étude du
tossage. Comme les hypothèses sur lesquelles elle repose s’avèrent restrictives,
Wagner développa une théorie plus évoluée [129] mais l’expression de la pression
est singulière sur le bord du dièdre. Watanabe améliora cette théorie pour con-
trecarrer ce problème [131] et développa l’expression de la pression en prenant en
compte le jet : il raccorda la solution de Wagner avec la solution du jet sur une
plaque. Cointe [37] donna, indépendamment, une théorie similaire. Zhao et Faltin-
sen, dans l’article [137], reprirent aussi l’expression de la pression de Wagner pour
la raccorder asymptotiquement avec une solution externe pour le jet. La méthode
des similitudes repose sur la théorie de Dobrovolskaya développée dans [42]. Cette
théorie qui fut aussi l’objet des articles [100] et [54], considère que l’écoulement
est autosimilaire réduisant ainsi l’étude à la résolution d’une équation intégrale
non-linéaire unidimensionnelle.

1.4 Formulation ALE des équations de Navier-Stokes

En interaction fluide-structure, les équation de Navier-Stokes sont résolues dans

un domaine mobile ou déformable. Parmi les approches utilisées, la méthode
Lagrangian-Eulérian Arbitraire (ALE) suscite beaucoup d’intérêt [117]. L’intérêt
d’utiliser une formulation ALE [44] pour décrire les phénomènes est de combiner
les avantages de la description Lagrangienne et ceux de la description Eulérienne.
Soit X les coordonnées Lagrangiennes d’un point dans le repère matériel, x ses
coordonnées Eulériennes dans le repère spatial et χ ses coordonnées dans un repère
ALE. La relation de passage [19]. L’ALE traite le maillage comme un repère qui
se déplace avec une vitesse arbitraire w :

f (x, t) = fb(χ, t) (1.1)

 ∂ fb ∂ fb
La dérivée mixte est donné par : |χ = ∂f | + w.∇x f , cela nous permet
∂t x
∂t ∂t
finalement d’écrire la relation de passage, pour la quantité f entre le repère matériel
et le repère ALE, soit :

df ∂f ∂ fb
= |x + u∇x f = |χ + (u − w)∇x f (1.2)
dt ∂t ∂t
∂x ∂x
où u(x, t) = |
∂t x
est la vitesse matérielle et w = |
∂t χ
est la vitesse arbitraire.
Cette propriété est très utile dans le développement d’une description cinématique
dans le repère ALE.
La loi de conservation de la quantité de mouvement exprimée en ALE :

∂u 1
|ξ + (u.5)u = − 5 p + η∆u + f (1.3)
∂t ρ

La relation 1.2 appliquée à u donne :

· ∂u ∂u (ξ, t)
u= |x + (u.5)u = |ξ + ((u − w) .5) u (1.4)
∂t ∂t

Enfin, la relation fondamentale qui permet de passer d’une loi exprimée en

variable Eulérienne à une loi équivalente exprimée en variables mixtes (ALE) :

∂u (ξ, t) 1
|ξ + ((u − w) .5) u = − 5 p + η∆u + f (1.5)
∂t ρ

On obtient donc l’écriture des équations de Navier-Stokes dans la description

ALE en incompressible :

 
∂u
ρ + ((u − w) .5) u = − 5 p + µ∆u + f (1.6)
∂t

Cette formulation assigne à chaque point de l’espace une vitesse de grille w, qui
décrit un mouvement arbitraire du maillage et une vitesse matérielle u.
 13

CHAPITRE 2

RÉSOLUTION DES ÉQUATIONS D’ÉCOULEMENT DE FLUIDES

PAR ELÉMENTS FINIS

2.1 Introduction : Méthode des éléments finis

La méthode des élément finis est fondée sur deux idées principales : une méthode
d’interpolation et une méthode variationnelle.
Etant donnés un domaine continu Ω et une fonction u solution de L(u) = 0 ,

• a) la méthode d’interpolation fournira une fonction proche de u construite à

partir d’un nombre finis de réels

• b) la méthode variationnelle remplacera L(u) = 0 par un nombre fini d’équations.

Si l’on suppose la fonction approchée connue en un certain nombre de points de

l’élément, elle est évaluée en n’importe quel point de l’élément par interpolation.
La restriction de la fonction approchée à un élément K (triangle) du maillage,
notée φK
h , est de la forme :

XNk
φK
h (x, y) = φi Ni (x, y) (2.1)
i=1

où φi désigne la valeur de la fonction approchée φh au point xi de coordonnées

(xi , yi ), Ni la fonction d’interpolation asociée à xi , Nk le nombre de fonctions
utilisées pour l’élément K.
Les coefficients φi de cette combinaison linéaire sont les inconnues - ou degré
de liberté - du problème discret. Les fonctions d’interpolation sont les fonctions de
base d’un espace vectoriel de dimension finie Nh dans lequel on cherche à priori la
solution discrète φh .
 14

Afin de déterminer les inconnus du problème discret, une méthode de résidus

pondérés est utilisée : on minimise l’erreur commise est minimisée en remplaçant
dans les équations aux dérivées partielles la fonction exacte par la fonction ap-
prochée. Pour cela, on projette l’erreur sur un espace de fonction Φh à préciser.
La méthode de Galerkin est un cas particulier de la méthode des résidus pondérés,
l’espace de projection est le même que l’espace d’approximation de la solution du
problème, autrement dit Φh = Nh .

2.2 Eléments d’analyse fonctionelle

• Espace L1 (Ω)

On note L1 (Ω) l’espace vectoriel des classes de fonctions (on identifie les fonc-
tions qui sont égales presque partout) intégrables pour la mesure de Lebesgue dx.
C’ est un espace de Banach (normé complet) lorsque l’on prend pour norme de
f ∈ L1 (Ω) :

Z
kf kL1 = |f (x)| dx (2.2)
Ω

• Espaces Lp (Ω)

Soit 1 ≤ p ≤ ∞. On définit

Lp (Ω) = f ∈ F (Ω; R) : f est mesurable et |f |p ∈ L1 (Ω)


(2.3)

c’ est un espace de Banach pour la norme

Z 1/p
p
kf kLp = |f (x)| dx (2.4)
Ω
 15

c’ est un espace de Banach pour la norme

kf kL∞ = inf { C : |f (x)| ≤ C p.p. sur Ω} (2.5)

• Espaces de Sobolev

On appelle espace de Sobolev sur Ω , l’espace fonctionnel

W k;,p (Ω) = {u ∈ Lp (Ω)∀α ∈ Z, |α| ≤ k : Dα u ∈ Lp (Ω)} (2.6)

où Dα u est la dérivée de u au sens des distributions.

On note H k (Ω) = W k;2 (Ω), en particulier

H 1 (Ω) = . u ∈ L2 (Ω).∀i ∈ {1; ....; d} : uxi ∈ L2 (Ω) .


(2.7)

Soit l’application H 1 (Ω) × H 1 (Ω).−→ R+ , et soit

Z
(u, v) −→ (uv + ∇u.∇v) dx (2.8)
Ω

est un produit scalaire et munit H 1 (Ω) d’une structure d’espace de Hilbert. On

note la norme associée par :

Z
kukH 1 (Ω) = (u2 + |∇u|2 )dx (2.9)
Ω

• Inégalité de Poincaré

Soit Ω un ouvert borné de Rd , alors il existe une constante C(Ω) telle que

Z Z
∀u ∈ H01 (Ω) : 2
u dx ≤ C(Ω) |∇u|2 dx (2.10)
Ω Ω
 16

Soit Ω un ouvert borné de Rd , alors

Z 1/2 Z 1/2
2 2 2
kukH 1 (Ω) = (u + |∇u| )dx et kukH 1 (Ω) = |∇u| dx (2.11)
0
Ω Ω

sont des normes équivalentes sur H01 (Ω) .

2.3 Mise en oeuve de la méthode des éléments finis

2.3.1 Formulation variationnelle forte

La formulation forte d’un problème aux conditions aux limites s’écrit alors
comme suit :
Trouver u(x) ∈ espace des solutions US tel que

L(u) − F = 0 ∀x ∈ Ω B(u) = u0 ∀x ∈Γu0 (2.12)

On note que l’ espace U s dans lequel nous recherchons la fonction u , doit être
défini judicieusement.

2.3.2 Formulation variationelle faible

La formulation variationnelle faible peut être obtenue par les étapes suivantes :

1. approcher la solution u, partout dans le domaine Ω, par u

e telque :

n
X
u
e = Ni u
ei (2.13)
i=1

où Ni =Ni (x, y) désignent les fonctions de forme et u

e i des valeurs nodales
indéterminées.

2. introduire (2.13) dans (2.12) il en résulte une erreur ou un résidu  tel que :
 17

 = L(e
u) − F =
6 0 (2.14)

où  = 0, u
e = u c.à.d lorsqu’une solution exacte existe.

3. obtenir une valeur nulle de  en introduisant une fonction de point arbitraire

W tel que :
Z
u) − F) dΩ = 0, k = 1, 2...N
Wk (L(e (2.15)
Ω

où N Nombre d’inconnues indépendants du système (2.12).

En formulation de Galerkine les fonctions tests sont assimilées aux fonctions de
forme. On pose

Z n
!
X
Wk = Nk =⇒ Nk .L Ni u
e i dΩ = 0 (2.16)
Ω i=1

On dit que u ∈ H01 (Ω) est une solution faible de (2.18) si :

Z
∀v ∈ H01 (Ω) : a(e
u, v) = f vdΩ (2.17)
Ω

a(: , :) est la forme bilinéaire associée à L Soit Ω un ouvert borné de Rd et ∂Ω

est régulier. On considère alors le problème :

Lu = f sur Ω et u = 0 sur ∂Ω (2.18)

2.3.3 Conditions aux limites

Trouver la solution qui satisfait un système d’équations aux dérivées partielles

et les conditions aux limites (conditions intiales et conditions frontières) forme un
problème bien posé si les trois conditions suivantes sont satisfaites.
 18

La définition correcte des conditions aux limites est un des éléments essentiels
pour obtenir un problème bien posé, on rappel les différents types des conditions
aux limites communément rencontrées sont :

condition de Dirichlet : u=f

condition de Neumann u,n = f
condition de Cauchy u,n = f et u = g

2.4 Les équation du mouvement des fluides incompressible

2.4.1 Équation de continuité

Dans un milieu continu en mouvement l’équation de continuité doit être vérifiée

pour tout volume Ω, :
∂ρ ∂ ∂
+ (ρu) + (ρv) = 0 (2.19)
∂t ∂x ∂y

Dans le cas d’un fluide incompressible, le champs de vitesse vérifie nécessairement

:
∂u ∂v
∇.u = + =0 (2.20)
∂x ∂y

2.4.2 Conservation de la quantité de mouvement

On considère un écoulement d’un fluide viscqueux incompressible, irrotationnel,

soumis à des forces volumiques F et de surface σ. L’équation de quantité de
mouvement s’écrit :

∂u ∂u2
 2
∂ u ∂2u

∂(uv) 1 ∂p Fx
+ +v + −η 2
+ 2 − =0 (2.21)
∂t ∂x ∂y ρ ∂x ∂x ∂y ρ
∂(uv) ∂v 2 1 ∂p
 2
∂ v ∂2v

∂v Fx
+u + + −η 2
+ 2 − =0 (2.22)
∂t ∂x ∂y ρ ∂y ∂x ∂y ρ
∂u ∂v
∇.u = + =0
∂x ∂y
 19

2.5 Ecriture variationnelle faible dans le cas de fluides incompressibles

En choisissant les fonctions test sous la forme :

 
Ψ
 u 
  
 

 Ψu   
Ψ= , Ψu = Ψ (2.23)
 Ψ   v 
p 
 

Ψw
 

l’écriture faible s’écrit :

W = WN −S + WCont (2.24)

où l’indice N − S désigne les termes de conservation de la masse pour les équations
de Navier-Stokes et l’indice Cont désigne termes de continuité.
De l’équation 2.21 on obtient :

∂u ∂u2
Z    
∂(uv) Fx 1 ∂Ψu ∂Ψu ∂u ∂Ψu ∂u
WN −S = Ψu + +v − − p+ µ + µ dΩ+
∂t ∂x ∂y ρ ρ ∂x ∂x ∂x ∂y ∂y
Ω
∂v ∂v 2
Z    
∂(uv) Fy 1 ∂Ψv ∂Ψv ∂v ∂Ψv ∂v
Ψv + +v − − p+ µ + µ dΩ+
∂t ∂y ∂y ρ ρ ∂y ∂x ∂x ∂y ∂y
Ω
I I I   I  
∂u ∂u ∂v ∂v
Ψv pny ds + Ψu pnx ds − Ψu nx + ny ds − Ψv nx + ny ds
∂x ∂y ∂x ∂y
∂Ω ∂Ω ∂Ω ∂Ω
Z  
∂u ∂v
WCont = − Ψp + dΩ avec u = us sur Su (2.25)
∂x ∂y
Ω

n = (nx , ny ) étant la normale à la frontière dirigée vers l’extérieur de la frontière

S = ∂Ω.
Conditions aux limites associées sont : à la paroie Paroi Condition d’adhérence
∂u
u = 0, à l’Entrée u = u et à la Sortie p = p̄, .
∂η
 20

2.5.1 Discrétisation par éléments finis

1. Le domaine de calcul Ω est représenté par un ensemble de sous domaines

élémentaires Ωe , donc

N elt
X N elt
X
Ω= Ωie =⇒ W = Wei (2.26)
i=1 i=1e

N elt =Nombre total d’élements ,Wei : formulation faible élémentaire

2. Approximation nodale des variables :

Pour chaque élément on a :


u = hN i {un } , v = hN i {vn } , P = N {Pn } (2.27)

La forme conservative des équations de Navier-Stokes est discrétisée en utilisant

la même approximation pour u, v et u2 , uv, v 2 , c’est à dire:

u = hNi i {ui } v 2 = hNi i {vi2 } vδv = hNi vi i {δvi }

v = hNi i {vi } uv = hNi i {ui vi } uδv = hNi ui i {δvi } (2.28)

u2 = hNi i {u2i } uδu = hNi ui i {δui } vδu = hNi vi i {δui }

3. Approximation nodale des fonctions test (Formulation Galerkine)

Dans la formulation Galerkine, la fonction de forme est identique pour les fonc-
tions test et les variables.


Ψu = hN i {δΨun } , Ψv = hN i {δΨvn } , Ψp = N {δΨPn } (2.29)


• hN i sont des fonctions d’interpolation pour les vitesses et N sont des fonc-
tions d’interpolation pour la pression.
 21

• {un } , {vn } et {Pn } sont les variables nodales associées à chaque élément.

4. Formulation matricielle élémentaire

L’introduction des fonctions solution (u, v, P ) et des fonction tests (Ψu , Ψv , Ψp )

dans la formulation faible élementaire, on obtient :

 
•
e e e e
W = hδΨi [M] {U } + [K] {U } − {F } =0 (2.30)
•
=⇒ [M]e {U } + [K]e {U } − {F }e = 0 (2.31)

• [M]e [ndle, ndle] : Matrice masse élémentaire

• [K]e [ndle, ndle] : Matrice de rigidité élémentaire

• {F }e [ndle] : Vecteur force élémentaire

La forme matricielle générale du problème est donnée par :

N elt  
X •
W = Wei =⇒ W = hδΨi [M]{U } + [K]{U } − {F } = 0 (2.32)
i=1e

Le système à résoudre est alors :

•
[M]{U } + [K]{U } = {F } (2.33)

{U (t0 )} = {U0 }

{U } : Vecteur des variables globales

[M] [ndlt, ndlt] : Matrice masse globale
[K] [ndlt, ndlt] : Matrice de rigidité globale
{F } [ndlt] : Vecteur global des sollicitations extérieures.
 22

2.6 Méthodes de stabilité numérique

Pour mettre en évidence les différentes méthodes de stabilité on considère le

cas simple d’un écoulement incompressible stationnaire de fluides newtoniens. Les
composantes des contraintes normales et tangantielles s’expriment en fonction de
pression et des taux de déformation :

∂u
σx = −p + 2µ
∂x
∂v
σy = −p + 2µ (2.34)
∂x
τxy = τyx = 0

Les équations linéarisées de mouvement du fluide vérifient :


 5·σ(u, p) + f = 0
, ∀x ∈ Ω (2.35)
 5·u=0

Avec les conditions aux limites suivantes :

σ(u, p).n = g, ∀x ∈ Γn ; u = u0 , ∀x ∈ ΓD (2.36)

[Link] Formulation faible mixte

b ·(5·σ(u, p) + f >= 0, ∀b
<u u ∈ Û (2.37)

En effectuant une intégration par parties et en utilisant le théorème de la di-

vergence, on obtient :
 
u, u) − b(p, u
 < a(b u) ∀b
b , ) >= f (b u ∈ Û
,
 p5·u = 0 ∀b
p ∈ Pb

a(b
u, u) =< (b
u), b(p, u
b , ) = p5.b u) = hb
u, f (b u.f i + hhb
[Link] (2.38)

[Link] Formulation discrète

La Formulation discrète consiste à écrire une approximation de u et p sous la

forme suivante

n
X n
X
u (x)= ui φui (x),
h
p (x)= pi φpi (x)
h
(2.39)
i=1 i=1

où ui et pi sont des valeurs nodales inconnues, tandis que φui (x) et φpi (x) sont
des fonctions de forme spécifiées. On a volontairement utilisé des fonctions de
forme de type différent pour la pression et les déplacements.
La formulation discrète s’écrit :

T rouver uh (x) et ph (x) ∈ Û et c

P tels que :

 < a(b uh , uh ) − b(ph , u
b h ) >= f (b
uh ) ∀b
uh ∈ Û h



 ph , uh ) ∀b
b(b ph ∈ Pbh

(2.40)
 24

CHAPITRE 3

MODÈLE FLUIDE : TRAITEMENT NUMÉRIQUE ET

VALIDATION

3.1 Résolution numérique

Le code que nous avons développé pour les écoulements de fluides incompress-
ibles est conçu suivant l’architecture du code Reflex [10], développé à l’Université
de Technologie de Compiègne. Le maillage, le post-traitement et les graphes ont
été effectués et exécutés à l’aide du logiciel GID.

3.1.1 Algorithmes de résolution : Discrétisation Temporelle

La Discrétisation Temporelle consiste en l’utilisation du schéma d’Euler (1er

Ordre). Ce schéma consiste à discrétiser (2.33) par la méthode de différences finis
décentrées en introduisant le paramètre d’implicitation 0 ≤ α ≤ 1.

• 1
{U t+α∆t } = ({Ut+∆t } − {Ut })
∆t
Ut+α∆t = α {Ut+∆t } + (1 − α) {Ut } (3.1)

• Lorsque α = 0, ce schéma est explicite et son domaine de stabilité est re-

streint. Il est donc nécessaire de choisir un petit ∆t

• Lorsque α = 0.5, ce schéma est est appelé Crank-Nocholson. Ce schéma de

2ème ordre assure une meilleur stabilité et une une excellente précision.

• Schéma implicite
 25

soit s = α∆t,ona :

•
[M]{U }s+1 + [K]{U }s+1 = {F }s+1
• {U }s+1 − {U }s • •
{U }s+1 = = (1 − α){U }s + α{U }s+1 (3.2)
∆ts+1
• {U }s+1 − {U }s (1 − α) •
=⇒ {U }s+1 = − {U }s
α∆ts+1 α

•
On introduisant l’expression de {U }s+1 dans l’équation :

•
[M]{U }s+1 + [K]{U }s+1 = {F }s+1 (3.3)

On obtient :

•
[M]{U }s+1 − [M]{U }s + α∆ts+1 [K]{U }s+1 = [M]∆ts+1 (1 − α){U }s + α∆ts+1 {F }s+1

([M]s+A + α∆ts+1 [K]s+A ){U }s+1 = ([M]s − (1 − α)[K]s ∆ts+1 ){U }s

| {z } | {z }
e s+1
[K] e s
[K]

+ ((1 − α){F }s + α{F }s+1 ) ∆ts+1

| {z }
e s+1
[F]

e s+1 {U }s+1 = [K]

[K] e s {U }s + [F]
e s+1 (3.4)

Le schéma devient implicite car le calcul de la matrice K nécéssite le cal-

cul préalable de la matrice Masse au pas de temps précédent. Ce qui a pour
conséquence d’augmenter le temps de calcul.
La discrétisation du terme temporel, pour un pas de temps ∆t, est faite par le
biais du schéma d’Euler explicite telle que :

 
· {Ut+∆t } − {Ut }
U = (3.5)
t ∆t
 26

Ce schéma est inconditionnellement stable et plus particulièrement autorise des

pas de temps plus grands. Nous obtenons finalement :

([M ] + ∆t [K]) {Ut+∆t } = ∆t {Ft+∆t } + [M ] {Ut }

Pour la partie fluide, nous avons utilisé essentiellement des élément triangu-
laires. La discrétisation des formulations faibles issus des équation de Navient-
Stokes conduit à un système global non linéaire, mixte en (u,p). Ce système est
résolu à l’aide d’une méthode itérative de Newton Raphson basée sur une approche
asymptotique [62].

3.1.2 Linéarisation : Méthode itérative de Newton-Raphson

Plusieurs méthodes peuvent être utiliséés pour traiter les termes non-linéaires
de l’équation (3.4). Dans les présentes applications nous avons utilisés la méthode
itérative de Newton-Raphson [62]. Elle consiste à supposer dans un premier temps
que le résiduel R est nul, soit :

R(U ) = [K({U })] {U } − {F } = {0} (3.6)

Ensuite on développe cette dernière équation en série de Taylor en ne retenant

que les termes d’ordre 1. Il vient :

∂ {R}r
R(U ) = {R}r + ({U }r+1 − {U }r ) + O({U }r+1 − {U }r )2
∂ {U }
= {R}r + K T δ {U } + O(δ {U })2 = {0}
 
(3.7)
 T  ∂ {R}r
K = à {U } = {U }r : Matrice Tangente
∂ {U }
δ {U } = {U }r+1 − {U }r : Solution incrémentale
 27

Enfin, de l’équation (3.7) on déduit que :

 −1 −1
δ {U } = K T {R}r = K T {U }r ({F } − [K({U }r )] {U }r )

(3.8)

La solution à l’itération (r+1)th est obtenue en utilisant :

{U }r+1 = {U }r + δ {U } (3.9)

3.1.3 Stabilité de la solution : Condition de Ladysenskaya-Brezzi-Babuska

(LBB)

Considérons la formulation discrète mixte (2.40) et analysons la stabilité de ce

problème. Une condition naturelle de stabilité numérique devrait s’énoncer comme
suit.
Il existe une constante C telle que si une solution (uh ; ph ) ∈U h ×P h satisfait
(2.40), alors :
u + ph ≤ 1 kf k
h
1 0 −1 (3.10)
C

où les normes k.k1 et k.k0 sont les normes des espace de Sobolev H1 et H0
= L2 . La norme d’une forme b et la norme k.k1 du vecteur u = (u; v; w) sont
respectivement définies par les expressions

|f (u)|
q
u = kuk2 + kvk2 + kwk2
h
kf k−1 = sup h
, 1 1 1 1 (3.11)
u∈U h ku k1

Intuitivement, on souhaite qu’une petite erreur sur les données (ici, f) n’ait
qu’une conséquence finie sur la norme de la solution du problème mixte.
Il est possible d’obtenir une inégalité (3.10), si la condition suivante dite de
Ladysenskaya-Brezzi-Babuska est satisfaite :
Il existe une constante CLBB > 0 telle que :
 28

(3.12)
 h h 
b(p , u )
inf sup ≥ CLBB
ph ∈P h uh ∈U h kuh k1 kph k0

Pour simplifier la démonstration mathématique, on supposera que : Uh = Û h

b h . Cela conduit à prendre alors u
et Ph = P b h = uh et pbh = ph dans la formulation
discrète mixte (2.40) , on peut écrire :


 a(uh , uh ) + b(ph , uh ) = f (uh ) ∀uh ∈ U h
(3.13)
 b(ph , uh ) ∀ph ∈ P h

L’obtention de l’inégalité de stabilité à partir de la condition LBB se réalise en

deux parties.
- En soustrayant les deux lignes de (3.13), on obtient :

a(uh , uh ) = f (uh )
2
a uh 1 ≤ a(uh , uh ) = f (uh )
u ≤ 1 kf k
h
≤ kf k−1 uh 1 ⇒ 1 −1 (3.14)
a

- En ne considérant que la première ligne de l’équation (3.13) , on écrit :

−b(ph , uh ) = −a(uh , uh ) + f (uh ) (3.15)

En vertu de la condition LBB (3.12), on obtient après un calcul intermédiaire :

 29

 h h     
b(p , u ) a(uh , uh ) f (uh )
CLBB ≤ sup ≤ sup + sup (3.16)
uh ∈U h kuh k1
uh ∈U h kuh k1 h
uh ∈U h ku k1
| {z } | {z }
c =kf k−1
≤ kf k−1
a
h 1 c
⇒ p 0 ≤ ( + 1) kf k−1 (3.17)
CLBB a

3.1.4 Méthode SUPG (Streamline Upwind Petrov-Galerkin)

Cette méthode est utilisée dans les présentes applications. Elle consiste à
décentrer les fonctions de poids suivant les lignes de courant de l’écoulement. Pour
le noeud i d’un élément, la fonction de poids s’exprime :

 
αLe ∂Ni ∂Ni
νi = Ni + Vx + Vy (3.18)
2 |V| ∂x ∂y

où Le est la longueur de l’élément suivant la direction de la vitesse et mesurée

en passant par le point d’intégration de Gauss [74].
Contrairement à la méthode du produit, un seul facteur a est calculé pour tout
l’élément. Il fait intervenir Le et la norme de la vitesse :

1
α = coth Pe − (3.19)
Pe

Cette méthode a d’abord été introduite par l’intermédiaire d’une diffusion artifi-
cielle anisotrope orientée dans le sens de l’écoulement [74]. Elle présente l’avantage
de s’appliquer indifféremment à des éléments carrés et triangles. La méthode SUPG
réduit l’influence des noeuds situés perpendiculairement à l’écoulement par rapport
à la méthode produit. Hughes et Brooks [68] montrent que la méthode produit tend
à lisser exagérément les résultats dans la direction perpendiculaire à l’écoulement.
Ils parlent de ”crosswind diffusion”.
 30

Le décentrement des termes convectifs est obtenu en éléments finis par la

méthode SUPG (Streamline Upwind Petrov Galerkin). Pour cela on considère:

φi = Ni + γu. 5 Ni (3.20)

où Ni est la fonction de base choisi, γ un paramètre de stabilisation et u le

champ de vitesse associé au terme convectif.
Pour un élément en amont du noeud i, u. 5 Ni est positif. La fonction de
pondération associée sur cet élément au noeud i devient donc plus importante.

3.2 Validation du module fluide

La phase de validation du module fluide est effectuée suivant trois tests académiques,
et concernent : Poiseuille, écoulement d’un fluide autour d’un cylindre fixe, écoulement
dans une cavité carrée, écoulement dans une conduite à section variable (2 tubes de
diamètres différents), écoulement à travers une marche bas, étranglement du canal
et enfin l’écoulement d’un fluide dans une conduite où l’on a placé : un obstacle
vertical ou une bosse circulaire. Le test de Poiseuille permet d’obtenir une solution
stationnaire. Tandis que les autres tests sont effectués pour le cas transitoire.
Afin d’aider à identifier les caractéristiques physiques des écoulements les équations
de Navier-Stokes sont adimensionalisées, en remplaçant la densité ρ (constante) par
1 et la viscosité dynamique µ par 1/Re. Le nombre de Reynolds Re représente le
rapport des forces d’inertie aux forces de viscosité et est défini par :

Re = U L/η (3.21)

où U est une vitesse caractéristique de l’écoulement étudié (par exemple la

vitesse du mouvement d’un corps ou la vitesse moyenne dans une section déterminée
du canal, etc.), η = µ/ρ est la viscosité cinématique du fluide et L est une longueur
 31

caractéristique (par exemple le diamètre ou la longueur d’un corps, la largeur

intérieure d’un canal, etc...).

3.2.1 Cas test de Poiseuille

Figure 3.1: Domaine de calcul pour le problème de Poiseuille et Maillage (2000

éléments).

Il s’agit de l’écoulement d’un fluide dans un canal 1m × 10m avec un profil

de vitesse parabolique imposée en entrée et d’une viscosité η = 0.002m2 /s. Ce
problème possède une solution analytique suivante :

4y −4x
u=( (H − y), 0), et p = + constante. (3.22)
H2 H2

où H est la hauteur du canal.

Le maillage et les différentes conditions aux limites du domaine sont présentées
sur la figure 3.1.
L’analyse des résultats des tests numériques montre que les champs de vitesses
Figure 3.2 ont bien un aspect parabolique qui confirme l’allure de la solution analy-
tique Figure 3.22, ces résultats sont obtenus au bout de 5P as. Quant à la pression
Figure 3.3 on constate que son allure perturbée n’est pas linéaire contrairement
à la solution analytique (Figure 3.4). Pour pallier à ce problème et améliorer la
stabilité de la solution concernant la pression on a utilisé la méthode de stabité
SUPG (le principe de la méthode est donné dans la section Figure 3.1.4).
 32

Figure 3.2: Isovaleurs de vitesse u dans le canal.

Figure 3.3: Isovaleurs de pression correspondantes.

 33

Figure 3.4: Pression en fonction de x au milieu du canal, comparer avec la pression

analytique et stabiliser par SUPG.

3.2.2 Cas test d’écoulement autour d’un cylindre

Ce problème est caractérisé par un écoulement instationnaire et le développement

d’un détachement tourbillonnaire à l’arrière du cylindre, connue sous le nom d’allées
de Von Karman. Ce problème a été étudié numériquement et expérimentalement
par divers auteurs [27; 15]. Il est aussi considéré comme un test de référence (bench-
mark) pour évaluer la performance des méthodes numériques pour la résolution
des équations de Navier-Stokes instationnaire. En effet, selon les expériences
numériques réalisées par Buffat en 1991, il s’avère que les schémas de discrétisation
en temps trop diffusifs (par exemple les schémas Upwind d’ordre 1) ne sont pas en
mesure de prédire correctement l’évolution temporelle de la solution.
 34

• Description du problème et procédures de calcul

La figure 3.5 présente le domaine de calcul et les conditions aux limites. Le

domaine de calcul a pour dimensions 13m × 7m unités de longueur adimensionnelle
Figure 3.6 avec un cylindre de diamètre égal à 1.

Figure 3.5: Domaine de calcul et conditions aux limites.

Figure 3.6: Maillage du domaine, 11000 noeuds et 23000 éléments.

Un écoulement uniforme (u = 1, v = 0) est imposé à l’entrée et aux frontières

supérieure et inférieure du domaine. Le fluide adhère à la paroi du cylindre. A la
sortie pas de condition. Le maillage est construit en tenant compte de la solution
 35

du problème. Il est suffisamment raffiné dans le sillage du cylindre, afin d’avoir un

champ de vitesse et une pression correcte figures (3.8, 3.9 et 3.10) et de capturer
la formation des tourbillons figure 3.7, comparés avec les résultats [48] figure 3.11.
Le comportement de l’écoulement est étudié en fonction du nombre de Reynolds
Re = DU/η où D est le diamètre du cylindre, U est la vitesse uniforme en entrée
et η est la viscosité cinématique du fluide tel que : η = 0.01m2 /s.

Figure 3.7: Lignes de courant pour Re = 100.

Figure 3.8: Isovaleurs de la vitesse u.

 36

Figure 3.9: Isovaleurs de la vitesse v.

Figure 3.10: Isovaleurs de la pression.

Figure 3.11: Lignes de courant pour Re = 100. Résultats Filippone [48].

 37

Figure 3.12: Allure de vitesse et de la pression en amont de l’obstacle (couleur

jaune et rouge) et en aval de l’obstacle (couleur verte).

Pour étudier la sensibilité de notre modèle au nombre de reynolds, nous avons

choisis le même cas de cylindre avec un maillage homogène Figure 3.13, les champs
de vitesses obtenus pour Re = 100, Re = 400, Re = 1000 et Re = 3000 sont donnés
sur les figures 3.14, 3.15, 3.16 et 3.17. Ce cas test nous a permis d’évaluer la
stabilité de notre élément.

Figure 3.13: Maillage homogène, 2689 noeuds.

 38

Figure 3.14: Champ de vitesses, Re = 100.

£
Figure 3.15: Champ de vitesses, Re = 400.

Figure 3.16: Champ de vitesses, Re = 1000.

 39

Figure 3.17: Champ de vitesses, Re = 2500.

Les figures 3.18 et 3.19 présentent les isovaleurs vitesses en fonction de plusieurs
steps (60, 65, 70 et 75).

Figure 3.18: Isovaleurs vitesses, step = 60, step = 65.

Figure 3.19: Isovaleurs vitesses, step = 70, step = 75.

 40

3.2.3 Cas test d’une Cavité carrée

Dans ce test, consiste en une cavité carrée de dimension 1 × 1 et de longueur

infinie (Figure 3.20 à gauche). La paroi supérieure se déplace à une vitesse uniforme
entraı̂nant ainsi le fluide. Celui-ci forme un large tourbillon dont le centre dépend
du nombre de reynolds. C’est un cas test classique utilisé en simulation numérique
des écoulement laminaire de fluide incompressible. Les résultats seront confrontés
aux résultats de [67] et [19].
Les conditions aux limites en vitesse-pression, sont les suivantes :

u(x, y) = (0, 0) si (x, y) ∈ Γ0 (3.23)

u(x, y) = (1, 0) si (x, y) ∈ Γ1

Figure 3.20: Géométrie de la cavité et le maillage (30 × 30).

Nous représentons sur les figures suivantes les champs des composantes de la
vitesse suivant les axes x et y figure 3.21 pour Re = 1 et figure 3.22 pour Re = 100.
 41

Figure 3.21: A gauche : vitesse u. A droite : vitesse v pour Re = 1.

Figure 3.22: A gauche : vitesse u. A droite : vitesse v pour Re = 100.

Les variations des vitesses u, v en fonction de y, x respectivement pour plusieurs

coupes Figure 3.23.
 42

Figure 3.23: A gauche : vitesse u(y). A droite : vitesse v(x) pour Re = 100.

Sur la figure 3.24 à droite, nous avons représentés les profils de la composante
de la vitesse u en fonction de y pour x=0.5 et ce pour les nombres de Reynolds de
100 et 400. Les résultats des simulations sont en bon accord avec les résultats de
[67] et [19].

Figure 3.24: Vitesse en fonction de nombre de Reynolds.

 43

3.2.4 Cas test d’écoulement autour d’un obstacle

Les problèmes d’écoulement autour d’un obstacle font l’objet de plusieurs ap-
plications. Des exemples de tels écoulements se retrouvent dans des problèmes
d’échangeurs de chaleur, d’écoulement du vent autour de maison, écoulement dans
les artères en présence de plaques, etc. Parmi cette classe d’écoulement, le problème
de l’écoulement autour d’un obstacle rectangulaire entre deux plaques est con-
sidéré dans cette étude. Le domaine de calcul est composé d’un canal avec un
rétrécissement (cas d’une marche) suivi d’un canal avec un élargissement (cas d’une
marche inverse) figure 3.25.
Ce cas test a été choisi pour les raisons suivantes : d’abord, la géometrie est
simple et l’écoulement est caractérisé par la présence de plusieurs zones de re-
circulation tel que les zones fermées et des allées [2]; et finalement, à travers les
expériences numériques réalisées par [86], il s’avère que ce problème est un test
plus contraignant que les problème standard de la marche inversée. En effet, con-
trairement au problème de la marche inversée, les erreurs numériques provenant de
l’avant de l’obstacle peuvent influencer de façon très significative la prédiction de
l’écoulement dans tout le domaine.

• Description du problème et procédures de calcul

Le domaine de calcul est constitué d’un canal à l’intérieur du quel se trouve un

obstacle de hauteur H et de largeur h, les conditions aux limites sont illustrées sur
la figure 3.25. Le fluide rentre dans le canal avec un profil de vitesse parabolique

3
u(y) = ( y(2 − y), 0) pour 0 ≤ y ≤ 2, (3.24)
2

adhéré aux parois solides (condition de non-glissement). Le nombre de Reynolds

b = 2 umax =0.1 m/s) et la hauteur H
est basé sur la vitesse moyenne à l’entrée (U 3

de l’obstacle. Le maillage uilisé (4000 noeuds et 8000 éléments) est bien raffiné
 44

autour de l’obstacle afin de capturer les caractéristiques de base de l’écoulement

(vitesses, pression et lignes de courant). Les figures 3.26, 3.27, 3.28, 3.29, 3.30 et
3.31 représentent respectivement la pression, la vitesse u et les lignes de courant
obtenues par notre modèle et comparé aux résultats obtenues par TDYN.

Figure 3.25: Maillage de domaine (4000 noeuds et 8000 éléments) et conditions

aux limites.

Figure 3.26: Pression obtenue par notre modèle, Re = 1.

Figure 3.27: Pression obtenue par TDYN, Re = 1.

 45

Figure 3.28: Vitesse u obtenue par notre modèle, Re = 10.

Figure 3.29: Vitesse u obtenue par TDYN, Re = 10.

Figure 3.30: Les lignes de courant obtenues par notre modèle Re = 10.

Figure 3.31: Les lignes de courant obtenues par TDYN Re = 10.

 46

Dans ce test, on cherche à prédire les caractéristiques de l’écoulement (zone

de recirculation et la composante horizontale de la vitesse à différentes stations en
aval de l’obstacle. Notre solution est comparée aux résultats numériques fournies
par N’Dri [96], Armaly [3] et Cruchaga [35].
Notre cas tests est composé de deux marche, nous avons choisis après la marche
pour évaluer la longueur de la zone de recirculation. Pour cela nous considérons le
domaine et les conditions aux limites donnés sur la figure 3.32

Figure 3.32: Domaine et conditions aux limites.

La longueur de la zone de recirculation est estimé selon les auteurs entre cinq
et six fois la hauteur h de la marche. Les résulats obtenus pour Re = 100, 250, 500
sont représentés sur les figures 3.33, 3.34 et 3.35. Pour Re = 250 la longueur
obtenue (5.66 × h) par notre modèle est résonable par rapport aux résultats de
Armaly [3] et Cruchaga [35].

Figure 3.33: Les lignes de courant, Re = 100.

 47

Figure 3.34: Les lignes de courant, Re = 250.

Figure 3.35: Les lignes de courant, Re = 500.

Dans le canal avec l’obstacle vertical on observe la présence de zone de recir-

culation le long de paroi solide. Sur la figure 3.37 la distribution de vitesse est
comparée aux résultats numériques et expérimentaux [96].
 Figure 3.36: Coupe de la composante horizontale de la vitesse, nos résultats.

Figure 3.37: Coupe de la composante horizontale de la vitesse, résultats Kouakou

[96].
 49

CHAPITRE 4

MODÈLE SOLIDE : TRAITEMENT NUMÉRIQUE ET

VALIDATION

4.1 Introduction

Dans ce module, il a été considéré le cas de structures minces (coque). Toute

fois, une coque est une structure tridimensionnelle caractérisée par une épaisseur
très petite par rapport aux autres dimensions. Ce type de structures apparaı̂t
fréquemment dans des constructions courantes (ponts, toits de bâtiments,...), dans
la conception industrielle (turbines, pièces de mécanique, carrosserie de voiture),
et même en biomécanique (artères, bronches,...).
Afin de simuler numériquement le comportement des structures minces, il est
important de formuler des méthodes d’éléments finis de coques qui soient ro-
bustes vis-à-vis du verrouillage numérique, mais aussi consistantes. On trouve
fréquemment des structures minces de géométrie complexe qui nécessitent l’utilisation
de maillages surfaciques non-structurés qui incluent forcément des éléments finis
triangulaires [88].
Il est bien connu que les techniques d’élements finis standards marchent bien
dans des situations à membrane dominante, mais en général celles-ci donnent des
solutions approchées trop rigides pour les problèmes à flexion dominante, prin-
cipalement lorsque l’épaisseur est très petite, phénomène connu sous le nom de
verrouillage numérique [32]. Le verrouillage numérique est l’une des difficultés
majeures dans la formulation d’éléments finis robustes pour l’analyse numérique
des structures minces, mais ce phénomène concerne aussi les formulations incom-
pressibles de la mécanique des milieux continus en particulier [9], [25], [10] et [69].
Des nombreux travaux ont été réalisés dans la recherche de remèdes contre cette
 50

pathologie. Les méthodes mixtes dont les fondements théoriques ont été établis par
[6] et [25] constituent un outil efficace pour la conception et l’analyse rigoureuse
des méthodes numériques qui échappent au verrouillage, en particulier dans la
mécanique des fluides incompressibles et les modèles de plaques [25], [33] et [112].

4.2 Synthèse bibliographique

L’obtention des modèles de plaques et coques a été un sujet de recherche pri-

mordial en mécanique [28], [53], [91], [108], [124], [126] et [134]. L’idée de base dans
la théorie de ces modèles est d’utiliser des hypothèses et simplifications physique-
ment justifiables à travers l’épaisseur. Le but visé est d’obtenir la déformation
d’une structure mince tridimensionnelle à partir d’un problème formulé sur sa
surface moyenne. Par exemple, on néglige l’influence de la composante normale
des contraintes, les effets des déformations à travers l’épaisseur, ou on utilise sou-
vent des hypothèses cinématiques du type Reissner-Mindlin qui supposent que les
composantes tangentielles u, v du déplacement sont linéaires en z, tandis que la
composante transversale est constante [77], [87], [92] et [109]. Mais on peut ausssi
penser à des hypothèses cinématiques d’ordre plus élévé, en particulier [31].
Les problèmes de robustesse des méthodes d’éléments finis pour les coques, et
notamment le phénomène de verrouillage numérique de cisaillement et de mem-
brane qui intervient lorsque la structure considérée est mince, sont désormais bien
compris, en particulier [32] et ses références. Même si l’on ne connaı̂t pas, pour
l’instant, d’élément dont on ait pu démontrer qu’il serait parfaitement robuste, on
dispose de méthodologies détaillées et rigoureuses pour évaluer la robustesse des
éléments de coques au moyen de cas-tests [14], [32]. De fait, certains éléments ex-
istants comme l’élément quadrangulaire Q4γ24 ”Elément de Coque” [90] montrent
un très bon comportement à travers ces tests. Cet élément est développé dans la
présente version de REFLEX , utilisée dans le cadre des présents travaux.
 51

4.3 Description des éléments finis utilisés dans le modèle structure

4.3.1 Eléments de coques de type Kirchhoff discret

Les premiers éléments de plaques et de coques développés étaient basés en

grandes parties sur l’hypothèse de Kirchhoff-Love. Cette hypothèse, qui consiste à
annuler l’energie due aux déformations de cisaillement transverse, exige une conti-
nuité C 1 (pour un modèle déplacement compatible). Ce modèle, qui ne s’applique
que pour des structures minces, a été utilisé par [135]. La continuité C 1 , qui exige
la continuité de w et des dérivées premières, est très difficile à assurer [140]. Pour
éviter ces difficultés, plusieurs approches ont été proposées. Parmi ces approches,
on trouve la classe des éléments de type Kirchhoff discret. Cette approche permet
de satisfaire la continuité inter-éléments des variables essentielles (u, v, w, βx , βy )
et permet encore d’obtenir des résultats numériques très satisfaisants. Les hy-
pothèses de Kirchhoff sous forme discrète furent introduites par Wempner et al.
[133] et Dhatt [39]. En 1978 de nombreux chercheurs, développeurs de codes et
d’industriels ont porté un intérêt marqué pour les éléments triangulaires DKT et
quadrilatéraux DKQ. Ces éléments ont été utilisés pour l’analyse des plaques et
des coques minces dans les domaines linéaires et non linéaires , constituées de
matériaux isotropes ou composites et soumises à des sollicitations statiques et dy-
namiques Fafard et al. [46]. Les éléments de type Kirchhoff discret ont été proposés
comme une alternative pour définir des éléments à faible nombre de noeuds et de
ddl sans utiliser la théorie continue de Kirchhoff. Dans ce cadre, nous présentons
deux éléments finis de coques minces , l’un de type quadrilatéral nommé DKQ20 et
l’autre triangulaire noté DKT 15. Ces derniers sont basés sur une approche tridi-
mensionnelle dégénérée qui tient compte de l’hypothèse de Kirchhoff sous forme
discrète, et permettent de simuler le comportement des coques minces élastiques
de formes quelconques.
 52

[Link] Description géométrique de DKQ20 et DKT 15

L’interpolation géométrique d’un vecteur position d’un point q s’exprime en

fonction des coordonnées ζ et η de l’élément de référence sous la forme suivante :

1
Xq = Xp + ζhn (4.1)
2

X X hi
Xq (ξ, η, ζ) = Ni (ξ, η) Xi + ζ Ni (ξ, η) ni (4.2)
i=1,nd i=1,nd
2

nd : nombre de noeud (3 pour DKT 15 et 4 pour DKQ20)

Ni : fonctions de formes linéaire en ζ et η :
hN1 N2 N3 i = h1 − ζ − η ζ ηi pour DKT 15
hN1 N2 N3 N4 i = 14 h(1 − ±ζ) (1 − ±η)i pour DKQ20

Figure 4.1: Géométrie et variables nodales de l’élément DKQ2.

 53

Figure 4.2: Géométrie et variables nodales de l’élément DKT 15.

[Link] Champs de déplacements

Le champ de déplacements d’un point quelconque de la coque s’écrit avec

l’hypothèse de sections droites de Hencky/Mindlin/Reissner :

h
Uq (ξ, η, ζ) = Up (ξ, η) + ζ β (ξ, η) (4.3)
2

Up : vecteur déplacement du point p (ζ = 0)

β (ξ, η) : vecteur rotation orthogonal à la normale, s’écrit également :

β = θ ∧ n = (θ1 t1 + θ2 t2 ) ∧ n = θ2 t1 − θ1 t2 (4.4)

L’interpolation des déplacements est faite à partir des vecteurs déplacements

aux noeuds et des rotations autour des vecteurs tangents t1i et t2i :

2nd
!
X ζ X X
Uq = Ni Upi + Ni hi (θ2i t1i − θ1i t2i ) + Pk tsk hk αk (4.5)
i=1,nd
2 i=1,nd k=nd+1

Pk sont des fonctions linéaires et quadratiques compatibles associées respective-

 54

ment aux éléments DKT 15 et DKQ20, représentées sur le tableau ci-dessous.

Eléments Pk (k = nd, 2nd)

P4 = 4ξ (1 − ξ − η)
DKT 15 (k = 4 à 6) P5 = 4ξη
P6 = 4η (1 − ξ − η)
P5 = 12 (1 − ξ 2 ) (1 + η)
P6 = 21 (1 − η 2 ) (1 + ξ)
DKG20 (k = 5 à 8)
P7 = 21 (1 − ξ 2 ) (1 − η)
P8 = 21 (1 − η 2 ) (1 − ξ)

[Link] Champs de déformations

Les déformations en un point de coordonnées ξ, η, ζ s’écrivent :

{s } = {e} + ζ {χ} (4.6)

Pour DKQ20 et DKT 15 les expressions des composantes curvilignes de déformations

dans la base [Q] sont :

• déformations de membrane :

hei = hex ey exy i (4.7)

ex = t1 up,x ; ey = t2 up,y ; exy = t2 up,x + t1 up,y

• déformations de flexion (courbures) :

hχi = hχx χy χxy i (4.8)

h


χx = t1 (up,ξ bc11 + up,η bc21 ) + t1 2
β ,x
h


χy = t2 (up,ξ bc12 + up,η bc22 ) + t2 2
β ,y
h h



χxy = t2 (up,ξ bc11 + up,η bc21 ) + t1 (up,ξ bc12 + up,η bc22 ) + t2 2
β ,x + t1 2
β ,y
 55

En tenant compte des hypothèses de Kirchhoff, les déformations de cisaillement

transversal constantes en ζ sont négligeables. Les expressions γxz = 0 et γyz = 0
peuvent être utilisées pour relier β aux gradients de up .

[Link] Matrices de rigidité

En introduisant les expressions des déformations en fonction des variables nodales,

on peut obtenir les matrices de rigidité en intégrant explicitement suivant ζ tel que
(−1  ζ  1) :

m f
W (up , β) = Wint − Wext avec : Wint = Wint + Wint (4.9)

Sachant que :


 W m = R e he∗ i [Hm ] {e} dA = hu∗ i [km ] {un }
int A n
(4.10)
 Wf = R 4
hχ∗ i [H ] {χ} dA = hu∗ i [k ] {u }
int Ae h2 f n f n

avec :


 [km ] = R e [Bm ]T [Hm ] [Bm ] dA
A
(4.11)
 [k ] = R 4 [B ]T [H ] [B ] dA
f Ae h2 f f f

√
où dA = adξdη
La matrice de rigidité élémentaire (DKT 15 ou DKQ20) s’exprime par :


 [k e ]
e (15×15) pour DKT 15
[k ] = [km ] + [kf ] ⇒ (4.12)
 [k e ]
(20×20) pour DKQ20
 56

Pour DKT 15 :
* [km ] représente la rigidité de membrane.
* [kf ] représente la rigidité de flexion qui est obtenue en considérant un schéma
d’intégration numérique de Hammer à trois points suivant ξ, η.
* h représente l ’épaisseur de la coque.

Pour DKQ20 :
* [km ] représente la rigidité de membrane. Cette matrice est relative à Q4.
* [kf ] représente la rigidité de flexion qui est obtenue en considérant un schéma
d’intégration numérique de Gauss à deux points suivant ξ, η.
Nous remarquons que les termes bcij dans l’expression des courbures de flexion
(χx , χy et χxy ) interviennent que les éléments soient coplanaires ou non.

4.3.2 Eléments de coques de type Reissner-Mindlin

La théorie de Reissner-Mindlin consiste à tenir compte du cisaillement transver-

◦
sal. Cette théorie permet de considérer une continuité C des déplacements (et
rotations) alors qu’une continuité C 1 est nécessaire pour un élément conforme basé
sur la théorie de Kirchhoff. Elle offre aussi un domaine d’application plus large
◦
(coques minces et épaisses). Plusieurs éléments finis de coque de continuité C ont
été construits selon cette théorie, et sont actuellement les plus utilisés même dans
le cas où la structure modélisée vérifie bien l’hypothèse de Kirchhoff. Parmi les
◦
éléments à continuité C , la formulation introduite par Ahmad [1], [122] permet
de formuler une famille d’éléments de coque isoparamétrique de façon commode et
efficace qui assure la continuité géométrique. En se basant sur cette cinématique,
nous avons développé deux éléments finis de coque avec cisaillement transverse,
l’un de type quadrilatéral nommé Q4γ20 et l’autre triangulaire noté T 3γ15. Ces
derniers restent dans le cadre d’une formulation en déplacement et s’appuient sur
des techniques dites ”d’interpolation des distorsions transverses”.
 57

[Link] Elément Q4γ20

Concerne la formulation de l’élément quadrilatéral à quatre noeuds non coplanaires

ayant comme degrés de liberté (ddl) : trois composantes de déplacement et deux
rotations. La formulation est inspirée de l’élément Q4γ24 développé par Batoz
et Dhatt [11]. Cet élément est représentée par la surface moyenne A et la fibre
épaisseur hn. Le vecteur position d’un point quelconque s’écrit en fonction des
coordonnées isoparamétriques [12] :

1
xq (ξ, η, ζ) = xp (ξ, η) + ζhn (4.13)
2

xp est le vecteur position du point p de la surface moyenne, défini par :

X
xp (ξ, η) = Ni (ξ, η) xi
i=1,4

Figure 4.3: Elément de coque courbe Q4γ20.

 58

[Link] Elément T 3γ15

Nous présentons brièvement la définition de l’élément triangulaire à trois noeuds,

simple et efficace mais qui n’est pas parfait. L’élément T 3γ15 présente une pseudo-
normale continue et donc une pseudo-courbure. Il est ainsi bien adapté aux géométries
courbes, et il a comme ddl : trois composantes de déplacements et deux com-
posantes de rotations. L’élément T 3γ15 est inspiré de l’élément développé par
Boisse et al. [20], lui même inspiré de l’élément de Hughes et al. [71]. T 3γ15 est
complètement définie par sa surface moyenne, sa fibre épaisseur hn. Le vecteur
position d’un point quelconque s’écrit :

1
xq (ξ, η, ζ) = xp (ξ, η) + ζhn (4.14)
2

4.3.3 Elément à deux dimention en coordonnées cartisiennes, Elément

triangulaire à trois noeuds (T3)

Figure 4.4: Elément de référence à gauche et l’élément réel à droite

Les fonctions d’interpolation sont :

hN i = h1 − ξ − η ξ ηi
 59

La matrice jacobienne inverse [j] étant :

   
∂ξ ∂η
y31 −y21
[j] =  ∂x ∂x = 1   (4.15)
∂ξ ∂η 2A −x x21
∂y ∂y 31

J = det[J] = 2A = y31 x21 − x31 y21

Les déformations dans le plan x, y sont :

 


 u,x 


 
{ε} = v,y = [B] {un } (4.16)

 

 
u,y + v,x
 

avec :

hun i = hu1 v1 : u2 v2 : u3 v3 i

 
y23 0 y31 0 y12 0
1  
[B] = (4.17)
 
 0 x32 0 x13 0 x21
2A 


x32 y23 x13 y31 x21 y12

Le domaine d’intégration est :

h
Z Z 1 Z 1−ξ Z
2
(...)dv = ( (...)dz)2Adηdξ (4.18)
ve 0 0 −h
2

La matrice de rigidité est :

Z 1 Z 1−ξ
T
[k] = [B] ( [Hm ] dηdξ) [B] 2A (4.19)
0 0

avec :
 60

Z h
2
[Hm ] = ([H])dz (4.20)
−h
2

où h est l’épaisseur dans la direction z et [H] dépend des hypothèses (contraintes
ou déformations planes) et du matériau . Pour un matériau isotrope :

 
H1 H2 0
 
[H] =  H2 H1 0 (4.21)
 

 
0 0 G

E(1 − aυ)
H1 = (4.22)
(1 + υ) (1 − v − aυ)
υH1
H2 = (4.23)
1 − aυ

Avec : E module d’élasticité ; υ coefficient de Poisson.

E
G= (4.24)
2 (1 + υ)

a = 0 en contraintes planes (CP) ; a = 1 en déformation planes (DP).

La matrice [k] est donnée explicitement dans le tableau 2.6.4 (pour[H] définie
par 2.6.35 et h constant).
La matrice[m] associée aux variables nodales u1, u2, u3 (idem pour v1, v2, v3)
est :

Z 1 Z 1−ξ
[m] = ρm {N } hN i 2Adηdξ
0 0

avec :
Z h
2
ρm = ([H])dz (4.25)
−h
2

où ρ est la masse volumique en un point x, y, z.

Pour ρm constant :
 61

 
2 1 1
1  
[m] = ρm A  1 2 1  (4.26)
 
12  
1 1 2

La matrice masse diagonale est simplement :

 
1 0 0
A 
[m] = ρm (4.27)
 
 0 1 0
3


0 0 1

Le vecteur sollicitation dû à des forces de volume fx (x, y) (suivant x) et fy (x, y)

(suivant y) est :

Z 1 Z 1−ξ
hfnv i = hN1 fx N1 fy N2 fx N2 fy N3 fx N3 fy i h.2A.dη.dξ (4.28)
0 0

Pour fx et fy constants :

Ah
hfnv i = hfx fy fx fy fx fy i (4.29)
3

Le vecteur dû aux contraintes initiales [σ0 ]et thermiques [σth ] est :

h
!
Z 1 Z 1−ξ Z
2
{fσ } = − [B]T ({σ0 } + {σth })dz 2A.dη.dξ (4.30)
0 0 −h
2

hσ0 i = hσ0x σ0y σ0xy i (4.31)

hσth i = hσthx σthy σthxy i (4.32)

Les contraintes thermiques dépendent du matériau et des hypothè ses de con-

 62

traintes ou déformations planes. Pour un matériau homogène isotrope :

EαT
hσth i = − h1 1 0i (4.33)
1 − aυ

avec : a = 1 en (CP) et a = 2 en (DP).

Les contraintes σx , σy , σxy en un point x, y, z de l’élément sont :

{σ} = [H] [B] {un } + {σ0 } (4.34)

En général, les contraintes ainsi calculées sont discontinues d’un élément à

l’autre. On affecte les contraintes au centre de gravité ({σ} étant constant si{σ0 }
est constant ). Un champ de contraintes continues est obtenu par une technique
de lissage (moyennes avec ou sans pondération, moindres carrés, ...).

4.4 Résultats numériques : tests de validation

4.4.1 Modèles éléments finis

Parmi les objectifs visés dans les présents travaux de thèse concerne la formu-
lation d’éléments finis de coque avec continuité complète de la géométrie avec une
intégration explicite dans l’épaisseur. En s’inspirant de la formulation d’éléments
dits isoparamétriques dégénérés décrite en [12], nous avons développé des éléments
de type triangulaires et quadrilatéraux à cinq degrés de liberté par noeud. Certains
éléments prennent en compte les effets de cisaillement transversal (T3γ 15, Q4γ
20) ; d’autres sont obtenus en considérant les hypothèses de Kirchhoff sous forme
discrète (DKT15, DKQ20).
Les principales caractéristiques des éléments utilisés sont :
- la variation linéaire des composantes de déformations en ζ , permettant une
intégration explicite dans l’épaisseur,
 63

- la prise en compte du gauchissement de la surface moyenne dans le cas des

éléments quadrilatéraux DKQ20 et Q4g 20,
- la définition de cinq degrés de liberté par noeud qui sont les trois composantes
de translation U, V, W et deux composantes de rotations θ1 , θ2 .

4.4.2 Panneau cylindrique soumis à son poids propre.

Le problème représenté à la figure 4.5 est un problème qui a été le plus étudié
pour vérifier les formulations et comparer les performances d’éléments plats ou
courbés, triangulaires, cylindriques ou isoparamétriques. Il s’agit d’un panneau
cylindrique, modilisé par Naceur. H. Professeur à l’UTC, dont les bords rectilignes
sont libres et dont les côtés courbés reposent sur des diaphragmes rigides dans leur
plan. La coque est soumise à une charge verticale de gravité (poids propre) de fz =
−0.625 × 104 P a . Ce problème décrit dans plusieurs références [1, 12, 45, ...], sert à
tester l’aptitude d’un élément à simuler des états de déformations et de contraintes
complexes avec une certaine domination de la partie membrane. Les rapports
L/h = 200 et R/h = 100 sont importants et les déformations de cisaillement
transverse sont négligeables. Les déformations de membrane sont importantes par
rapport à celles de flexion, d’où un grand risque de blocage en membrane et en
cisaillement.
Dans la littérature, deux solutions sont proposées. L’une est basée sur la théorie
des coques surbaissées [118] ; l’autre est obtenue en considérant une théorie des
coques profondes. La différence entre les deux théories étant de l’ordre de 3%:
Valeurs de référence (théorie des coques ”profondes”) :
WB = −3.61cm ; WC = 0.541cm
Solution analytique (théorie des coques surbaissées) :
WB = −3.703cm ; WC = 0.525cm
Les données et les conditions aux limites dans ce cas sont :
 64

L = 6m ; R = 3m ; h = 0.03m ; E = 3 × 1010 P a ; φ = 40◦

U = W = φ2 sur AD ; U = φ2 = 0 sur CD ; V = φ1 = 0 sur CB

Figure 4.5: Panneau cylindrique soumis à son poids propre.

Nous avons considéré le quart de la structure à cause de la double symétrie en

considérant des maillages réguliers avec N = 2, 4, 6, 8, 10, 12, 16 et 20 éléments suiv-
ant AB et BC. L’orientation des éléments triangulaires est indiquée sur la figure
4.5. Le maillage de type A (Figure 4.6) est utilisé pour évaluer les performances
des éléments triangulaires.
 65

Figure 4.6: Type de maillage utilisé.

Figure 4.7: Convergence du déplacement W au bord libre B. Eléments à 3 noeuds.

 66

Figure 4.8: Convergence du déplacement W au centre C. Eléments à 3 noeuds.

Les courbes de convergence en terme de DDL total des déplacements verticaux

au point B (milieu du bord libre), et au point C (centre de la structure) sont
présentées sur les figures 4.7, 4.8 pour les éléments à trois noeuds et sur les figures
4.9, 4.10 pour les éléments à quatre noeuds.
 67

Figure 4.9: Convergence du déplacement W au bord libre B. Eléments à 4 noeuds.

Figure 4.10: Convergence du déplacement W au centre C. Eléments à 4 noeuds.

 68

Le déplacement vertical W suivant C −B et le déplacement tangentiel V suivant

D−A sont présentés sur les figures 4.11 et 4.12 en considérant un maillage de 20×20
avec les éléments DKT 15, T 3γ15, Q4γ20 et DKQ20.

Figure 4.11: Déplacement V le long de DA.

Figure 4.12: Déplacement W le long de CB.

 69

4.4.3 Poutre 2D encastré-libre, soumise à un effort concentré :

Le but de ce test est de valider le modèle E.F utilisé (Q4 mixte-hybride). Pour
cela, on considère une poutre 2D tels que les déplacements suivant pour x = 0 sont
nuls (Figure 4.13).

Figure 4.13: Domaine et maillage (707 noeuds, 600 éléments).

Les données dans ce cas sont les suivantes :

E = 2, 1.1011 P a ν = 0.3 L = 1m h = 0.03m
Les conditions aux limites sont :
En x = 0, u = v = 0 ; En x = L, F = −1000N
L’épaisseur de la poutre est égale à 0.03m et le moment est 6.75.10−8 m4 .
Le domaine déformé (×0.005) et les déplacements suivant l’axe y sont données
sur les figures 4.14 et 4.15.

Figure 4.14: Domaine déformé, f est la flèche).

 70

Figure 4.15: Vecteurs déplacements vy en mm.

La flèche obtenue par notre modèle est égale à f = −22.487mm, cette valeur
est comparée à la solution calculé, ci-dessous, par la résistance des matérieux.
On sait que : EIy ” = −M tel que : M (x) = −F (L − x), on obtient :

EIy ” = −F (L − x) (4.35)

2
Une première intégration donne : EIy 0 = F Lx−F x2 +c1 , suivie d’une deuxième
intégration qui donne l’évolution de la flèche en fonction de x et des constantes c1
et c2 . soit :

x2 x3
EIy = F L −F + c1 x + c2 (4.36)
2 6

Les coditions aux limites sont :


 y 0 (0) = 0
(4.37)
 y (0) = 0

On ontient : c1 = c2 = 0 et donc :

F 2
y= x (3L − x) (4.38)
6EI
 71

Pour ymax , la flèche maximale est :

F L3
f= (4.39)
3EI
F L3
Dans notre cas test, le flèche au point x = L est : f = 3EI
= −23.516mm.
Nous avons comparés l’évolution de cette flèche avec la solution analytique
F.x2
y= 6EI
(3L − x), la figure 4.16 représente la solution calculé par notre modèle en
noir et la solution analytique en rouge. Une concordance parfaite entre les deux
résultats.

Figure 4.16: L’évolution de la solution y = f (x).

 72

4.4.4 Poutre 2D en flexion encastré-encastré.

Soit une poutre droite de longueur L = 5m, de section constante 0.3m ×0.005m
et de moment d’inertie I, qui supporte en son milieu une charge F = 10000N
(Figure 4.17).

Figure 4.17: Domaine et maillage (707 noeuds et 600 éléments).

• Effort tranchant (T) :

L
si 0 ≺ x ≺ 2
⇒ T = − F2
L
si 2
≺ x ≺ L ⇒ T = − F2 + F = F
2

• Moment fléchissant (M) :

L
si 0 ≺ x ≺ 2
⇒ M = − F2x
L Fx
si 2
≺x≺L⇒M = 2
− F (x − L2 ) = F
2
(L − x)
Le moment fléchissant est au maximum au milieu de la poutre (x = L/2) est :
FL
Mmax = 4

• Calcul des flèches :

On sait que EIY ” = −Mf

Si 0 ≺ x ≺ L2 , M = − F2x , donc EIY ” = − F2x
Une première intégration donne : EIY 0 = − F4 x2 + c1
 73

Quand la flèche est à son maximum, x = 21 , Y 0 = 0, alors :

F F L2
EIY 0 = − x2 +
4 16
F 3 F L2
Une deuxième intégration donne : EIY = − 12 x + 16
x + c2
La poutre est toujours en appui, pour x = 0, Y = 0 donc c2 = 0.
L’équation de la déformée est :

1 F 3 F L2
Y = ( x + x) (4.40)
EI 12 16

Pour Ymax :

1 F L 3 F L2 L F L3
Y =f = (− + )= (4.41)
EI 12 8 16 2 48EI

Avec : E = 2.1010 P a le module de Young et I = 4.457.10−10 m4 le moment

quadratique.
On obtient :

10000 × 53
f= = 2.2735.10−3 m (4.42)
48 × 2.1010 × 4.457.10−10 m4

La figure 4.18 présente les vecteurs déplacements. La valeur maximale de la

flèche obtenue par notre modèle est :

f = 2.3585.10−3 m (4.43)

Sur la figure 4.19, l’évolution de la solution obtenue (en noir) et la solution

analytique de l’équation 4.40 (en rouge).
 74

Figure 4.18: Vecteurs déplacements et domaine déformé (×200).

Figure 4.19: L’évolution de la solution calculée par notre modèle et la solution

analytique.
 75

4.4.5 Déformation d’une arche elliptique

Ce cas test [41] concerne l’étude de la répartition des déplacements dans une
arche elliptique de pont (R1 × R2 = 6m × 3m). Les seuls déplacements possibles
sont suivant l’axe y. Les paramètres choisis dans cet exemple sont : E = 2.1010 P a
et ν = 0.3. L’arche a pour hauteur H = 3.5m et pour longueur L = 6.5m. Le
domaines, les conditions aux limites et le maillage choisi pour ce cas sont présentés
sur les figures 4.20 et 4.21.

Figure 4.20: Géométrie et conditions aux limites d’après G. Dhatt [41].

 76

Figure 4.21: Maillage : 181 noeus et 150 éléments.

Les vecteurs déplacements, suivant l’axe y, obtenus par notre modèle sont
présentés sur le figure 4.22.

Figure 4.22: Vecteurs déplacements (mm) suivant l’axe y et maillage déformé

(facteur= 200).

Nous avons comparés les résultats obtenus pour les troix noeuds de l’extrémité
gauche de l’arche avec cels obtenus par G. Dhatt [41], le tableau ci-dessous montre
un accord raisonable entre les deux résultats.
 77

Noeuds Déplacement-Reflex Dé[Link][41]

1 -0.0011542 -0.0011428
2 -0.0011503 -0.0011456
3 -0.0011561 -0.0011440

Figure 4.23: Les déplacements suivant l’axe y pour N1 , N2 et N3 .

De ces résultats on constate que la déformation maximale est située à l’extrémité

opposée du point d’encastrement de la structure. Nous avons donc comparé la
déformation simulée avec les résultats donnés dans [41], au niveau où la déformation
est maximale. L’erreur entre les deux résultats est de l’ordre de (1 − 0.99) = 1%
soit :

|∆X|1 = |Xf − Xi | = 0.11542.E −2 ; nos simulations

|∆X|2 = |Xf − Xi | = 0.11428.E −2 ; d’après [41] (4.44)
|∆X|1
rapport = |∆X|2
= 0.9901
 CHAPITRE 5

COUPLAGE FLUIDE-STRUCTURE

5.1 Formulations eulériennes, lagrangiennes et mixtes

5.1.1 Formulation eulérienne

La description eulérienne étudie le passage d’une certaine quantité de matière à

travers une région donnée de l’espace. Le système de référence associé est fixe et est
appelé Système de Référence Spatial (SRS). Dans cette description, l’observateur se
place en un point fixe et étudie comment varient les caractéristiques de la particule
passant par ce point à l’instant considéré. Cette formulation est généralement celle
qui est adoptée en mécanique des fluides.
Dans le cadre de la méthode des éléments finis, la formulation eulérienne im-
pose un maillage fixe dans l’espace. Ce maillage définit un volume de contrôle
invariable au cours du temps et une particule matérielle donnée peut migrer d’un
élément fini vers un autre. Certaines particules matérielles peuvent même quitter
le domaine défini initialement et, dans ce cas, toute information ultérieure relative
à ces particules est perdue.
L’utilisation de formulations eulériennes fait apparaı̂tre des difficultés lors de
l’étude de problèmes à frontières variables. Les frontières du corps étudié ne
coı̈ncident généralement ni partout ni à chaque instant avec le bord d’un élément.
Il en résulte de grandes difficultés pour prendre en compte les conditions aux lim-
ites matérielles ainsi qu’une définition peu précise des frontières matérielles lorsque
celles-ci évoluent au cours du temps.
En revanche, les formulations eulériennes autorisent d’importantes distorsions
de la matière. Ceci explique pourquoi elles sont largement utilisées en mécanique
 79

des fluides ou dans le cas de la mise en forme des matériaux pour l’étude des
processus stationnaires.

5.1.2 Formulation lagrangienne

La description lagrangienne étudie le mouvement du corps considéré comme

une certaine quantité de matière se déplaçant dans l’espace à partir d’une référence
connue. Le système de référence, appelé Système de Référence Matériel (SRM), est
donc attaché à la matière. Dans ce cas, l’observateur peut être considéré comme
attaché à une particule matérielle donnée et il étudie l’évolution temporelle des
propriétés de cette particule au cours de son mouvement. Cette formulation est
classiquement adoptée en mécanique des solides.
Dans le cadre de la méthode des éléments finis, la formulation lagrangienne
impose à chaque nœud de rester en coı̈ncidence avec la même particule matérielle
tout au long du processus de déformation. Cette propriété est intéressante d’un
point de vue numérique puisque, lors de l’intégration des lois de comportement, les
termes convectifs n’apparaissent pas dans les équations. D’autre part, un élément
donné contient toujours le même ensemble de particules matérielles et est donc
entraı̂né (convecté) par la déformation. En particulier, un nœud situé à la frontière
du corps reste sur cette frontière au cours de la simulation et, dans ce cas, le suivi
d’une frontière mobile ne présente pas de difficultés.
En revanche, en présence de grandes déformations, le maillage peut être forte-
ment distordu et entraı̂ner une imprécision des résultats voire l’impossibilité de
conduire la simulation à terme. Des opérations de réactualisation de maillage
(remaillage) sont alors indispensables, opérations coûteuses en temps homme et
machine.
 80

5.1.3 Formulation ALE

[Link] introduction et revue bibliographique

La formulation ALE est particulièrement bien adaptée à l’étude des problèmes

présentant de grandes déformations et plus généralement à l’étude des problèmes
couplés.
Elle permet de combiner les avantages de la représentation lagrangienne (bonne
définition des frontières variables, facilité d’imposer des conditions aux limites
matérielles) et ceux de la formulation eulérienne (possibilité de prise en compte de
grandes distorsions). Le principe de cette formulation mixte est le développement
de maillages se déplaçant à taux indépendant (bien que non nul) de celui de la
matière. En conséquence, la matière se déplace par rapport au maillage, comme
dans une formulation eulérienne pure, mais la forme extérieure des éléments est
contrôlée par les conditions aux limites du problème, comme dans une formulation
lagrangienne.
Hirt et al. [64] utilisent une technique ALE associée à des discrétisations par
différences finies pour l’étude de l’effet de chocs dans des tuyaux, d’écoulements
autour d’obstacles, d’écoulements instables ou d’écoulements dans l’atmosphère
soumise à une explosion. Hughes et al. [70] utilisent une discrétisation éléments
finis pour le calcul d’écoulements à surface libre. Les degrés de liberté normaux
à la surface libre sont lagrangiens et ceux tangents sont eulériens. Ailleurs, ils
sont purement eulériens. Ce choix permet de gérer le problème de frontières mo-
biles. Brooks et Hughes [27] développent une formulation éléments finis qui permet
le calcul d’écoulements multidimensionnels autour d’un cylindre en résolvant les
équations de Navier Stokes dans la cas incompressible.
Dans le domaine nucléaire, Donea et al. [43] et Kennedy et Belytschko [75]
adaptent la formulation ALE aux problèmes couplés. La structure est traitée
 81

de manière lagrangienne tandis que le fluide est traité de façon eulérienne dans
l’écoulement et de façon lagrangienne à l’interface avec la structure.
Nomura [94] et Nomura et Hughes [95] simulent des écoulements dans lesquels
un solide, considéré comme indéformable, est en mouvement. D’autres, dans le
même esprit, considèrent des structures mobiles avec des modèles fluides plus ou
moins complexes (modèle turbulent par exemple). Boschitsch et Quackenbush [23]
utilisent la méthode des volumes finis pour le fluide alors que Blom et Leyland [21]
utilisent des maillages éléments finis. Undreiner et Frei [125] présentent les résultats
d’une étude couplée fluide-structure portant sur la simulation du comportement
dynamique de la zone fluide entre un flasque de roue et un corps de pompe.

Figure 5.1: Evolution du maillage suivant la formulation utilisée.

La figure 5.1 illustre les difficultés associées aux formulations eulérienne et la-
grangienne ainsi que l’intérêt des formulations de type ALE.

[Link] Mise en œuvre des formulations ALE

Ce paragraphe présente les étapes principales de la mise en œuvre des formu-

lations ALE. Pour plus de détails, il convient de se reporter, par exemple, aux
travaux de Ponthot [104] et de Souli [121] et [120].

a) Relation fondamentale
 82

Pour la formulation ALE, trois configurations sont distinguées :

- La configuration initiale C0 où chaque particule est associée à un nœud.
- La configuration Ct . C’est la configuration après évolution du domaine (déforma-
tions, mouvement) pour laquelle l’équilibre est recherché.
- La configuration de référence Cref . C’est dans cette configuration que les
nœuds sont positionnés et associés à des particules.
Pour la formulation lagrangienne, Cref et C0 sont identiques alors que pour la
formulation eulérienne, ce sont Cref et Ct qui sont identiques. En ce qui concerne
la formulation ALE, Cref est choisie de la façon la plus judicieuse possible. Dans ce
cas, lors du processus d’évolution, la particule est transportée par une application
alors que le nœud est transporté par une application ϕ. Dans la configuration Ct ,
le nœud n’est alors plus associé à la même particule que dans C0 .

Configurations Coordonnées
C0 X = hX1 , X2 , X3 it
Ct x = hx1 , x2 , x3 it
Cref χ = hχ1 , χ2 , χ3 it

Lors de la mise en équation, on est amené à calculer la dérivée par rapport au

temps de champs tels que ceux relatifs à la vitesse, à la masse volumique. Le choix
du référentiel conditionne alors la forme des expressions des dérivées.
Soit f une fonction représentant un champ. Lorsque le système de référence
est le système de référence spatial (formulation eulérienne), la fonction est notée
f (x, t) . Si le système de référence est matériel (formulation lagrangienne), la
fonction est notée fe(X, t). Enfin, pour un système de grille (formulation ALE), la
fonction est notée fb(χ, t). Il vient donc :

f (x, t) = fe(X, t) = fb(χ, t) (5.1)

 83

En calculant les dérivées de chaque fonction par rapport au temps, on obtient :

- Dans le cas d’une formulation eulérienne :

df ∂f ∂f ∂xi ∂f
= + = + ∇f.u (5.2)
dt ∂t ∂xi ∂t ∂t

- Dans le cas d’une formulation lagrangienne :

dfe ∂ fe ∂ fe ∂Xi ∂ fe
= + = (5.3)
dt ∂t ∂Xi ∂t ∂t

- Dans le cas d’une formulation ALE :

dfb ∂ fb ∂ fb ∂χi ∂ fb ∂ fb
= + = + ∇fb.(u − w) = + ∇fb.c (5.4)
dt ∂t ∂χi ∂t ∂t ∂t

où :
∇f : Gradient de la fonction f
∇fb : Gradient de la fonction fb
u : Vitesse de la particule par rapport au système de référence spatial
w : Vitesse du maillage par rapport au système de référence spatial
c : Vitesse de la particule par rapport au système de référence de grille (vitesse
convective)
L’expression (ALE) permet de retrouver le cas eulérien (b
u) et le cas lagrangien
(u = w).
En combinant les équations Lag, Eul et ALE, on obtient la relation suivante,
entre les expressions des dérivées dans les différents systèmes de référence :

∂f ∂ fe ∂ fb
+ ∇f.u = = + ∇fb.c (5.5)
∂t ∂t ∂t

b) Principales étapes de la mise en œuvre de la formulation ALE

 84

La mise en équations de type ALE conduit à un système non carré. En effet,

les déplacements de la matière (n inconnues) et les déplacements du maillage (n
inconnues) sont à déterminer avec seulement n équations. Pour résoudre le système,
il est donc nécessaire de réduire le nombre d’inconnues ou d’augmenter le nombre
d’équations.
L’augmentation du nombre d’équations conduit à établir n équations supplémentaires
pour déterminer le déplacement des nœuds du maillage. Ces équations supplémentaires
peuvent être obtenues en connaissant à priori la vitesse des nœuds du maillage. Par
exemple, cette technique est quelquefois employée en mécanique des fluides, dans le
cas où la vitesse du maillage est donnée par des équations de type lignes de courant
sous forme discrétisées ([70] par exemple). Hu et Liu [65] étendent cette technique
à la mécanique des solides.
La réduction du nombre d’inconnues est retenue par de nombreux auteurs.
Cette réduction est réalisée en partitionnant l’opérateur eulérien-lagrangien en :

• Une phase lagrangienne,

• Une phase eulérienne.

La phase lagrangienne met en jeu un calcul classique où le système d’équations

ne comporte pas de termes convectifs (système carré n x n). Cette phase permet
de résoudre les équations régissant le comportement du domaine étudié tout en
respectant les conditions aux limites. Le maillage utilisé est matériel. La phase
eulérienne est ensuite basée sur une opération de remaillage, le maillage est alors
arbitraire, puis de transfert des données du maillage matériel au maillage arbitraire.
Les techniques classiques de remaillage, appliquées régulièrement au cours d’une
simulation, sont basées sur des critères géométriques, éviter les distorsions par
exemple, ou sur des critères de qualité de la solution.
 85

L’étape de transfert des données est nécessaire après une phase de remaillage.
Les données utiles sont transférées du maillage matériel au maillage arbitraire.
Pour cela, plusieurs méthodes sont également disponibles dans la littérature. Dans
le cadre de la méthode des éléments finis, la valeur des champs connus aux nœuds
du maillage (déplacements, vitesses, forces externes) est aisément calculable en
tout point de l’espace en utilisant les fonctions de forme des éléments. Pour les
champs connus uniquement aux points de Gauss, le lissage de ces champs est une
étape préliminaire au transfert de données. Martinet [83] propose des méthodes
globale et locale des moindres carrés. En fonction des valeurs du champ connu aux
points de Gauss et des fonctions de forme, il est possible d’obtenir les valeurs aux
nœuds du maillage. Une fois les champs à transférer continus, le transfert effectif
est réalisé. La méthode globale des moindres carrés est ici aussi souvent utilisée.
L’idée est de minimiser l’écart entre les deux champs (celui du maillage matériel
f M et celui du maillage arbitraire f A ) :

Z

2
err = fM − fA d∆ (5.6)
∆A

où ∆A représente le nouveau maillage couvrant le domaine ∆.

Le fait que la fréquence de remaillage soit élevée et que le remaillage soit à
topologie constante permet de minimiser l’étape de transfert. En effet, il est pos-
sible de considérer qu’à un point du maillage matériel (nœud ou point de Gauss)
ne correspond qu’un seul point du maillage arbitraire. De plus, ces deux points
restent à des distances assez faibles et ceci permet d’utiliser un développement au
premier ordre. Pour un champ f donné, on a alors :

∂f
f A = f A + dα (5.7)
∂xα
 86

∂f
où dα et ∂xα
sont respectivement les composantes du vecteur joignant les points
des deux maillages et les composantes du gradient du champ f ; f A et f A sont les
valeurs du champ f sur le maillage arbitraire et sur le maillage matériel.

5.2 Procédure du couplage

L’action de couplage fluide/structure doit nécessairement commencer par une

phase de mise en place d’algorithme de couplage entre trois opérateurs principaux :
un solveur de l’écoulement du fluide, un solveur pour la structure et une procédure
de mise à jour du maillage fluide.

5.2.1 Conditions de couplage

Le couplage entre le fluide et la structure intervient au niveau de l’interface entre

les deux domaines. Les conditions du couplage sont donc imposées sur la frontière
Γc . Soit Γc la frontière d’interface entre le fluide et la structure tel que Γc =
T
∂Ωf ∂Ωs . Les conditions de couplage consistent en une opération de transfert de
champs entre le fluide et la structure, et concerne le transfert des contraintes de
pression au solide, qui en se déformant transfert les champs de déplacements au
fluide.
Ces conditions doivent vérifier :

• la continuité des vitesses (condition cinématique) : vf = vs

• la continuité des efforts σf .nΓc = σs .nΓc

Les étapes de transfert des données se déroulent de la manière suivante (Figure

5.2) :
1. le code fluide transmet les pressions, nodales de l’interface, au code structure.
2. le code structure calcule la solution et remet à jour les coordonnées des
noeuds
 87

3. la structure transmet les nouvelles positions des noeuds de l’interface

4. le code fluide mis à jour le maillage, en tenant compte des nouvelles positions
des noeuds, et calcule la solution (pressions), telle que :

X = X0 + ∆X, Y = Y0 + ∆Y (5.8)

X0 et Y0 désignent respectivement les positions initiales suivant x et y ; ∆X et ∆Y

désignent l’incréments des déplacement suivant x et y.
5. retour à l’étape 1.

5.2.2 Algorithme de couplage

La procédure de couplage consiste en la mise à jour régulière de transfert de

champs entre le code fluide et le code structure. Cette remise à jour des variables
entre codes est faite par l’intermédiaire d’un coupleur en fortran mise en place dans
le cadre de cette thèse. L’algorithme de ce coupleur est présenté sur les figures (5.2
et 5.3).

Figure 5.2: Transfert des données entre les deux codes.

 88

Figure 5.3: Algorithme d’interaction fluide-structure.

5.3 Résultats de la procédure de couplage

La procédure de couplage utilisée est basée sur un algorithme explicite. A

chaque itération tn , après convergence du code fluide, on obtient la pression cor-
respondante Pn , qui sera transmise au code structure. Pn sera considérée comme
une sollicitation extérieure aux frontières de la structure.
Le code de calcul des structures, développé à l’UTC de compiègne, est station-
naire et non-linéaire. De ce fait, la convergence ne concerne que la matrice de
rigidité KT , en utilisant la méthode de Newton-Raphson, et non la matrice masse.
Après convergence du code structure, on extrait les nouvelles positions xn et yn ,
puis on actualise la géométrie pour un calcul du code fluide.
 89

5.3.1 Cas test d’un canal plan en présence d’un obstacle vertical

Ce cas test concerne l’écoulement d’un fluide stationnaire (η = 1), visqueux et

incompressible dans un canal de dimensions (L × H = 30m × 3m), à l’intérieur
duquel est placée verticalement une plaque flexible de dimension (l × h = 1m ×
0.1m), (Figure 5.4). La plaque, considérée, comme obstacle déformable, est situéé
à 13m de l’entrée du canal. Les conditions aux limites sont : (u = 10m/s, v = 0)
à l’entrée; (p = 1bar) à la sortie et (u = 0, v = 0) sur les autres parois.

Figure 5.4: Domaine et conditions aux limites.

Figure 5.5: Maillage : 793 noeuds et 1400 éléments.

Figure 5.6: Zoom sur l’interface fluide-structure.

 90

Figure 5.7: Champ de pression le long de la conduite.

Figure 5.8: Déformation de la structure, E = 2.104 M P a.

Le champ de pression est représentée dans (Figure 5.7), la déformation de la

structure dans (Figure 5.8), les lignes de courant avec déformation de la structure
dans (Figure 5.9) et la déformation de la structure pour un module de young
E = 104 M P a.
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.
C’est un couplage dit faible avec des effets rétroactifs du fluide sur la structure,
qui, en se déformant agit sur la nature de l’écoulement du fluide. Les figure (5.10)
montrent que la déformation de la structure dépend de sa rigidité (Module de
Young E) et de la nature de l’écoulement qui exerce sur elle des forces de pression.
 91

Figure 5.9: Les lignes de courant, Re = 1.

Figure 5.10: Module de Young= 104 M P a.

Toutefois, ces forces de pression doivent être calculées judicieusement à partir

du modèle fluide . Ainsi, le profil de l’écoulement a été vérifié en plusieurs point du
canal. On constate qu’à la sortie du canal, ce profil est parfaitement parabolique,
tandis qu’au niveau de l’obstacle une dissymétrie du profil est constatée.
 92

5.3.2 Écoulement dans un canal plan en présence de deux obstacles

Dans ce cas test, on considère l’écoulement fluide (l’eau de mer, η = 1, 9.10−3 )

dans un canal où est immergée deux structures rigides, considérées comme des
obstacles à l’écoulement Figure 5.11. Ces deux obstacles sont fixés sur les parois
intérieurs de ce canal .
Le canal est de longueur L = 30m et de largeur H = 3m, et les deux obstacles
sont identiques, avec une hauteur l = 1m et une épaisseur h = 0.1m .

Figure 5.11: Domaine de calcul et conditions aux limites.

La figure 5.12 représente l’évolution de la pression le long du canal. Nous

avons constatés une depression importante au niveau du deuxième obstacle, ce qui
explique la différence de déplacements des deux obstacle.
 93

Figure 5.12: Pression le long du canal : rouge à travers le premier obstacle et jaune
à travers le deuxième obstacle.

Les figures (5.13, 5.14, 5.15 et 5.16) montrent les déformations des obstacles et
la pression dans le canal en fonctions du module de Young.
En effet des figures (5.13 et 5.15) on constate deux régimes de baisse de pression
dans le canal. Une première depression après le premier obstacle et une deuxième
 94

depression plus importante après le second obstacle. Ceci quelque soit la rigidité
de la structure (quelques soit la valeur du module de Young).

Figure 5.13: Champs des pressions, Modude de Young E = 4.104 M P a.

Figure 5.14: Déformation des structures avec un modude de Young E = 4.104 M P a.

Figure 5.15: Champs des pressions, Modude de Young E = 104 M P a.

 95

Figure 5.16: Déformation des structures avec un modude de Young E = 104 M P a.

Sur la figure (5.17) sont représentées les lignes de courant de l’écoulement,

en fonction de la déformation des deux structures, qui correspond au module de
Young E = 4.104 M P a. De cette figure on constate la présence de structures
tourbillonaires en aval de chaque structure (Figure 5.17).
La taille de ces tourbillons est proportionnelle à la valuer de la viscosité du
fluide et de la longueur de l’obstacle ainsi que de sa déformation.

Figure 5.17: Lignes de courant, Modude de Young E = 4.104 M P a et Re  1.

 96

5.3.3 Conduite en forme de T, avec une paroi flexible

On considère l’écoulement d’un fluide incompressible dans une conduite en

forme de T, composée d’une part d’un tronçon de canal vertical à parois rigides et
d’autre part d’un canal horizontal dont une partie de la paroi supérieure est flexi-
ble. La géométrie et les conditions aux limites sont données dans la figure 5.18. La
partie flexible est fixée à ces extrémités à la paroi rigide. Les déformations de la
partie flexible en fonction du module de Young sont données sur les figures (5.23,
5.24 et 5.25).

Figure 5.18: Domaine et conditions aux limites.

Les champs de pression et des vitesses dans la conduite, sont présentés sur les
figures (5.19, 5.20 et 5.21), pour la valeur du module de Young E = 104 M P a.

Figure 5.19: Champ de pression (Pa), ( E = 104 M P a).

 97

Figure 5.20: Vitesse u(m/s) de l’écoulement (E = 104 M P a).

Figure 5.21: Vitesse v(m/s) de l’écoulement ( E = 104 M P a).

Figure 5.22: Lignes de courants ( E = 104 M P a).

De ces dernières figures on constate que l’écoulement vertical, issu de la branche

verticale, perturbe l’écoulement (vitesse u) dans la branche horizontale. Cet écoulement
 98

vertical a pour effet d’amplifier la vitesse v au niveau de l’intersection des deux

canaux, ce qui a pour conséquences de provoquer la déformation de la paroi supérieure,
notamment au niveau de la partie flexible (5.23, 5.24 et 5.25).
En aval de l’intersection entre les deux branches, l’écoulement devient uniforme
(Figure 5.22).

Figure 5.23: Déformation de la paroi ( E = 104 M P a).

Figure 5.24: Déformation de la paroi ( E = 75.103 M P a).

Figure 5.25: Déformation du paroi ( E = 50.103 M P a).

Les vecteurs déplacements correspondants à ces trois cas de figures sont représentés
sur les figures 5.26, 5.26 et 5.26.
 99

Figure 5.26: Vecteurs déplacements (E = 104 M P a).

Figure 5.27: Vecteurs déplacements (E = 75.103 M P a).

Figure 5.28: Vecteurs déplacements (E = 50.103 M P a).

Pour se rendre compte de l’influence de ce courant vertical dans cette conduite,

on a effectué des coupes longitudinale suivant le canal horizontal (Figure 5.29), afin
de mettre en évidence l’ordre de grandeur des vitesses au niveau de chaque coupe.
 100

Figure 5.29: Positions des coupes longitudinales.

Figure 5.30: Evolution de la vitesse u(x), en fonction de la position de la coupe.

Figure 5.31: Evolution de la vitesse v(x), en fonction de la position de la coupe.

Des figures (5.30 et 5.31), on constate que le maximum des vitesses se trouve
au niveau de l’intersection entre les deux canaux, avec une perturbation notable
de la vitesse v au niveau immédiat de la jonction. Cette perturbation, amplifie
cette vitesse vertical v de l’écoulement, ce qui provoque la déformation de la paroi
flexible (Figure 5.23 et Figure 5.25).
 102

CONCLUSION

Le modèle d’interaction fluide-structure développé dans cette thèse fait appel

à des méthodes développées à partir des éléments finis. La méthode arbitraire
lagrangienne-eulérienne, les équations de Navier-Stokes et de la dynamique d’un
solide élastique nous ont aidé à établir ce modèle avec les dernières approches issues
des méthodes des éléments finis.
Les géométries du fluide et du solide sont envisagé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 mouvement propre sans
passer par un remaillage systématique.
Le solide est considéré suivant une formulation lagrangienne où l’on donne au
maillage qui lui est associé la possibilité de se déplacer en fonction de sa position
actuelle. Cette approche nous a donc permis de tester de manière satisfaisante le
couplage en étudiant des mouvements de structures immergées dans un canal dans
laquelle s’écoule un fluide.
Nous avons résolu les équations de solide élastique avec succès et nous avons ap-
profondi cette résolution pour de grands déplacements dans le cadre de l’interaction
fluide-structure. Nous avons ainsi pu modéliser les déplacements d’un obstacle à
l’intérieur d’un canal plan et obtenir plusieurs résultats intéressants suivant un
écoulement de Poiseuille constant imposé à l’entrée du canal.
Plusieurs perspectives de développement restent encore envisageables :

• L’automatisation du maillage et du remaillage doivent être développés afin

de pouvoir plus facilement aborder des géométries plus complexes par des
maillages structurés ou non structurés.

• C’est à partir de cette amélioration que pourront être développées les tech-
 niques utilisant par exemple la non conformité géométrique ou d’autres types
d’éléments, les tétraèdres afin de rendre plus souple encore la génération de
maillages et d’aborder plus facilement des géométries plus complexes en 2D
et en 3D.

• Le problème lié aux grands déplacements est que nous pouvons être confrontés
à des mouvements non linéaires. Ceci peut fortement diminuer l’éventail de
simulations réalisables car les théories numériques sont construites sur la
résolution de problèmes linéaires ou linéarisés.

Un travail sur la génération de maillage est un point clé de la résolution

de ce type de problème. La déformation de maillage au niveau de la struc-
ture ne permet que de partiellement aborder les grands déplacements. Une
manière plus générale pour décrire le cas de grands déplacements même avec
une technique de remaillage reste manquante et les efforts mathématiques et
numériques allant dans ce sens n’ont souvent été abordés que pour de petits
déplacements ou en considérant un mouvement de corps rigide connu a priori.
 104

BIBLIOGRAPHIE

[1] Ahmad S., Irons B.M. and Zienkiewicz O.C., Analysis of thick and thin shell
structures by curved finite elements, I.J.N.M.E. , Vol. 2, p 419-451, 1970.

[2] Ait Moudid L., Ouahsine A., Taki et M., Korotenko K., Estimation de la
diffusion effective dans un écoulement instationnaire , Revue européenne de
génie civil, Vol. 9 N◦ 7-8/2005, p 993-1009.

[3] Armaly B. F. , Durst F. , Pereira J. C. F., Schönung B., Experimental and

theoretical investigation of backward-facing step flow. J. Fluid Mech. 127, 473-
496, 1983.

[4] Asuaje M., Méthodologie et optimisation dans la conception et l’analyse des

performances des turbomachines à fluide incompressible, Thèse de doctorat,
Ecole Nationale Supérieure d’Arts et Métiers Centre de Paris, 2003.

[5] Auricchio F. and Lovadina C., On the enhanced strain technique for elasticity
problems. Computers and Structures, 81, 777-787, 2003.

[6] Babuska I., The finite element method with Lagrangian multipliers. Numer.
Math., 20, 179-192, 1973.

[7] Baiocchi C. and Lovadina C., A shell classification by interpolation. Math.

Models Methods Appl. Sci., 12(10), 1359-1380, 2002.

[8] Bathe K.J., Chapelle D. , and Lee P.S., A shell problem ’highly sensitive’ to
thickness changes. Internat. J. Numer. Methods Engrg., 57, 1039-1052, 2003.

[9] Bathe K.J., Finite Element Procedures. Prentice Hall, Englewood Cliffs, 1996.

[10] Batoz J.L. and Dhatt G., Modélisation des Structures par Eléments Finis.
Hermes, 1990.
 105

[11] Batoz J.L. and Dhatt G., Modélisation des structures par éléments finis, vol.
1, éditions Hermès, Paris, 1992.

[12] Batoz J.L. and Dhatt G., Modélisation des structures par éléments finis, Vol.
3 : Coques. Hermès Editeur, Paris 1992.

[13] Barrenechea G. and Valentin F.. An unusual stabilized finite element method
for a generalized Stokes problem. Numer. Math., (92), 653-677, 2002

[14] Bathe K.J., Finite Element Procedures, Prentice-Hall, Englewood Cliffs, 1995

[15] Behr M., Liou J., Shih R., et Tezdwar T., Vorticitystreamfunction formulation
of the unsteady incompressible flow past a cylinder : sensibility of the com-
puted flow field to the location of the outflow boundary, International Joumal
for Numerical Methods in Fluids, 12, 323-342, 1991.

[16] Bendjeddou Z., Méthodologie pour la simulation numérique des vibrations

induites par écoulements dans les faisceaux de tubes, Thèse de doctorat, Uni-
versité des siences et technologies de Lille, 2005.

[17] Benjamin AS. , Denny E ., On the convergence ofnumerical solutions for 2-D
fiows in a cavity at large Re, J. Comput. Phys. 33, 340-358, 1979.

[18] Benabid F., Modélisation des écoulements turbullents à faibles Reynilds par la
méthode des éléments finis. Thèse de doctorat, Institut national des sciences
appliquées de rouen 2002.

[19] Benaouicha M., Contribution au Développement d’Algorithmes de Couplage

en Interaction Fluide-Structure, Thèse de doctorat, Université de La Rochelle,
2007.

[20] Boisse P., Daniel J.L. and Gelin J.C., A C◦ three-node shell for non-linear
structural analysis. Computers and Structures, Vol. 58, N◦ 2, p. 249-261, 1996.
 106

[21] Blom F. and Leyland P., Analysis of fluid-structure interaction on moving

airfoils by means of an improved ALE method. AIAA paper, V. 1770, 11p,
1997.

[22] Bodard N., Minteraction fluide-structure par la méthode des éléments spec-
traux, Thèse de doctorat, École Polytechnique Fédérale de Lausanne, 2006.

[23] Boschitsch A. H. and Quackenbush T. R., High accuracy computation of fluid-

structure interaction in transonic cascades. 31st Aerospace Sciences Meeting
and exhibit, AIAA 93-0485. 14p, 1993.

[24] Brezzi F., On the existence, uniqueness and approximation of saddle-point

problems arising from lagrangian multipliers. R.A.I.R.O., Anal. Numér., 8,
129-151, 1974.

[25] Brezzi F. and Fortin M., Mixed and Hybrid Finite Element Methods. Springer-
Verlag, New York, 1991.

[26] Brezis H., Analyse fonctionnelle, Masson

[27] Brooks A. N. and Hughes J. R., Streamline upwind/Petrov-Galerkin formula-

tions for convection dominated flows with particular emphasis on the incom-
pressible Navier-Stokes equations. Computer methods in applied mechanics
and engineering, vol 32, p 199-259, 1982.

[28] Calladine C.R., Theory of Shell Structures, Cambridge University Press, Cam-
bridge, 1983.

[29] Cescutti J. P., Wey E. and Chenot J. L., Finite element calculation for hot
forging with continuous remeshing. in J.L. Chenot, E. Onate, Modelling of
metal forming processes : Proc. of the Euromech 233 Colloquium, Sophia
Antipolis, France, Dordrecht : Kluwer Academic Publishers, p 207-216, 1988.
 107

[30] Chapelle D., Ferent A., and Bathe K.J., 3d-shell finite elements and their
underlying model. Math. Models Methods Appl. Sci., 14(1), 105-142, 2004.

[31] Chapelle D., Ferent A., and Bathe K.J., 3d-shell finite elements and their
underlying model. Math. Models Methods Appl. Sci., 14(1), 105-142, 2004.

[32] Chapelle D. and Bathe K.J., The Finite Element Analysis of Shells - Funda-
mentals. Springer-Verlag, 2003.

[33] Chinosi C. and Lovadina C., Numerical analysis of some mixed finite element
methods for Reissner-Mindlin plates. Comput. Mech., 16, 36-44, 1995.

[34] Ciarlet P. The Finite Element method for elliptic problems, North Holland
1978

[35] Cruchaga M. A., A study of the backwards-facing step problem using a gen-
eralised streamline formulation. Communications in Numerical Methods in
Engineering, 14, 697-708, 1998.

[36] Conca C., Duran M., J. Planchard, A quadratic eigenvalue problem involving
Stokes equations, Comput. Methods Appl. Mech. Engrg. 100 (3) 295-313, 1992.

[37] Cointe, R. Two-dimensional water-solid impact. ASME [Link] Mech. Arc.

Eng., 111, 1989.

[38] Cointe, R., Geyer, P., King, B., Motin, B. and Tramoni, M. Nonlinear and
linear motions of a rectangular barge in a perfect fluid. In Proceedings 18th
Symposium on Naval Hydrodynamics (pp. 85-99). Ann Arbor, MI, 1990.

[39] Dhatt G., An efficient triangular shell element, A.I.A.A. Journal, Vol. 8, N◦ 11,
p. 2100-2102, 1970.
 108

[40] Dhatt G. , finite element modelling of fluids, Lecture, von Karman Institue,
Mars 1992.

[41] Dhatt G. et Touzot G., Une présentation de la méthode des éléments finis,
Maloine S.A, Editeur, Paris France, 1984.

[42] Dobrovol’skaya, Z. On some problems of similarity flow of fluid with a free

surface. JFM, 36, 1969.

[43] Donea J., Giuliani S. and Halleux J. P., An arbitrary lagrangian-eulerian finite
element method for transient dynamic fluid-structure interactions. Computer
method in applied mechanics and engineering, vol 33, p 689-723, 1982.

[44] Donea J., Huerta A., Ponthot J.-Ph., and Rodriguez-Ferran A. Arbitrary
lagrangian-eulerian methods. Encyclopedia of Computational Mechanics, 1 :
Fundamentals, 2004.

[45] Domissy E., Formulation et évaluation d’éléments finis volumiques modifi2s

pour l’analyse linéaire et non linéaire des coques, Thèse de Doctorat, Université
de Technologie de Compiègne, 1997.

[46] Fafard M., Dhatt G. et Batoz J.L., A new discrete Kirchhoff plate/shell ele-
ment with updated procedures, Computers and Structures, Vol. 31, p 591-606,
1989.

[47] Faltinsen, O. Numerical solution of transient nonlinear free-surface motion

outside or inside moving bodies. In Wehausen, J. and [Link] (Eds.), Pro-
ceedings 2nd Conf. on Num. Ship. Hydro (pp. 347-357). [Link]: Uni-
versity Extension Publications, 1977.

[48] Filippone A., site web : [Link]

 109

[49] Formaggia L., Gerbeau J.F., Nobile F., and Quarteroni A., On the coupling
of 3D and 1D Navier-Stokes equations for flow problems in compliant vessels.
Comp. Meth. Appl. Mech. Engrg., 191(6-7), 561-582, 2001.

[50] Fortin A., Jardak M., Gervais J.J., and Pierre R. . Localization of Hopf bifur-
cations in fluid flow problems. [Link] Fluids, 24(11), 1185-1210,
1997.

[51] Fortin M., Glowinski R. Méthode de lagrangien augmenté. Application à

la résolution numérique des problèmes aux limites. Collection méthodes
mathématiques de l’informatique, Dunod, Paris, 1982.

[52] Fourestey G., Simulation numérique et contrôle optimal d’interaction fluidein-

compressible/structure par une méthode de Lagrange-Galerkin d’ordre 2. Ap-
plication aux ouvrages d’art., Thèse de doctorat, Ecole Nationale des Ponts
et Chaussées, 2005.

[53] Flügge W. , Stresses in Shells, Springer-Verlag, New York, Heidelberg, 2nd

edition, 1973.

[54] Garabedian, P. Oblique water entry of wedge. Comm. Pure and Appl. Math.,
6, 1984.

[55] Gervais J.J., Lemelin D. , Pierre R., Some experiments with stability analysis
of discrete incompressible flows in the lid-driven cavity, Int. J. Numer. Meth.
Fluids 24, 477-492, 1997.

[56] Ghia U., Ghia K., Shin C.T., Higher-Re solutions for incompressible pou using
the Nauier-Stokes equations and a multigrid method, J. Comput. Phys. 48,
387-41 1, 1982.
 110

[57] Gingold, R. and Monaghan, J. Smoothed particle hydrodynamics - theory

and application to non-spherical stars. Royal Astronomical Society, Monthly
Notices, 181, 375-389, 1977.

[58] Glockner S. Aquilon, Notice théorique et numérique. Technical report, 2001.

[59] Godlewski E. and Raviart P.A., Hyperbolic systems of conservation laws,

SMAI, Mathematiques Applications 3/4, Ellipses-Edition Marketing, Paris,
1991..

[60] Grandmont C., Guimet V. et Maday Y.. Existence de solutions d’un problème
de couplage fluide-structure bidimesionnnel instationnaire. Math. Prob. in
Mech., 326(1), 525-530, 1998.

[61] Gustafson K. and Halasi K.,Cavity flow dynamics at higher reynolds number
and higher aspect ratio, J. Comput. Phys. 70, 271-283. 1987.

[62] Hadji S. et Dhatt G. , Asymptotic-Newton method for solving incompressible

flows, Int. [Link]. Methods fluids, 44, 861-878, 1997.

[63] Havu V. and Pitkäranta J., Analysis of a bilinear finite element for shallow
shells. Approximation of inextensional deformations. Math. Comp., 71, 923-
943, 2002.

[64] Hirt C. W., Amsden A. and Cook J. L., An arbitrary lagrangian eulerian
computing method for all flow speeds. Journal of Computational Physics, vol
14, p 227-253, 1972.

[65] Hu Y. K. and Liu W. K., ALE finite element formulation for ring rolling
analysis. Int. journal for numerical methods in engineering, vol 33, p 1217-
1236, 1992.
 111

[66] Hu Y. K. and Liu W. K., ALE hydrodynamic lubrication finite element method
with application to strip rolling. International journal for numerical methods
in engineering, vol 36, p 855-880, 1993.

[67] Hubert G., Modélisation d’écoulement de fluides incompressible par la

méthode des éléments finis. Thèse de doctorat, Université de Compiègne 1984.

[68] Hughes T.J.R. and Brooks A.N.. A multidimensional upwind scheme with no
crosswind diffusion. Finite Element Methods for Convection Dominated Flows,
ed. T.J.R. HUGHES, AMD, 34, 19-35 ASME, New York, 1979.

[69] Hughes T.J.R., The Finite Element Method. Prentice-Hall, 1987.

[70] Hughes J. R., Liu W. K. and Zimmerman T. K., Lagrangian Eulerian finite
element formulation for incompressible viscous flows. Computer methods in
applied mechanics and engineering, vol 29, p 329-349, 1981.

[71] Hughes T.J.R. and Taylor R.L., The linear triangle bending elements, In The
Mathematics of Finite Element and Application IV, MAFELAP 1981, London,
Academic Press, p. 127-142, 1982.

[72] Hunter S. C., Mechanics of continuous Halsted Press, Sheffield, 1976.

[73] Jeng Y. H. and Payne U. J., An Adaptive tvd Limiter, J. Comput. Phys. 118,
229-241, 1995.

[74] Kelly D.W., Nakasawa S., Zienkiewicz O.C., and HeinrichJ.C., A note on up-
winding and anisotropic balancing dissipation in finite element approximations
to convective diffusion problems. International Journal of Numerical Methods
in Engineering, 15, 1705-1711, 1980.
 112

[75] Kennedy J. M. and Belytschko T., Theory and application of a finite ele-
ment method for arbitrary lagrangian eulerian fluids and structures. Nuclear
engineering and design, vol 68, p 129-146, 1981.

[76] Kreiss H.O., Numerical methods for solving time independent problems for
partial differential équations, Sém. de Math. Sup. Univ. Montreal 1978.

[77] KoiterW. T., On the nonlinear theory of thin elastic shells. Proc. Kon. Ned.
Akad. Wetensch., B69, 1-54, 1965.

[78] Longuet-Higgins, M. and Cokelet, E. The deformation of steep surface waves

on water: I. a numerical method of computation. Proc. R. Soc. Lond, 350,
1-26, 1976.

[79] Lucy, L. A numerical approach to the testing of the fission hypothesis. Astro-
nomical Journal, 82, 1013-1024,.1977.

[80] Luo X.Y., Pedley T.J., Numerical simulation of steady flow in a 2-D Collapsible
channel, J. Fluids Struct. 9, 149-197, 1995.

[81] Majda A., Osher S., A systematic approach for correcting nonlinear instabili-
ties, Num. Math.,30, pp.429-452, 1978.

[82] Marcotte J. , Méthodes itératives pour la résolution, par éléments finis, du

problème de Stokes non linéaire, Thèse de doctorat, Université de Montréal,
2000.

[83] Martinet F., Etude et mise en oeuvre d’une formulation ALE Application au
contact lubrifié en mise en forme. Thèse de doctorat : Université Aix-Marseille
II 204 p, 1998.
 113

[84] Mendes P.A. et Branco F.A.. Analysis of fluid-structure interaction by an Ar-

bitrary Lagrangian-Eulerian finite element formulation. Int. J. Numer. Meth.,
30, 897-919, 1999.

[85] Monaghan, J. and Gingold, R. Shock simulation with the particle method sph.
J. Comp. Phys., 52, 374-389, 1983.

[86] Moondra S. and Upadhyay A. and Sushanta K. Mitra, Investigation of Fluid

Structure Interaction of an Elastic Membrane in a Microchannel. Excerpt from
the Proceedings of the COMSOL Users Conference. Bangalore, 2006.

[87] Mindlin R.D., Influence of rotary inertia and shear on flexural motion of
isotropic elastic plates. J. Appl. Mech., 18, 31-38, 1951.

[88] Naceur H., Contribution à l’optimisation de forme de structures minces en

présence de non linéarités géométriques et matérielles. , Thèse de doctorat,
Université de Compiègne, Compiègne, 1998.

[89] Nichols, B. and Hirt, C. Volume of fluid (vof) method for the dynamics of free
boundaries. Journal of Computational Physics, 39(201), 1981.

[90] Hammadi F., Formulation et evaluation d’élements finis à continuite CR de

la geometrie pour l’analyse lineaire et non lineaire des coques = formulation
and evaluation of finite element models with CR continuity of the geometry
for linear and non linear analysis of shells . , Thèse de doctorat, Université de
Compiègne, Compiègne, 1998.

[91] Naghdi P.M., Foundations of elastic shell theory, In Progress in Solid Mechan-
ics, Vol. 4, 1-90. North-Holland, Amsterdam, 1963.

[92] Novozhilov V.V., Thin Shell Theory. Wolters-Noordhoff Publishing, Gronin-

gen, 2nd edition, 1970.
 114

[93] Nomura T. et Hughes T.J.R., An Arbitrary Lagrangian-Eulerian finite element

method for interaction of fluid and a rigid body. Comput. Meth. Appl. Mech.
Engrg., 95, 115-138, 1992.

[94] Nomura T., ALE finite element computations of fluid-structures interaction

problems. Computer methods in applied mechanics and engineering, vol 112,
p 291-308, 1994.

[95] Nomura T. and Hughes T. J. R., An arbitrary lagrangian eulerian finite ele-
ment method for interaction of fluid and a rigid body. Computer methods in
applied mechanics and engineering, vol 95, p 115-138, 1992.

[96] N’dri K. D., Formulation éléments finis espace-temps pour les équations de
Navier-Stokes, Thèse de doctorat, Université de Montréal, 2001.

[97] Paidoussis M.P., Dynamics of flexible slender cylinder in axial flow, Part I,.
Theory,. J. of Fluid Mechanics,. 26, p. 717-736, 1966.

[98] Pian T.H.H. and Sumihara K., Rational approach for assumed stress finite
elements, International Journal of Numerical Methods in Engineering, vol. 20,
pp. 1685-1695, 1984.

[99] Piperno S., Simulation numérique de phénomènes d’interaction fluide-

structure, Thèse, ENPC, jun 1995

[100] Pierson, J. The penetration of fluid surface by a wedge. Stevens Inst. Tech.,
Expt. Towing Tank Report, 387, 1950.

[101] Piperno S., Simulation numérique de phénomènes d’interaction fluide-

structure, Thèse de doctorat, Ecole Nationale des Ponts et Chaussées, 1995.

[102] Piperno S. , Numerical simulation of aeroelastic instabilities of elementary

bridgedecks, technical report, INRIA, RR-3549, 1998.
 115

[103] Pironneau O., On the transport diffusion algorithm and its applications to
the Navier- Stokes equations, Numer. Math. Vol. 38, pp. 309-332, 1982.

[104] Ponthot J. P. and Belytschko T., Arbitrary lagrangian eulerian formulation

for element free Galerkin method. Computer methods in applied mechanics
and engineering, vol 152, p 19-46, 1997.

[105] Pitkäranta J., Leino Y., Ovaskainen O., and Piila J., Shell deformation states
and the finite element method : a benchmark study of cylindrical shells. Com-
put. Methods Appl. Mech. Engrg., 128, 81-121, 1995.

[106] Raviart P.A., Les méthodes d’éléments finis en mécanique des Coll. DER-
EDF Eyrolles 1981.

[107] Raviart P.A.,Thomas J.M., Introduction a l’analyse des équations aux

dérivées partielles, Masson, 1983.

[108] Reissner E., Stress strain relations in the theory of thin elastic shells. J. Math.
Phys., 31, 109-119, 1952.

[109] Reissner E., The effect of transverse shear deformation on the bending of
elastic plates. J. Appl. Mech., 67, A69-A77, 1945.

[110] Renou J.Y., Une méthode eulérienne de calcul numérique de forces fluide-
élastiques . Thèse de. doctorat : Université Paris VI, 1998.

[111] Rifai M.S. and Simo J.C., A class of mixed assumed strain methods and
the method of incompatible modes. Internat. J. Numer. Methods Engrg., 29,
1595-1638, 1990.

[112] Roberts J. E. and Thomas J.M., Mixed and hybrid methods. In Handbook
of numerical analysis, Vol. II, Handb. Numer. Anal., II, pages 523-639. North-
Holland, 1991.
 116

[113] Rosenfeld M. , Kwak D. and Vinokur M.. A Solution Method for The Un-
steady And Incompressible Navier-Stokes Equations In Generalized Coordi-
nate Systems, ALAA Paper, 88-718, 1988.

[114] Rusanov V., On difference schemes of Third order of accuracy for non linear
Hyperbolic Systems J. Comp. Phys., 5, 507-516, 1970.

[115] Rubin K. and Hosla, Navier-Stokes calculations with a coupled strongly irn-
pliczt method : finite difference solutions, Computers and Fluids 4, 163-180,
1981.

[116] Sanchez-Hubert J. and Sanchez-Palencia E., Coques Elastiques Minces - Pro-

priétés Asymptotiques. Masson, Paris, 1997.

[117] Sarrate J., Huerta A., and Donea. J. Arbitrary lagrangian-eulerian formula-
tion for fluid-rigid body interaction. Computer Methods in Applied Mechanics
and Engineering, 190 :3171-3188, 2001.

[118] Scordelis A.C. and Lo K.S., Computer analysis of cylindrical shells, J. Amer.
Concrete Institute, Vol. 61, p. 539-561, 1969.

[119] Simo J.C. ans Reddy B.D., Stability and convergence of a class of enhanced
strain methods. Siam. J. Numer. Anal., 6, 1705-1728, 1995.

[120] Souli,M., Ouahsine, A., Lewin,L., ALE formulation for fluid-structure inter-
action problems. Comput. Methods Appl. Mech Engrg,190,659-675, 2000.

[121] Souli, M. and Zolesio, J.P., Arbitrary Lagrangian-Eulerian and free surface
methods in fluid mechanics, Comput. Methods Appl. Mech. Eng., 191, 2001.

[122] Stanley G., Park K.C. and Hughes T.R.J., Continuum-based resultant shell
elements, In Hughes T.J.R and Hinton E. (Eds), Finite Element Methods for
 117

Plates and Shell Structures, Vol 1, Pineridge Press, Swansea, UK., p. 1-45,
1986.

[123] Thomée V. , Stability of difference schemes in the maximum norm.

[Link] équations,1,273-292, 1965.

[124] Timoshenko S. and Woinowsky-Krieger S., Theory of Plates and Shells.

McGraw-Hill, New York, 1959.

[125] Undreiner S. and Frei A., Simulation du comportement dynamique de la zone

lame fluide entre le flasque de la roue et le corps de pompe. Colloque de la
société hydrotechnique de France, Chatou, 19-20 Novembre, p 185-191, 1997.

[126] Valid R., The Nonlinear Theory of Shells through Variational Principles. John
Wiley and Sons, Chichester, 1995.

[127] Vinje, T. and Brevig, P. Nonlinear ship motions. In Proceedings 3rd Inter-
national Symp. Num. Ship Hydro. Paris, FRANCE, 1981.

[128] Von Karman, T. The impact on seaplane floats during landing. NACA TN,
321, 1929.

[129] Wagner, H. über stoss-und gleitvorgänge an der oberfläche von flüssigkeiten.

ZAMM, 12, 1932.

[130] Walsh C., Flutter in One-Dimensional Collapsible Tubes, J. Fluids Struct.,

9, pp. 393-408, 1995.

[131] Watanabe, Y. On the properties of hydrodynamic impact pressure by

matched asymptotic expansion technique. T. West-Japan Soc. Nav. Arch.,
32, 1966.
 118

[132] Wei G. and Kirby J.T.,Time-dependent numerical code for extended Boussi-
nesq équations, Journal of Waterway, Port, Coastal and Ocean Engineering
121, 251-261, 1995.

[133] Wempner G.A., Oden J.T. et Kross D.A., Finite element analysis of thin
shells, Engng. Mech. Div., Proc. ASCE, Vol. 94, EM6, p. 1273-1294, 1968.

[134] Wunderlich W., On a consistent shell theory in mixed tensor formulation. In

Proc. 3rd IUTAM Symposium on Shell Theory : Theory of Shells, Amsterdam.
North Holland, 1980.

[135] Yang H.T.Y, Saigal et Liaw., Advances of thin shell finite elements and some
applications - version 1, Computers and Structures, Vol. 35, p. 481-504, 1990.

[136] Youngs, D. Time dependent multi-material flow with large fluid distortion.
Numerical Methods for Fluid Dynamics (pp. 273-285), 1982.

[137] Zhao, R. and Faltinsen, O. Slamming loads on high-speed vessels. In 19th

ONR Conference. Korea, 1992.

[138] Zhao, R. and Faltinsen, O. Water entry of two-dimensional bodies. JFM, 246,
1993.

[139] Zhong, Z. H. Finite Element Procedures for Contact-Impact Problems. Ox-

ford: Oxford Science Publications, 1993.

[140] Zienkiewicz O.C et Taylor R.L., La méthode des éléments finis. Formulation
de base et problèmes linéaires, Edition française AFNOR, Paris 1991.