UNIVERSITÉ DU QUÉBEC À CHICOUTIMI

THÈSE PRÉSENTÉE À
L'UNIVERSITÉ DU QUÉBEC À CHICOUTIMI
COMME EXIGENCE PARTIELLE
AU DOCTORAT EN INGÉNIERIE

PAR
SÉBASTIEN PERRON

RÉSOLUTION NUMÉRIQUE
D'ÉCOULEMENTS 3 DIMENSIONS
AVEC UNE NOUVELLE MÉTHODE DE VOLUMES FINIS
POUR MAILLAGES NON STRUCTURÉS

15 Novembre 2001
 UIUQAC
bibliothèque
Paul-Emile-Bouletj

Mise en garde/Advice

Afin de rendre accessible au plus Motivated by a desire to make the

grand nombre le résultat des results of its graduate students'
travaux de recherche menés par ses research accessible to all, and in
étudiants gradués et dans l'esprit des accordance with the rules
règles qui régissent le dépôt et la governing the acceptation and
diffusion des mémoires et thèses diffusion of dissertations and
produits dans cette Institution, theses in this Institution, the
l'Université du Québec à Université du Québec à
Chicoutimi (UQAC) est fière de Chicoutimi (UQAC) is proud to
rendre accessible une version make a complete version of this
complète et gratuite de cette œuvre. work available at no cost to the
reader.

L'auteur conserve néanmoins la The author retains ownership of the

propriété du droit d'auteur qui copyright of this dissertation or
protège ce mémoire ou cette thèse. thesis. Neither the dissertation or
Ni le mémoire ou la thèse ni des thesis, nor substantial extracts from
extraits substantiels de ceux-ci ne it, may be printed or otherwise
peuvent être imprimés ou autrement reproduced without the author's
reproduits sans son autorisation. permission.
 Résumé

L'objectif de cette thèse de doctorat est de développer un nouveau schéma numérique

de résolution des équations 3D de Navier-Stokes. Les écoulements considérés sont incom-
pressibles et les propriétés physiques telle la viscosité peuvent être variables. De plus, des
quantités scalaires telle l'enthalpie ou une fraction volumique peuvent être transportées par
l'écoulement. Ce schéma est basé sur une discrétisation spatiale de type volumes finis clas-
sique (une seule inconnue par volume de contrôle) qui permet l'utilisation de maillages non
structurés. Les équations du modèle mathématique sont découplées et une méthode de pro-
jection de type pas de temps fractionnaire est utilisée pour calculer un champ de vitesse
satisfaisant l'équation de conservation de la masse. La pression et la vélocité sont couplées à
l'aide d'une ré-interpolation de la vélocité sur les faces des volumes de contrôle, l'utilisation
d'un maillage décalé ("staggered grid") n'est donc pas nécessaire.
On présente une revue détaillée du modèle mathématique décrivant les écoulements
thermiques et le modèle de turbulence k — e de Launder et Spalding. De plus, on fait un rappel
de quelques algorithmes de résolution numérique des systèmes d'équations linéaires et des
concepts fondamentaux qui leur sont associés.
On discute d'une librairie d'objets et de programmes utilitaires qui ont été développés
pour encourager l'utilisation et le développement du schéma par d'autres chercheurs.
Des résultats numériques obtenus avec le schéma sont présentés : écoulement entre
deux plans parallèles, écoulement engendré par le déplacement d'une paroi, écoulement de
Boussinesq dans une cavité carrée, écoulement turbulent au-dessus d'une marche, écoule-
ments autour d'un disque et d'un cylindre et écoulement thermique dans une conduite cylin-
drique.
 Remerciements

Je tiens à remercier M. Sylvain Boivin, mon directeur, pour m'avoir si bien dirigé.
De plus, sa patience et son aide lors de mes démarches pour l'obtention d'un emploi ont été
grandement appréciées.
Je remercie également M. Guy Simard du G.R.I.P.S pour ses précieux conseils et de
m'avoir donné la chance d'établir mes premiers contacts avec l'industrie.
Je tiens aussi à souligner l'aide que m'a apporté M. Jean-Marc Hérard d'Électricité
de France dans mon travail de recherche.
Finalement, je remercie M. Vinko Potocnik du CRDA pour m'avoir fait confiance et
ainsi permis d'acquérir de l'expérience avec l'industrie privée.
 Table des matières

Remerciements i

Résumé ii

Nomenclature ix

Introduction 1

1 Modèle mathématique, écoulement laminaire 11

1.1 Notions fondamentales 11
1.2 Conservation de la masse 14
1.3 Conservation de la quantité de mouvement 16
1.4 Équation de conservation de l'énergie 22
1.5 Modèle complet 25
1.6 Nombres sans dimension 26

2 Modèle mathématique, écoulement turbulent 31

2.1 Modélisation de la turbulence 32

3 Définitions et résultats essentiels sur les méthodes de volumes finis classiques 41

3.1 Discrétisation temporelle 44
3.2 Discrétisation spatiale 45
3.3 Résolution de l'équation de conservation de la masse, méthodes de projection 56
 IV

4 Une nouvelle méthode de volumes finis pour des maillages non structurés 62
4.1 Discrétisation temporelle 63
4.2 Discrétisation spatiale 64
4.3 Approximation discrète sur les interfaces 68
4.4 Calcul des flux d'interface 69
4.5 Opérateur d'approximation du gradient Grad((j)) 70
4.6 Équations discrètes, opérateur de convection-diffusion CD (• ••) 72
4.7 Projection 73
4.8 Discrétisation des lois de paroi pour les écoulements turbulents 76
4.9 Algorithme de résolution 80
4.10 Remarques et résultats théoriques 80

5 Résolution numérique des systèmes d'équations linéaires 83

5.1 Caractérisation de la matrice des coefficients et détermination des valeurs
propres 84
5.2 Méthodes directes de résolution des équations linéaires 87
5.3 Méthodes itératives de résolution numérique des systèmes d'équations linéaires 92
5.4 Conditionnement 97
5.5 Résultats comparatifs et choix de l'algorithme de résolution 102

6 Mise-en-oeuvre 105
6.1 Langage de programmation et programmation orientée objet 105
6.2 Description des objets et de la librairie 106
6.3 Programmes pour le pré et le post-traitement 116
6.4 Installation et disponibilité 117

7 Résultats numériques 118

7.1 Écoulement entre deux plans parallèles 119
7.2 Écoulement engendré par le déplacement d'une paroi 122
7.3 Convection naturelle dans une cavité carrée 127
7.4 Écoulement turbulent au-dessus d'une marche 131
7.5 Écoulement autour d'un cylindre 136
7.6 Écoulement thermique 3D dans un cylindre 146
8 Extension de la méthode aux écoulements avec densité variable 153
8.1 Discrétisation temporelle 154
8.2 Calcul de la densité d'interface 155
8.3 Projection 155
8.4 Choix de l'algorithme de résolution des équations linéaires 157
8.5 Respect du principe du maximum 158

Conclusion 159
 Table des figures

3.1 Volume de contrôle pour un maillage Voronoï 48

3.2 Extrapolation du gradient de pression 58
3.3 Illustration du problème de l'échiquier pour la vitesse 59
3.4 Maillage décalé structuré 60
3.5 Élément d'un maillage non structuré "décalé" 60

4.1 Position des I.D.O.s 64

4.2 Configurations des I.D.O.s 65
4.3 Illustration du passage des éléments géométriques aux volumes de contrôle . 67
4.4 Problème d'interpolation à l'interface 68
4.5 Algorithme de calcul de la vitesse de frottement 78
4.6 Volume de contrôle chevauchant deux milieux hétérogènes 82

5.1 Algorithme de Gauss 87

5.2 Factorisation LU, décomposition de Crout 89
5.3 Factorisation de Choleski 91
5.4 Algorithme du Gradient conjugué 94
5.5 Algorithme Orthomin2 95
5.6 Algorithme d'Arnoldi 96
5.7 Algorithme GMRES 97

6.1 Format AMDBA d'un maillage 114

6.2 Diagramme de relation entre les classes 115

7.1 Géométrie, écoulement entre deux plans parallèles 118

 Vil

7.2 Convergence, écoulement entre deux plans parallèles 119

7.3 Vitesse horizontale v, écoulement entre deux plans parallèles 121
7.4 Pression, écoulement entre deux plans parallèles 122
7.5 Coupes, écoulement engendré par le déplacement d'une paroi 124
7.6 Vélocité, écoulement engendré par le déplacement d'une paroi 125
7.7 Pression, écoulement engendré par le déplacement d'une paroi 126
7.8 Convergence, écoulement engendré par le déplacement d'une paroi 126
7.9 Convergence, convection naturelle dans une cavité carrée 129
7.10 Vélocité, convection naturelle dans une cavité carrée 130
7.11 Pression, convection naturelle dans une cavité carrée 130
7.12 Température, convection naturelle dans une cavité carrée 131
7.13 Géométrie, écoulement turbulent au-dessus d'une marche 131
7.14 Profils de vitesse, écoulement au-dessus d'une marche 134
7.15 Convergence, écoulement au-dessus d'une marche 135
7.16 Géométrie, écoulement 2D autour d'un cylindre 136
7.17 Coupes du maillage, écoulement 2D permanent autour d'un cylindre 138
7.18 Convergence, écoulement permanent 2D autour d'un cylindre 140
7.19 Géométrie, écoulement 3D autour d'un cylindre 142
7.20 Convergence, écoulement 3D permanent autour d'un cylindre 145
7.21 Vitesse transversale, écoulement thermique 3D dans un cylindre 149
7.22 Profil de vitesse à la sortie, écoulement thermique 3D dans un cylindre . . . . 149
7.23 Profil de température, écoulement thermique 3D dans un cylindre 150
7.24 Gradient de pression, écoulement thermique 3D dans un cylindre 151
7.25 Convergence, écoulement thermique 3D dans un cylindre 152
 Liste des tableaux

3.1 Schémas convectifs 53

5.1 Temps de résolution pour différents algorithmes de résolution, problème de

l'écoulement 3D thermique dans une conduite 103

7.1 Résultats comparatifs, écoulement engendré par le déplacement d'une paroi . 124
7.2 Résultats comparatifs, convection naturelle dans une cavité carrée 128
7.3 Coefficients, écoulement 2D permanent autour d'un cylindre 140
7.4 Coefficients, écoulement 2D périodique autour d'un cylindre 142
7.5 Coefficients et statistiques, écoulement 3D permanent autour d'un cylindre . . 146
 Nomenclature

Caractères usuels
dS : Élément de surface 2

dV : Élément de volume m3
D : Longueur de référence m
cp : Chaleur massique à pression constante J-kg-1
cv : Chaleur massique à volume constant J-kg-1
e • Énergie interne par unité de masse J• kg'1
£K Ensemble des interfaces du volume K
c-ext : Ensemble des interfaces sur le bord du domaine
£int Ensemble des interfaces à l'intérieur du domaine
F Vecteur des flux par unité de surface
F Flux total
I Vecteur des forces de volume N -m-3
h Enthalpie massique J• kg'1
hc Coefficient de transfert convectif W-m-' •K-1
hr Coefficient de transfert radiatif W -m-' •H'1
k Énergie cinétique turbulente m2 • s~2
m(..-) Mesure d'un hyperplan
M Nombre de Mach
n Vecteur normal unitaire.
P Pression N-m-2
P Pression cinétique m2 • s2
 Caractères usuels
-2
ç" : Vecteur des flux de chaleur par unité de surface W-m
r : Constante des gaz parfaits par unité de masse J• kg~l K -i
s : Fonction de production par unité de volume
t : Coordonnée de temps s
t(n) : Vecteur des contraintes normales N-m~2
T : Température K
u* : Vitesse de frottement à la paroi m • s~l
u+ : Vitesse de frottement sans dimension
uT : Vitesse tangente à la paroi m
v : Norme du vecteur vitesse m
Vi : Composante du vecteur vitesse m
v : Vecteur vitesse (vélocité) m
x, y, z : Coordonnées spatiales m
x_ : Vecteur position m
XK '• Position d'un point associé à un volume ou un élément K
y : Distance fictive entre la paroi réelle et la frontière de calcul m
y+ : Distance fictive sans dimension entre la paroi réelle
et la frontière de calcul (valeur sans dimension)
 XI

Caractères grecs
a Paramètre de diffusion
/3 Coefficient d'expansion volumique K-1
Delta de Kronecker
e Émissivité thermique
e Taux de dissipation de l'énergie turbulente m • S ô
À Conductivité thermique W - m"1 • K-1
Viscosité dynamique moléculaire kg- m-1- s'1
Viscosité dynamique turbulente kg- m'1 • s'1
v Viscosité cinématique moléculaire 2
s~l
1
vt Viscosité cinématique turbulente m2 • s"
Volume de contrôle arbitraire m3
Volume de contrôle matériel dépendant du temps m3
4> Variable scalaire
p Masse volumique kg- m~3
Constante de Boltzman 5.67 x 10""8 • m'2 • K~4
OK,L Interface séparant les volumes (éléments) K etL
g_ Tenseur des contraintes visqueuses N • m -2
Transmittivité de l'interface aK r
 Xll

Principales notations
Tenseur d'ordre 1 (ou vecteur)
Tenseur d'ordre 2 (ou matrice)
ð Valeur moyenne de la variable 4>
Fluctuation de la variable <j>
® Produit dichotomique entre deux vecteurs
V• Divergence d'un vecteur
Gradient d'une fonction scalaire
 Introduction

Cette thèse de doctorat porte sur le développement d'une nouvelle méthode numérique
de résolution d'écoulements fluides. Les écoulements considérés sont incompressibles et les
propriétés physiques telles la viscosité ou la densité peuvent être variables. Ces problèmes
sont caractérisés par la résolution d'un système d'équations de transport du type

[vj)(t>-aV(f>}-ndS = fs(---)dV,
n ot Jan Jn
où:

• 0. est un volume de contrôle arbitraire, dCl est sa surface ;

• n est un vecteur unitaire normal à la surface d£l ;

• p est la masse volumique qui n'est pas nécessairement constante ;

• <f> est une variable scalaire qui dépend de la position x_ et du temps t ;

• / [up^] ' 21 dS est le flux convectif, y_ étant la vélocité ;

Jan

I — a V0-n dS est le flux diffusif, a étant un paramètre de diffusion telle la viscosité ;

JdQ

s (• • • ) dV est un terme source, s étant une fonction de production ou de destruction

En général, les problèmes d'écoulements fluides sont très complexes et les équations du
modèle mathématique doivent être résolues avec des schémas numériques. Ces schémas
sont basés sur des méthodes de discrétisation des équations de transport qui peuvent être
regroupées en quatre familles :

1. les méthodes de différences finies ;

2. les méthodes d'éléments finis ;
3. les méthodes de volumes finis de type éléments finis (ou volumes finis de type Ga-
lerkin);

4. les méthodes de volumes finis classiques.

La discrétisation par la méthode des différences finies (voir [60], [34] et [38]) est basée
sur une approximation des opérateurs de dérivées. Ces opérateurs sont approximés à l'aide de
développements par série de Taylor ou bien par la dérivation d'un polynôme d'approximation.
En général, les schémas basés sur une telle discrétisation sont peu flexibles et il est difficile
de tenir compte de propriétés physiques (telle la conductivite) qui peuvent être discontinues
lorsque le milieu est hétérogène. De plus, le maillage est généralement structuré, orthogonal
et le pas d'espace constant pour une direction donnée : la résolution d'écoulements fluides
dans des geometries complexes est alors problématique.
Les schémas de la famille des éléments finis (voir [15],[20],[34],[43] et [61]) sont
basés sur une formulation variationnelle du problème continu et du problème discret. Cette
formulation est obtenue en multipliant les équations aux dérivées partielles par une "fonction
test". Les inconnues sont ensuite approximées par une combinaison linéaire de "fonctions
formes" (ou fonctions de base). Les équations discrètes sont obtenues en intégrant l'équation
résultante sur tout le domaine. Les schémas de cette famille permettent l'utilisation de mail-
lages non structurés et sont particulièrement bien adaptés aux geometries complexes. Les flux
sont conservés globalement. Localement, la conservation des flux n'est toutefois pas assurée.
Les schémas de cette famille sont particulièrement bien convenir à l'adaptation de maillages,
l'adaptation du maillage permettant dans de nombreux cas une amélioration significative de
la qualité des solutions.
Les méthodes de volumes finis de type éléments finis sont elles aussi très bien adap-
tées aux geometries complexes (voir [36],[6] et [28]) : le maillage étant composé d'éléments
géométriques pouvant avoir des formes différentes (quadrangles, triangles, tétraèdres, hexaè-
dres, ...). Sur les éléments, les inconnues sont approximées avec des fonctions polynômiales.
L'équation discrète pour une inconnue est obtenue en intégrant les équations de transport
sur des parties d'éléments associés à l'inconnue considérée. Tout comme pour les méthodes
d'éléments finis, ces méthodes de volumes finis sont compatibles avec l'adaptation de mail-
lages.
En ce qui concerne les méthodes de volumes finis classiques (voir [21],[34],[36] et
[59]), le domaine peut être discrétisé avec des éléments convexes divers (quadrangles, tétraè-
dres...). Les inconnues sont approximées avec des fonctions constantes par volume de con-
trôle et une position de référence est associée à chacun des volumes. Les équations discrètes
sont obtenues en intégrant les équations de transport sur les volumes de contrôle, les flux entre
les volumes étant approximés avec des formules aux différences. Dans la pratique, pour les
problèmes de convection-diffusion, les maillages sont essentiellement structurés ou obtenus
avec une technique Voronoï. Pour toutes les méthodes de volumes finis (incluant les vol-
umes finis de type éléments finis), les flux sont localement conservés lorsque les équations de
transport sont écrites sous forme conservative et il est possible de tenir compte de propriétés
physiques discontinues. Les schémas de cette famille ne se sont toutefois pas compatibles
avec l'adaptation de maillages : la qualité de la solution peut être dégradée par retirement
des éléments. Dans ces conditions, compte tenue du faible degré des fonctions d'approxima-
tion, les dérivées secondes ne seront pas être bien évaluées. Une propriété qualitative propre
à ces schémas est le respect du principe du maximum discret lorsqu'il s'applique, et ce, pour
des maillages généraux.
Finalement, on peut mentionner les méthodes de volumes finis où les fonctions d'ap-
proximations ne sont pas nécessairement les mêmes dans tous les termes des équations de
transport.. Pour ces méthodes, les inconnues sont constantes par volume de contrôle dans le
terme transitoire. Pour évaluer les flux, des polynones d'approximation d'un degré supérieur
peuvent être utilisés afin d'obtenir une approximation consistante des flux lorsque le maillage
n'est pas orthogonal (voir [72]). Toutefois, en augmentant l'ordre des fonctions d'approxima-
tion, le principe du maximum discret ne sera pas nécessairement respecté.
Lorsque l'écoulement est incompressible, trois catégories de difficultés sont inhéren-
tes à toutes les méthodes de discrétisation :

1. la détermination de la pression ;
2. le couplage vitesse-pression ;
3. la discrétisation du terme de convection.
Avec des conditions de bord adéquates, les équations de Navier-Stokes et l'équation de la
quantité de mouvement peuvent être résolues directement (sans découplage des équations) et
la matrice représentant le système d'équations est alors singulière. Pour faciliter la résolution,
une dérivée temporelle de la pression peut être ajoutée à l'équation de continuité (méthode de
compressibilite artificielle). Chorin[24], en 1967, a été le précurseur de cette approche qui est
mal-adaptée aux problèmes transitoires : une boucle devant être effectuée à chacun des pas
de temps pour que la solution soit satisfaisante [23]. Les méthodes les plus populaires sont
basées sur un découplage des équations et sont souvent appelées méthodes de correction de
pression. En 1965, Harlow et Welsh [37] ont proposé la première méthode où les équations
sont découplées, soit le schéma "Marker-and-Cell" (MAC). Leur méthode était explicite et in-
efficace pour la résolution d'écoulements permanents. En 1972, Patankar et Spalding[58] ont
proposé l'algorithme "SIMPLE" ("Semi-Implicit Method for Pressure Linked Equations")
où l'équation de correction de pression est déterminée à l'aide des coefficients de l'équation
discrétisée de la quantité de mouvement. Pour être vraiment efficace, cet algorithme néces-
site l'optimisation de paramètres de relaxation. Une version améliorée et largement utilisée
encore aujourd'hui a ensuite été proposée, il s'agit de l'algorithme SIMPLER ("SIMPLE
Revised")[59]. Ce dernier algorithme peut, lui aussi, diverger lorsque des paramètres de re-
laxation ne sont pas adéquatement déterminés. En 1984, Doormal et Raithby[33] ont introduit
l'algorithme SIMPLEC ("SIMPLE Consistent") qui ne nécessite aucun paramètre de relax-
ation.
Parallèlement à la famille SIMPLE, les méthodes à pas fractionnaires se sont dévelop-
pées à partir de 1968 avec le schéma proposé par Chorin[25], soit le schéma de projection-1.
Ce schéma est basé sur une séparation complète des opérateurs : le terme de pression est
omis des équations de conservation de la quantité de mouvement et une deuxième équa-
tion transitoire est construite en utilisant l'équation de continuité et le terme de pression. Un
champ de vitesse satisfaisant l'équation de continuité est déterminé avec cette dernière équa-
tion. La séparation complète des opérateurs n'est pas sans conséquence, la solution finale
étant fortement dépendante du pas de temps (même lorsque l'écoulement est permanent). En
1981, Leveque[54] a démontré que pour résoudre des équations évolutives avec une sépa-
ration complète des opérateurs, il fallait modifier les conditions de bord. En 1985, à partir
de ce résultat, Kim et Moin[46] ont proposé un schéma à pas fractionnaires avec lequel il
est possible d'obtenir de bons résultats avec une discrétisation temporelle explicite du terme
de convection. D'autres chercheurs ont proposé de conserver le gradient de pression dans
les équations de conservation de la quantité de mouvement et d'introduire une correction de
pression (voir [67], [68], [69] et [55]). Avec ces schémas, il est possible d'obtenir de bons
résultats lorsque la discrétisation temporelle est implicite. Les schémas de cette famille sont
souvent appelés schémas de projection-2 ou projection-3. A priori, contrairement aux algo-
rithmes SIMPLE et SIMPLER, ces méthodes de projection ne nécessitent aucun paramètre
de relaxation.
Finalement, on peut mentionner qu'au lieu de procéder avec des méthodes de pro-
jection (méthodes à pas fractionnaires et algorithmes de la famille SIMPLE), il est aussi
possible de résoudre les équations de conservation de la quantité de mouvement et de la
masse à l'aide d'un changement de variable. Ce changement de variable est effectué en
calculant le rotationnel des équations de la quantité de mouvement et en introduisant les
fonctions de vorticité(u;) et de courant(^)[44]. Bien que l'on puisse, avec ces méthodes, ré-
soudre très rapidement les écoulements qui ne dépendent pas d'autres quantités scalaires,
il y a des désavantages. Entre autre chose, la formulation des conditions de bord est prob-
lématique et il est difficile d'adéquatement tenir compte de propriétés physiques variables.
Finalement, l'extension 3D est peu pratique, le nombre de variables à déterminer passant de
4 à 6 (ipx, tpy, ijjz, u>x, ujy, uz).

Pour bien résoudre les écoulements incompressibles, les chercheurs ont constaté qu'il
fallait introduire des espaces d'approximation différents pour la vélocité et la pression. Cette
construction est nécessaire car elle permet de faire un bon couplage entre la pression et la
vitesse et d'éviter des oscillations non réalistes de la fonction de pression (ces oscillations
sont souvent appelées faux modes de pression). En 1965, Harlow et Welsch[37] ont proposé
d'utiliser des maillages décalés ("staggered grid") où la pression et la vélocité ne reposent pas
aux mêmes positions. Durant de nombreuses années, cette formulation a pratiquement été la
seule avenue possible pour les méthodes de volumes finis et de différences finis. En 1983,
Rhie et Chow[26] ont montré que les faux modes de pression pouvaient être évités à l'aide
d'une ré-interpolation de la vélocité. Ainsi, le même maillage peut alors être utilisé pour la
vitesse et la pression. Il a été démontré que la ré-interpolation de Rhie-Chow et les maillages
décalés sont équivalents[48]. En ce qui concerne les méthodes de la famille des éléments
finis, des espaces d'approximation différents pour la vélocité et la pression sont directement
introduits à l'aide d'éléments dits "admissibles". La combinaison d'espaces est admissible
si et seulement si la condition Brezzi-Babuska[34] est satisfaite. Les principales familles
d'éléments admissibles sont les éléments de Taylor-Hood[42] et de Crouzeix-Raviart[29].
Outre la problématique liée à l'incompressibilité, la résolution numérique d'équations
de convection-diffusion introduit des difficultés qui lui sont propres. En effet, lorsque le rap-
port entre les forces d'inertie et les forces visqueuses ou de diffusion (nombres de Reynolds
et de Peclet) est grand, les solutions obtenues avec une approximation du deuxième ordre
(différence centrée) du terme de convection présentent généralement des oscillations qui sont
physiquement irréalistes. La convergence des schémas est alors dépendante du pas d'espace,
une solution réaliste ne pouvant être obtenue qu'avec un maillage comportant un très grand
nombre d'éléments. Dans la pratique, particulièrement en 3D, cette solution est extrêmement
coûteuse en temps de calcul et en espace mémoire et ne peut pas être envisagée. L'approxima-
tion (ou discrétisation) en amont du terme de convection s'est très tôt avérée une solution à ce
problème. Les schémas amont n'ont pas tendance à produire de solutions irréalistes lorsque le
nombre de Peclet (ou le nombre de Reynolds) est grand. Pour de petits nombres de Peclet ou
de Reynolds, ils sont toutefois moins précis que les schémas utilisant une différence centrée
du terme de convection.
En 1962, Spalding[70] proposa le schéma hybride qui fut par la suite perfectionné
par d'autres chercheurs ([7] et [59]). Le schéma hybride tente de combiner les avantages
des schémas amont et centrés : le terme de convection est approximé avec une discrétisation
amont dans les régions où la convection prédomine et une différence centrée dans les autres
régions. Ce schéma donne de bons résultats lorsque l'écoulement est permanent et orienté
avec le maillage. De plus, les termes sources doivent être constants et négligeables (voir [30],
[62] et [53]). D'autres schémas n'ayant pas les désavantages du schéma hybride ont aussi été
proposés pour les maillages structurés (Leonard[52] et Raithby[63]). Le schéma de Leonard
résoud avec succès les problèmes liés à la présence du terme source et de la fausse diffusion
causée par le non-alignement de l'écoulement avec le maillage. Toutefois, il n'assure pas une
solution monotone ou satisfaisant le principe du maximum.
 Un effort considérable a été fait pour résoudre le problème de la fausse diffusion
causée par un écoulement oblique aux lignes directrices du maillage1. En 1979, Van Leer[51]
a introduit les schémas d'approximation basés sur une reconstruction de la solution par vol-
ume de contrôle, les flux d'interface étant ajustés en fonction de paramètres d'interpolation.
Ces schémas, qui sont souvent appelés schémas MUSCL ("Monotone Upstream-centered
Schemes for Conservation Laws"), sont particulièrement bien adaptés aux méthodes de vol-
umes finis classiques. Ils ne sont toutefois pas adaptés à une discrétisation temporelle im-
plicite et le pas de temps doit respecter des contraintes de type CFL : le pas de temps est alors
limité par un quotient qui dépend de la vitesse et de la taille des mailles. De plus, cette plus
grande précision est atteinte au détriment de la stabilité, la solution n'étant pas nécessaire-
ment bornée. Afin d'éviter les oscillations et d'obtenir une solution monotone, de nombreux
schémas d'ordre 2 ont par la suite été construits avec l'approche "TVD" ('Total Variation
Diminishing"). Cette approche est assez générale : l'idée centrale étant de limité les flux pour
que qu'il n'y ait pas d'accroissement des variations. En 1987, Yee[75], proposa une approche
générale pour la construction de schéma dits "TVD". Les solution calculées avec ces sché-
mas ne présenteront pas d'oscillations causées par un trop grand pas d'espace. De plus, ces
schémas sont bien adaptés à la résolution d'écoulements transitoires où les équations de con-
servation n'ont pas de terme source. En effet, la discrétisation temporelle est souvent explicite
et les limiteurs de flux ne tiennent pas compte des termes sources.
En 1980, Baliga et Patankar ([5], [4] et [5]) ont proposé d'approximer les scalaires et
les composantes du vecteur vitesse avec une fonction solution d'un problème de convection
diffusion mono-dimensionnel. Cette approche a été utilisée avec succès avec des méthodes
de volumes finis de type éléments finis, de volumes finis classiques et d'éléments finis. L'a-
vantage de cette approche est d'introduire beaucoup moins de dissipation numérique que
les schémas amont. Il est toutefois possible d'obtenir des solutions qui ne respectent pas le
principe du maximum.
Pour les méthodes d'éléments finis, la construction de schémas dont la stabilité n'est
pas conditionnelle à la taille des mailles ni au pas de temps est encore aujourd'hui un prob-
lème ouvert. Des schémas plus stables peuvent être construits en introduisant implicitement
'Précisons que ces lignes directrices peuvent ne pas existées lorsque le maillage est non structuré (non
orienté avec un système de coordonnées orthogonales).
 8

de la dissipation numérique (voir [11], [12] et [14]) ou en modifiant les fonctions tests. Ces
schémas numériques sont souvent appelées méthodes d'éléments finis Petrov-Galerkin dif-
fusion ligne de courant (SUPG) ou méthodes de diffusion ligne de courant (méthodes SD),
elles ont été introduites par Brooks et al.[35] et elles sont discutées en détail dans l'ouvrage
collectif [28].
Il n'existe donc pas de solution idéale au problème de la discrétisation du terme
de convection. Les méthodes les plus stables sont théoriquement moins précises et souvent
conçues en étudiant des problèmes mono-dimensionnels. Les méthodes les plus précises sont
construites au détriment de la stabilité où ne peuvent être utilisées que pour certains prob-
lèmes en particulier. Finalement, on tient à préciser que la fausse diffusion causée par un
écoulement oblique aux lignes directrices du maillage est atténuée lorsque le maillage est
non structuré (on rappelle que ces maillages ne possèdent pas nécessairement de lignes di-
rectrices).

À la lecture de ces résultats, il s'avère que la discrétisation spatiale par la méthode

des volumes finis classiques soit la seule qui, théoriquement pour des maillages et des prob-
lèmes de convection-diffusion généraux, permette de construire un schéma qui respecte les
propriétés suivantes :

1. la conservation locale des flux ;

2. le respect du principe du maximum discret lorsqu'il s'applique.

Ces deux propriétés sont particulièrement importantes pour la résolution de problèmes indus-
triels où l'écoulement peut dépendre des quantités scalaires transportées. En effet, lorsque
ces propriétés sont respectées, les quantités scalaires transportées par l'écoulement seront
toujours en accord avec la physique du problème, et ce, sans qu'il soit nécessaire d'utiliser
un maillage possédant un très grand nombre d'éléments. Malheureusement, pour utiliser
une approximation constante par élément (méthodes de volumes finis classiques) lorsque
la géométrie est complexe, une ou l'autre de ces techniques de construction de maillage doit
être envisagées :

1. une décomposition par bloc du domaine et l'utilisation d'un maillage structuré par bloc
avec changement local de coordonnées ;
2. un maillage construit avec une technique Voronoï.
En ce qui concerne les maillages structurés par bloc, il y a des problèmes liés à la conserva-
tion des flux entre les blocs. De plus, lorsque le domaine est relativement complexe, il peut
être difficile (si non impossible) de le décomposer en blocs. Les maillages de type Voronoï
sont quant à eux difficiles (ou impossibles) à générer lorsque le domaine est tridimensionnel
et qu'il comporte des obstacles. Finalement, dans les deux cas, il peut être impossible de con-
struire un maillage qui soit bien adapté à la solution. On tient à préciser que ces contraintes
sur le maillage sont essentiellement liées à une approximation consistante des flux diffusifs
entre les cellules. Lorsqu'il n'y a pas de diffusion et qu'aucune dérivée seconde n'a à être
approximée, le maillage n'a pas à être orthogonal ou bien Voronoï (ce sujet sera discuté dans
le chapitre 3).
Le schéma numérique présenté dans cette thèse est une nouvelle méthode de volumes
finis classiques. Ce schéma est inspiré des résultats théoriques publiés par Eymard et al.[36]
et d'une autre méthode de volumes finis classique 2D pour maillages non structurés proposée
par Boivin et al. ([9], [10] et [19]). Les principales caractéristiques de ce nouveau schéma
sont les suivantes :

1. toutes les variables scalaires et les composantes du vecteur vitesse sont approximées
avec des fonctions constantes par volume de contrôle ;

2. les bilans sont calculés sur des volumes de contrôle arbitraires composés d'un assem-
blage de tétraèdres ;

3. les flux sont localement conservés ;

4. le principe du maximum discret est localement respecté lorsqu'il s'applique.

L'idée principale de cette méthode de volumes finis est d'associer à tous les éléments géomé-
triques de la discrétisation (triangles en 2D ou tétraèdres en 3D) une position de référence
pour le calcul des flux. Les éléments géométriques sont ensuite combinés pour former les vol-
umes de contrôle. Les fonctions d'interpolation pour toutes les variables du système d'équa-
tions sont constantes par volume de contrôle. Ce schéma numérique est convergent, stable et
permet de respecter des lois importantes des modèles continus : notamment la conservation
locale des flux et le principe du maximum lorsqu'il s'applique.
 10

Cette thèse est organisée comme suit :

• Le chapitre 1 présente un rappel de la dérivation des équations de conservation qui

déterminent les écoulements incompressibles.

• Dans le chapitre 2, on présente une description du modèle de turbulence qui a été utilisé
pour valider le schéma.

• Les résultats théoriques qui ont été utilisés pour développer le schéma sont présentés
dans le chapitre 3.

• Le schéma numérique est présenté dans le chapitre 4. La discrétisation temporelle,

la discrétisation spatiale, la formation des volumes de contrôle, le calcul des flux et
l'approximation des gradients sur les volumes de contrôle y sont discutés.

• Dans le chapitre 5, on présente les algorithmes de résolution numérique de systèmes

d'équations linéaires qui ont été mis en oeuvre avec notre schéma.

• Suite aux résultats obtenus avec notre schéma, une librairie informatique destinée a être
utilisée par d'autres chercheurs a été développée. Cette librairie est présentée dans le
chapitre 6.

• Les principaux résultats numériques obtenus lors de validation du schéma sont présen-
tés dans le chapitre 7.
• Finalement, dans le chapitre 8, on propose des modifications au schéma pour les prob-
lèmes où la densité est variable.
 Chapitre 1

Modèle mathématique, écoulement laminaire

Introduction

De nombreux problèmes d'ingénierie sont liés au transport d'une ou plusieurs quan-

tités scalaires par un fluide newtonien. Lorsque le fluide est homogène, le modèle mathéma-
tique décrivant l'écoulement est composé des équations de conservation de la masse, de la
quantité de mouvement et de l'équation de conservation de l'énergie. Dans ce chapitre, on
présente ces équations ainsi qu'un rappel de leur développement.

1.1 Notions fondamentales

Toutes les équations de conservation sont dérivées en établissant un bilan sur un mi-
lieu isolé séparé du milieu extérieur par une surface réelle ou imaginaire. Ce milieu isolé
est souvent appelé volume de contrôle. Ce volume de contrôle peut être fixe (volume de
contrôle fixe), avoir un déplacement arbitraire (volume de contrôle arbitraire) ou bien se
déplacer avec l'écoulement (volume de contrôle matériel).
Soit / un flux quelconque (flux de masse, d'énergie,...) ou une force s'exerçant sur
la surface dQ d'un volume de contrôle arbitraire fi. Le bilan surfacique de ce flux peut-être
calculé à l'aide d'une intégrale volumique de la divergence de ce flux à l'intérieur du volume.
Ce passage d'une intégrale de volume à une intégrale de surface (ou l'inverse) est effectué à
l'aide du théorème de Gauss-Ostrogradsky (ou théorème de la divergence).
 12

Théorème 1.1 Théorème de Gauss-Ostrogradsky(Joyal et Mercier[47])

Considérons le volume f2 borné de H 3 dont le bord dQ est une réunion finie de sur-
faces fermées lisses par morceaux et orientables. Si (p(x) : R 3 —» H 3 est de classe C1 sur
Q, U dQ, alors1
f V-0dV= [ $-ndA,
où:

• n est le vecteur unitaire normal pointant vers l'extérieur de dVt ;

• dA est un élément de surface.

Considérons une particule et une quantité scalaire 0. Lorsque cette particule est fixe,
le taux de variation de <f> par rapport au temps est exprimé à l'aide de la dérivée partielle

~dt'

Lorsque (j> est une variable d'un écoulement, cette dérivée est souvent appelée dérivée locale.
En générale, la quantité 4> dépend de la position dans l'espace. Donc, lorsque cette particule
se déplace, la variation de (f) par rapport au temps doit aussi tenir compte de la variation de cj)
par rapport à l'espace. La dérivée convective

(v • V)<f>

permet de calculer cette variation. Pour tenir compte à la fois de la variation locale et de la
variation spatiale, le taux de variation de 4> par rapport au temps est calculé à l'aide de la
dérivée matérielle (ou dérivée particulaire)

Cette dérivée exprime le taux de variation de <j> par rapport au temps lorsque l'on suit une
particule se déplaçant à une vitesse v_.
Considérons une quantité scalaire <>
/ et un volume de contrôle isolé Q,. Le taux de
variation de 4> dans ce volume de contrôle est égal à la somme du bilan de la variation locale
'On rappelle qu'une fonction de classe C 1 est une fonction au moins une fois derivable et ses dérivées sont
continues.
 13

de 4> dans Q, et du flux de <f> traversant la surface dQ,. Ce taux de variation peut être calculé à
l'aide du théorème de transport (ou règle de Leibnitz).

Théorème 1.2 Théorème de transport(Candel[16J)

Soit Q,m(t) un volume de contrôle matériel quelconque limité par une surface dQ.m{t)
et 4>{x_, t) une fonction définie sur Q. Le taux de variation dans le temps de l'intégrale

<t>dV
• "

est donné par l'expression

<f>dV = / — dV + / (j>v_-ndA (1.1)

où il est la vitesse locale du fluide et n est la normale extérieure de <

Remarques.
1. Le théorème de transport peut aussi être appliqué (sous une forme différente) à un
volume de contrôle arbitraire ou fixe.
2. Lorsque le volume matériel Œm(t) et le volume fixe Q, coïncident à l'instant t, le taux de
variation pour un volume matériel peut être exprimé à l'aide de l'expression suivante :

(f,dV = [ ^rdV+ [ cj>v-ndA.

dtja

3. Pour alléger la notation, on laissera tomber le sous-indice m(t) pour indiquer si le

volume est matériel ou fixe.
Le théorème de la divergence permet de transformer l'intégrale de surface de l'équation (1.1)
en une intégrale de volume :

4- [ 4>dV = f ^ d V + f (v • V(t> + <t> V • v) dV.

dt Jçi Jo. at Jn

Cette équation fait apparaître deux termes dont la signification physique est importante :

1. IV' V(j) dV détermine le taux de variation de 0 causé par la variation spatiale de 0

Jn
(les particules se déplaçant) ;
 14

2. / <f> V • v dV est le taux de variation déterminé par la dilatation (ou la compression)

du volume. Cette dilatation (ou compression) du volume est le taux de changement du
volume par rapport au temps. Il est égal à la divergence de la vélocité V • v.

1.2 Conservation de la masse

Considérons un volume matériel Q, . La masse m contenue dans ce volume peut être

calculée à l'aide de l'expression
m — j pdV
Jn
où p est la masse volumique (ou densité) qui peut dépendre du temps et de l'espace. Si ce
volume ne contient ni source ni puit, la masse contenue dans ce volume est constante. Le taux
de variation de la masse contenue dans ce volume doit donc être nul :

[ [[
dt JQ Jn dt Jaçî
Lorsque l'on applique le théorème de la divergence à l'intégrale de surface on obtient :

dm d

Puisque le volume de contrôle est arbitraire, seul l'intégrant peut être considéré :

§£ + v.(«o = o.
Souvent, cette équation est appelée équation de continuité.

Remarque. Bien que seul le volume matériel soit considéré, l'équation de conservation
pour tout autre volume de contrôle est la même, seul le domaine d'intégration est modifié.
Lorsque l'écoulement est incompressible, la pression et la température ne font varier
que très faiblement la densité du fluide. Si le fluide est homogène, la densité peut alors être
considérée comme étant constante. Pour ces cas, le respect de l'équation

V-v = 0
 15

en tout point impliquera le respect de la conservation de la masse sur tout volume de contrôle.
D'un point de vue physique, il ne peut donc pas y avoir de compression ou de dilatation d'un
volume de contrôle matériel lorsque l'écoulement est incompressible. Toutefois, ce volume
peut être déformé puisque le cisaillement

n'est pas nécessairement nul.

Il est important de préciser que la distinction compressible/incompressible s'applique
à l'écoulement plutôt qu'au fluide. En effet, la densité d'un fluide est en général dépendante
de la température et de la pression. Toutefois, la compressibilité de l'écoulement peut être
caractérisée par le nombre de Mach
M = V-
c
où:

• c = —rT est la vitesse locale de perturbation acoustique,

cv
• Cp est la chaleur spécifique à pression constante,

• cv est la chaleur spécifique à volume constant,

• r est la constante des gaz parfaits par unité de masse.

Ap
En effet, le rapport — est de l'ordre du nombre de Mach. Lorsque M < 1 l'écoulement
P
peut être considéré incompressible.

Forme conservative d'une équation de conservation

Toutes les équations de bilan sont des équations de conservation et elles sont conserva-
tives. Il existe plusieurs façons de représenter ces équations. Pour s'assurer de la conservation
locale des flux entre les volumes de contrôle, les équations de bilan doivent être écrites sous
forme conservative. Soit (j> une quantité scalaire quelconque, une équation de conservation
pour (j) est sous forme conservative lorsqu'elle peut s'écrire :
 16

où s est un terme source et / représente un mécanisme de transport telle la convection ou la

diffusion. Cette formulation des équations de bilan convient bien à la résolution numérique.
Entre autre, elle permet :

1. d'éliminer les "faux" puits ou sources de masse ;

2. d'abaisser l'ordre des dérivées à l'aide du théorème de la divergence ;

3. d'obtenir plus facilement une formulation stable des équations de convection-diffusion

à l'aide d'une dérivée amont (section 3.2.3).

1.3 Conservation de la quantité de mouvement

La quantité de mouvement contenue dans un volume de contrôle matériel est :

L vpdV,

où v (x_, t) et p (x, t) peuvent dépendre de la position et du temps. La variation de quantité de

mouvement est déterminée par le principe fondamental de la dynamique. Cette variation est
égale à la somme des forces s'éxerçant sur le volume :

Y [ vpdV = [ £dV+ [ t(n)dS

dt JdQ Jçi Jdn

où :

• I fdV est l'ensemble des forces volumiques (telle la gravité) s'éxerçant sur Q ;
Jn ~
• / t(n) dA représente les forces appliquées sur la surface dSl, t(n) étant le vecteur
Jan
des contraintes.

1.3.1 Vecteur des contraintes

Lorsque le fluide est au repos, la contrainte agissant sur un volume est normale à sa
surface :
t(n) = —pn,
 17

p étant la pression. Lorsque le fluide est en mouvement, en plus de subir la contrainte de

pression, le volume doit opposer une résistance à la déformation et il y a friction entre les par-
ticules. Cette résistance est déterminée par le tenseur des contraintes visqueuses a. Ainsi,
pour un fluide en mouvement, t(n) est formé de deux composantes :

t(n) — ~PR + H ' R •

La forme générale du tenseur des contraintes visqueuses peut être établi à partir des hy-
pothèses suivantes (Landau et Libshitz[49]) :

• Le volume ne peut se déformer que lorsque les particules se déplacent avec des vitesses
différentes. Ainsi ce tenseur ne doit dépendre que de la variation spatiale de la vitesse.
Lorsque celle-ci est nulle ce tenseur doit être nul.

• Le fluide est newtonien, la relation entre la résistance du volume à la déformation et la

variation de vitesse est linéaire.

• Lorsque la rotation du fluide est uniforme, il n'y a pas de friction entre les particules.
Dans ce cas, ce tenseur doit être nul.

Une forme générale qui permet de respecter toutes ces hypothèses est :

a = ad + WV • v_,

où :

1. a et b sont des constantes qui sont associées à la viscosité du fluide

2. d est le tenseur des taux de déformation :

1 ( dui d%
dij = « 1 â— + â"
2 \ OX • £73

3. S est le delta de Kronecker :

0,
 18

Expérimentalement, on peut déterminer que a — 2/x où \i est la viscosité dynamique du

fluide. Pour b la situation est plus complexe et controversée (White[73j) :

• Stokes a déterminé que b devait être posé égale à — | / i .

• Karim et Rosenhead (1952) ont obtenu expérimentalement des résultats qui contredis-
ent la valeur déterminée par Stokes : pour la majorité des liquides, cette valeur serait
positive et beaucoup plus grande que //.
• Truesdell (1954) a mis en doute les résultats des expériences de Karim et Rosehead.

D'après White, cette seconde viscosité ne devrait être considérée que pour les problèmes
suivants :

1. absorption d'ondes par des liquides ;

2. les chocs.

Une des formes les plus courantes du tenseur des contraintes visqueuses est obtenue par
l'introduction de la viscosité volumique (ou seconde viscosité) C = b — \y, :

g = 2/i U - -£V • v J + C£V • v.

Lorsque le gaz est mono-atomique, le coefficient de viscosité Ç est nul (en général, on n'en
tient pas compte). De plus cette formulation est plus "réaliste", les coefficients de viscosité
étant toujours positifs. Lorsque le fluide est incompressible, la divergence du vecteur vitesse
est nulle et a peut être simplifiée :

g = l* (Vu + Vu*) .

Lorsque l'écoulement est incompressible, l'équation de conservation de la quantité de mou-

vement est :

dV +
/ s* / [v®(pv)-iJ.(Vv + 'Vvt)]-ndS = / f_dV - / pndS.

Généralement, cette équation est simplifiée. Lorsque la viscosité est constante le terme
 19

est nul. Lorsque l'écoulement est incompressible, bien que la viscosité puisse être variable, il
est fréquent de considérer que ce terme n'influence pas significativement l'écoulement et de
ne pas le conserver(Chabard[20], Wilcox[74])2 :

r <H0_dV+ f [ y ® ^ -^Vv)-ndS = [ idV- [ pndS. (1.2)

Pour regrouper tous les termes sous une intégrale volumique on a qu'à appliquer le théorème
de divergence. Puisque notre volume d'intégration est quelconque seul l'intégrant de l'équa-
tion résultante peut être considéré :

^ v)-fj,Vv} = /-Vp. (1.3)

1.3.2 Conditions de bord pour l'équation de conservation de la quantité

de mouvement

La frontière dQ peut être séparée en N éléments de frontière. En général, on dis-

tinguera trois types de frontières :

1. Entrée : la vitesse normale (ou le débit massique) à l'entrée est imposée :

v-n = Vb-

2. Sortie : la vitesse de sortie est inconnue, la contrainte normale est imposée :

—pÔ + /j, (Vv) = g_-n.

Souvent, la sortie est libre et la contrainte normale est imposée égale à zéro à l'aide
d'une isobare et d'une dérivée normale nulles.

3. Paroi : la vitesse normale et la vitesse tangente sont nulles :

v • n - 0, v • 7 j = vTl , v • r2 = vT2,

Hi e* I2 s o n t d es vecteurs orthogonaux tangents à la surface.

2
Pour être tout à fait rigoureux, il faudrait conserver le terme JgQ (/xVt;*) • ndS lorsque la viscosité est
variable.
 20

1.3.3 Convection naturelle, approximation de Boussinesq

Dans de nombreux procédés thermiques, l'écoulement peut être induit par des forces
de poussée liées à des fluctuations de température. Lorsqu'il y a de la convection naturelle, le
fluide est mis en mouvement par des gradients de masse volumique résultants des gradients de
température. Ces gradients de masse volumique entraînent des variations de force de volume
induites par le champs gravitationnel :

/ = pg,

g étant l'accélération gravitationnelle. Si l'on veut considérer la densité constante, cette force
peut être prise en compte via l'approximation de Boussinesq.
La pression est décomposée en deux parties : ph pour la pression hydrostatique et
pd la pression dynamique causée par le mouvement du fluide. La pression dynamique est
inconnue et déterminée par les équations de conservation de la quantité de mouvement et
de la masse. Lorsque le milieu ambiant est calme, le gradient de pression hydrostatique est
connu et proportionnel au champs gravitationnel :

= pog,

Po étant la masse volumique associée à l'état de référence. La force induite par la pression
hydrostatique peut ainsi être incorporée aux forces volumiques :

/ = (p~Po)g-

Puisque l'écoulement est incompressible, la masse volumique du fluide ne varie qu'en fonc-
tion de la température. En procédant à une expansion en série de Taylor de p on obtient :

P = *>+($) ( r " T o ) + ( 0 ) (T-To)2... (1.4)

Cette expression peut être simplifiée en introduisant le coefficient d'expansion volumique

0:
 21

B-

En substituant ce coefficient dans l'équation (1.4) on a

-^=^ 1 -» 0. La force volumique

/ p
s'exerçant sur le fluide est alors approximée par :

/ S -g0p(T - To) .

Toutefois, pour pouvoir utiliser cette approximation, deux autres conditions doivent être re-
spectées. Premièrement, le rapport entre la variation de densité et la densité doit être petit :

Deuxièmement, dans l'équation de continuité, la variation de densité doit être négligeable par
rapport à la variation de vitesse :

Avx &Vy kvz\ /A/) Ap Ap\

A
+ + )+v- M - , - ^ , 7 - = 0
~ \Ax Ay AzJ
1 (Avx Avy AvA j_l (Ap Ap A/A
Ay AzJ \\v\\p'\Ax'Ay'Az)

où A(f> est la variation de </>. Dans les cas où

Au Ap
ll^ll P
la densité peut être considérée constante dans l'équation de conservation de la masse puisque

(
Ap Ap A p \
——, ——, —— j est négligeable par rapport à celle de
{Avx | Avy | A^
P
\Ax Ay Az
 22

1.4 Équation de conservation de l'énergie

L'équation fondamentale pour établir un bilan d'énergie est l'équation de conservation

de l'énergie totale. Elle tient compte de l'énergie interne, de l'énergie cinétique et lorsque
nécessaire de l'énergie potentielle. Lorsque l'écoulement est compressible et que la densité
dépend à la fois de la pression et de la température, cette équation est essentielle pour avoir
un système complet.

1.4.1 Équations de conservation de l'énergie totale, cinétique et interne

L'équation de l'énergie totale est déterminée par l'application du premier principe de

la thermodynamique : la variation d'énergie d'un système est égale au travail effectué sur
ce système et à la chaleur apportée à ce système. L'application de ce principe à un volume
matériel £2 donne :

^- f p(e + ]-v2) dV = [ {s + pl-v) dV+ [ t_(n) • vdS - [ ç" • ndS,

dt Jn \ l ) Jn Jan Jan

où:

• — / p ( e + -v2 ) dV est la variation d'énergie interne et cinétique ;

dt Jn \ 2 /

• / / • y. dV est le travail des forces de volume ;

./fi ~

• / t{ll) -y. dS représente le travail effectué par les contraintes sur la surface du volume ;

• — / q" -ndA représente la chaleur qui quitte le volume Q, par sa surface ;

Jan ~
• s est une source ou un puits de chaleur.

En appliquant le théorème de transport au membre de gauche de cette dernière équation et en

explicitant le tenseur des contraintes on obtient :

/ J dV+ p (e + -v2 )v • n dA

= / s + pf-vdV+ [(-pâ->r vÇVy))-v]-ndA- (f'-ndA.

Jn ~ Jan Jan
 23

Ensuite, le théorème de Gauss-Ostrogradsky permet de transformer les intégrales de surfaces

en intégrales de volume. Le volume de contrôle étant quelconque, seul l'intégrant peut être
considéré :

dp(e + W) „ \ ( 1 2\

= s + p / • v + V • [ (-p£ + // (Vv)) -v] - V • g".

Cette équation établi le bilan d'énergie totale. Pour obtenir l'équation de bilan de l'énergie
interne, l'énergie cinétique doit être soustraite de cette équation. L'énergie cinétique est cal-
culée en effectuant le produit scalaire entre l'équation de bilan de la quantité de mouvement
écrite sous forme locale et le vecteur vitesse :

+P
=~

=
dt \ 2PV<2 I + V
' ( 2Pv2-) +
~ ' V P ~ - ' V ' ^ V ~ ^ ~ P- ' ~ = ° '

Cette équation détermine la variation d'énergie cinétique par rapport au temps. Lorsque sous-
traite de l'équation de conservation de l'énergie totale, on obtient l'équation de conservation
de l'énergie interne :

• Ve = s — pV • v_ — V • q + V • (g_v) •

1.4.2 Équation de bilan de l'énergie thermique

Lorsque l'on ne s'intéresse qu'à l'énergie thermique, l'énergie interne doit être rem-
placé par l'enthalpie massique h = e + P - :
P
dh _, _ , , dp N ,
P^r + pvVh = -VV • (vp) £ - V • q + V • (av)
dt + pvVh = (vp) - ut
•£ + s.
ut - \— '
De plus, il peut être démontré (Chabard [20]) que la dérivée particulaire de l'enthalpie mas-
sique peut être écrite sous cette forme :

v
dt dt pK H
' dt
 24

Cette équation permet d'écrire une équation de bilan thermique en fonction de la température :

pcp— + pcpv • VT = pT(-£ + v-Vp) - V -q + V • (gv) + s .

Pour de nombreux cas, les simplifications suivantes sont appropriées :

• lorsque le nombre de Mach est faible, la dissipation visqueuse est faible :

V- (gv) ~ 0 ;

• les effets de dilatation thermique sont négligeables :

dp \
— + v • Vp) ~ 0 ;
ot )
• le flux local de chaleur par diffusion respecte la loi de Fourier :

A est la conductivité thermique qui peut être variable.

En définitive, l'équation de bilan thermique locale est :

dT
pCp-^ + pCpV-VT = V-(A;VT) + s.

Il est important de remarquer que cette équation n'est pas écrite sous forme conservative. Elle
n'est donc pas appropriée pour la résolution numérique par volumes finis. Toutefois, lorsque
la chaleur massique ne présente que de faibles variations spatiales et que l'écoulement est
incompressible, cette équation peut être exprimée sous forme conservative :

f d(pCpT) dV+ f [(vpcJr)-n-\VT\-ndS = [sdV,

Jn ot JaQ Jn

= s.
 25

Conditions de bord de l'équation de l'énergie

La frontière du domaine peut être séparée en n éléments de frontière. Sur chacun de

ces éléments de frontière, la condition de bord peut être du type :

1. Dirichlet, température Tt imposée

T = Tb;

2. Neuman, flux diffusif q'i imposé

3. Mixte (condition de Robin), flux convectif imposé

hc{T-TO0)=<£,

où hc est le coefficient de transfert convectif, T^une température de référence ;

4. Radiatif, le flux radiatif est imposé

hT = eaB(T2 + T

où hr est le coefficient de transfert radiatif, aB la constante de Bolzman et e l'émis-

sivité thermique.

1.5 Modèle complet

Lorsque les hypothèses suivantes sont respectées :

1. le nombre de Mach est petit : M < 1 ,

2. la viscosité est constante, ou bien lorsqu'elle est variable, les termes -— ( u—^-
ox \ dx
peuvent être négligés,
3. la chaleur massique ne présente que de faibles variations spatiales,

l'écoulement incompressible d'un fluide newtonien homogène peut être modélisé avec les
équations de
 26

• conservation de la masse
——h V • (pv_) = 0, (1.5)

• conservation de la quantité de mouvement

d(pv)
o 4- V • \pv <g> v — p5 — /J,VV\ = /, (1.6)

• conservation de l'énergie

-^- h V • [pcpTv_ — ÀVT] = s. (1-7)

1.6 Nombres sans dimension

Dans bien des cas, les principales caractéristiques d'un écoulement peuvent être déter-
minées à l'aide de nombres sans dimension. Dans cette section, on présente les nombres sans
dimension qui sont utilisés dans ce document. La majorité de ces nombres apparaît naturelle-
ment lors de la construction des équations sans dimension [20]. Ces dernières équations ne
seront toutefois pas présentées.

Nombre de Reynolds

La quantité
pvD

où

• D est une longueur caractéristique qui dépend de la géométrie du problème,

• v est une vitesse de référence telle la vitesse maximale (ou moyenne) d'entrée dans une
conduite,

est appelée nombre de Reynolds. Il compare les forces d'inertie (termes convectifs) aux
forces de viscosité. Il caractérise les écoulements forcés dans une conduite ou autour d'un
obstacle. Ce nombre est souvent utilisé pour corréler des résultats expérimentaux à des ré-
sultats numériques. De plus, le passage du régime laminaire au régime turbulent est souvent
 27

caractérisé par ce nombre. Par exemple, pour une conduite cylindrique, le passage du régime
laminaire au régime turbulent a lieu dans l'intervalle 2100 < Re < 2500.

Nombre de Mach

Le nombre de Mach au carré

c2
ou

• c = —rT est la vitesse locale de perturbation acoustique,

cv
• cp est la chaleur spécifique à pression constante,
• cv est la chaleur spécifique à pression constante,
• r est la constante des gaz parfaits par unité de masse,
mesure le rapport entre l'énergie cinétique du gaz et son énergie interne. Pour les écoule-
ments à faible nombre de Mach, ce rapport est faible et seul l'énergie thermique (interne)
peut être considérée pour ces écoulements. De plus, le nombre de Mach mesure les effets de
compressibilité. Pour un faible nombre de Mach, les effets de compressibilité peuvent être
négligés et l'écoulement considéré incompressible. On peut aussi s'en servir pour qualifier la
vitesse d'écoulement :

• 0 < M < 1, vitesse subsonique;

• M > 1, vitesse supersonique.

La structure des ondes de choc autour d'un objet étant dépendantes du nombre de Mach, il
est souvent utilisé pour la description des principales caractéristiques des écoulements com-
pressibles à haute vitesse.

Nombre de Prandtl

Dans les écoulements thermiques, le rapport entre les transferts de quantité de mou-
vement par les forces visqueuses sur le transfert de chaleur par conductivité thermique est
appelé nombre de Prandtl :
 28

Ce nombre mesure l'efficacité des transferts de quantité de mouvement et d'énergie par dif-
fusion. Pour un gaz, il est de l'ordre de 1 et les transferts d'énergie thermique et de quantité
de mouvement par diffusion sont comparables. Pour un métal liquide, Pr <C 1 et le transfert
d'énergie par diffusion excède grandement le transfert de quantité de mouvement. Finale-
ment, le développement de la couche laminaire thermique est caractérisé par le nombre de
Prandtl puisque

épaisseur de la couche limite de la vélocité

épaisseur de la couche limite thermique

n est un nombre entier positif dépendant de la nature du fluide (Incropera[45]).

Nombre de Peclet

Similaire au nombre de Reynolds, le nombre de Peclet mesure le rapport entre le

transfert d'énergie thermique par convection et le transfert d'énergie par diffusion :

Les nombres de Reynolds et de Prandtl lui étant préférés, il est rarement utilisé pour la carac-
térisation des écoulements. Toutefois, ce nombre est souvent employé pour établir une con-
dition locale de stabilité lors de la discrétisation de l'équation de l'énergie avec des méthodes
de volumes finis, d'éléments finis ou de différences finies (section 3.2.3)3.

Nombre de Grashof

Le nombre de Grashof
3
Cr=zP>gD /3AT

où

• g est l'accélération gravitationnelle ;

• /? le coefficient d'expansion thermique ;
3
Lorsque l'on établi cette condition de stabilité, c'est le nombre de Peclet local qui est utilisé. Son analogue
pour l'équation de conservation de la quantité de mouvement est le nombre de Reynolds local.
 29

• AT un écart de température de référence (entre deux parois isothermes par exemple)

est employé pour comparer les forces de flottabilité causées par la différence de masse vo-
lumique entre le fluide chaud et le fluide froid et les forces visqueuses qui s'exercent sur un
volume élémentaire de fluide. Il apparaît dans les problèmes de convection mixte ou de con-
vection naturelle. Lorsque ce nombre est faible, la convection naturelle peut être négligée.
Ce nombre est particulièrement important pour les problèmes de convection libre. En effet,
pour ces problèmes, le nombre de Reynolds ne peut être calculé (il n'existe pas de vitesse
de référence) et le passage du régime laminaire au régime turbulent est alors fortement lié au
nombre de Grashof.

Nombre de Rayleigh

Un autre nombre utilisé pour l'analyse par similitude des écoulements engendrés par
la convection naturelle ou la convection mixte est le nombre de Rayleigh :

Ce nombre est similaire au nombre de Grashof. Il compare les forces de flottabilité dues à la
différence de masse volumique entre le fluide chaud et le fluide froid au produit des forces
visqueuses et thermiques qui s'exercent sur un volume élémentaire de fluide.
Pour faire une analyse par similitude d'un écoulement induit par la convection na-
turelle, une des deux combinaisons suivantes de nombres sans dimension est suffisante :

1. PretGr;
2. PretRa.

Nombre de Strouhal

Lorsque l'écoulement doit contourner un obstacle, il peut devenir instationnaire et

périodique. Alors, la vitesse dans le sillage de l'obstacle fluctue de façon quasi-sinusoïdale et
la fréquence fondamentale du passage des tourbillons est corrélée avec une distance D et une
vitesse caractéristique v via le nombre de Strouhal :
 30

Coefficient de portance

La portance est la force perpendiculaire à un corps en immersion dans un écoulement.

C'est la composante verticale de cette force qui assure la sustentation d'un avion . Le co-
efficient de portance permet de procéder à l'analyse par similitude de la portance de corps
placés dans un écoulement. Ce nombre est défini à l'aide du rapport

où

• v est la vitesse moyenne de l'écoulement autour du corps, telle la vitesse moyenne à

l'entrée d'une conduite,

• FL est la force de portance autour du corps,

• A est l'aire du corps.

Pour une aile d'avion ou une voilure, ce coefficient dépend entre autre de l'angle d'inci-
dence et du nombre de Reynolds, ce coefficient étant inversement proportionnel au nombre
de Reynolds.

Coefficient de traînée

La traînée est la force qui s'oppose à l'avancement d'un corps dans un fluide (ou si
l'on préfère, la résistance qu'oppose un corps fixe à un écoulement). Pour faire une analyse
par similitude de la force de traînée FD tangente à l'écoulement, on utilise le coefficient de
traînée
r
pv2 A
Dans le cas d'un écoulement aérodynamique autour d'une aile d'avion, on calculera la com-
posante horizontale de la force de traînée. Dans Landau et Lifchitz[49], on retrouve une dis-
cussion détaillée des forces de portance et de traînée pour les écoulements aérodynamiques.
 Chapitre 2

Modèle mathématique, écoulement turbulent

Introduction

Lorsque les composantes du vecteur vitesse, la pression ou la température évoluent

dans le temps de façon chaotique et que l'écoulement est instable, le régime de l'écoulement
est qualifié de turbulent. Ces écoulements ont les caractéristiques suivantes :

• Les variables, telles la vitesse, la pression et la température varient de façon aléatoire

et stochastique. Les variations spatiales sont nécessairement tridimensionnelles.

• Le rapport entre les forces d'inertie et les forces visqueuses ou le rapport entre les forces
de flottabilité et les forces visqueuses sont très grands. On rappelle que ces rapports sont
mesurés à l'aide des nombres de Reynolds et de Grashof (ou de Rayleigh).

• Ils sont généralement très dissipatifs, la diffusion de la quantité de mouvement, de la

masse et de la température est beaucoup plus grande que pour un écoulement laminaire.
Beaucoup d'énergie est éliminée et convertie en chaleur par dissipation visqueuse, l'é-
coulement à besoin d'un apport constant d'énergie. En fait, la turbulence accélère les
processus de mélange.

• II sont caractérisés par une large gamme d'échelles. Les phénomènes à grande échelle
sont liés à la configuration moyenne de l'écoulement. Les phénomènes à petite échelle
ont, quant à eux, une dynamique propre.

• Des tourbillons de différentes échelles de grandeur sont présents et ceux-ci sont con-
tinuellement créés et détruits.
 32

A priori, Les équations 3-D de Navier-Stokes sont valides à toutes les échelles de continuum,
elles permettent de décrire mathématiquement les écoulements laminaires comme turbulents.
Pour pouvoir simuler directement un écoulement turbulent, la dimension des mailles du mail-
lage doit être inférieure aux plus petits tourbillons. La taille de ces tourbillons peut être ap-
proximée à l'aide de la loi de Kolmogorov et elle est de l'ordre de vzl^. Pour un maillage
structuré, le nombre de points de la discrétisation est de l'ordre de O(Re9^). Pour un nombre
de Reynolds de 104, le nombre de points serait donc de l'ordre de 109. Lorsque l'on considère
que beaucoup de problèmes d'aérodynamique sont caractérisés par des nombres de Reynolds
de l'ordre de 106, il est impensable de pouvoir résoudre directement ce type de problèmes
avec les ordinateurs présentement disponibles.

2.1 Modélisation de la turbulence

Lors de la résolution de problèmes d'écoulements, l'ingénieur s'intéressse souvent au

comportement macroscopique de l'écoulement en non pas aux comportements à plus petite
échelle. Donc, ce sont les valeurs moyennes de la vitesse, de la pression, de la température
ou d'autres variables qui l'intéressent. Puisque l'interprétation des phénomènes observés lors
de la simulation d'un écoulement turbulent ne demande pas de connaître les fluctuations à
petite échelle, mais bien des valeurs moyennes de l'écoulement, une approche statistique de
la turbulence est souvent utilisée.

Décomposition de Reynolds

Considérons une variable 4- Celle-ci peut être décomposée en une partie représentant
sa valeur moyenne (4) et une autre partie représentant sa fluctuation 4' :
 33

Soient <j> et 7 deux variables quelconque, l'opérateur de moyenne () possède les propriétés
suivantes :

• (f)' = 0 ;
• 07' = 0 ;

7 = 0 + 7;

Pour un écoulement turbulent, trois types de moyennes sont utilisées (Wilcox[74])'.

Soit 4>(x_, t) une variable de l'écoulement, T un intervalle de temps et Q un domaine quel-
conque :

1. Pour un problème où la turbulence n'est pas homogène et l'écoulement est stationnaire,

la moyenne temporelle est
_ -i pt+T
4>(x) = lim — / (/>(x_, t) dt.
T->oo 1 Jt

2. Pour un problème où la turbulence est homogène, la moyenne spatiale est

4>{t) — lim — / (j>(x_, t) dV.

n-Hx il J

3. Pour un problème où l'écoulement est transitoire et la turbulence non-homogène, la

moyenne d'ensemble est

N
~$(x,t) = lim —Y] (l>i(x,t),
N-^oo TV r—'

iV étant un nombre de réalisations de l'expérience.

Les équations de conservation pour l'écoulement moyen sont obtenues en posant

'Lorsque l'écoulement est compressible, la moyenne de Favre est plus appropriée. Il s'agit d'une moyenne

pondérée avec la densité :<f>=- lim — / p(x,t)<j>(x,t)dt..

p r-+oo T Jt
 34

Ces valeurs sont ensuite substituées dans les équations de Navier-Stokes auxquelles on ap-
plique un opérateur de moyenne :

1. conservation de la masse :

2. conservation de la quantité de mouvement :

d (pi)) _ r _ _ — . ,_ _ fs -,
£ + V • [ pw (8> u + pv' g) u' - pô - /i (Vu + V i ) J = / ,

3. conservation de l'énergie :

^ ^ + V.[pc
L pTv + pcpW-\VT] J = s
at

Ces équations décrivent l'écoulement moyen et sont exactes puisqu'aucune simplification

ou hypothèse n'a été faite. Ce système n'est cependant pas fermé, puisque les corrélations
pv' (^) v! e* Tv! sont inconnues. Le terme y[ (££) v! est le tenseur de Reynolds ou tenseur des
contraintes turbulentes.
Le problème fondamental de la modélisation de la turbulence réside dans la modéli-
sation des ces contraintes. Pour obtenir des équations de transport pour ces 6 contraintes (le
tenseur étant symétrique), on définit l'opérateur N(vi) (Wilcox[74]) :

( | ) ] - fi = 0

Ensuite, on forme la moyenne temporelle

) + VùSf{vj) = 0.

Après simplification, on obtient bien 6 nouvelles équations de transport. Toutefois, il y aura

22 nouvelles corrélations inconnues. Le système d'équations reste indéterminé, le nombre
d'inconnues étant supérieur au nombre d'équations. Pour fermer le système d'équations, une
des approches les plus utilisée est basée sur l'hypothèse de Boussinesq.
 35

2.1.1 Hypothèse de Boussinesq

En 1877 Boussinesq a émis l'hypothèse que le tenseur des contraintes turbulentes est
déterminé par les taux de déformation de l'écoulement moyen :

2
-pv' <g> v' = ut (Vw + Vu*) - -pk§_, (2.1)

Ht étant la viscosité turbulente et k = \ ^ vfy l'énergie cinétique des fluctuations turbu-

lentes (énergie turbulente). Le terme -pkô_ permet de tenir compte des contraintes normales.
ô
En effet, la somme des termes de la diagonale principale est égale à :

Uz Uz Z t\> 7 e \J ,
I

Par analogie, les corrélations Tv' sont déterminées à l'aide de la conductivite turbulente Xt
et du gradient de température :
-pcpTv! = XtVT. (2.2)

En substituant les équations (2.1) et (2.2) dans les équations de conservation des valeurs
moyennes de l'écoulement, on a :

d(pv) r _

d(pcPT)
= s.

2
Le terme -pk peut être omis. En effet, le gradient de ce terme peut être absorbé par le
3
_ 2
gradient de pression (voir Chabard[20]). La pression statique est alors remplacée par p+-pk.
o
Ce système d'équations n'est toutefois pas fermé. En effet, les valeurs de la viscosité
turbulente et de la conductivite turbulente sont inconnues. Pour déterminer ces valeurs, un
modèle de turbulence doit être utilisé. Pour ce travail de recherche, on a choisi le modèle
k - e de Launder et Spalding[50]. Il s'agit d'un modèle général qui donne des résultats
physiquement réalistes pour de nombreux cas. Toutefois, ce modèle est déficient pour certains
types d'écoulements (Mohammadi[57], Wilcox[74], Speziale[71], Chang[22]) :

• Pour un écoulement turbulent dans une conduite carrée, le modèle k — e ne permettra

 36

pas de reproduire les zones de recirculation près des coins de la conduite, les variables
k et e prédites par ce modèle étant indépendantes de l'état de rotation du fluide.

• Ce modèle a tendance à surestimer le niveau de turbulence lorsque l'écoulement est

perpendiculaire à un mur. Ainsi, pour tous les problèmes de contournement d'obstacles
et de jets, les niveaux de turbulence prédits par ce modèle seront trop élevés.

• Pour les problèmes à bas nombre de Reynolds, les niveaux de turbulence ne seront pas
prédits adéquatement.
• Lorsque le problème est fortement influencé par un gradient de pression inversé (tel le
problème de l'écoulement au-dessus d'une marche), ce modèle sous-estimera toujours
la longueur de la zone de recirculation.

Pour tous ces problèmes, des variantes du modèle k — e ont été proposées par différents
chercheurs et permettent d'obtenir de meilleurs résultats. Finalement, ce modèle a été utilisé,
sous sa forme originale ou sous d'autres variantes, dans plusieurs codes industriels et par
d'autres chercheurs. Ces faiblesses sont bien connues et pour de nombreux cas, il peut être
adapté au problème considéré. Ces qualités en font un bon modèle de turbulence pour la
validation d'un nouveau schéma numérique.

2.1.2 Équations de bilan du modèle k — e

L'étude détaillée du modèle k - e ne fait pas partie des objectifs de cette thèse, la
dérivation des équations de bilan pour A; et e ne sera donc pas présentée. Pour une présentation
détaillée de la dérivation des équations de ce modèle, on propose le livre de Mohammadi et
Pironneau[57].
Contrairement aux autres équations de bilan, ces équations ne sont pas basées sur la
modélisation d'un phénomène physique. Elle sont obtenues à l'aide de manipulations mathé-
matiques et la fermeture de ces équations est basée sur de nombreuses hypothèses simplifica-
trices. Il est toutefois important de spécifier que cette "modélisation des mathématiques" est
un problème commun à tous les modèles statistiques de turbulence (Wilcox[74]).
 37

Les variables k et e sont déterminées à l'aide des équations de transport suivantes

[73]:

d(pk) , _ f . ( H,
+ V • Ipvk - [n + — VA; = P-pe + G,
dt

= £[C(P.
où:

• P = fj,t \u + y.1 est un terme de production de l'énergie cinétique turbulente par ci-
saillement ;
• G = /3 (Xt/Cp) q-VTest un terme source dû aux forces de gravité, ce terme est présent
dans les problèmes de convection naturelle ;

• CM = 0.09, C^T = 0.09, aK = 1.0, ae = 1.3, ,CeX = 1.44, ,C e 2 = 1.92, ,C e3 =

(G > 0 : 0, G < 1 : 1), sont des constantes qui ont été déterminées expérimentale-
ment.

La viscosité turbulente et la conductive turbulente sont calculées à l'aide des équations :

k2
Ht = pC»—
k2
Xt = pCpCfj,T— •

2.1.3 Conditions aux limites

Lorsque le domaine de calcul est borné par des parois, les conditions limites pour la
vitesse sont en général des conditions de vitesse nulle et de forts gradients de vitesse sont
présents près des parois. Il y a donc une couche limite où les phénomènes visqueux sont
prépondérants et le modèle k — e est inapproprié pour cette région (en général, le modèle
k — € est mal adapté aux régions où le nombre de Reynolds est petit). Pour bien modéliser la
couche limite, il faudrait un très grand nombre d'éléments près de la frontière. Pour éviter la
discrétisation de la couche limite, et ainsi réduire les coûts de calcul, la paroi est modélisée
par son influence sur l'écoulement du fluide : la frontière de calcul est déplacée.
 38

On considère généralement que la couche limite peut être découpée en trois couches
(Chabard[20], Cardot[17]) :

1. une sous-couche visqueuse ;

2. un zone de transition ;

3. une sous-couche externe où l'écoulement turbulent est régi par les grosses structures,
les effets visqueux y étant négligeables.

Dans les applications, la frontière de calcul est déplacée dans la sous-couche externe.
Pour ce faire, on définit les quantités suivantes :

• y la distance fictive entre la paroi réelle et la frontière de calcul ;

• uT la vitesse tangente à la paroi ;

1/2
duT
• u' = ' -• la vitesse de frottement au bord ;
y=0j
Vil
+
• y = — la distance sans dimension entre la paroi physique et la paroi fictive ;
• u+ = / I j la vitesse de frottement sans dimension (/ est une fonction qui sera
précisée plus tard).

Pour la vélocité, les conditions de bord sur une paroi sont les suivantes :

• v_- n = 0;

0" + Vt)-Q^{v-u) = 0;
V
= -r— ;
k\
= n x ri
 39

Pour l'énergie turbulente et les taux de dissipation les conditions de bord sont les suivantes :

1. Entrée:

• énergie et taux de dissipation imposés.

2. Sortie :

• V/c • n = 0 ;

• Ve • n = 0.

3. Parois :

lu*'3
• e = .. , K = 0.41 est la constante de Karman.
VKy

Calcul de la vitesse de frottement

du
Dans la sous-couche visqueuse, la contrainte est proportionnelle à a-—
dy

De plus, le profil est linéaire :

u+ = — = y+ , 0 < y+ < 5.
u*

Dans la zone de transition, u+ vérifie une loi du type

u+ = f(y+) 5<2/+<30.

Finalement, dans la zone externe, la vitesse est déterminée par une loi logarithmique

u+ = — ln(y+) + Ci, 30 < y+ < 100 .

K

Pour déterminer u+, on utilisera la loi de Reichardt :

y+\ y+ ,
~TT ~ TT e x P(-°-
 40

Cette loi permet de raccorder les différentes couches. Lors de la résolution numérique du
système d'équations, la valeurs de la distance fictive entre la paroi réelle à la frontière de
calcul y et la vitesse de frottement u* doivent être déterminées pour chacun des éléments de
frontière. L'algorithme utilisé pour déterminer ces valeurs sera présenté dans la section 4.8.
 Chapitre 3

Définitions et résultats essentiels sur les méthodes de

volumes finis classiques

Introduction

Pour mettre en oeuvre un schéma numérique de résolution d'écoulements incompress-

ibles, de nombreux aspects autres que physique doivent être étudiés : la discrétisation tem-
porelle, la discrétisation du terme de convection, la stabilité, la convergence, le choix des
espaces d'approximation. Dans ce chapitre, tous ces aspects seront discutés en détails.
Pour simplifier la présentation, la densité est considérée constante. De plus, les équa-
tions de conservation sont présentées dans leur forme locale :

- ^ + V • £ = s(---), (3.1)

où:

1. (f> est la quantité scalaire conservée ;

2- E. — [n<i> — aV</>] est le vecteur des flux convectif et diffusif ;

3. v est le vecteur vitesse ;
4. a est un paramètre physique de diffusion telle la viscosité ou la conductivité ;

5. s est un terme source.

 42

Dans ce chapitre, on fera souvent référence à trois définitions fondamentales de l'analyse

numérique :

1. la consistance ;

2. la stabilité;

3. la convergence.

Afin éviter les erreurs d'interprétation, on présente un bref rappel de ces définitions.
Pour simplifier la présentation de ces définitions, on suppose que le problème continu est
indépendant du temps.

Consistance

Une approximation discrète d'un problème continu est dite consistante lorsqu'on peut
montrer que l'erreur de discrétisation tend vers zéro lorsque le pas d'espace h tend vers zéro.

Définition 3.1 (Consistance (au sens des différences finies)) (Andersont et [Link]])
Soit P une équation aux dérivées partielles et Ph le problème discrétisé, nous dirons
que Ph est consistant par rapport à P si :

lim Ph (x) - P (x) = 0, Viefi,,

/i->0

Çl étant le domaine de résolution.

Par la suite, la différence Ph (x^ — P (ZJ) entre le problème continu et le problème

discret sera appelée résidu (noté r j .

Convergence

Un schéma convergent est un schéma pour lequel la solution numérique tend vers la
solution exacte lorsque le pas d'espace tend vers zéro.
 43

Définition 3.2 (Convergence) (Cuvelier et al.[34])

Soit <j) (XJ) la solution exacte à une position x_{ et fa la solution du problème discrétisé
en ce point, le problème discrétisé Ph est dit convergent lorsque

lim \l (x)-£\ = 0.

Par la suite, la différence (j) ( x j — <j>. entre la solution exacte et la solution discrète sera
appelée erreur d'approximation (notée et).
On serait tenté de croire qu'une approximation consistante du problème continu im-
plique la convergence, ce n'est toutefois pas le cas. Pour qu'un schéma soit convergent, il doit
aussi être stable.

Stabilité

Définition 3.3 (Stabilité) (Cuvelier et al. [34])

Soit ri le résidu au point x_{ et ej l'erreur d'approximation en ce point, un schéma est
stable en norme infinie (\- • • l^) lorsque

IÊL $ c Id.» >

C étant une constante indépendante du pas d'espace.

Souvent, la stabilité d'un schéma dépend du pas d'espace ou du pas de temps et la

convergence du schéma vers la solution exacte est alors dépendante du pas de temps ou du
pas d'espace, peu importe si la discrétisation est consistante ou non. On peut aussi construire
des schémas convergents (vers la bonne solution) dont la discrétisation spatiale n'est pas con-
sistante au sens des différences finies avec la forme locale des équations de transport. C'est le
cas des méthodes de volumes finis classiques lorsque le pas d'espace n'est pas constant (voir
[36]). Essentiellement, si on peut démontrer la convergence d'un schéma de volumes finis et
que l'on respecte les contraintes de stabilité, on est théoriquement assuré qu'un raffinement
adéquat du maillage donnera un approximation plus précise de la bonne solution. Finalement,
si la discrétisation est consistante et que les contraintes de stabilité sont repectées, le schéma
est nécessairement convergent [1].
 44

3.1 Discrétisation temporelle

L'équation transitoire (3.1) est résolue en introduisant une série croissante (£ n ) ngN où
en général to = 0, et la différence ôtn+\ = tn+i — tn est appelée pas de temps. La dérivée
temporelle est approximée en introduisant un polynôme de Lagrange de degré m

et en le dérivant par rapport au temps. Pour un schéma à deux niveaux, tel le schéma d'Euler,
l'approximation est la suivante :

~dt t=t +i St
n

Pour un schéma à trois niveaux (tel le schéma de Gear) et un pas de temps constant, la dérivée
temporelle est approximée à l'aide de l'expression

dt t=t 2ôt
n+1

Le choix dépend de la précision et de la stabilité recherchées. Lorsque le membre de droite

de l'équation (3.1) est évalué au temps t = tn+\ (schéma implicite), que la discrétisation
spatiale n'est pas considérée et que le terme source est indépendant de </>, la stabilité de ces
deux méthodes est indépendante du pas de temps. Lorsque le champ convecteur est connu,
le schéma à trois niveaux est consistant d'ordre 2, le schéma à deux niveaux est quant à lui
consistant d'ordre 1 (Pozrikidis[60]).
La discrétisation temporelle peut aussi être basée sur l'intégration en temps des flux
et des termes sources. Dans ce cas-ci, l'équation (3.1) est intégrée sur un pas de temps
générique :
rt=t
rt=tn+1
<t>(x,tn+l)-<f>(x,tn) = /
Jt=tn
 45

Lorsque les flux et le terme source sont interpolés linéairement en temps, on a :

= y { [ - V • F(x, t n + l ) + s ( t n + 1 , • • • ) ] + [ - V • F{x, t n ) + s ( t n , • • • ) ] }

Cette discrétisation temporelle est souvent appelée discrétisation de Crank-Nickolson ou

méthode du trapèze. Ce schéma est consistant et d'ordre deux. Cette discrétisation n'intro-
duit pas de dissipation et elle est plus précise que le schéma de Gear. Puisque la méthode du
trapèze n'introduit pas de dissipation, elle ne peut pas atténuer les petites oscillations et la
stabilité dépend du pas de temps utilisé.

3.2 Discrétisation spatiale

Pour les méthodes de volumes finis classiques, la discrétisation spatiale est obtenue
en introduisant un maillage T^ du domaine Q. Ce maillage est constitué d'un ensemble d'élé-
ments géométriques. Généralement, les bilans sont calculés en intégrant les équations de con-
servation sur ces mêmes éléments. Les éléments géométriques sont alors appelés volumes de
contrôle. En tout temps t = tn, une seule inconnue §\ est associée à un volume de contrôle
K.

3.2.1 Équations discrètes

L'équation discrète associée à une inconnue <%+1 est obtenue en intégrant (3.1) sur un
volume K générique. Lorsque la discrétisation temporelle des flux est implicite, on a :

dt . = -X>-"+m
où:

• — dépend de l'approximation de la dérivée temporelle;

dt ..
t—In+l
• SK est l'ensemble des faces du volume K ;

est l'approximation du flux total à l'interface a au temps t = tn+i ;

 46

• m(-) est la mesure (volume, aire) d'un hyperplan ;

• s^- est une valeur moyenne du terme source par unité de volume.

La construction de l'équation discrète se termine par l'approximation du flux Fa sur chacune

des interfaces1. Cette approximation du flux n'est pas arbitraire. En effet, pour s'assurer de
la convergence d'un schéma vers la solution exacte, on doit construire un approximation du
flux du volume K vers le volume L,

FK,L « / £(£,<)• nKL dS,

JK/L

qui respecte les propriétés suivantes (Eymard et al.[36]) :

1. les flux doivent être conservés :

2. l'approximation doit être consistante :

FK,L -» / E(x, t) • nKL,dS

JK/L
IK/L

lorsque la taille des volumes de contrôle tend vers zéro.

Bien sur, pour qu'il y ait convergence, les conditions de stabilité doivent être respectées.

3.2.2 Maillages admissibles

Pour que la discrétisation des flux d'interface soit consistante, le maillage doit pos-
séder certaines propriétés que voici :

Définition 3.4 Maillage admissible, domaines tridimensionnels(Eymard et al. [36])

Soit Q un domaine polygonal ouvert et borné de H 3 . Un maillage Th admissible pour
Çl est constitué

• d'un ensemble de volumes de contrôle, lesquels sont des polyèdres convexes ouverts
sous-ensemble de H 3 dont le volume n'est pas dégénéré (ou nul),
'Pour alléger la présentation, on omettra le sur-indice de temps.
 47

• d'un ensemble d'hyperplans £ de R 3 (ce sont les faces des volumes de contrôle) dont
l'aire n'est pas dégénérée (ou nulle),
• d'un ensemble de points V.

Ces ensembles doivent satisfaire les contraintes suivantes :

(i) La fermeture de l'union de tous les volumes de contrôle est la fermeture du domaine
(noté H ).
(ii) Soit £ = OKÇT £K Pour tous
volume de contrôle K G fl, il existe un sous-ensemble
£K de £ tel quedK = K\K = (Jff€£jf ô.
(Hi) Pour tous (K, L) e T%, K ^ L, la mesure de Kf^LestO ou bien Kf]L = â. Pour
tous a G £, cette mesure sera notée K\L.
est
(iv) L'ensemble V = XKerh tel Que pour tous K G Th, XK G K.
(v) Si a = K \ L, on assume que XK ^ XL.
(vi) La droite T>K,L qui passe par les points XK et XL est orthogonale à OK,L-

(vii) Pour tous a G d£l, soit T>K,C la droite qui passe par XK est orthogonale à O~K, si
XK ^o alors VK,c Ç\o^QetXa = £K,a f l ° •

Remarques
1. On rappelle que la fermeture d'un ensemble Çl est l'ensemble des points adhérents à Çl
(si Q, = [0,1), alors H = [0,1]).
2. La condition XK G K (pour tous K G Th) peut être relaxée (item (iv)).
Un point associé à un volume de contrôle peut être à l'extérieur de celui-ci.
3. La contrainte a G £ext, ^K,<r f] G ¥" ° P e u t être relaxée (item (vii)).
Sur la frontière, le point d'intersection Xa de la droite passant par XK et orthogonale
à aK peut ne pas être à l'intérieur de a. En effet, cette contrainte n'est un problème
que lorsque la condition de bord est de Dirichlet. Cette contrainte peut être négligée en
imposant 4>K = 9{XK) lorsqu'elle n'est pas respectée, g étant la condition de bord.
4. La contrainte (v) implique que la transmittivité

(XL — XK) • RK,L

 48

Dual, triangulation

de Delaunay

FIG. 3 . 1 - Volume de contrôle pour un maillage Voronoï

de l'interface OK,L est définie pour toutes les interfaces du domaine.

5. Ensemble, les contraintes (iv) et (v) assurent la positivité de la transmittivité.
6. Dans le cas ou la contrainte (v) n'est pas respectée à l'intérieur du domaine, un nou-
veau volume de contrôle peut-être construit en réunissant les deux volumes considérés.
Ainsi, la contrainte (v) sera respectée.
Plus-tard, lorsque des résultats théoriques seront présentés, trois types de maillages seront
considérés :

1. Maillages de type M\ : la quantité TK/L est strictement positive et la contrainte XK e K

est respectée sur tout le domaine.

2. Maillages de type M2 : la quantité TK/L est strictement positive sur tout le domaine, la
contrainte XK G K n'est pas respectée.

3. Maillages de type M3 : la quantité TK/L est négative pour au moins une interface.

En 3D, une façon (qui est la seule à notre connaissance) pour directement construire des mail-
lages admissibles est d'utiliser une technique Voronoï (figure 3.1). Un diagramme Voronoï est
une partition du domaine en cellules. Une cellule K est associée à un point xK (appelé généra-
teur) et est constitué des points plus proches de xK que de n'importe quel autre générateur.
 49

Définition 3.5 (maillage Voronoï) (Eymard et al. [36])

Soit Q, un sous-espace polygonale ouvert borné de R n . Soit V une famille de points
de Cl. Les volumes sont définis par rapport aux positions x_ des points de V :

Kx = {y G Cl | \x - y] < \z - y] , Vz G V, z # x) •

La figure 3.1 illustre un partie d'un maillage 2D Voronoï. Les volumes de contrôle
sont délimités par des lignes continues. Sur cette figure, on présente le dual du diagramme
Voronoï, ce dual est appelé triangulation de Delaunay.

3.2.3 Approximation des flux

Flux convectif

La discrétisation du terme de convection

f (#(£)) • n(x) dS
JK,L

doit être consistante et ne doit pas favoriser l'apparition d'oscillations qui ne sont pas physi-
quement réalistes. Considérons une discrétisation uniforme du domaine Q. de dimension 1. Il
peut être démontré (Cuvelier [34]) que l'approximation du flux convectif par unité de surface
avec la différence centrée
<j>L + <l>K ifv \
v « vcj)(Xa)
conduit à une solution fortement oscillatoire lorsque le nombre de Reynolds local (ou le
nombre Peclet local) est supérieur à deux :

h étant la taille de la maille.

Pour assurer la stabilité du schéma lorsque le nombre de Peclet local est supérieur à
deux, il est essentiel d'utiliser une approximation amont du terme de convection. Considérons
les volumes de contrôle K et L, le flux convectif total du volume K au volume L peut être
 50

approximé avec le schéma amont

/ (vé(x)) • n(x) dS
JK/L

où :

v<T=y_o- nKL est la vitesse normale à l'interface a, y^ est la vélocité à l'interface ;

, si vff > 0
, autrement

Cette expression est une approximation du premier ordre du flux convectif à l'interface. En
effet, soit :

m x = f XK, ^-(XL-XK) >0

° [XL, v, • (XL - XK) <0 '
• x_crK L = Xo + rAx le point d'intersection de la droite joignant les points XL et XK et
de l'interface OK,L,

• (j) une fonction définie et derivable sur l'intervalle [XK, X L ) .

Alors, il existe une position £ comprise entre XK et XL tel que

va (/)(Xo + rAx_) = va 4>(XQ) + va 4>' (Ç(r)) \\Ax}\ .

Ainsi, l'expression
va <j){X0) « va<f>(Xa)

est bien une approximation d'ordre 1 du flux convectif à l'interface.

Flux diffusif

Considérons les volumes de contrôle adjacents K et L d'interface OK,L • Une approx-

imation consistante du flux diffusif total du volume K au volume L est :

f
I aV(f)(xj • n[x) dS.
JK/L
 51

Elle est consistante puisque l'approximation du flux par unité de surface

4>{xL) -

est d'ordre un. En effet, considérons une droite passant par les points XK et XL, cette droite
est décrite par le paramétrage

x(r)=XK + r(XL-XK).

Soit <j) une fonction définie et trois fois derivable sur l'intervalle [XK, XL]. Alors, il existe f
compris entre XK et XL tel que

P-\\X{T)-XK\

La dérivée de cette expression par rapport à r est :

= <t>(XL) - 4>{XK) + R'(r)

où
è (f(r))
K
£(r)
R'(r) = * -^-}\\x{r)~XK\\-\\x{r)-XL\

est la dérivée du reste. De plus, lorsque 0 < r < 1

è'
9
R'(r) =

< ||A£, — A/cll I 2 J ^ 2!

 52

où C < 1 est une constante. Étant donné que

(XL-XK)=nKtL-\\XL-XK\\,

l'expression de l'approximation de la dérivée peut être simplifiée

- XK\\

Dans le cas où l'interface est à mi-chemin entre XK et XL cette approximation est d'ordre
deux et la démonstration est similaire. Toutefois, au lieu de paramétrer avec un polynôme de
degré un, un polynôme de degré deux passant par les points XK , Xa et XL est utilisé.

Important. Dans certains codes de volumes finis classiques, les positions XK et XL re-
posent au centre de gravité de leur volume de contrôle respectif. Souvent, l'approximation

f
JK
K/L

est utilisée pour évaluer les flux d'interface. Cette approximation n'est valide que pour un
maillage rectangulaire et parallèle au système de coordonnées de référence. Pour tout autre
maillage, cette approximation des flux n'est pas consistante puisque la droite reliant les points
XK et Xi n'est pas perpendiculaire à l'interface &K,L Un schéma basé sur cette technique
de discrétisation ne sera pas nécessairement convergent vers la bonne solution. En pratique,
pour un tel schéma, un raffinement adéquat du maillage n'augmentera pas nécessairement la
précision de la solution numérique.

Anti-diffusion et diffusion numérique

On peut obtenir des solutions qui sont physiquement réalistes avec une discrétisation
en amont du terme de convection. Cependant, il y a un problème bien connu inhérent à cette
discrétisation : elle est trop diffusive. Il peut être démontré que la viscosité artificielle in-
vh
troduite par un tel schéma est de l'ordre de —- lorsque la taille h des mailles est constante
(Cuvelier et al[34]). Afin d'atténuer la diffusion introduite par cette discrétisation, différents
schémas antidiffusifs ont été développés. On peut construire ces schémas à partir de la solu-
 53

Schéma (t><r = p
f $x±*Lt Rei < 2 f T, Rei < 2
Hybride
\ (pa,+, autrement \ 0, autrement

Loi de puissance T • max [0, (1 - [Link])5)

P Rei
Exponentiel
eRe> - 1
TAB. 3.1 - Schémas convectifs

tion du problème de convection-diffusion unidimensionnel

— (u(j)-a— j = 0, (j)(x0) = (j>o, XQ<X<XX (3.2)

où u est un champ convecteur connu et constant. La solution à ce problème est donnée par
l'expression

(3.3)
- <t>o
L'étude de cette solution permet la construction des schémas présentés dans le tableau 3.1
(Patankar [59]). Le schéma hybride est le plus simple des trois et permet une approximation
d'ordre deux du terme de convection dans les régions où la convection n'est pas dominante.
Toutefois, il est délicat à mettre en oeuvre, l'approximation du terme de convection avec une
différence centrée pouvant nuire à la stabilité2. Les schémas "loi de puissance" et "exponen-
tiel" donnent essentiellement les mêmes résultats. Ils permettent de conserver la stabilité du
schéma amont et d'atténuer efficacement la diffusion numérique. Toutefois, le schéma "loi
de puissance" est moins coûteux en temps de calcul, il nécessite une seule évaluation de la
fonction puissance comparativement à deux évaluations de la fonction exponentielle pour le
schéma "exponentiel". À notre connaissance, les schémas hybrides et "loi de puissance"
n'ont été utilisés que dans le cadre de méthodes de volumes finis classiques sur mail-
lages structurés. Le schéma "exponentiel" a, quant à lui, été utilisé avec succès avec des
méthodes de volumes finis classiques sur maillages structurés, d'éléments finis et de volumes
finis de type Galerkin sur maillages non structurés triangulaires (Baliga et al. [3], [4] et [5]).
2
En fait, lors du développement du schéma, le schéma hybride n'a pas donner de très bons résultats. Il n'a
donc pas été conservé.
 54

3.2.4 Quelques résultats théoriques

Gallouet et al.([36] et [39]) ont démontré de nombreux résultats théoriques sur la

convergence et la stabilité des méthodes de volumes finis classiques. Dans cette section, on
présente les principaux résultats qui nous ont inspiré pour le développement de notre schéma.
Soit cj> une quantité scalaire obéissant à une loi de conservation, on considère le
schéma de volumes finis suivant :

K,* = m (aK,L) v^lf - TKtLa {(j>nL+1 - ffi1) , a Gint£i

FK,<T - m {a) 1V02J1 - Taa (g (Xff) - ^ + 1 ) , a G £ext, Condition de Dirichlet

FK,O = m (cr) va<Pa^+ TT / â'fe) ^ a e

£exu Condition de Neumann

F<x,,ff = hcm(a) {<j>a - (j)^ , a € £exU Condition de Robin

où:

• — est une approximation discrète de la dérivée temporelle ;

• £K est l'ensemble des interfaces du volume K, £int l'ensemble des interfaces à l'in-
térieur du domaine, £ext l'ensemble des interfaces sur le bord du domaine ;

• va est la vitesse normale à l'interface a ;

• 4%^+ est la valeur en amont de l'interface ;

K
• rKL = j — ' est la transmittivité de l'interface oK L ;
(XL - XK) • nKL
• XL, et XK sont les positions de calcul des flux associées aux volumes de contrôle K
et L, nKL est le vecteur unitaire normal à l'interface OK,L et pointant vers l'extérieur
deK;
 55

• / £ est une valeur moyenne du terme source sur le volume K au temps tn ;

• g(x.) est une fonction qui définit la valeur au bord (condition de Dirichlet), ou le flux
au bord (condition de Neumann) ;

• hc est un coefficient de transfert convectif et (j)^ une valeur de référence.

Pour un maillage admissible, lorsque l'écoulement est permanent et que la condition de bord
est de Dirichlet, de Neumann ou de Robin, Eymard et al. ont démontré que :

1. ce schéma converge vers une solution unique et que l'erreur d'approximation de dépend
que de la taille du maillage3dim (7^) :

où C > 0 est une constante qui est indépendante du pas d'espace.

2. le principe du maximum est respecté (condition de Dirichlet).

Pour une condition de bord de type Dirichlet, lorsque l'écoulement est transitoire et que
la discrétisation temporelle est à deux niveaux (le schéma d'Euler), Eymard et al.[36] ont
démontré que l'erreur est estimée par l'expression

(<f>(XK) - <f>K)2Tn(K) <C(h + St),

où :

• h est la taille du maillage et 6t le pas de temps ;

• C > 0, est une constante indépendante du pas de temps et du pas d'espace.

Remarque Dans certains cas, l'ordre de l'erreur d'approximation observée expérimentale-

ment est souvent supérieur. Dans le cas de ce schéma de volumes finis, pour des problèmes
de diffusion (elliptiques) en 2D lorsque le maillage est non structuré et composé de triangles,
l'ordre de l'erreur d'approximation observée est égale à 2 (voir Cayré [18]).
3
On rappelle que la taille du maillage est le diamètre de la plus grande maille.
 56

3.3 Résolution de l'équation de conservation de la masse, méthodes de

projection

Considérons un écoulement incompressible évolutif. Cet écoulement est déterminé

par les équations de conservation de la quantité de mouvement et de la masse :

dv
-Ï= + X? -[v®v-i;Vv} + VP = / (3.4)
V-v = 0 (3.5)

P étant la pression cinétique. Le système étant fortement couplé et non-linéaire, la résolu-

tion numérique de ces équations est problématique. Bien que ce système d'équations puisse
être résolu en considérant l'ensemble des équations, il est en général préférable de découpler
les équations. Lorsque le problème considéré est instationnaire, le découplage des équations
réduit significativement le temps nécessaire à la résolution (Turek[66]). Toutefois, pour cer-
tains problèmes stationnaires où le découplage n'est pas nécessaire, cette approche peut être
beaucoup plus onéreuse en temps de calcul.

3.3.1 Projection

Les schémas les plus utilisés pour résoudre le système d'équations (3.4-3.5) sont
des méthodes dites de projection. Ce sont les schémas développés à partir de l'algorithme
"Simple" de Spalding et Patankar[59] et du schéma de projection initialement proposé par
Chorin[25].
L'algorithme proposé par Chorin est basé sur un découplage des opérateurs. Le dé-
couplage des opérateurs est une procédure qui permet de séparer un pas de temps en une
suite de pas de temps dans lesquels les opérateurs sont pris en compte séparément. Ainsi, le
système d'équations évolutif

dv . . . .
-5= = Am + A2v + A3 + ... + An
ot

(les Ai... An sont des opérateurs qui modélisent les forces s'exerçant sur le fluide) peut être
 57

résolu en plusieurs étapes :

dv
dt
dv
MV
tt = ~
dv_ _
= An
m -
La dernière solution est celle prévalant pour tout le pas de temps. Dans le schéma proposé
par Chorin, le système d'équations (3.4-3.5) est résolu en deux étapes :
1. Convection-diffusion :
vn+l/2 _ ^n
+ V • [v®v - vVv] = f.
It
2. Projection :
- --VF , v . -w.

La vélocité u n+1 / 2 est une vélocité intermédiaire qui ne satisfait pas nécessairement la con-
trainte d'incompressibilité. De plus, la discrétisation temporelle des termes de diffusion et de
convection dépend du schéma employé. Pour obtenir des résultats réalistes avec ce type de
schéma, le terme de convection doit absolument être discrétisé explicitement4. De plus, les
conditions de bord doivent être modifiées pour tenir compte de la séparation complète des
opérateurs (Kim et al. [46], Leveque et al.[54]). Pour éviter la limitation du pas de temps et
stabiliser le schéma, la discrétisation suivante est préférable :

1. Convection-diffusion :

— i Y7 , F,, (+ + \ <> ,,"+1/2 ,.V7..n+1/2] , r r r>n /o £\

4
Lorsque le terme de convection est discrétisé implicitement, on peut démontrer que la solution finale d'un
écoulement permanent est fortement liée à la vélocité et au pas temps.
 58

-• /3 = 2, V P n + 2

-• p = §, V F n + 2

ypn+l

VFn

tn+1 tn+2

FIG. 3.2 -Extrapolation du gradient de pression

2. Projection :

vn+l _ vn+l/2
- V (SPn+1), V (SPn+1) = V P n + 1 (3.7)
P Tt
V • t; n+1 = 0 (3.8)

où j3 est un paramètre de relaxation. Il permet de faire une interpolation ou une extrapolation

de la valeur du gradient de pression au pas de temps suivant. Trois cas peuvent être considérés
(figure 3.2) :

1. P < 1, interpolation, le gradient de pression est une combinaison linéaire des gradients
de pression aux temps t = tn et t = tn+i.

2. p = 1, le gradient de pression est celui du temps précédent (t = tn+i) ;

3. P > 1, le gradient de pression est une extrapolation (sur-relaxation).

Théoriquement (Shen[67]), le choix P = 2 conduit à une discrétisation temporelle

d'ordre deux. Toutefois, dans la pratique, la valeur optimale peut dépendre du problème con-
sidéré.
Pour respecter l'équation de conservation de la masse. On procède ainsi : vn+1 est
 59

(2,0) (0,0) (2,0) (0,0) (2,0)

> O <» O i> O <• 0 <> o «
(1.0) (1,0) (1,0) (1,0) (1,0) (1,0)

FlG. 3.3 - Illustration du problème de l'échiquier pour la vitesse

isolée dans l'équation (3.7) :

vn+1 = -jV {6Pn+1) +vn+1'2 (3.9)

En substituant cette expression dans l'équation de continuité on obtient une équation où seule
la variation de pression est inconnue :

(3.10)
St

La pression est modifiée en fonction de la solution à cette équation et la vitesse est ensuite
mis à jour en appliquant l'équation (3.9).

3.3.2 Maillages décalés

Un autre problème associé à la résolution des équations de Navier-Stokes pour un

écoulement incompressible est lié aux espaces d'approximation de la vitesse et de la pres-
sion. Si ces espaces d'approximation ne sont pas choisis judicieusement, la pression ne sera
pas définie de façon unique à une constante près. Dans certains cas, la solution numérique
présentera d'importantes oscillations qui ne sont pas physiquement réalistes (faux modes de
pression). Ce problème est causé par les dérivées premières

VP et Vv

qui sont présentes dans les équations de conservation de la quantité de mouvement et de

la masse. Si les espaces d'approximation ne sont pas choisis adéquatement, il est possible
de construire une fonction P (ou v) non constante telle que V P = 0 (ou Vv. = 0). Cette
construction est illustrée sur la figure 3.3. Sur cette figure, bien que la composante horizontale
de la vélocité varie entre 0 et 2 sur tous les éléments, la valeur interpolée aux interfaces est
 60

FIG. 3.4 - Maillage décalé structuré

FlG. 3.5 - Élément d'un maillage non structuré "décalé"

constante et égale à 1. L'approximation de la dérivée pour cet exemple est nulle sur chacun
des éléments, bien que la fonction ne soit pas constante.
Pour résoudre ce problème, on peut utiliser des maillages différents pour la pression
et la vélocité. Ces maillages sont souvent appelés maillages décalés. Pour ces maillages, les
degrés de liberté pour la vélocité et la pression reposent en des positions différentes. Pour un
maillage structuré, la construction est illustrée sur la figure 3.4. Les degrés de liberté pour la
pression reposent au centre des éléments. Pour la vélocité, les degrés de liberté sont au centre
des interfaces. Ainsi, le volume de contrôle pour les composantes de la vélocité est "décalé"
par rapport à celui de la pression.
Pour un maillage non structuré constitué de triangles ou de tétraèdres, Boivin et al[13]
ont introduit une méthode pour éviter les faux modes de pression. Des espaces d'approxima-
tion différents pour la pression et la vitesse sont introduits en interpolant la vélocité sur les
interfaces du maillage. Ces valeurs sont ensuite utilisées pour effectuer la projection. Boivin
 61

et al. ont démontré que la pression calculée après les étapes de ré-interpolation et de projec-
tion est définie de façon unique à une constante près. Un élément pour ce maillage "décalé"
est illustré sur la figure 3.5. On tient à préciser que la démonstration de Boivin et al. s'applique
au schéma qui est présenté dans cette thèse.
 Chapitre 4

Une nouvelle méthode de volumes finis pour des maillages

non structurés

Notre nouveau schéma numérique est inspiré d'une méthode de volumes finis de ré-
solution d'écoulements incompressibles 2D sur maillages non structurés proposé par Boivin
et al.([13], [9], [10], et [19]). Pour ce schéma, les volumes de contrôle sont les triangles qui
forment le maillage (les éléments géométriques). Les positions de référence des volumes de
contrôle ne sont pas au barycentre des cellules, mais placées à l'intersection des bissecteurs
orthogonaux. Notre schéma est quant à lui plus général, le domaine peut être tridimensionnel
et les volumes de contrôle ne sont pas nécessairement les éléments géométriques mais des
assemblages d'éléments.
Les principales caractéristiques de notre schéma sont :

1. la méthode de discrétisation spatiale est de type volumes finis classiques ;

2. ce schéma permet l'utilisation de maillages non structurés généraux ;

3. pour toutes les variables du système d'équations les fonctions d'interpolation sont des
constantes par volume de contrôle ;

4. la discrétisation temporelle est semi-implicite à pas fractionnaires.

Bien que le schéma puisse être utilisé lorsque la densité est variable, celle-ci est considérée
constante pour l'instant.
 63

4.1 Discrétisation temporelle

La discrétisation temporelle est semi-implicite et basée sur la méthode de projection

décrite dans la section 3.3.1. Soit,

-z , cas stationnaire
—. àt
dt
autrement
l'approximation discrète de la dérivée temporelle et

. . I -vn vn~l vélocité

v(t = tn+1)= < 2~ 2"
I vn+1 autrement
l'approximation du champs convecteur au temps t = tn+i, la discrétisation temporelle des
systèmes d'équations est la suivante :

1. Extrapolation du champs de vitesse :

dv
+ V • [v (t = tn+1) ® vn+1/2 - uVvn+1/2} + V F " = /" . (4.1)
~dt t=tn+i

2. Projection :

(4.2)

V-t; n + 1 = 0. (4.3)

3. Convection et diffusion des autres variables scalaires telle la température :

-$• + V • \v (tn+i) </>n+1 - aV(j)n+1} = sn . (4.4)

at t 1

Remarque. Le problème est linéarisé à cette étape, le champs convecteur étant considéré
connu pour le temps t = tn+i. De plus, lorsque l'on ne tient pas compte de l'erreur intro-
duite par la discrétisation spatiale, l'utilisation du schéma de Gear pour l'approximation de
la dérivée temporelle nous donne un schéma qui est d'ordre 2 en temps pour la vélocité et les
variables scalaires (Shen[55]).
 64

FIG. 4.1 - Position des LD.O.s

4.2 Discrétisation spatiale

4.2.1 Fonctions d'approximation

Toutes les variables du système d'équations sont approximées avec des fonctions con-
stantes par volume de contrôle. Il n'y a qu'un seul degré de liberté (ddl) par volume contrôle.

4.2.2 Eléments géométriques et volumes de contrôle

Éléments géométriques

Un ensemble d'éléments géométriques est utilisé pour discrétiser le domaine Q. Cet

ensemble d'éléments est constitué exclusivement de tétraèdres. Une position de référence
XK est associée à tout tétraèdre K de la discrétisation. Cette position est le centre de la
sphère circonscrite au tétraèdre K. Ainsi, la droite reliant les positions XK et XL de deux
tétraèdres adjacents est toujours perpendiculaire à l'interface OK,L (Eymard[36]). Par la suite,
cette position sera appelée intersection des droites orthogonales (I.D.O.). La position de
l'I.D.O. est illustrée sur lafigure4.1.
 65

X,

FIG. 4.2 - Configurations des I.D.O.s

 66

Volumes de contrôle

Les volumes de contrôle sont des assemblages d'éléments (de tétraèdres) sur lesquels
on réalise les bilans. À priori, les volumes de contrôle correspondent aux tétraèdres. Cepen-
dant, dans certains cas (nous verrons pourquoi plus-tard), des volumes de contrôle sont con-
stitués en réunissant plusieurs tétraèdres. Même pour ces cas, le nombre de degrés de liberté
associés à un volume de contrôle demeure égale à 1, toutes les variables étant constantes par
volume de contrôle.
La construction des volumes de contrôle est déterminée par la valeur de la transmit-
tivité de l'interface aKtL séparant deux tétraèdres K et L :

TKIL =
(XL - XK) • nK,L '

Les possibilités sont les suivantes (ces configurations sont les mêmes en 2D, elles sont illus-
trées sur la figure 4.2 de la page 65 pour deux triangles voisins K et L) :

1. Configurations M\ et M2 :
TK,L > 0.

Les volumes de contrôle sont les éléments géométriques.

2. Configuration M 3 :
TK,L < 0,

XK et XL ne reposent pas dans leurs tétraèdres associés et sont inversés. Les éléments
K et L sont combinés en un seul volume de contrôle. Ce volume de contrôle possède
au moins deux I.D.O : XK et XL.

3. Configuration M4 :
\TK,L\ - K » ,

au moins deux I.D.O. occupent la même position (à epsilon près). Les éléments K et
L sont combinés en un seul volume de contrôle. Ce volume de contrôle ne possède
qu'un seul I.D.O. (tous les I.D.O. de l'ensemble des tétraèdres qui forment le volume
occupent la même position).
 67

1
X
la« (X,,

2 /

/ 3 •

\
X4b
\
4
X•4a .
\

\
Maillage initial et Volumes de contrôle et
points de calcul des flux points de calcul des flux

FiG. 4.3 - Illustration du passage des éléments géométriques aux volumes de contrôle

Exemple de construction des volumes de contrôle

La figure 4.3 illustre le passage des éléments géométriques aux volumes de contrôle.
En premier lieu le maillage est constitué de 6 éléments. À chacun des éléments, on a associé
une position de calcul des flux (un I.D.O.). Les points de calcul des flux des éléments 1 et 2
occupent la même position. Ces éléments sont combinés pour former le volume 1. Pour les
éléments 3 et 4, aucun traitement n'est nécessaire, les éléments deviennent les volumes de
contrôle 2 et 3. Pour les éléments 5 et 6, il y a inversion, ces éléments sont combinés pour
former le volume 4. Ce volume possède deux positions de calcul des flux, les points Xj a et
Xa,. Le calcul du flux entre les volumes 3 et 4 est effectué à partir des points Xia et X3.

Remarque. Bien que des positions de calcul des flux puissent être à l'extérieur de leur
volume de contrôle, c'est la solution associée à tout le volume de contrôle qui est utilisée
pour calculer les flux.
 68

FIG. 4.4 - Problème d'interpolation à l'interface

4.3 Approximation discrète sur les interfaces

Les variables et paramètres du système d'équations sont interpolées sur les interfaces à
l'aide d'une moyenne géométrique. Soit <j) une variable ou un paramètre telle la conductivité,
l'approximation discrète sur l'interface GK,L est :

m(K)(j)K + m(L)(j)L
<t>K,L = m{K) + m(L)

L'approximation géométrique permet de respecter localement le minimum et le maximum de

la solution. Mais, elle est moins précise que l'interpolation linéaire

4>K,L = (1 - tK,L)<i>K

, _ {XK - XKL - XK)

lK L
' " (xL-xK)-(xL-xK) •
En effet, cette dernière interpolation préserve l'ordre d'approximation aux positions XK et
XL lorsque 0 < tx,L < 1- Malheureusement, en 3D, il est fréquent que tKtL » 1 ou ÎK,L <C 0
(figure 4.4). En ce qui concerne l'interpolation linéaire, le minimum et le maximum de la
solution approximative ne sont pas préservés localement. Lors du développement du schéma,
on a constaté que l'interpolation linéaire de la solution aux interfaces pouvait faire diverger
le schéma : l'approximation géométrique a donc été privilégiée.
 69

4.4 Calcul des flux d'interface

4.4.1 Opérateur de convection C(<j)o,+)

L'approximation discrète du flux convectif entre les volumes K et L est calculée à

l'aide du schéma amont :

f
/ {voix)) ' Q(ZÙ dS ~ i
JK/L
où:

va = v_ (t = tn+i) • nKL est une approximation de la vitesse normale à l'interface a au

temps t = tn+\ ;

, autrement

4.4.2 Opérateur de diffusion D(<j)K, <I>L)

Considérons les volumes de contrôle adjacents K et L d'interface OK<L . Le flux dif-

fusif total du volume K au volume L est approxime à l'aide de l'expression

/ aV(j){x) n(x_) dS « aa TK,L {4>L ~ 4>K) ,

JK/L

où:

• aa est le coefficient de diffusivité à l'interface a ;

• T~K,L est la transmit!vitté de l'interface OK,L-

Calcul du coefficient de diffusivité

Pour un cas générique, le coefficient de diffusivité A est un tenseur d'ordre 2 qui

dépend de la direction et de la position. Pour avoir une bonne approximation du flux entre les
point XK et Xi, le coefficient de diffusivité est approxime comme suit :
 70

1. Calcul de la valeur moyenne pour chacun des volumes de contrôle :

AK = — ^ f A(x) dV = A(XK).

2. Calcul de la valeur absolue de la projection dans la direction normale à l'interface :

OLK,a = \AK -Ha] •

3. Interpolation à l'aide d'une moyenne géométrique sur l'interface ;

m(K)aKtC + m(L)aLi(T

( . .

m(K) + m(L)
T— O t Cint

Ot-K,a O e £ext
Anti-diffusion
On utilise une dérivée amont pour discrétiser le terme de convection. On rappelle
que cette dérivée introduit implicitement de la diffusion numérique. Pour les problèmes qui
demandent une plus grande précision (tels ceux où l'on veut calculer les coefficients de por-
tance et de traînée), le coefficient de diffusion moléculaire est ajusté avec le schéma "loi de
puissance" :

Remarque. Cet ajustement du coefficient de diffusion moléculaire permet de conserver

les propriétés du schéma. Le schéma reste consistant puisque h —> 0 => Rei -> 0 et le
coefficient de diffusion tend vers sa valeur réelle. De plus, la stabilité du schéma est conservée
puisque le terme de diffusion numérique introduit par la discrétisation amont du terme de
convection n'est pas modifié.

4.5 Opérateur d'approximation du gradient Grad{<j))

A priori, puisque nos fonctions d'approximation sont constantes par volume de con-
trôle, on ne peut pas calculer explicitement le gradient d'une fonction sur un volume. Pour
 71

calculer un gradient, on procède avec l'approximation suivante.

Le gradient Vçi> d'une fonction scalaire cj) sur un volume de contrôle K est calculé à
partir du gradient normal aux interfaces, ce dernier étant connu. On suppose que la projection
du gradient pour le volume de contrôle K sur l'interface OK,L est égale au gradient normal à
l'interface :

(4.5)
(4.6)

En appliquant l'équation (4.5) à toutes les interfaces de la frontière du volume K, on obtient

le système d'équations

o\

où, de façon plus compacte

En général, ce système d'équations est incompatible. La "meilleure" solution pour ce système

d'équations est calculée avec la méthode des moindres carrés : c'est le système d'équations

N*NV<I>K = A^'

qui est résolu.

 72

4.6 Équations discrètes, opérateur de convection-diffusion CD(- • • )

Pour tout volume de contrôle K, les équations de convection-diffusion peuvent être

écrites sous la forme :

f dA dV +
JK dt t=tn+l

dK
•ndA = [ (4.7)
K

L'évaluation des intégrales de surface à l'aide des opérateurs de convection et de diffusion

ainsi que l'évaluation des intégrales volumiques nous donne :

t=tn+i

Ho) va (t = tn+1) C ï - ct <j>nK+l)} = m{K)sl (4.8)

Les termes sources pour les différentes équations sont les suivants1 :

1. Composantes du vecteur vitesse :

2. Énergie cinétique turbulente :

m(K)sK = m(K)
* J

3. Énergie de dissipation turbulente :

m(K)st = m(K)% C,

'Bien sur, ces termes sources peuvent être modifiés en fonction du problème considéré.
 73

4.7 Projection

Après l'étape de convection-diffusion, l'équation de conservation de la masse n'est

pas satisfaite. Lorsque la densité est constante, la vélocité n'est pas nécessairement à diver-
gence nulle. Une étape de projection est donc nécessaire. Celle-ci se déroule en deux étapes :

1. Une étape d'extension (de ré-interpolation) sur les interfaces qui permet d'éviter les
faux modes de pression.
2. Une étape de projection pour corriger la pression et la vélocité.

O p é r a t e u r d ' e x t e n s i o n E(- • • ) .

Pour pouvoir procéder à la projection de la vélocité dans l'espace des fonctions solé-
noïdales, la vélocité doit être correctement définie sur les interfaces. Puisque le nombre de
faces pour chacun des volumes de contrôle est variable, cette vélocité d'interface ne peut être
définie à l'aide de la même fonction d'interpolation pour tous les volumes de contrôle. La
vélocité d'interface ne peut donc pas être calculée explicitement.
La vitesse d'interface est calculée implicitement. Au même titre que la vélocité sur les
volumes de contrôle, l'évolution de la vélocité sur les interfaces doit satisfaire les équations
de Navier-Stokes. Soit v_KL la vélocité à l'interface OK,L-, son évolution du temps t = tn au
temps t = tn+i/2 doit satisfaire l'équation

dy_
= - V • [v {t = tn+1) <g> vn+l'2 - vVvn+1'2} + V P n + f1 (4.9)
~dt

En intégrant le membre de gauche de l'équation (4.9) sur les volumes de contrôle associés à
l'interface OK,L, on a :

dv_ m(K) dy_ m(L) dv_

dt a, t=tn+i m{K) + m(L) dt K,*=Un+1 L,t=tn+1

OU
dv
dv
dt L,t=t
dt
désignent les approximations discrètes de la dérivée temporelle sur les volumes K et L. La
vélocité d'interface est calculée à partir de l'expression
 74

= E(vK,vL)

n m(K)+m(L)
( n+1/2
\V-K

pour un écoulement permanent. Lorsque l'écoulement est transitoire, la vélocité d'interface

est déterminée par l'équation suivante :

= E(vK,vL)

^ m(K)+m(L)

Opérateur de projection P(- • • ).

L'opérateur de projection est utilisé pour calculer un champ convecteur à divergence

nulle et évaluer le champ de pression. Cet opérateur est défini à partir des équations

-^ = -V{5Pn+l) (4.10)

V.y_n+1 = 0 (4.11)

On rappelle que fi est un paramètre qui permet de faire de l'interpolation ou de l'extrapolation

(section 3.3.1).
L'opérateur de projection est appliqué en quatre étapes :

1. Correction du champ de pression.

2. Correction du champ convecteur.

3. Mise-à-jour de la vélocité sur les volumes de contrôle.

4. Mise-à-jour du gradient de pression sur les volumes de contrôle.

 75

Correction de la pression. Pour obtenir une équation dont la seule inconnue est la correc-
tion de pression, on considère l'équation (4.10) que l'on substitue dans l'équation (4.11) :

Ensuite, cette équation est intégrée sur un volume générique K et le théorème de Gauss est
appliqué pour obtenir l'intégrale de surface

/ V (ôPn+1) -ndS= [ ^vn+xl2 -ndS (4.12)

dK JôK à*

L'opérateur de diffusion permet d'obtenir l'équation discrète

rff (ÔPL - 5PK) = £ ^m(a)

Deux types de conditions de bord peuvent être imposées :

1. vitesse normale v^ imposée :

V (ôPn+1) • n = Y (vn+1/2 - vb) ;

2. pression imposée :
ÔP = Pn+1 - Pn .

Correction du champ convecteur. La résolution de l'équation (4.12) sur tout le domaine

permet de projeter le champ de vitesse dans l'espace des fonctions solénoïdales. Après cette
résolution, la vitesse d'interface est corrigée à l'aide de l'équation (4.10) :

St rn

Mise-à-jour de la vélocité sur les volumes de contrôle. Lorsque la vélocité est un champ
solénoïdal et que les volumes de contrôle sont des tétraèdres, la vélocité sur les volumes de
contrôle doit satisfaire l'équation

,.n+l . _ n+1
 76

sur toutes les interfaces (Cayre [19]). Toutefois, puisque dans notre cas les volumes de con-
trôle ne sont pas nécessairement des tétraèdres, cette équation ne peut être satisfaite pour
toutes les interfaces. Pour tous les volumes de contrôle, c'est l'équation (4.10) qui est utilisée.
La correction de vélocité est calculée en résolvant le système d'équations

-wA (yôPn+1) 0"!

ou :
• SvnK+1 = vn+1 - vn+1/2 ;

• -64 (V6Pn+1) = - ^ - y - r (SPL - ÔPK) = 6vZ±\ , est la correction de vitesse

sur l'interface OK,L \

CQO, paroi imperméable

1 autrement
• ÉQO est un grand nombre.

Ce système d'équations est généralement sur contraint (incompatible) et la solution est cal-
culée à l'aide d'une minimisation par la technique des moindres carrés.
Le paramètre Wa est utilisé pour s'assurer que pour tous les volumes de contrôle
adjacents aux parois, la vélocité au bord du domaine respecte la physique du problème. Les
faux puits ou les fausses sources de masse près des parois imperméables sont ainsi évités.

Mise-à-jour du gradient de pression sur les volumes de contrôle. Le gradient de pression

est approximé à l'aide de l'opérateur d'approximation discrète du gradient (page 70) :

VP = Grad(P).

4.8 Discrétisation des lois de paroi pour les écoulements turbulents

Dans la section 2.1.3 (page 37) les conditions de bord pour les parois imperméables
ont été présentées. On rappelle que ces conditions de bord modelisent l'influence de la paroi
sur l'écoulement. Ces conditions de bord sont les suivantes :
 77

V. ' UL = 0)

(A* + fh) fa-(v-l.2) = 0,

= T-r, r 2 = n X T! ,

• e= \u*f
.. , ii' = 0.41 est la constante de Karman,
vKy
• u*est la vitesse de frottement, elle satisfait l'équation

/( y +) = fai = 2.5 In (1 + 0Ay+)

+ 7.8 ^1 - exp ( - ^ Q - ^ e x p (-0.33y+)) (4.13)

y+ = —

L'équation (4.13) doit être résolue sur toutes les faces qui discrétisent les parois imper-
méables. Pour résoudre cette équation, on emploie un algorithme qui a été utilisé avec succès
par Cardot[17]. Cet algorithme est basé sur une méthode de sur-relaxation et il est présenté à
la page 78.
 78

Soit a l'interface du volume K sur la frontière du domaine, y = dfcj(T/10 la distance initiale

entre la paroi physique et la frontière de calcul et vn = v_n • T_\ la vitesse tangente à la paroi.
ÏVTI1
1. Poser «* = -
f
Pourk = 1, 2 , . . .
o Pouri = 1, 2 . . .5

Poser u* =

„. lut5 — u*A
4I
o Si ' < e

• alors u* = ut, et fin.

• sinon, poser

MQ = (1 — 9) • u[ -
+ 6(1 - 9)292 • u* + 4(1 - 9)93 • u\ + 94 • u*5

et recommencer, 9 étant un paramètre de relaxation dont la valeur est

posée à 0.5.
, y • u*
2. A convergence, poser y"1" = ,

• y+ > 11.63, y+ appartient à la sous-couche logarithmique et la solution est

valide.
• si y+ < 11.63, poser y = —'• et recommencer la résolution pour u*, aller en
u*
1.

FlG. 4.5 - Algorithme de calcul de la vitesse de frottement

 79

Intégrales de bord pour les parois imperméables

En 3D, la vélocité est égale à la somme de trois vecteurs orthogonaux :

v. — (y. • n j i i + (v • 1.2)12 + {v. • R)R >

où Tj, r 2 sont des vecteurs tangents à la paroi et n la normale à la paroi. Puisque

V. • R — V. • Z2 = 0 ,

l'égalité suivante est vérifiée :

Ô r. . . . / \ 1 ^ ((^ ' Z l ) l l i )
V • R = / ^ e ^ l(î> • ZiJlIi + (V • I 2 J I 2 + (^ • ^ ) ^ J = A*e ^ = -p

Lorsque l'écoulement est turbulent, la dérivée normale à la paroi est calculée avec cette
expression.
Pour les parois imperméables, les intégrales de bord sont les suivantes :

1. quantité de mouvement : / /ieV;u • ndA = —m(a)p(u*)2Zi ',

Jcr

2. énergie turbulente : / ( \x H—- I VA; • ndA = ra ( n -\—- 1 I A;^ —(u*f

3. dissipation turbulente : / I u + — ) Ve • ndA = ra [ u-\—-

 80

4.9 Algorithme de résolution

En résumé, lorsque la densité est constante, la résolution d'un problème d'écoulement

se déroule comme suit :

• Étant donnée une solution initiale vT1, P"1, 4>° :

o Appliquer l'opérateur de projection : (v°, P°, VP°) = P (•••)•

• Étant donnée la solution y_n, pn, 4>n :

o pour un écoulement turbulent, calculer les coefficients de diffusion turbulente et

la vitesse de frottement sur la paroi ;

o appliquer l'opérateur de convection-diffusion à la vélocité : vn+1/2 = CD (•••)»

o appliquer l'opérateur de projection : (vn+1, Pn+1, VPn+1) = P (•••);

o appliquer l'opérateur de convection-diffusion aux champs scalaires telle l'énergie

cinétique turbulente et la température : 4>n+l — CD{- • • ).

4.10 Remarques et résultats théoriques

Notre schéma permet la résolution d'écoulements incompressibles 3D dans des geo-

metries complexes. Les bilans sont calculés sur des volumes de contrôle arbitraires. Lorsque
les volumes de contrôle ne possèdent qu'une seule position de calcul des flux, les résultats
théoriques présentés dans la section 3.2.4 s'appliquent à notre schéma :

1. le schéma est convergent et l'erreur d'approximation est O(h), h étant le diamètre de

la plus grande maille du maillage ;

2. lorsqu'il s'applique, le principe du maximum est respecté pour toutes les variables
scalaires.
 81

4.10.1 Maillage

On rappelle que pour démontrer la convergence du schéma, le maillage doit être "ad-
missible" :

1. une seule position de référence est associée à chacun des volumes de contrôle ;

2. la droite qui relie les positions de calcul des flux de deux volumes de contrôles adjacents
est perpendiculaire à l'interface qui les sépare ;

3. les transmittivités de toutes les interfaces du maillage doivent être positives.

On rappelle que notre domaine de calcul est en premier lieu discrétisé avec des tétraè-
dres. Ensuite, nos volumes de contrôle sont construits en suivant les conditions présentées
dans la section 4.2.2 à la page 66. Dans tous les cas, notre approximation des flux convectifs
et diffusifs entre deux volumes de contrôle est consistante. Toutefois les résultats théoriques
pourraient ne pas s'appliquer dans les situations suivantes :

1. un volume de contrôle possède plusieurs positions de calcul des flux diffusifs qui n'oc-
cupent pas la même position ;

2. un volume de contrôle chevauche deux milieux hétérogènes ;

3. un point de calcul des flux est à l'extérieur du domaine.

Lorsqu'un volume de contrôle a plusieurs positions de calcul des flux qui occupent des po-
sitions différentes, la valeur de l'inconnue associée à ce volume est supposée égale en au
moins deux positions différentes. Localement, l'ordre d'approximation des flux peut ne pas
être suffisant et on ne peut pas évaluer l'erreur à convergence. Dans Eymard[36], il est toute-
fois mentionné que le schéma peut continuer à converger vers la bonne solution lorsque

m
52 (a) „ -,
m
(a)
où a E S est l'ensemble des interfaces où l'approximation du flux diffusif n'est pas consis-
tante et a E Cl est l'ensemble des interfaces du domaine.
Lorsque les éléments qui séparent deux milieux hétérogènes A et B (une phase liquide
et une phase solide par exemple) sont séparés par une interface à transmittivité négative, ils
 82

Liquide

Solide

FIG. 4.6 - Volume de contrôle chevauchant deux milieux hétérogènes

sont combinés en un seul volume de contrôle qui chevauche deux milieux dont les propriétés
physiques ne sont pas les mêmes (figure 4.6).
Lorsqu'une interface de bord à une transmittivité négative et que la condition de bord
est de Dirichlet ou de Robin, on doit imposer la solution au volume associé à l'interface. En
effet, si on ne procède pas ainsi, le schéma peut diverger (lors d'expérimentations, c'est ce
qui a été observé). Pour une condition de Dirichlet, on impose la valeur à l'interface et, pour
la condition de Robin, la valeur à l'infini. Une solution physiquement réaliste peut alors être
obtenue. Néanmoins, localement, on ne peut pas évaluer l'erreur à convergence. On tient à
préciser que de très bons résultats ont été obtenus avec ce traitement (le problème de l'écoule-
ment engendré par le déplacement d'une paroi du chapitre 7 par exemple).

4.10.2 Modèle de turbulence

Pour la mise-en-oeuvre du schéma on a utilisé le modèle de turbulence k — t. On

tient à préciser que notre schéma n'est aucunement dépendant de ce modèle de turbulence.
Pour des problèmes où le modèle k — t n'est pas adéquat, il serait possible d'utiliser d'autres
modèles de turbulence.
 Chapitre 5
Résolution numérique des systèmes d'équations linéaires

Introduction

La discrétisation temporelle priviligiée par notre schéma implique la résolution d'un

ensemble de systèmes d'équations linéraires1

Ax = b

où:

• A est la matrice des coefficients ;

• x est le vecteur des inconnues ;
• b est le vecteur des termes constants.

Dans ce chapitre, on présente les algorithmes qui ont été mis en oeuvre avec notre schéma
de résolution d'écoulements incompressibles. De plus, un rappel sur la classification des sys-
tèmes d'équations linéaires et le conditionnement de ces systèmes est présenté. Les références
pour ce chapitre sont les suivantes :

1. caractérisation de la matrice des coefficients : Atkinson et al.[2], Ciarlet[27] et Green-

baum[41] ;
2. méthodes directes de résolution : Atkinson et al.[2], Cuvelier[34] et Fortin[38] ;
3. méthodes itératives de résolution : Cuvelier[34] et Greenbaum[41] ;
4. algorithmes de préconditionnement : Ciarlet[27] et Greenbaum[41].
'Afin d'alléger la notation pour ce chapitre, les matrices seront simplement notées à l'aide d'une lettre
majuscule.
 84

5.1 Caractérisation de la matrice des coefficients et détermination des

valeurs propres

5.1.1 Caractérisation de la matrice des coefficients

II est bien connu que le comportement des méthodes de résolution numérique des sys-
tèmes d'équations linéaires dépendent du type de matrice considérée. Ces matrices peuvent
être qualifiées comme suit :

1. matrice symétrique,
aij = aji, l<i,j<n;

2. matrice symétrique à diagonale strictement dominante,

Kl >

3. matrice symétrique à diagonale faiblement dominante,

4. matrice symétrique définie positive,

x*Ax > 0, Va; 6 R n , xl ï (0,0,... 0) ;

5. matrice symétrique semi-définie positive,

xlAx > 0, V£ G R n , x* # (0,0,... 0).

Une matrice symétrique à diagonale positive dominante est toujours définie positive. De plus,
une matrice symétrique à diagonale positive faiblement dominante est au moins semi-définie
positive et peut être définie positive.
 85

5.1.2 Détermination des valeurs propres

Les méthodes de résolution numérique des systèmes d'équations linéaires sont forte-
ment dépendantes des valeurs propres de la matrice des coefficients. Malheureusement, dans
la pratique, il est aussi coûteux de calculer ces valeurs propres que de résoudre le système
d'équations. On s'intéressera plutôt au domaine de définition des valeurs propres.
Pour toutes les matrices carrées, le théorème de Gerschgôrin permet de connaître le
domaine de définition des valeurs propres en fonction des coefficients de la matrice. Soit A
une matrice carrée et
n
Ri{A) = ^ \ciij\, 0 < i < n,

il peut être démontré que toutes les valeurs propres de A sont situées à l'intérieur de l'union
des disques
n
\J{ze:C\z-aii\<Ri(A)} , 0 < i < n.

5.1.3 Matrices des coefficients associées à notre schéma

Pour notre schéma, la matrice des coefficients pour les équations de convection-
diffusion est une matrice non-symétrique à diagonale positive dominante. La matrice pour
la correction de pression est, quant à elle, symétrique à diagonale positive faiblement domi-
nante.
Dans tous les cas, la matrice des coefficients est une matrice creuse (une matrice dont
le nombre de coefficients nuls est largement supérieur au nombre de coefficients non-nuls).
En effet, pour tout volume K, seules les inconnues associées aux volumes adjacents à K
pourront générer des termes non-nuls dans la matrice.
 86

Équations discrètes de convection-diffusion

On rappelle que toutes les équations discrètes de convection-diffusion peuvent être

écrites sous la forme

UV
t=tn+l
+l n 1
53 [m{a)va{t = tn+l)^-aara{(j)l -4> + )} = m{K)s% (5.1)

ou :
1. £K est l'ensemble des interfaces de la frontière du volume K ;
2
- T,veeK m{a)va = 0 ;

v
3.«s=! **• ; - \ •
autrement
Puisque les flux sont toujours conservés, l'inéquation

5 3 a°T° ^ 5 Z a°T<J
est toujours satisfaite, £int étant l'ensemble des interfaces qui ne sont pas sur la frontière du
domaine. Pour les termes de convection, l'inéquation

rn(a)va> JZ M*K| (5.3)

0 O€£K,

est elle aussi toujours satisfaite :

m(a)vff = 0 =>

Puisque les inéquations (5.2) et (5.3) sont satisfaites et que le terme transitoire est toujours
présent, la diagonale principale est positive et strictement dominante. La matrice est donc
toujours inversible pour toute valeur de ôt. De plus, la matrice inverse est à coefficient tous
positifs. Cette propriétés garantit le principe de positivité discrète pour <\>\ en tout temps et
sur tout le domaine.
 87

Soit A une matrice carrée inversible d'ordre n.

1. Construction de la matrice triangulaire supérieure U
Pouri = 1,.. .n — 1
• Pour j — i + 1,.. .n
o Poser mji = a^/aa.
o Pour k = i + 1 , . . . n.
• Poser ajk = ajk - mji • aik.
o Poser bj = bj — m J2 • 6t-.
2. Résolution, une remontée triangulaire
Poser xn = bn/ann. Pour i = n — l,n — 2,..., 1

• Poser Xi = (bi - J2"=i+i aHxi) faa-

FiG. 5.1 - Algorithme de Gauss

Pour nos systèmes d'équations de convection-diffusion, le théorème de Gershgôrin

permet d'affirmer que :

1. la partie réelle des valeurs propres sera toujours positive ;

2. une valeur propre nulle n'est pas une valeur propre admissible.

Matrice de correction de pression

Pour la matrice de correction de pression, seul les termes de diffusion sont présents.
L'inéquation (5.2) est satisfaite et la matrice est symétrique à diagonale faiblement domi-
nante. Cette matrice est au moins semi-définie positive et peut être définie positive.

5.2 Méthodes directes de résolution des équations linéaires

Les méthodes directes de résolution numérique des systèmes d'équations linéaires

permettent de calculer directement la solution du système d'équations. Parmi ces méthodes,
on retrouve les algorithmes de Gauss, la factorisation LU et la factorisation de Choleski
(figures 5.1 à 5.3).
 88

5.2.1 Algorithme de Gauss

L'algorithme de Gauss consiste à ramener la matrice A à une matrice U triangulaire

supérieure et de procéder à la résolution du système linéaire équivalent. Lorsque A est in-
versible la matrice U peut toujours être calculée. Pour obtenir ce résultat, il est toutefois
nécessaire d'adopter une stratégie de traitement des pivots nuls :

1. pivot partiel : lors de l'élimination de la colonne A;, si akk = 0, on permute la ligne k

avec la ligne i(i > k) dont l'élément \a,ik\ est le plus grand ;

2. pivot total : lors de l'élimination de la colonne k, si akk = 0, on choisit l'élément

\a,ij\ (i > k), le plus grand, et on permute la colonne j avec la colonne k (ce qui
modifie l'ordre des inconnues), puis on permute la ligne k avec la ligne i.

De plus, pour éviter les erreurs d'arrondi, une stratégie de permutation des lignes doit être
définie. Par exemple, avant l'élimination de la colonne k, la ligne k sera permutée avec la ligne
i (i > k) dont l'élément |a,-fc| est le plus grand. Toutefois, pour une matrice A à diagonale
strictement dominante, aucune stratégie de pivot n'est nécessaire. Pour ces matrices, il n'y a
pas d'accroissement des erreurs d'arrondi ni de pivots nuls.
Cet algorithme nécessite le stockage de tous les éléments de la matrice des coeffi-
2
cients . De plus le nombre d'opérations arithmétiques nécessaire pour toute la résolution est :

4n 3 + 9n2 - In
6 '

5.2.2 Factorisation LU

Pour cette méthode, la matrice A est factorisée (décomposée) en un produit d'une

matrice triangulaire inférieure L et une matrice triangulaire supérieure U :

A = LU.
2
Pour l'algorithme de Gauss, la factorisation LU et la décomposition de Choleski, le coût du stockage et le
nombre d'opérations arithémiques peuvent être significativement réduits en réduisant la largeur de bande de la
matrice.
 89

Soit A une matrice carrée inversible d'ordre n.

1. Factorisation LU, décomposition de Crout
• Première colonne de L :
o Pour i = 1 , . . . n , poser lu = an
• Première ligne de U :
o Pour i — 2 , . . . n , poser un = au
• Pouri = 2 , . . .n — 1

o Calcul du pivot : lu — au —

o Pour j = i + 1,... n.

• Calcul de la ie colonne de L : Iji = aji — Y^

• Calcul de la ie ligne de U : Uij = I a^ — ^ kkUkj ) /h

• Poser lnn = ann - J2lZi inkUkn-

2. Descente et remontée triangulaires
• Descente triangulaire pour résoudre Ly_ = b
o Poser y i = bi/ln.
l
( ^
o Pouri = 2,...n,yi= I bt - 2_JkkVk
\ k=\

• Remontée triangulaire pour résoudre Ux. = y

u ruser xn — yn.
n
o Pour i = n — 1, n — 2,... ,1, Xi = yi —
k-i+l

FIG. 5.2 - Factorisation LC7, décomposition de Crout

 90

La solution du système Ax = b est ensuite obtenue en résolvant les systèmes

Ly = b et Ux_ = y.

Les solutions de ces deux systèmes sont directement calculées à l'aide d'une descente et
d'une remontée triangulaire.
Puisque les matrices L et U ne sont pas uniquement définies, la diagonale principale
de la matrice L (décomposition de Doolittle) ou de la matrice U (décomposition de Crout)
est imposée égale à 1. Tout comme pour la méthode de Gauss, il faut adopter une stratégie
de traitement des pivots nuls. Dans ce cas-ci, la stratégie est plus complexe, l'historique des
permutations devant être conservé. Cependant, lorsque l'algorithme de Gauss ne nécessite
pas de pivots (lorsque A est à diagonale strictement dominante), la matrice A pourra toujours
être factorisée sans permutation.

Tout comme pour l'algorithme de Gauss, cet algorithme nécessite l'entreposage de

tous les coefficients (d'une matrice pleine). Le nombre d'opérations arithmétiques nécessaire
est
4n 3 - 3n2 -n , ,,
# pOur j a decomposition,
0
• 2n2 — n pour la remontée et la descente triangulaires.

4^ | Qn yn
Ainsi, opérations arithmétiques sont nécessaire pour décomposer et résoudre
6
un système, le même que pour une élimination de Gauss. Cet algorithme est très avantageux
lorsque plusieurs systèmes d'équations ont la même matrice des coefficients. En effet, la
résolution de tous les systèmes d'équations ne nécessite qu'une factorisation et pour chacun
des systèmes, une remontée et une descente triangulaires sont requises.
5.2.3 Factorisation de Choleski

Lorsque A est symétrique définie positive, il peut être démontré qu'il existe une ma-
trice triangulaire inférieure C telle que
 91

Soit A une matrice carrée symétrique inversible d'ordre n.

1. Factorisation de Choleski CC% — A(C est une matrice triangulaire inférieure).
• Poser en = y/âîï
• Première colonne de C :
o Pour i = 2 , . . . n , poser Cu = au/eu
• Pour k = 2 , . . . n
o Calcul du pivot :
A-l
c
Ckk = Ckk - T. kjCki
\

o Pouri = k + 1,... n.
• Calcul de la ie ligne de C :

Cik = ( CLik - 2_^ CijCkj ) /Ckk

3=1

2. Descente et remontée triangulaires

• Descente triangulaire pour résoudre Cy = b
o Poser y i = bijc\\.
o Pouri = 2,.. .n,

\ k=i
l
• Remontée triangulaire pour résoudre C x_ = y
o Poser xn = yn.
o Pouri = n — 1, n — 2 , . . . , 1,

k=i+\

FiG. 5.3 - Factorisation de Choleski

 92

Cette factorisation est appelée factorisation (ou décomposition) de Choleski. Elle peut être
obtenue de la factorisation LU en imposant lu = un. Cette factorisation a l'avantage d'exiger
un moins grand nombre d'opérations arithmétiques :

2n3 + 15n2 + n
6

(à peu près la moitié des opérations requises pour la méthode de Gauss et la factorisation
LU).

5.2.4 Application des méthodes directes

Bien qu'on puisse toujours calculer la solution avec ces algorithmes, ils ne sont pas
adaptés à nos systèmes d'équations. A priori, lorsqu'aucun algorithme de réduction de la
largeur de bande n'est utilisé, tous ces algorithmes nécessitent le stockage de tous les coeffi-
cients de la matrice et un très grand nombre d'opérations arithmétiques.
En format double précision, le stockage de tous les coefficients d'une matrice requiert
8 • n 2 mots mémoire. Par exemple, pour un petit système avec 10000 variables, l'espace
nécessaire au stockage de la matrice est de 800 mégaoctets.
Pour une station de travail équipé d'un microprocesseur Pentium III possédant une
vitesse d'horloge de 733Mhz et capable d'effectuer 60 millions d'opérations arithmétiques
par secondes (60 mégaflops), le temps nécessaire pour résoudre un seul système d'équations
ou pour faire une factorisation LU est de 3 heures. Ce temps passe à 1 heure 30 minutes pour
une factorisation de Choleski. Puisque la simulation d'un problème d'écoulement 3D requiert
la résolution d'au moins 4 systèmes d'équations par pas de temps, il est indiscutable que ces
méthodes ne sont pas adaptées à notre schéma.

5.3 Méthodes itératives de résolution numérique des systèmes d'équa-

tions linéaires

Les méthodes itératives permettent de résoudre les systèmes d'équations linéaires à

l'aide d'une série d'approximations de la solution. À chaque itération, une nouvelle approxi-
 93

mation de la solution

est calculée, c^ étant une matrice ou un scalaire etp, la direction de descente. Les algorithmes
les plus simples sont :

1. la méthode de Jacobi : a^ = D~l où D est la matrice des éléments de la diagonale

principale de A et pk = b — Ax_k ;
2. la méthode de Gauss-Siedel : a* = L'x où L est la matrice des éléments de la partie
triangulaire inférieure de A et pk = b — Ax_k >
3. la méthode de relaxation (SOR) : ak = w~ l D — L où u est un paramètre de relaxation,
D est la matrice des éléments de la diagonale principale, L est la matrice des éléments
de la partie triangulaire inférieure de A (excluant la diagonale principale) et

Malheureusement, pour tous ces algorithmes, la convergence est linéaire. Il n'y a donc pas
d'accélération de la convergence lorsque le résidu r^ = b — Ax_k diminue. Pour ces méthodes,
le coefficient (ou la matrice) <XK et la direction de descente^sont statiques. Aucun effort n'est
fait pour les optimiser. Des méthodes itératives plus performantes peuvent être construites en
procédant à des modifications appropriées du coefficient a# et de la direction de descente
p . Ces méthodes itératives sont basées sur la minimisation de l'erreur \\x_ — a^H (x. étant la
solution) ou du résidu ||7^. — 2jt_i||- La prochaine direction de descente est calculée à l'aide
d'une minimisation sur un espace vectoriel

Ce sous-espace de dimension k — i +1 est appelé espace de Krilov. Les vecteurs qui forment
la base de cet espace sont orthogonaux. En général, les algorithmes qui convergent le plus
rapidement sont ceux dont l'espace de Krilov est le plus grand. On présentera trois méthodes
itératives de ce genre :

1. le gradient conjugué ;
2. Orthomin2;
3. GMRES.
 94

Soit A une matrice carrée inversible d'ordre n, x^ une approximation de la solution, calculer

LQ — b — AXQ et poser po = r^.

Pourk = 1, 2 . . . ,
£fc 1
• Calculer ak-\ = ~ ' f*"A
k-x • lk~i
Poser xk = xk_x + ak-iP_k_v Lk = I f c -i - aife_i

Calculer bk = ~fc ' Lk .

• Poser p^ = r^. + 6*-i

FlG. 5.4 - Algorithme du Gradient conjugué

5.3.1 Algorithme du gradient conjugué

Lorsque la matrice des coefficients est symétrique définie positive, la norme ^Ae^. de
l'erreur à l'itération k peut être minimisée sur tout l'espace

V = ÇQ + span {AÇQ, A2e0,...,Ak~1e0} = e0 + j ^ , p ^ . . . ,p fc _ l

Pour ce faire, on choisit

comme prochaine direction de descente. L'algorithme du gradient conjuque est présenté sur
la figure 5.4 (page 94).

Remrques.
1. La norme du résidu \\r_k\\ n'est pas monotone décroissante.

2. La convergence est assurée en au plus N itérations pour r € HN.

 95

Soit A une matrice carrée inversible d'ordre n, x$ une approximation de la solution, calculer

£o = b — AXJQ et poser p0 = ro-

Pourk = 1, 2 . . . ,
r_k_x • Ap
• Calculer ak-\ = — ~fc~1 .

• Poserx_k = Xu_, + ak-ip ,,Lk— L k - \ ~ ak-iAp .

—K — 1 —/C — 1

Ar_k • Ap
• Calculer bk = — .

• Poserp^ = Ifc - h-iPf.^-

FlG. 5.5 - Algorithme Orthomin2

5.3.2 Algorithme Orthomin2

Soit A une matrice carrée inversible, lorsque

Ar_k • Ap
p, . = rk — T=—p.
-K+i "• j±p . Ap -k
—A:—1 — k — 1
est choisi comme prochaine direction de descente, la norme du résidu peut-être minimisée
sur l'espace

• V = io + s P a n {A1JQI -42Io>• • • > AkLo] = Io + span jAp^, Apv ..., Apk_x \ lorsque
A est symétrique,

• v = 2^_i + span | virA;-1, AP_k_l} autrement.

L'algorithme Orthomin2 est présenté sur la figure 5.5. Il peut être démontré (Greenbaum[41])
que cet algorithme ne peut pas échouer lorsque toutes les valeurs propres sont différentes de
zéro.
On rencontre moins souvent cet algorithme car ces itérations sont plus coûteuses que
celles du gradient conjugué. Cet algorithme requiert à la fois un produit matrice vecteur ainsi
qu'un produit scalaire supplémentaires à chacune des itérations. Il est toutefois plus poly-
valent que l'algorithme précédent, la matrice des coefficients n'ayant pas à être symétrique
 96

Soit q., j = 1 . . .nun ensemble de vecteurs unitaires de dimension net A une matrice carrée

d'ordre n. Pour j = 1 , 2 . . . ,

• Pouri = 1 , . . . , j

o Poser hij = q~j • qit qj+1 = qj+1 - h^q.

• Poser hj+hj = \\qj+i\\, qj+1 = gj+1/hj+1<j.

FiG. 5.6 - Algorithme d'Arnoldi

définie positive.

5.3.3 Algorithme GMRES

Pour des matrices non-symétriques, on peut vouloir minimiser la norme sur un espace
plus grand en construisant explicitement la base de l'espace de Krilov. L'algorithme GMRES
de Saad et Schultz[64] (figure 5.7 ) permet d'augmenter la taille de l'espace de Krilov pour
les matrices qui ne sont pas symétriques. La base de l'espace est construite à l'aide d'une
procédure pour orthogonaliser tel l'algorithme d'Arnoldi (figure 5.6). Pour cet algorithme,
le résidu est minimisé sur tout l'espace

V = u + s P a n {^Io> ^2Io> • • • » ^*Io} = Io + s P a n {AP& AEv • • • » AEk-i

Pour qu'il soit efficace, l'algorithme GMRES doit converger assez rapidement. En
effet, la construction de la base de Krilov est assez coûteuse et demande le stockage d'un
vecteur pour chacune des itérations de GMRES. Pour diminuer les coûts liés à la construction
de la base de l'espace de Krilov, l'algorithme GMRES est redémarré après un nombre m
prédéterminé d'itérations. Le nombre optimal d'itérations ne peut toutefois pas être évalué et
il dépend du problème considéré. Ce nombre ne peut être déterminé que par la pratique.
Cet algorithme est le seul des algorithmes présentés dans cette section qui puisse
 97

Soit XQ une approximation de la solution, calculer r^ = b — Ax_0 et poser qx = TQ/ ||ro||.

Pourk = 1, 2 . . . ,

• Calculer qUl,ethik, i = l,k + len utilisant l'algorithme d'Arnoldi. Entreposer ces

vecteurs dans la matrice Q.

• Calculer Xj. = Xn + Qy , où yh minimise ||r o || ex — H y . H étant la matrice des

coefficients de l'algorithme d'Arnoldi et ex le vecteur unitaire ( 1 , 0 , . . . ) .

FlG. 5.7 - Algorithme GMRES

efficacement résoudre les problèmes où la diagonale principale n'est pas au moins faiblement
dominante.

5.4 Conditionnement

5.4.1 Notions de base et opérateurs de conditionnement

Pour tous les algorithmes présentés dans la section précédente, le nombre d'itérations
nécessaire pour converger ne dépend pas seulement de la dimension de l'espace de Krilov
mais aussi du conditionnement de la matrice des coefficients. Le conditionnement d'une
matrice mesure la sensibilité de la solution du système d'équations

Ax = b

vis-à-vis des variations des coefficients de A et b. Une matrice (ou le système linéaire associé)
est bien conditionnée lorsque cette sensibilité est petite. Soit A une matrice inversible, le
conditionnement de A est mesuré à l'aide du produit

= \\A\\i\\A-1\\i>l,
 98

où j | . . -H; est une norme matricielle :

1. \\A\\i =

2. ||i4||2 = y/p{AlA), p{AtA) étant la plus grande valeur propre de A*A ;

3. ||A||iB/ = m a x ^ | a y | .
3

Plus ce nombre est petit, meilleur est le conditionnement de la matrice. Pour une ma-
trice normale (donc symétrique), le nombre de conditionnement cond2(^4) peut être déterminé
à l'aide du rapport entre la plus grande et la plus petite de ses valeurs propres :

Malheureusement, dans la pratique, on ne connaît pas les valeurs propres de la matrice, ni

son inverse. A priori, il est donc difficile de connaître le conditionnement d'une matrice.
Il est toutefois possible d'améliorer le conditionnement d'un système d'équations en
appliquant des opérateurs de préconditionnement. Symboliquement, ces opérateurs sont ap-
pliqués comme suit :
1 l
y = M;lb, y = Mdx

Mg et Md étant des matrices d'ordre n tels que M~lAM^1 « / , / étant la matrice identitée.
La matrice Mg est un préconditionneur à gauche et Md un préconditionneur à droite. Soit E
une matrice de correction et D la diagonale principale de A, les méthodes de précondition-
nement les plus connues sont :

1. la méthode de factorisation incomplète ILU : Mg = L, Md = Û, LU = A - Ë

2. la méthode de factorisation de Choleski incomplète IC (A doit être symétrique) : Mg =

C, Md = C\ CCi = A-É;

3. la méthode de Jacobi :

(a) matrice symétrique : Mg = y/~D, Md =

(b) matrice non-symétrique :Mg — D, Md — I.

 99

Remarque. Lors de la mise-en-oeuvre, les matrices de la factorisation ne sont pas inver-

sées. Pour tous les algorithmes itératifs (Orthomin2, gradient conjugué, GMRES...) plusieurs
étapes sont ajoutées :

1. Avant de démarrer la méthode itérative :

(a) on effectue le changement de variable x^ = M^XQ , x^ étant la solution initiale ;

(b) le système d'équations Mgb = b est résolu.

2. Tous les produits matrice vecteur Av_ sont remplacés par :

(a) la résolution du système M^ = v;

(b) le calcul du produit f>2 = Avx ;
(c) la résolution du système Mgv = v2.

3. À la fin de la résolution, on effectue un changement de variables en résolvant le système

d'équations Mdx_ = y, y étant la solution de l'algorithme itératif.

On rappel que les matrices Mg et Md sont des matrices triangulaires, la résolution des sys-
tèmes d'équations à l'intérieur des algorithmes itératifs se résume au calcul explicite de la
solution.
En ce qui concerne les méthodes de factorisation incomplète, les algorithmes des fig-
ures 5.2 et 5.3 sont utilisés avec des stratégies de remplissage. La stratégie la plus simple
consiste à ne conserver que les coefficients associés à des entrées non-nulles de la matrice
A (factorisations ILU(0) et 7(7(0)). Malheureusement, les résultats obtenus avec ce type de
factorisation sont quelque peu aléatoires :

1. L'apparition de pivots nuls peut obliger l'utilisation d'une stratégie de traitement des
pivots nuls (même si théoriquement, la matrice des coefficients pourrait être décom-
posée sans stratégie de pivot).
2. La matrice du produit MgM^ peut dans bien des cas être significativement différente
de A. Alors, le système d'équations

M;lAMîly = M;lb, y = Mdx

peut être moins bien conditionné que le système initial et plus difficile à résoudre.
 100

Pour tous les opérateurs de préconditionnement basés sur une factorisation incomplète de la
matrice des coefficients, les problèmes sont essentiellement les mêmes :

1. le contrôle de l'espace mémoire utilisé pour la factorisation ;

2. la mise-en-oeuvre difficile des algorithmes les plus performants ;

3. il est difficile de prédire si le préconditionnement sera réellement efficace.

Une revue de la littérature sur ce sujet est présenté par Saad[65].

5.4.2 Équation de Poisson

Considérons l'équation de Poisson

aAcj) = s

où a est un paramètre qui n'est pas nécessairement constant, <j> une variable quelconque et s
un terme source. Les propriétés des systèmes d'équations discrètes résultant de la discrétisa-
tion par des méthodes de différences finies ont été étudiées en détail (Greenbaum[41]). Pour
les problèmes où

• le domaine est carré (ou cubique) de longueur unitaire,

• le maillage est structuré et le pas d'espace constant,

• a est constant,

la matrice des coefficients est constante le long des diagonales, symétrique et tridiagonale par
bloc (TST par bloc). Les valeurs propres pour une telle matrice sont connues explicitement.
Il peut être démontré que pour cette catégorie de matrices le nombre de conditionnement est

cond(A) = 1

h étant le pas d'espace de la discrétisation (Greenbaum[41]). Ainsi, le conditionnement de

la matrice est inversement proportionnel au pas d'espace. Lorsque le pas d'espace n'est pas
constant ou lorsque le maillage n'est pas structuré, on peut présumer que le conditionnement
est similaire.
 101

On termine cette section en précisant qu'il existe des résoluteurs et des précondion-
neurs spécialement adaptés à ce problème et d'autres plus généraux (Greenbaum[41]) :

1. matrice TST par bloc : résoluteurs directs d'équations de Poisson ("fast Poisson solver-
s'');
2. matrice symétrique définie positive : méthode SOR avec paramètre w optimal déter-
miné par la plus grande valeur propre (j3 < 1)

u =

3. problèmes généraux : méthodes multi-grilles.

5.4.3 Conditionnement utilisé dans notre logiciel

Pour notre logiciel, on a opté pour la technique de conditionnement de Jacobi. Cette

méthode ne nécessite le stockage que d'un seul vecteur (la diagonale principale de A ) et
sa mise-en-oeuvre est très simple. Pour toutes les résolutions de systèmes d'équations, les
étapes suivantes sont ajoutées à l'algorithme itératif de résolution :

1. matrice symétrique :

(a) avant la résolution :

i. mise à l'échelle des colonnes : â^ = a^jy/âj], 1 <i <n, 1 < j < n, ;

ii. mise à l'échelle des lignes : ô^ = â^/^/âû, 1 < j <n, 1 < i < n ;

iii. mise à l'échelle des variables :yi = Xi- y/â^, 1< i < n;

(b) après la résolution :

• mise à l'échelle des variables :x{ = yi/^/âH, 1 < i < n.

2. matrice non-symétrique :

• Avant la résolution, mise à l'échelle des lignes :

âij = ay/an, 1 < j <n, 1 <i <n.

 102

5.5 Résultats comparatifs et choix de l'algorithme de résolution

Pour comparer la performance des algorithmes de résolution et mesurer l'efficacité

de notre préconditionneur, un problème d'écoulement thermique en régime permanent3 a
été résolu avec les algorithmes du gradient conjugué, Orthomin2 et GMRES. L'algorithme
GMRES a été testé avec deux configurations :

1. 30 itérations et 5 pas (4 redémarrages) ;

2. 50 itérations et 3 pas (2 redémarrages).

Critères de convergence et d'arrêt

Soit k le nombre d'itérations et rj. le résidu à l'itération k, les critères de convergence

et d'arrêt pour chacun des algorithmes sont les suivants :

1. équations de la conservation de la quantité de mouvement :

(a) orthomin2 :

||£/t||L2 < 1 x lu"7 ou ||l A IL 2 /||lolL 2 <lxl0- 3 ou fc>40;

(b) GMRES :

| | r j i 2 < 1 x HT7 ou |tall L 2 /||rol| L 2 <lxlO- 3 ou A;>30;

2. équations de correction de pression :

(a) gradient conjugué et Orthomin2 :

U r J ^ < 1 x 10" 7 ou ||lfc||L2/j|r0||i2<lxl0-3 ou A;>300;

(b) GMRES :

\\rk\\L2 < 1 x lu" 7 ou | | r j | i 2 / ||ro|| L2 < 1 x 10" 3 .

3
Ce problème est présenté dans le chapitre 7 à la page 146.
 103

Algorithme Temps CPU Pas de temps Norme infinie

(secondes) pour convergence du bilan massique
à convergence
Avec
préconditionnement
Gradient conjugué 492.85 59 3.04 • 10- 9
et Orthomin2
Orthomin2 548.09 53 1.15-10- 9
GMRES(3Ox5) 1000.55 70* 4.42 • 10~7
GMRES(50x3) 1254.79 70* 3.64 • 10" 7
Sans
préconditionnement
Gradient conjugué 571.33 70* 3.33 • 10" 7
et Orthomin2
Orthomin2 703.16 70* 2.32 • 10- 8
GMRES(30x5) 963.72 70* 1.85 • 10~7
GMRES(5Ox3) 1203.28 70* 1.48 • 10- 7
* indique que le critère de convergence \wk+1 — wk\oo < l.e — 4 n'a pas été satisfait après 70
pas de temps.

TAB. 5.1 - Temps de résolution pour différents algorithmes de résolution, problème de l'é-
coulement 3D thermique dans une conduite

Pour toutes les simulations, le critère d'arrêt est

n étant le nombre de pas de temps.

Résultats

Pour tous les algorithmes, deux résolutions ont été effectuées : une avec précondition-
nement et une autre sans préconditionnement. Le tableau 5.1 présente le temps CPU néces-
saire pour la résolution, le nombre de pas de temps pour satisfaire le critère de convergence
et la norme infinie du bilan massique à convergence. Cette dernière information nous permet
déjuger de la qualité de la solution à convergence.
Les résultats montrent que le préconditionnement a facilité la résolution des systèmes
d'équations. Lorsqu'il n'y a pas de conditionnement, aucun des algorithmes n'a permis de
 104

résoudre tout le problème d'écoulement en moins de 70 pas de temps.

Avec et sans conditionnement, c'est l'algorithme Orthomin2 qui a permis d'obtenir
"la meilleure" solution (dans les deux cas, la norme du bilan massique est inférieure à celle
des autres algorithmes). De plus, c'est avec cet algorithme que la solution a été obtenue le
plus rapidement (en nombre de pas de temps).
En ce qui concerne la vitesse d'exécution, c'est l'association gradient conjugué et
Orthomin2 qui est la moins coûteuse en temps CPU (10% à 20%).
L'algorithme GMRES a été le moins efficace des trois algorithmes. Les temps de
résolution sont significativement plus élevés et l'algorithme n'a permis dans aucun cas de
satisfaire le critère de convergence pour tout le problème.
En ce qui concerne l'espace mémoire utilisé par ces algorithmes, on peut mentionner
qu'aucun des algorithmes ne nécessite la construction explicite du système d'équations. Les
vecteurs nécessaires à la résolution sont au nombre de 3 pour le gradient conjugué, 4 pour
Orthomin2 et de la dimension de l'espace de Krilov plus 2 pour GMRES.
C'est l'algorithme Orthomin2 avec préconditionnement qui a été utilisé pour résoudre
tous les systèmes d'équations linéaires des problèmes présentés dans le chapitre 7,
 Chapitre 6
Mise-en-oeuvre

Introduction

L'objectif de ce chapitre est de présenter une librairie informatique qui a été dévelop-
pée dans le cadre de cette thèse de doctorat. Un première librairie de routines a été mise-en-
oeuvre pour faciliter le développement de notre schéma numérique. Bien que les routines de
cette première librairie soient efficaces, elles ne sont pas adéquates pour le développement du
schéma par d'autres chercheurs. Suite aux résultats obtenus avec notre schéma, on a décidé
de mettre en oeuvre une librairie informatique qui soit plus appropriée à l'exploitation du
schéma par d'autres chercheurs.
Notre librairie permet la résolution de problèmes 2D et 3D en régime permanent ou
transitoire. Les écoulements peuvent être laminaires ou turbulents. Pour les écoulements tur-
bulents, seul le modèle de turbulence k — e a été mis en oeuvre. Il est toutefois possible
d'ajouter d'autres modèles de turbulence. Des problèmes avec densité constante ou variable
peuvent être résolus avec cette librairie. De plus, le nombre de variables scalaires transportées
par l'écoulement est variable. Ainsi, avec cette librairie, des écoulements thermiques avec
plusieurs composants chimiques peuvent être résolus.

6.1 Langage de programmation et programmation orientée objet

La librairie informatique a été mise en oeuvre avec le langage de programmation C++.

Ce choix est en partie justifié parce que ce langage supporte la programmation procédurale
(ou programmation structurée) et la programmation orientée objet ("POO"). La pro-
grammation procédurale consiste à résoudre un problème en construisant un ensemble de
 106

fonctions (les algorithmes) et par la suite de déterminer la méthode de stockage des don-
nées. Les langages de programmation Fortran, C, Pascal sont des langages de programmation
structurée. La POO inverse cet ordre. En premier lieu, on cherche à connaître les données
nécessaires pour la représentation de notre problème. Ensuite, on associe des fonctions pour
la manipulation de ces données. En POO, on combine les données (les attributs) et les fonc-
tions (les méthodes) pour former des objets. Les principaux langages POO sont le JAVA, le
Smalltalk et le C++. Le C++ est le seul de ces langages qui supporte à la fois une approche
procédurale et une approche POO. Les langages Smalltalk et JAVA sont des langages pure-
ment POO : toutes les données devant être représentées par un objet (le programme principal
est lui-même un objet). Deux concepts sont fondamentaux à la POO : 1'encapsulation et
l'héritage.
L'encapsulation (ou isolement des données) consiste à regrouper les attributs et in-
terfaces (les interfaces sont des méthodes) entre l'objet et l'utilisateur dans un même con-
tenant et de dissimuler l'implémentation de ces interfaces aux utilisateurs des objets. Lorsque
l'encapsulation est respectée, on ne peut modifier ou accéder aux attributs d'un objet qu'en
passant par une interface visible à l'utilisateur. Pour les gros projets, le respect de l'encap-
sulatation permet de plus facilement déboguer les programmes : les attributs d'un objet ne
pouvant être modifés que par un nombre déterminé de fonctions.

L'héritage permet d'étendre la fonctionnalité d'un objet A à un (ou plusieurs) objet B.

L'objet B possède toute la fonctionnalité de l'objet A (méthodes et attributs) et de nouvelles
propriétés. Ainsi, avec l'héritage, on peut réutiliser des méthodes déjà existantes avec de
nouveaux objets plus évolués.
Finalement, en C++, il est possible de créer des fonctions ou des objets génériques
qui sont conçus pour manipuler des données hétérogènes. On dit que ces fonctions ou objets
sont paramétrés. En C++ les fonctions ou objets sont paramétrés à l'aide de template.

6.2 Description des objets et de la librairie

Dans cette section, on présente une description des objets construits pour supporter la
librairie. On tient à préciser qu'en C++ les objets sont construits via les classes.
 107

6.2.1 Classe VecteurBase

C'est la classe principale utilisée pour l'allocation dynamique de vecteurs d'objets

(Boivin[8]). La classe est assez générale et est destinée a être spécialisée par héritage. Ses
principales caractéristiques sont :

• elle est entièrement paramétrée (template) ;

• les opérations sont optimisées pour les processeurs RISC ;

• l'accès au coefficient i est effectué via l'opérateur [i] ;

• l'accès aux éléments peut être fait avec validation des indices.

6.2.2 Classe Vecteur

C'est une classe pour la manipulation de vecteurs numériques, elle est dérivée de la
classe VecteurBase (Boivin[8]). Ses principales caractéristiques sont :

• elle est paramétrée (template) ;

• les opérations sont optimisées pour les processeurs RISC ;

• l'accès au coefficient / est effectué via l'opérateur fi] ;

• elle permet d'effectuer les principales opérations arithmétiques vectorielles : produit

scalaire, addition (soustraction), multiplication (division) par un scalaire, recherche de
la valeur maximale (minimale), calcul de la moyenne, calcule de normes discrètes,
multiplication vectorielle pour les vecteurs de dimension 2 et 3.

6.2.3 Classe Matrice

Les matrices pleines sont manipulées avec la classe Matrice (Boivin[8]). Ses princi-
pales caractéristiques sont :

• elle est entièrement paramétrée (template) ;

• les opérations sont optimisées pour les processeurs RISC ;

• les données de la matrice sont stockées via un objet Vecteur ;
 108

• l'accès au coefficient sur la ligne i et la colonne j de la matrice est effectué via l'opéra-
teur (i,j) ;
• elle permet d'effectuer les principales opérations mathématiques associées aux matri-
ces numériques pleines : le produit matriciel, l'addition (la soustraction) matricielle, le
produit matrice-vecteur, le produit matrice-scalaire ;
• le calcul de déterminants pour les matrices carrées d'ordre 4 ou moins ;
• la résolution directe par la méthode de Cramer de matrices carrées d'ordre 4 ou moins ;
• la résolution directe par la méthode de Gauss.

Les 3 dernières caractéristiques ont été ajoutées par l'auteur de cette thèse.

6.2.4 Classe Element

C'est une classe générique pour la représentation d'éléments géométriques ou de vol-
umes de contrôle . Un Element possède une mesure (son aire, volume) et un Vecteur d'entiers
pour l'entreposage de données liées à la connectivité du maillage (les sommets pour les élé-
ments, les interfaces pour les volumes de contrôle).

6.2.5 Classe Face

C'est une classe qui est utilisée pour la représentation des faces d'un maillage, elle est
dérivée de la classe Element. Elle est spécialisée à l'aide

• d'un Vecteur de nombres réels pour représenter le vecteur normal,

• de deux entiers pour identifier les numéros des éléments qui lui sont associés.

6.2.6 Classe Interface

Cette classe est utilisée pour la représentation des interfaces entre les volumes de
contrôle, elle est dérivée de la classe Face. Elle est spécialisée à l'aide

• des numéros de référence des positions de calcul des flux qui lui sont associées ;
• d'une valeur réelle pour la transmittivité ;
• d'une valeur réelle pour l'interpolation.
 109

Remarque. Pour les classes Element, Face et InterFace, la recherche d'erreurs de program-
mation est facilitée par :

1. un mode développement pour contrôler les accès aux Vecteurs ;

2. l'emploie des macros _LINE_ et _FELE_ pour afficher la ligne et le fichier lorsqu'une
opération non-permise est effectuée (ces macros sont standards aux langages C et C++).

6.2.7 Classe EquLin

La représentation et la manipulation des systèmes d'équations linéaires dont la matrice

des coefficients est creuse est effectuée via la classe EquLin. Seules les entrées non-nulles de
la matrice des coefficients et le vecteur des termes constants sont stockés. Quatre attributs
sont utilisés pour représenter les systèmes d'équations :

1. un entier pour le nombre de lignes de la matrice ;

2. une liste de Vecteurs d'entiers pour la structure topologique ;
3. une liste de Vecteur de nombres réels de format double précision pour les coefficients ;
4. un Vecteur de nombres réels en format double précision pour représenter le vecteur des
termes constants.

Soit i et y des entiers, les principales opérations possibles avec cette classe sont :

• l'accès direct aux attributs via des interfaces qui vérifient la validité des accès en mode
développement ;
• l'accès direct à la structure topologique via l'opérateur LI(i,j) ;
• l'accès direct aux coefficients via l'opérateur A(i,j) ;
• l'accès direct aux termes constants via l'opérateur B(i) ;
• la multiplication matrice-vecteur ;
• la multiplication matrice-transposée-vecteur ;
• la multiplication d'une ligne par un vecteur ;
• la mise à l'échelle des lignes ou des colonnes ;
• le calcul du résidu ;
• le calcul de la norme matricielle.
 110

6.2.8 Classes VFNS et USAGER

Une instance de la classe VFNS est un algorithme de résolution d'un problème d'é-
coulement. Cette classe est utilisée via la classe USAGER (ou une autre classe), VFNS étant
une classe template. C'est la classe USAGER qui sert d'interface entre l'usager qui doit ré-
soudre un problème d'écoulement et la classe VFNS. La classe USAGER est constituée d'un
attribut qui est une instance de la classe VFNS et d'une série de méthodes dont la programma-
tion est effectuée par l'usager. Ces méthodes sont utilisées pour la description des conditions
de bord, l'initialisation des paramètres qui décrivent l'écoulement (turbulent ou laminaire,
transitoire ou permanent...), le calcul des termes sources et des paramètres de diffusion.
Cette classe comporte de nombreuses méthodes. Ici, seules les méthodes publiques1
sont décrites.

Méthode Resoud()

Cette méthode est appelée pour débuter la résolution.

Méthode SauvegardeSolO

Cette méthode sauvegarde la solution. La solution peut être sauvegardée en format

binaire ou ASCII.

6.2.9 Méthodes de la classe USAGER

La modélisation des différents problèmes d'écoulements est supportée via la program-

mation du corps de certaines des méthodes de la classe USAGER.

Méthode InitParametresO

Les principaux paramètres nécessaires à la résolution sont initialises via cette méth-
ode :
'On peut distinguer deux types d'attributs : les attributs publiques et les attributs privées. Les attributs
publiques sont visibles aux clients de la classe. Les attributs privés ne sont visibles qu'aux objets de la classe.
Cette distinction est nécessaire pour 1'encapsulation.
 Ill

1. écoulement laminaire ou turbulent ;

2. écoulement permanent ou transitoire ;

3. valeur du paramètre d'anti-diffusion ;

4. nombre de variables scalaires à résoudre (la température, une concentration,...) ;

5. dimension d'espace (2 ou 3) ;

6. le nom du fichier de maillage ;

7. le format des entrées-sorties (ASCII ou binaire).

Méthode ConditionFrontiere()

Le type et la valeur des conditions de frontière sont déterminés par la méthode Con-
ditionFrontiere().

Méthode InitSOQ

La solution initiale, pour toutes les variables, est déterminée à l'aide de cette méthode.

Méthode ChoixDtQ

Avant chacune des itérations (ou pas de temps), le valeur du pas de temps est déter-
minée par un appel à cette méthode. L'usager peut donc modifier son pas de temps en fonction
du temps ou du nombre d'itérations.

Méthode DensiteQ

Avant la projection, la valeur de la densité peut être calculée en fonction de l'état de

l'écoulement (température, pression, concentrations volumiques). Cette méthode autorise la
résolution d'écoulements avec densité variable.

Méthode DiffusionO

Cette méthode est appelée pour déterminer le coefficient de diffusion avant le calcul
des équations discrètes.
 112

Méthode Source()

Les termes sources des équations de conservation sont déterminés avec cette méthode.

Méthode ParametresResolutionQ

Avant la résolution des systèmes d'équations linéaires, cette méthode est appelée pour
déterminer les paramètres requis par le résoluteur :

1. choix du résoluteur (gradient conjugué, Orthomin2, GMRES)

2. conditionnement du système d'équations linéaires ;

3. critère de convergence ;
4. nombre maximal d'itérations.

Méthode ArretQ

Cette méthode est appelée après chacun des pas de temps pour déterminer s'il faut
arrêter la résolution. Le critère d'arrêt peut être choisi en fonction du nombre de pas de
temps, de la norme résiduelle, de la variation des variables.

Méthode PostTraitement()

Après chacun des pas de temps, cette méthode est appelée pour permettre à l'usager
d'effectuer des opérations de post-traitement (le calcul des coefficients de portance ou de
traînée par exemple).

6.2.10 Classe Géométrie

La manipulation des éléments du maillage est effectuée avec les objets de la classe
Géométrie. Ses membres sont essentiellement des méthodes qui servent à la manipulation
des structures de données qui représentent un maillage.
 113

Méthode ConstruitFacesNSO

À partir de la liste des sommets des éléments du maillage, cette fonction construit la
liste des faces du maillage. Le maillage peut-être constitué de triangles en 2D ou de tétraèdres
en 3D.

Méthodes MesureTetraO et MesureTriangIe()

Calcule le volume (l'aire) d'un tétraèdre (d'un triangle).

Méthode CalculelDBOO

Détermine la position du point d'intersection des bissecteurs orthogonaux d'un trian-

gle en 2D ou 3D.

Méthode CalculeIDO()

Détermine le centre de la sphère circonscrite à un tétraèdre.

Méthode LectureAMDO

Lecture d'un maillage en format AMDBA (figure 6.1, page 114). Le maillage peut
être 2D ou 3D et écrit en format binaire ou ASCII.

Méthode ConstruitVCO

Cette fonction construit la liste des volumes de contrôle à partir des éléments géomé-
triques du maillage.

6.2.11 Diagramme de relation entre les classes

La figure 6.2 (page 115) présente un diagramme de relation entre les classes.
 114

Soit,
• ns le nombre de sommets (nombre entier),
• ne le nombre d'éléments (nombre entier),
• fr un entier qui peut indiquer le numéro de frontière,
• x, y et z des nombres réels qui indiquent les coordonnées d'une position,
• s i , . . . s 4 des nombres entiers qui indiquent une référence à un sommet,
le format AMDBA d'un maillage est :
1. Maillage 2D
ns ne
lxy fr
2 xy fr

ns x y fr
1 si s2 s3 fr
2 si s2 s3 fr

nt si s2 s3 fr

2. Maillage 3D
ns ne
lxy z fr
2xy z fr

ns x y z fr
1 si s2 s3 s4 fr
2 si s2 s3 s4 fr

nt si s2 s3 s4 fr

FiG. 6.1 - Format AMDBA d'un maillage

 115

InterFace
Face
jFace
Element iElement
1
InterFace
J Vecteur
VecteurBase"

f
,int ;
j double |
Vecteur
, double i T T^
71 ^"1
USAGERi CT^

[Matricfe]" VFNS

Usager

Géométrie EquLin

Légende
Classe A Classe A est Classe B
dérivée
de la classe B

Classe B possède
Classe A au moins une Classe B
instance
de la Classe A

Classe B possède
Classe A au moins une Classe B
référence
de la Classe A

Classe A est
Classe* Al' un template
de classe T

FiG. 6.2 - Diagramme de relation entre les classes

 116

6.3 Programmes pour le pré et le post-traitement

Une suite de programme a été crée avec la librairie. Ces programmes sont de petits
utilitaires qui sont utilisés pour le pré et le post-traitement.

6.3.1 Programme NopoToAMD

La plupart des maillages utilisés pour ce travail de recherche ont été créés avec la li-
brairie MODULEF. Cette librairie a été développée à l'INRIA en France et elle est disponible
gratuitement via internet (http ://[Link]/modulef/). Les maillages créés avec cette
librairie ne sont pas tous sauvegardés en format AMDBA. Le format standard pour cette li-
brairie est le format NOPO qui est binaire2. Les maillages 2D et 3D de format nopo peuvent
être traduits au format AMDBA avec le programme NopoToAMD3.

6.3.2 Programme VcToAMD

Une solution constante par volume de contrôle peut être convertie en une solution
constante par élément à l'aide du programme VcToAMD.

6.3.3 Programme AmdToGMV

Le programme AmdToGMV est utilisé pour convertir une solution constante par élé-
ment et un maillage AMDBA en un format qui puisse être interprété par le logiciel GMV.
GMV est un logiciel de visualisation de solutions 2D ou 3D développé au département de
physique appliqué du laboratoire de Los Alamos aux États-Unis. Ce logiciel est disponible
gratuitement pour de nombreuses architectures (SparcSolaris, HP-UX, Intel Linux...) à
l'adresse web suivante : http ://[Link]/XCM/gmv/[Link].
Comparativement au format ASCII, le format binaire ne peut pas être lu avec un éditeur de texte. Cependant,
les accès à ce type de fichier sont beaucoup plus rapides. De plus, ces fichiers occupent beaucoup moins d'espace
disque et il n'y a aucune de perte d'information. Malheureusement, les fichiers binaires ne peuvent pas être
transférés d'une architecture à une autre (d'un PC à une station Solaris par exemple).
3
Pour compiler le programme NopoToAMD, il est nécessaire de posséder la librairie MODULEF.
 117

6.4 Installation et disponibilité

La librairie informatique est disponible gratuitement à l'adresse web suivante :

http ://[Link]/~sboivin/sperron.
Sur une station Unix ce fichier est décompressé à l'aide de la commande :
gunzip [Link] : tar -[Link].
Lorsque l'archive est décompressée, le répertoire vfns est créé. Pour compiler toute la librairie
sur une station Linux, il suffit de se placer dans le répertoire vfns/bin et de faire la commande
make all. Pour les autres architectures Unix, il faut éditer le fichier Makefile et modifier les
paramètres de compilation (ils ont été optimisés pour une architecture Intel Linux). Avec cette
commande, toute la librairie est compilée ainsi que les programmes utilitaires et quelques
programmes de démonstration. En ce qui concerne la documentation, la thèse de doctorat
est disponible à la même adresse dans les formats PostScript compressé ([Link]) et PDF
([Link]).
 Chapitre 7

Résultats numériques

Introduction

Dans ce chapitre, on présente quelques uns des résultats numériques qui ont été obte-
nus avec notre schéma. Les problèmes résolus sont des écoulements 2D et 3D, stationnaires
ou instationnaires, en régime laminaire et turbulent. L'objectif était de valider le schéma
numérique avec des problèmes où les difficultés sont bien identifiées, et non pas de résoudre
des problèmes où les phénomènes sont nombreux et les difficultés difficilement identifiées.
On tient à préciser que bien que la majorité des écoulements présentés dans ce chapitre
soient bidimensionnels, ils ont tous été résolus dans un domaine à trois dimensions.

FiG. 7.1 - Géométrie, écoulement entre deux plans parallèles

 119

6 9 12 15
Pas de Temps

FIG. 7.2 - Convergence, écoulement entre deux plans parallèles

7.1 Écoulement entre deux plans parallèles

II s'agit d'un écoulement entre deux plans parallèles dont le profil de vitesse est
développé à l'entrée, La sortie est libre et les plans sont des parois imperméables. L'axe
x est un axe de symétrie. Le mai liage est constitué de parallélépipèdes droits (des macro-
éléments formés de 6 tétraèdres), 60 suivant l'axe z, 20 suivant l'axe y. La solution initiale et
les conditions de bord sont les suivantes :

Domaine :
[0.0,0.1] x [0.0, 1.0] x [0.0,0.2]

Conditions de bord :

Entrée :
 120

Sortie :

dv
Ri)
= 0, w(y = L,z)=Q, P(y = L,z) =
dn

Parois :

Solution initiale
v(x) = 0, P(x) = 0.

Nombres caractéristiques :

ne = = 100, p = 1.0, ÎWI; = 1, L = 1.0, // = 0.01

Pas de temps :
St = 0.1.

Critère de convergence :
\vn+l_vn <2xlO"5.
I oo

L'objectif, pour ce problème, est de vérifier si le schéma converge vers la solution

analytique pour toutes les variables du système d'équations. Pour ce faire, on a calculé pour
dp
les variables v, p et -^- les valeurs absolues des plus grandes différences entre les valeurs
dy
calculées numériquement et la solution analytique. La solution analytique pour ce problème
est la suivante :

v(y,z) = lOOz (h — z) .

dy

Les courbes de convergence sont présentées sur la figure 7.2 (page 119). Ces courbes mon-
trent que le schéma converge rapidement vers la solution. Après seulement 15 pas de temps,
 121

0.2 1 11 '1

0.18

0.16

0.14
N
0.12

0.1 Notre so
\
0.08
s olulion c Lit ion
xacic - -
/
0.06

0.04 y
0.02

0 L_ _ l1
0.2 0.4 0.6 0.8
Vitesse v

FIG. 7.3 - Vitesse horizontale v, écoulement entre deux plans parallèles

le système d'équations a satisfait le critère de convergence. À convergence, on remarque que

la plus grande différence entre la vitesse calculée numériquement et la solution analytique est
de l'ordre 1 x 10~3. À la page 121 sur la figure 7.3, notre solution numérique et la solution
analytique de la vitesse horizontale près de la sortie (y = 0.98) sont présentées. De plus,
sur la figure 7.4 (page 122) on montre la solution numérique et la solution analytique de la
pression au centre de la conduite (z = 0.1).
Ce problème montre que notre façon de mettre à jour la vitesse sur les volumes de con-
trôle ainsi que notre approximation du gradient de pression donnent des résultats numériques
qui sont consistants avec la solution analytique du problème.
 122

2 1

1.8 veil re S( 1 Li Li on
S(Jution ;xacte
1.6

1.4 "1
K
1.2 i

«V.D 11 s

£
0.8
.73]
0.6

0.4

0.2 J j
0
0 0.1 0.2 0,3 0.4 0.5 0.6 0.7 0.8 0.9 1
Position y

FIG. 7.4 - Pression, écoulement entre deux plans parallèles

7.2 Écoulement engendré par le déplacement d'une paroi

Le problème considéré est celui de l'écoulement dans une cavité carrée engendré par
le déplacement d'une paroi. L'objectif, pour ce problème, est de vérifier :

1. si le schéma anti-diffusif permet d'améliorer significativement la précision de la solu-

tion numérique ;
2. l'accord entre le champ de vitesse calculé lors de la projection et la solution des équa-
tions de la quantité de mouvement.

Les conditions de bord, la solution initiale et les nombres caractéristiques pour ce problème
sont les suivants :

Domaine :
[0.0,0.1] x [0.0, 1.0] x [0.0, 1.0].

Conditions de bord :

1 z = L
?
— Isn — — —
0 ailleurs
 123

Solution initiale :
v(x) = 0, P(x) = O.

Nombres Caractéristiques :

pVmaxL
Re = = 1000 , p=1.0, vmax = l, L = 1.0, // = 0.001
4
A

Pas de temps :
ôt = 0.1.

Critère de convergence :

vn+l _ n+1 _ JQ
-5
OO

Pour ce problème, on a utilisé plusieurs maillages constitués de prismes à base trian-

gulaire1 . La base des prismes repose dans le plan yoz .
Ce problème a été étudié par de nombreux chercheurs et nos résultats ont donc pu être
comparés. Dans le tableau 7.1, nos résultats sont comparés à ceux obtenus par Ghia et al.[40]
avec un schéma d'ordre 2 et 10000 inconnues. Les extrémas des composantes de la vélocité
sur les plans médians ainsi que la position où ils sont atteints y sont présentés. On tient à
préciser que pour la composante v la valeur maximale n'est pas présentée, celle-ci étant bien
sur atteinte lorsque z = 1.0 (sur la paroi supérieure).
On peut constater que l'on se rapproche de la solution de référence lorsque l'on aug-
mente le nombre d'inconnues. De plus, les résultats sont plus précis lorsque le schéma anti-
diffusif est utilisé. Pour bien juger de la précision de nos résultats, on montre, sur la figure
7.5 (page 124), les résultats sur les plans médians obtenus avec le maillage comportant 6178
volumes et les résultats de Ghia et al.
A des fins de comparaison, le champ de vitesse et les courbes de pression sont présen-
tés sur les figures 7.6 et 7.7 (pages 125 et 126).
1
On rappelle que lorsque le problème considéré se déroule dans un plan, les variables du système d'équations
sont indépendantes de la direction perpendiculaire à ce plan. Pour profiter de la symétrie et utiliser un maillage
non-structuré, on utilise des prismes à base triangulaire.
 124

Maillage Antidiffusion Min. v {z =) Min. w (y =) Max. w (y = )

2220 volumes non -0.299, (0.233) -0.432, (0.883) 0.286, (0.184)
2220 volumes oui -0.341, (0.234) -0.481, (0.883) 0.325, (0.184)
3982 volumes non -0.321, (0.233) -0.448, (0.885) 0.304, (0.184)
3982 volumes oui -0.359, (0.202) -0.492, (0.895) 0.347, (0.158)
6178 volumes non -0.351, (0.224) -0.481, (0.886) 0.337, (0.184)
6178 volumes oui [^0.387, (0.199) -0.519, (0.886) 0.373, (0.167)
Ghia et al. 10000 noeuds — -0.383, (0.172) -0.516, (0.906) 0.371, (0.156)

TAB. 7.1 - Résultats comparatifs, écoulement engendré par le déplacement d'une paroi

0.4 i
r- 0.3
y 0.2
0.8 - _f.

1
/
0,1
A..... \ r
0.6 - 0 i
r
d
S / \
•a -0.1 ""i

f
& +
0.4 - <m
Ghi KI al. -0.2
vfolre w lut ion
-0.3 — jliiaelal. +
L* solution
0.2 - .— -0.4
\ j
-0.5

-0.6 i
-Û.2 0 0.2 0.4 0.6 0.8 0.2 0.4 0.6 0.8
Vitesse v Position y

FlG. 7.5 - Coupes, écoulement engendré par le déplacement d'une paroi

 125

FIG. 7.6 - Vélocité, écoulement engendré par le déplacement d'une paroi

Lorsque l'on utilise une méthode à pas fractionnaires, il est important de vérifier si le
champ de vitesse calculé lors de l'étape de projection est bien une solution des équations de
Navier-Stokes. Lorsqu'un schéma est mal construit, un état dit stationnaire peut être atteint
bien que les équations de la quantité de mouvement ne soient pas bien résolues. Si c'était le
cas avec notre schéma, on constaterait que le résidu des équations de Navier-Stokes atteint
rapidement un point de stagnation bien que la variation des variables à chacun des pas de
temps tend vers zéro. Afin de vérifier s'il y a bien accord entre la solution projetée et les
équations de Navier-Stokes, on a tracé les courbes des résidus et celles des variation des
inconnues à chacun des pas de temps. La figure 7.8 (page 126) montre que les courbes des
variations et les courbes résiduelles sont en accord, la vélocité calculée lors de la projection
est bien une solution des équations de la quantité de mouvement.
 126

Cells

FIG. 7.7 - Pression, écoulement engendré par le déplacement d'une paroi

1 1 1 1 1 1 1 1 1 I i t

lo
-0.5 Sl0 , 2 Vitesse u
10
Sio 0 i°8io Ikil
-i - logio
L2 pressiutl P

-1.5 -1

-2 - -
-2 \ ^ ^ ^ ^ ^
-

-2.5
-3
-3

-3.5 -
" * • • - . . . ,

-4 -
-5 ' * • • ' - . ,

-4.5 . -

1 1 t 1 <:
-5
0 50 100 150 200 250 300 350 50 100 150 200 250 300 350
Pas de temps Pas de temps

FIG. 7.8 - Convergence, écoulement engendré par le déplacement d'une paroi

 127

7.3 Convection naturelle dans une cavité carrée

Le problème suivant est celui de l'écoulement thermique dans une cavité carrée. L'ob-
jectif, pour ce problème, est de vérifier le comportement du schéma lorsque l'écoulement
dépend non seulement des équations de conservation de la quantité de mouvement et de
la masse, mais d'une autre équation de conservation. Pour ce problème, toutes les équa-
tions sont fortement couplées, les solutions pour la température et la vélocité étant fortement
dépendantes l'une de l'autre.
L'écoulement est engendré par l'impulsion causée par le changement de densité, celle-
ci étant dépendante de la température. Toutefois, puisque la densité est considérée con-
stante, on modelise l'impulsion causée par le changement de densité avec l'approximation
de Boussinesq (chapitre 1).
Les données pour la mise-en-oeuvre du problème sont les suivantes :

Domaine :
[0.0,0.1] x [0.0, 1.0] x [0.0, 1.0].

Conditions de bord :

V. • n\aQ = 0,

= L,z) = l, ^ = 0,
z=O, z=L

Solution initiale :
v{x) = 0, p(x) = 0, T(x) = 0.

Nombres caractéristiques :

iiCv „ p2gL3PAT A^n „ _ fi

Pr = — - = 0.71, Gr = rl = 1408470, Ha = Pr • Gr = 10
k [i

Pas de temps :
6t = 0.001.
 128

(!/=)
iax ^77101(z=)
z 1/2 v = 1/2
Solution de référence 219.36 (y = 0.962) 64.63 (z = 0.150)
Nos résultats 220.35 (y = 0.965) 65.29 (z = 0.159)
TAB. 7.2 - Résultats comparatifs, convection naturelle dans une cavité carrée

Critère de convergence :

vn+1 - vn < 1 x 10~2 et wn+l _ < 1 X 10~2

Le maillage, pour ce problème, est non structuré et composé de prismes à base tri-
angulaire, le maillage comportant 6178 volumes de contrôle (c'est le même maillage qui a
été utilisé pour le problème de l'écoulement engendré par le déplacement d'une paroi). On
a considéré que l'état stationnaire avait été atteint lorsque la variation de vitesse est devenue
inférieure à 1 x 10~2 (soit « 0.01% de la vitesse maximale).
Dans le tableau 7.2 les valeurs maximales pour les composantes du vecteur vitesse
sur les plans médians sont comparées aux résultats obenus par d'autres chercheurs ([31] et
[32]). Ces valeurs sont en très bon accord avec les solution de référence. On tient à men-
tionner que la solution de références a été obtenue avec une méthode de discrétisation d'un
degré d'approximation supérieur et que cette dernière solution a été validée avec une solu-
tion expérimentale. On trouve sur la figure 7.9 (page 129), les courbes de convergence pour
la variation des variables et les résidus. Les courbes des variations sont en accord avec les
courbes résiduelles.
Finalement, pour des fins de comparaison, la vélocité, la pression et la température à
convergence sont présentées sur les figures 7.10, 7.11, et 7.12 des pages 130 et 131.
 129

20 40 60 80 100 120 140 160 180

Pas de Temps

4
i o l l r ; l l L 2 vitesse «
3 i o 1I T ;II L 3 viiessciu

lO H>N-|IL2 Passion P
2 i o I b i l l ^ a Tcni^rolureT

-1

-2

-3

-4

-5
20 40 60 80 100 120 140 160 180
Pas de Temps
FiG. 7.9 - Convergence, convection naturelle dans une cavité carrée
 130

FIG. 7.10 - Vélocité, convection naturelle dans une cavité carrée

Cllll
pression
a - ? ; » » os
|-a.l/«+3S
••!.M« -o

FlG. 7,11 - Pression, convection naturelle dans une cavité carrée

 131

CEIII

-0,6*1

-O.7S5
-0.741

-osm
-0-517

-O.W
-03SB
-nav

-0.Î49

FIG. 7.12 -Température, convection naturelle dans une cavité carrée

7.4 Écoulement turbulent au-dessus d'une marche

II s'agit d'un écoulement permanent turbulent au-dessus d'une marche (figure 7.13).
L'objectif pour ce problème est de vérifier le comportement du schéma lorsqu'un modèle de
turbulence est utilisé. Nos résultats seront comparés à des résultats expérimentaux, obtenus
par Makiola[56]. On tient à préciser que le modèle de turbulence qui a été mis en oeuvre,
le modèle k — e, n'est pas bien adapté à ce type d'écoulement (voir Speziale[71], Chang et
al.[22] et Wilcox[74]). En effet, avec ce modèle de turbulence, l'estimation de la longueur

2h

L=22h

FiG. 7.13 - Géométrie, écoulement turbulent au-dessus d'une marche

 132

de recirculation est toujours trop faible. L'objectif n'est donc pas d'obtenir une solution "ex-
acte", mais une solution physiquement réaliste et comparable à celle obtenue par d'autres
chercheurs. La solution initiale, les conditions de bord et les valeurs de référence pour ce
problème sont les suivantes :

Conditions de Bord :

Entrée :

v(x, y = 0, z) = Valeurs expérimentales, w(x, y = 0, z) = 0,

2 3

Sortie :

dv
= 0, w(x, y = 24h, z) = 0, P(x, y = 24h, z) = 0,
dn j/=24/i

= 0.
dn y=24h
an j/=24/i

Parois :

-T- (v-T1) = -p {U*f ,

on \v\
= nxz1, (// + fit) — {y • r2) = 0,
n

Solution initiale :

v(x) = 0, P(x) = 0,
= l.e-4, e(x) = L e - 4.
 133

Nombres caractéristiques :

pVmaxh
Re = = 15000, vmax = 1.0, h = 0.1, p = 1.0.

Pas de temps :
ôt = 0.02.

Critère de convergence :

vn+1 -vn < 1 x 10" 4 et wn+l -wn < 1 x 10~4 .

OO 00

Pour ce problème, tous nos volumes de contrôle sont des hexaèdres, le maillage étant
constitué de 12527 volumes de contrôle.
La figure 7.14 (page 134) montre le profil de la vitesse le long de la conduite. Bien
que nos résultats ne soient pas les mêmes que ceux obtenus expérimentalement, ceux-ci sont
réalistes. Soit x la longueur de la zone de recirculation, le rapport | est de 7.3, ce qui est de
11% inférieure à la valeur 8.2 prédite expérimentalement. Cette erreur est comparable à celle
d'autres chercheurs pour un problème du même type et le modèle k — e. En effet, l'erreur que
l'on retrouve dans la littérature pour des problèmes comparables varient de 14% à 31%.
Finalement, on trouve les courbes de convergence sur la figure 7.15.
 134

0.16 - O.'B X i.
•
V
0,1! . " ' " -
0 12

= 63 OUQQ

ce-1.

î
PB •

•'

(1
0 0.2 0.4 0.B O.S 1 0 O.S 0.4 0.6 O.S 1 •0.2 0 OS 0.1 O.« 0.8

0.3

. 5 16 3 IB
i1C ^

} •

O.'-B
/ •
«

3. se [Link]
«
A
0,04
* !
- . ' • '

Û M Q
-0.2 0 02 0.4 0.8 0.S 0 02 0.4 0.8 OS o 0.2 oj as

0.2

0 0.2 0.4 O.S 0 0.2 0.4 0r3

+ Résultats expérimentaux
Nos résultats

FiG. 7.14 - Profils de vitesse, écoulement au-dessus d'une marche

 135

-Pf
0
I;
-1 I*. • i r a
l-'-Jsf

-2
-3
^^Wrtfoiatai-
-4
-5
-6

-7
0 100 200 300 400 500
Pas de Temps

2 i 1 1 1—
lOBTnllrjllj.î vitHfBÏ

1 lûtfift 11 -rv II F i nresriou J3

login I k i l l r î JissipniiMiî
-
0 loiin \\TI\\ , s^ncrsiefc

Vite
-1
r \>
-2 " • - - • ^ ^ ~ ~ - ~ - ~ ^ ^ ^

X. -
-3 i \'i

-4
—
-5 -

i_ ' 1

0 100 200 300 400 500

Pas de Temps

FIG. 7.15 - Convergence, écoulement au-dessus d'une marche

 136

0.41

FlG. 7.16 - Géométrie, écoulement 2D autour d'un cylindre

7.5 Écoulement autour d'un cylindre

II s'agit d'une série de problèmes d'écoulements laminaires incompressibles autour

d'un obstacle. Le but est de vérifier si notre schéma peut adéquatement calculer les valeurs
des coefficients de portance et de traînée pour des cas transitoires et permanents en 2D et 3D.

7.5.1 Écoulement 2D permanent autour d'un cylindre

Les données pour la mise-en-oeuvre de ce problème sont les suivantes :

Domaine : Figure 7.16.

Conditions de bord :

Entrée :

_ 4 • vmoy -(H-z)
v(y — U, Z) — 2 ,

Sortie :

dv
= 0, w(y = 2.2, z) = 0, P(y = 2.2, z) = 0,
dn y=2.2

Parois :
V.' XL— 0 > V'Z=0.
 137

Solution initiale :
v(x) = 0, P(x) = 0.

Nombres caractéristiques :

pVmoy
Re = = 20, p=1.0, vmoy = 0.3, D = 0.1, /i = 0.001

Pas de temps :
ôt = 0.1.

Critère de convergence :

< 2 x 10~5 et wn+l -wn < 2 x 10" 5

oo

Ce problème a été résolu sur trois maillages différents. On pourra donc vérifier si le
maillage influence significativement la précision lors du calcul des coefficients de portance
et de traînée. Puisqu'il s'agit d'une géométrie relativement complexe, on a pu exploiter la
flexibilité du schéma et utiliser différents types de volumes de contrôle au sein du même
maillage :

• semi-structure à l'entrée et à la sortie ;

• semi-structure autour du cylindre ;

• non structuré dans la zone de contournement.

La figure 7.17 de la page 138 montre des coupes du maillage à l'entrée de la conduite et
autour du cylindre. On rappelle que les tétraèdres (les triangles) des parties semi-structurées
du maillage sont combinés pour former des hexaèdres. Cette configuration est préférable aux
triangles. En effet, ce maillage est mieux adapté à l'écoulement (le maillage ne nuit pas à
l'écoulement).
 138

/ / // /
//
t/ / / /V A/ / /
/
m m
i
1iI
s / / / / 7

f /

s/ /
> f
7
> !
s/ 7y
I1
~'-y /
r
//
/-/
/ /4 /-/ // /
/ /
7 j y

y
7 i i1
''
i
/

'7 y 7y ?7 7 7 7/ 77
ii
'/

i /
///> / '/
y J 77y yy/ / z )
y

m\/•>

/ 7 7/ / / / / / / /L/L/
i / / /

FiG. 7.17 - Coupes du maillage, écoulement 2D permanent autour d'un cylindre

 139

Pour ce problème, on n'a pas de solution de référence. Toutefois, on pourra faire des
comparaisons à l'aide des quantités suivantes :

1. le coefficient de traînée :

2-FD
CD =

2. le coefficient de portance :

2-FL
CL =

3. la différence de pression aux extrémités du disque :

A P = P(0, 0.15,0.2) - P(0, 0.25,0.2) ;

où:

• L = 0.10 est la profondeur, l'écoulement étant résolu en 3D ;

• r est un vecteur tangent à la surface du disque ;

• vT = y_ • r est la vitesse tangente à la surface du disque ;

• n est le vecteur normale unitaire pointant vers l'extérieur du disque ;

• dS est la surface du disque.

Les résultats obtenus pour les trois maillages sont présentés dans le tableau 7.3. Ces résultats
sont en accord avec ceux de Turek[66]. On constate toutefois que les maillages les plusfinsne
permettent pas d'améliorer significativement la précision des résultats. En effet, le meilleur
coefficient de portance est obtenu avec le deuxième maillage et le meilleur coefficient de
traînée avec le troisième maillage. On tient à préciser que les maillages les plus fins ont été
obtenus manuellement en ajoutant des noeuds lors de la génération du maillage, et non pas
en utilisant une méthode de raffinement automatique du maillage. Bien que cette façon de
 140

Volumes de contrôle CL cD AP
15352 0.00970 5.570 0.1165
23162 0,0110 5.567 0.1167
27085 0.0120 5.574 0.1169
Valeurs de référence 0.0104-0.0110 5.57-5,59 0.1172-0.1176
TAB. 7.3 - Coefficients, écoulement 2D permanent autour d'un cylindre

0 1 I I I T I \

-0.5
0 lo
-1 K i o l l r i l l L 2 vinsse tu

-1.5 k i l | I r 2 pression ,P
-1
-2
-2
-2.5
Wh_!^
-3 -3
-3.5
-4 -4
-4.5
-5
-5 ! îïï -
i i i i i I I I I I I I I
-5.5 -6
0 10 20 30 40 50 60 70 80 90 100 0 10 20 30 40 50 60 70 80 90 100
Pas de temps Pas de temps

FlG. 7.18-Convergence, écoulement permanent 2D autour d'un cylindre

procéder puisse donner de bons résultats lorsque la géométrie est simple (tel le problème
de l'écoulement engendré par le déplacement d'une paroi), elle paraît ne pas être efficace
lorsque la géométrie est plus complexe. Pour bien étudier la convergence du schéma lors du
raffinement du maillage, il faudrait utiliser un remailleur (déplacement et ajout de noeuds) ou
bien un algorithme de raffinement automatique (ou isotrope) du maillage.
On montre les courbes de convergence sur la figure 7.18. On constate que les varia-
tions des trois variables (particulièrement la pression) ne diminuent pas de façon asympto-
tique, mais cyclique. Ce comportement est toutefois en accord avec les courbes résiduelles.
Bien que les variations pour ces courbes soient beaucoup moins prononcées, elles aussi sont
cycliques.
 141

7.5.2 Écoulement 2D transitoire autour d'un cylindre

Ce problème est similaire à celui de la section précédente. Toutefois, la vitesse à

l'entrée est plus élevée et l'écoulement est périodique. En plus des valeurs comparatives du
problème précédent, le nombre de Strouhal St = ^ - (où / est la fréquence de séparation) est
calculé. Les périodes débutent lorsque le coefficient de portance atteint sa valeur maximale.
On présente également la différence de pression AP calculée au temps t0 + jrf, t0 étant
le début d'une période. Puisqu'il s'agit d'un écoulement périodique, on pourra vérifier si le
schéma prédit adéquatement l'évolution d'un écoulement transitoire laminaire.
Les données pour la mise-en-oeuvre du problème sont les mêmes que pour le prob-
lème précédent. Toutefois, la vitesse d'entrée est supérieure et il n'y a pas de critère de con-
vergence (l'écoulement étant périodique), mais plutôt un arrêt automatique après un nombre
pré-établi de pas de temps.

Nombres caractéristiques :

pVfn
Re = °yD = 100, 0=1.0, vmoy = 1.0, D = 0.1, n = 0.001.

Pas de temps :
St = {0.002, 0.0015, 0.001, 0.00075}

Les résultats obtenus sont présentés dans le tableau 7.4 à la page 142. Pour toutes les simula-
tions, la valeur du paramètre /5 de l'équation de projection

,, _ 7 ,n+l/2

0=—sf— = ~
est 3/2 pour une extrapolation du gradient de pression (voir figure 3.2, page 58). Le maillage
est similaire à celui utilisé pour le problème précédent. Il est toutefois plus fin, ce dernier
maillage étant composé de 38645 volumes de contrôle.
Nos résultats sont en accord avec ceux présentés par Turek[66]. Ce test montre que
la méthode de projection qui est privilégiée pour notre schéma (section 4.7), permet bel et
bien de simuler adéquatement des écoulements transitoires et non pas seulement de calculer
la solution d'écoulements permanents.
 142

ôt cD cL St = ^-
umoy
AP
t = t0 t = t0 t = to + J7
0.002 3.25 1.076 0.287 2.50
0.0015 3.24 1.043 0.287 2.50
0.001 3.23 1.001 0.288 2.50
0.00075 3.22 0.984 0.291 2.49
Valeurs de référence 3.22-3.24 0.99-1.01 0.295-0.305 2.46-2.50
TAB. 7.4 - Coefficients, écoulement 2D périodique autour d'un cylindre

0.41

0.41

FlG. 7.19 - Géométrie, écoulement 3D autour d'un cylindre

7.5.3 Écoulement 3D permanent autour d'un cylindre

II s'agit d'un écoulement 3D permanent autour d'un cylindre placé dans une conduite
carrée. Le domaine est montré sur la figure 7.19. Les données pour la mise-en-oeuvre sont
les suivantes :

Conditions de bord :

Entrée :

u(x,y = 0,z) = 0,
16 • Vmoy -X • Z- (H - X)(H - z)
v(x,y = 0,z)
w(x,y = 0, z) = 0.
 143

Sortie :

u(x,y = 2. 5,2) = o,
dv
dn y=25 o,
w(x,y = 2. 5,z) = o,
P(y = 2. 5, z) = 0.

Parois :
y • n = 0, w • r = 0.

Solution initiale :
«(£) = 0, P(z) = 0.

Nombres caractéristiques

= 0.1, ff = 0.41, /i = 0.005,

Pas de temps :
6t = 0.1.

Critère de convergence :

- U
-4
< 2 x 10"4, < 2 x 10"4 et < 2 x 10~4
 144

II s'agit du premier problème résolu où les phénomènes sont réellement tridimension-

nels. On peut donc vérifier si la méthode peut adéquatement prédire un écoulement 3D. Le
problème a été résolu avec trois maillages. Ces maillages sont similaires à ceux utilisés pour
résoudre l'écoulement 2D permanent autour d'un cylindre. Pour ces trois maillages, on a con-
struit un maillage 2D dans le plan yoz. Le maillage 3D a ensuite été obtenu en effectuant 25
translations du maillage 2D. Tout comme pour les deux problèmes précédents, les solutions
ont été validées à l'aide des quantités de références suivantes :

1. le coefficient de traînée :

V5 r\ TTTT '
P • Vmoy
moy D H - H
FD =
l
2. le coefficient de portance :

2-FL
p-v^-D-H-H'

3. la différence de pression pour deux points du cylindre :

AP = P (0.205, 0.45, 0.2) - P (0.205, 0.55, 0.2)

où:

• r est un vecteur tangent à la surface du cylindre ;

• vT est la vitesse tangente à la surface du cylindre ;
• n est le vecteur normale unitaire pointant vers l'extérieur du cylindre ;
• dS est la surface du cylindre.

Les résultats obtenus sont présentés dans le tableau 7.5 (page 146). Les quantités
calculées avec notre schéma sont en très bon accord avec celles présentées par Turek[66]. Tout
comme pour les autres problèmes, on montre les courbes de convergence pour le troisième
maillage sur la figure 7.20 de la page suivante.
 145

1 1 1 1 1 1

n +1 til
0 — _ •

1OE 10
*\ •-.• '. fi
p
-1 log 10 ' *L

-2

-3

-4

-5 -
i i i i
-6
10 20 30 40 50 60 70
Pas de Temps

loglO INIItS vilesseu

ïil
login l|Tj|li2 t

-4

-5
10 20 30 40 50 60 70
Pas de Temps
FIG. 7.20 - Convergence, écoulement 3D permanent autour d'un cylindre
 146

Maillage Inconnues Temps de Itérations Espace CL CD AP

résolution mémoire
(secondes) (Mégaoctets)
Ml 646600 7994 72 394 0.0115 6.201 0.170
M2 794000 9400 73 482 0.0105 6.187 0.170
M3 1069000 13514 72 607 0.00824 6.180 0.170
Valeurs 753664 à 88 à 0.008 à 6.05 à 0.165 à
de référence 12582912 3571 0.010 6.25 0.175

TAB. 7.5 - Coefficients et statistiques, écoulement 3D permanent autour d'un cylindre

7.6 Écoulement thermique 3D dans un cylindre

II s'agit d'un écoulement thermique forcé dans un cylindre. À l'entrée du cylindre la

vitesse et la température du fluide sont constantes. L'objectif est de vérifier si notre schéma
peut adéquatement simuler un problème de transport d'un scalaire passif. Pour ce problème,
la solution initiale, les conditions de bord et les valeurs de référence sont les suivantes :

Domaine :

0.05 • cos(0)
0.05 • sin(0) 0<^<TT/2, 0<Z<1.5.
z z
\

Conditions de bord :

Entrée :

u(x,y,z = 0) = 0, v(x,y,z = 0) = 0, w(x, y, z = 0) = 1,

T(x, y, z = 0) = 1.

Sortie :

u(x, y , z = 1.5) = 0 , v(x,y,z = , — = 0,

2=1.5

= 0.
dn 2=1.5
 147

Parois :
v-n = 0, V-T = 0, T = 0.

Solution initiale :
v(x) = 0, P{x) = 0 , T(x) = 0.

Nombres caractéristiques :

ReD — —T£!L— = 100, wmoy = 1.0, p = 1.0,

A4

Pas de temps :
5t = 0.1.

Wmoy est la vitesse moyenne transversale pour une section de la conduite, D est le diamètre
de la conduite. Pour profiter de la symétrie du problème, l'écoulement n'a été solutionné que
pour un quartier du cylindre. Le maillage est constitué de prismes à base triangulaire. La
base de ces prismes repose dans le plan x o y. Le maillage 3D a été obtenu en effectuant une
translation d'un maillage 2D. En tout, le maillage est composé de 16000 volumes de contrôle.
Il s'agit d'un des rares problèmes 3D dont on connaît la solution analytique pour une
partie du domaine. Pour un écoulement laminaire {Rep •< 2300), la longueur nécessaire pour
que l'écoulement se développe (c'est-à-dire pour que le profil de vitesse soit parabolique) est
donnée par l'expression
0.05ReD.
lam
Dans la région où l'écoulement est développé, la vitesse w et le gradient de pression peuvent
être calculées à l'aide des expressions

dP_ _
dz ~ (D/2)
 148

En ce qui concerne la température, la longueur nécessaire pour que l'écoulement soit ther-
miquement développé est donnée par

lam,T

Pour la température, on ne possède pas de solution algébrique. Toutefois, dans la

région où l'écoulement est thermiquement développé, la température axiale doit satisfaire
l'expression
d (Ts-T{x)
dz \TS-Tm(z),
où Tm(z) est la température moyenne axiale.
Les figures 7.21 à 7.24 (pages 149 à 151) illustrent les résultats obtenus.
La vitesse le long du cylindre est présentée sur la figure 7.21 pour des coupes suivant
4 positions radiales différentes à l'angle d — TT/4. Pour toutes ces coupes, la vitesse de-
vient constante et respecte la solution prédite théoriquement. En effet lorsque l'écoulement
est développé, le gradient ^ doit être nul. De plus, la distance pour que l'écoulement se
développe respecte la longueur prédite par la théorie :

( I ) t e m « 0.05 -100 -0.1 = 0.5.

La figure 7.22 compare la vitesse transversale à la solution analytique près de la sortie. La

vitesse transversale est fonction du rayon, l'angle 9 = 7r/4 est constant. On constate que la
solution obtenue respecte bien la solution analytique.
La figure 7.23 montre le graphe de la température adimensionnée

T. - T(x)
T =
Ts - Tm(z)

pour des coupes en 4 positions radiales différentes à l'angle 6 = TT/4. On constate que effec-
tivement, l'écoulement devient thermiquement développé. En effet, cette température devient
constante pour une position radiale donnée. De plus, la distance pour que l'écoulement soit
thermiquement développé respecte la valeur

0.05-100-2.0.1 = 1
 149

0 0.2 0.4 0.6 0.8 1.2 1.4

Position Ï

FlG. 7.21 - Vitesse transversale, écoulement thermique 3D dans un cylindre

2
1.8
1.6
1.4
1.2
i
0.8 Noire solution
Solution analytique
0.6
0.4
0.2
0
0,01 0.02 0.03 0.04
Rayon

FlG. 7.22 - Profil de vitesse à la sortie, écoulement thermique 3D dans un cylindre

 150

1.4

FiG. 7.23 - Profil de température, écoulement thermique 3D dans un cylindre

prédite par la théorie.

La figure 7.24 compare le gradient de pression calculé numériquement à ta solution
théorique
dP_ wmoy _ 3 2

dz
(D/2) 2
Ces coupes suivent les mêmes droites qui ont été utilisées sur les figures précédentes. Cette
figure montre que notre solution est en accord avec la solution prédite théoriquement. Le
gradient de pression dans la partie développée est constant et indépendant de la position. Fi-
nalement, comme pour les autres problèmes, on trouve les différentes courbes de convergence
sur la figure 7.25 à la page 152. Encore une fois, la vitesse calculée lors de la projection est
bien une solution des équations de Navier-Stokes. De plus, on peut constater que le système
a convergé assez rapidement vers la solution, la résolution ne nécessitant que 49 itérations.
 151

-3
-3.5

1 -4 i-O. ]
1-0.2
-4.5 r=0.3
solution analytique
1 -5

-8 -5.5
Gradient

-6

-6.5
-7

-7.5
0 0.2 0.4 0.6 0.8 1 1.2 1.4
Position z

FlG. 7.24 - Gradient de pression, écoulement thermique 3D dans un cylindre

 152

5 10 15 20 25 30 35 40 45
Pas de temps

10 15 20 25 30 35 40 45
Pas de temps

FIG. 7,25 - Convergence, écoulement thermique 3D dans un cylindre

 Chapitre 8

Extension de la méthode aux écoulements avec densité

variable

Introduction

Notre schéma peut aussi être utilisé pour la résolution d'écoulements incompressibles
où la densité est variable. Pour ces écoulements, la densité est souvent une fonction des
concentrations volumiques des constituants d'un mélange. Dans ce chapitre, on propose une
version plus générale du schéma qui supporte la résolution de cette catégorie d'écoulements.
Malheureusement, aucune validation de l'approche proposée n'a été effectuée.
Essentiellement, pour ces écoulements, la projection est effectuée après la résolution
de toutes les équations de transport et l'algorithme de résolution est :

• Étant donné une solution initiale v_~l, p~l, p°, <fP :

o Appliquer l'opérateur de projection : (y0, p°, Vp°) = P (•••)•

• Étant donnée", pn, pn, 4>n :

o pour un écoulement turbulent, calculer les coefficients de diffusion turbulente et

la vitesse de frottement sur la paroi ;

o appliquer l'opérateur de convection-diffusion aux champs scalaires telle l'énergie

cinétique turbulente, la température ou les concentrations volumiques :

• calculer la densité
 154

• appliquer l'opérateur de convection-diffusion à la vélocité : vn+1/2 = CD (•••);

• appliquer l'opérateur de projection : (va+1, p n + l , V p n + 1 ) = P (•••)•

Seules la discrétisation temporelle et la projection ont à être modifiées.

8.1 Discrétisation temporelle

Pour les écoulements où la densité est variable, l'approximation de la dérivée tem-

porelle est :

cas stationnaire
It
dt t=tn+l
autrement
2St
La densité n'est mise à jour qu'après la résolution des équations de conservation des variables
scalaires et la discrétisation temporelle de tout le système est :

1. Convection et diffusion des variables scalaires :

= sn.
dt

2. Mise à jour de la densité pour le temps t = tn+i : selon la physique du problème.

3. Extrapolation du champ de vitesse :

d(pv)
+ V • [v (t = tn+1) ® (pnvn+1/2) - n^vn+l12]+Vpn = P
dt t=tn+1

4. Projection :

«n+U.n+l _ nn.,n+l/2

St
pn+\ _ pn
+ V • (p n + V + 1 ) = 0.
It
 155

8.2 Calcul de la densité d'interface

Tout comme pour les autres propriétés physiques, la densité à l'interface est approx-
imée avec une moyenne géométrique. Soit l'interface OK,L séparant les volumes K et L, la
densité d'interface PK,z,est :

m(K)pK + m(L)pL
PK L
' m{K) + m(L) '

8.3 Projection

Après l'étape de convection-diffusion, l'équation de conservation de la masse n'est

pas nécessairement satisfaite. Pour qu'elle soit satisfaite, on effectue une projection en deux
étapes :

1. Une étape d'extension (de ré-interpolation) sur les interfaces qui permet d'éviter les
faux modes de pression.

2. Une étape de projection pour corriger la pression et la vélocité.

8.3.1 Opérateur d'extension E(- • • )

Soit v^ la vélocité à l'interface a, son évolution du temps t = tn au temps t = £n+i/2

doit satisfaire l'équation

d(pv)
= - V • [v {t = tn+1) (g) {pnvn+1'2) - //Vw n+1/2 ]+Vpn + f1
dt t=tn+1

En intégrant le membre de gauche de cette équation sur les volumes de contrôle associés à
l'interface a, on a :

d(pv) m(K) d {pv) m(L) d (pv)

dt m(K)+m(L) dt m(K) + m(L) dt L,t=tn+1
 156

La vélocité d'interface est calculée à partir des expressions :

1. écoulement permanent :

m(K) n+1/2
V.K
m(K) + m{L)
= E(vK, vL)=£ m(L) (/..n+1/2
n
€ tint
m{K) + m(L)
( n+1/2 \
[VK ~ V.K ) , £ £ext

2. écoulement transitoire :

n+l/2
v =_ E(vK, vL)

3pn

m{K)

n
v£112 - ApnvnK + pn-lv?Kl\ , a e£,ext

8.3.2 O p é r a t e u r d e p r o j e c t i o n P( • • • ) .

L'opérateur de projection est maintenant défini à partir des équations

n ii,. n+1/2
(8.1)
St
pn+\ _ pn
+ V • (p n + V + 1 ) = 0. (8.2)
It

Correction de la pression

Pour obtenir une équation dont la seule inconnue est la correction de pression, on
considère l'équation (8.1) que l'on substitue dans l'équation (8.2) :

St
 157

Ensuite, cette équation est intégrée sur un volume générique K et le théorème de Gauss est
appliqué pour obtenir l'intégrale de surface

/ V {ôpn+1) -ndS = - [ £ (pnvn+1/2) -ndS+ [ ^ ~ P dV

dK JdK àt JK àt
L'opérateur de diffusion permet d'obtenir l'équation discrète

~ SPK) = - £ Jt™

Deux types de conditions de bord peuvent être imposées :

1. flux de masse (débit massique) m imposé :

2. pression imposée :
Sp = pn+1 - pn.

Correction du champ convecteur

Après la résolution de l'équation de correction de pression, la vitesse d'interface est

corrigée à l'aide de l'équation

v"+1 =

Mise-à-jour de la vélocité sur les volumes de contrôle

La mise-à-jour de la vélocité est effectuée à l'aide de l'opérateur défini à la section

4.7.

8.4 Choix de l'algorithme de résolution des équations linéaires

Lorsque la densité varie avec le temps, la nature de la matrice des coefficients du

système d'équations discrètes est modifiée et elle n'est pas nécessairement à diagonale dom-
 158

inante. En effet, pour les écoulements où la densité est variable, l'inéquation

n'est pas nécessairement satisfaite. On rappelle que cette inéquation doit être satisfaite pour
que la matrice des coefficients soit toujours à diagonale dominante.
Dans ces conditions, les algorithmes du gradient conjugué et Orthomin2 ne sont pas
appropriés pour résoudre tous les systèmes d'équations linéaires. On recommande alors d'u-
tiliser les algorithmes Orthomin2 ou gradient conjugué pour la correction de pression et GM-
RES pour les équations de transport.

8.5 Respect du principe du maximum

Dans la section précédente, on a mentionné que l'inéquation

pa\m(o)va\ (8.3)
> Ver > 0 <J££K > "£•£int, Va < 0

n'est pas toujours satisfaite puisque

v-(p" + y + 1 ) = - p " pn
.

Malheureusement, pour que le principe du maximum soit respecté, l'inéquation (8.3) doit
être satisfaite. En effet l'inéquation V • y_ > 0 est une condition essentielle pour démontrer le
respect du principe du maximum lorsque l'écoulement est permanent et la densité constante.
En régime transitoire, lorsque la densité est variable, l'inéquation (8.3) n'est pas nécessaire-
ment satisfaite et on ne peut pas garantir le respect du principe du maximum. Toutefois, pour
un écoulement permanent convergé, cette inéquation est toujours satisfaite puisque la varia-
tion temporelle de la densité est nulle. Ainsi, à convergence, le principe du maximum sera
respecté lorsqu'il s'applique.
 Conclusion

Dans cette thèse de doctorat, on a présenté un nouveau schéma numérique de résolu-

tion d'écoulements incompressibles dont la discrétisation spatiale est de type volumes finis
classiques. Les éléments géométriques qui discrétisent le domaine sont des tétraèdres en 3D
ou des triangles en 2D. Les équations de bilan sont intégrées sur des volumes de contrôle com-
posés d'un ou plusieurs éléments géométriques. Les fonctions d'interpolation pour toutes les
variables du système d'équations sont constantes par volume de contrôle. Les équations de
conservation de la masse et de la quantité de mouvement sont découplées et résolues avec
une méthode à pas fractionnaires. Toutes les équations sont linéarisées et la discrétisation
temporelle est semi-implicite. Puisque ce schéma supporte l'utilisation de volumes de con-
trôle arbitraires, il se distingue des autres schémas de la famille des volumes finis classiques
qui nécessitent des maillages de type Voronoï ou structuré pour assurer un discrétisation con-
sistante des flux diffusifs.
À partir des résultats théoriques de Eymard et al. ([36] et [39]), on peut affirmer que
ce schéma est convergent et que le principe du maximum est localement respecté lorsqu'il
s'applique. De plus, puisqu'il s'agit d'une méthode de volumes finis, les flux sont locale-
ment conservés. Les principales qualités de ce schéma sont sa robustesse, sa flexibilité et sa
simplicité.
Le schéma a été utilisé pour résoudre des problèmes de type académique où les diffi-
cultés sont bien identifiées. Les résultats obtenus sont en accord avec les solutions analytiques
ou des résultats numériques obtenus par d'autres chercheurs lorsque les solutions analytiques
ne sont pas connues.
Afin que le schéma puisse être développé par d'autres chercheurs, une librairie infor-
matique a été mise en oeuvre et est disponible par internet à l'adresse web suivante :
http ://[Link]/~sboivin/sperron.
Avec cette librairie, il est facile de construire un programme pour résoudre des écoulements
laminaires ou turbulents, en 2D et 3D, où les propriétés physiques telles la viscosité ou la
densité peuvent être variables. De plus, la librairie supporte la résolution d'un nombre illim-
ité d'équations de transport. L'écoulement peut donc transporter plus d'une quantité scalaire.
 160

Finalement, des programmes de pré et post-traitement sont incorporés à la librairie. Ces pro-
grammes facilitent l'importation de maillages et l'exportation des résultats pour d'autres pro-
grammes de post-traitement.
Bien que le schéma numérique développé dans le cadre de cette thèse comporte de
grands avantages, il comporte aussi certains désavantages. Premièrement, la méthode à pas
fractionnaires peut ne pas être efficace (lente) pour résoudre des écoulements permanents où
le découplage des équations n'est pas nécessaire. Deuxièment, le gain en stabilité est fait au
détriment de la précision, le schéma étant d'ordre 1. Ce schéma étant d'un ordre peu élevé,
on peut présumer que la précision de la solution pourrait être plus fortement influencée par la
qualité du maillage que pour un schéma d'un ordre supérieur. Finalement, ce schéma est mal-
adapté aux maillages où de nombreux éléments sont très étirés. En effet pour ces maillages,
il existe de nombreuses interfaces à transmittivité négative : pour de tels maillages on ne peut
pas garantir la convergence du schéma vers la vraie solution du problème.
On croit qu'il serait possible d'utiliser une variante de ce schéma pour la résolution
d'écoulements compressibles à haute vitesse. Pour ce faire, il faudrait revoir les conditions
d'application des résultats théoriques publiés par Eymard et al.[36] et revoir les étapes de
l'algorithme de résolution. En effet, pour les écoulements compressibles, l'étape de projection
n'est pas nécessaire.
Pour cette thèse, seuls des problèmes de type académique ont été résolus. Mais,
compte tenu des solides fondements sur lesquels est basé notre schéma, des problèmes plus
complexes pourront fort possiblement être résolus avec celui-ci.
 Bibliographie

[1] D.A. Anderson, R.H. Pletcher, et J.C. Tannehill. Computational Fluid Mechanics and
Heat Transfer. Taylor and Francis, 1997.
[2] L.V. Atkinson, O.J. Harley, et J.D. Hudson. Numerical Methods with Fortran 77 : a
Practical Introduction. Addison and Wesley, 1990.
[3] B.R. Baliga et S.V. Patankar. A new finite element formulation for convection-diffusion
problems. Numerical Heat Transfer, Part B, 3 :393^K)9,1980.
[4] B.R. Baliga et S.V. Patankar. A control volume finite-element method for two-
dimensional fluid flow and heat transfer. Numerical Heat Transfer, Part B, 6 :245-261,
1983.
[5] B.R. Baliga et S.V. Patankar. Solution of some two-dimensional incompressible fluid
flow and heat transfer problems using a control volume finite-element method. Numer-
ical Heat Transfer, Part B, 6 :263-282,1983.
[6] F. Benkhaldoun, J.P. Chabard, et G. Pot. Projet N3S de mécanique des fluides. Résolu-
tion par volumes finis de l'étape de transport pour des problèmes d'écoulements turbu-
lents incompressibles. Rapport technique HE-41/89.25, Électricité de France, Direction
des Études et Recherche, 1988.
[7] F. G. Blottner. Numerical solution of convection-diffusion equations. Computers and
Fluids, 6:15-24,1978.
[8] S. Boivin. MEFOO : Un framework C++ pour la mise-en-oeuvre de la méthode des
éléments finis. Rapport technique HE-41/96/003/A, EDF, 1996.
[9] S. Boivin, F. Cayré, et J.M. Hérard. A finite volume method to solve the Navier-Stokes
equations for incompressible flows on unstructured meshes. International Journal of
Thermal Science, 39 :806-825,2000.
 162

[10] S. Boivin, F. Cayré, et J.M. Hérard. A finite volume scheme to compute incompressible
gas-solid two-phase flows. American Institute of Aeronautics and Astronautics, AIAA-
2000-2665, 2000.

[11] S. Boivin et M. Fortin. A new artificial viscosity method for compressible viscous flow
simulations by FEM. International Journal of Computational Fluid Dynamics, 1,1993.
[12] S. Boivin et M. Fortin. A nonisotropic artificial viscosity method : application to the
simulation of compressible viscous flows. International Journal of Computational Fluid
Dynamics, 7, 1996.

[13] S. Boivin et J.M. Hérard. Un schéma de volumes finis pour résoudre les équations de
Navier-Stokes sur une triangulation. Revue Européenne des Éléments Finis, 5 :461^490,
1996.

[14] Yves Bourgault. Méthodes d'Éléments Finis en Mécanique des Fluides, Conservation
et Autres Propriétés. Thèse de doctorat, Université Laval, 1996.
[15] D. S. Burnett. Finite Element Analysis. Addison and Wesley, 1988.

[16] S. Candel. Mécanique des fluides. Dunod, 1995.

[17] B. Cardot. Modélisation de la Turbulence par les Méthodes de Type k-epsilon et Ho-
mogénéisation. Thèse de doctorat, Université de Paris 6,1990.
[18] F. Cayré. Schémas volumes finis pour un problème elliptique sur un maillage triangu-
laire, étude numérique de convergence. Rapport technique HE-41/97/057/A, Électricité
de France, Direction des Études et Recherche, 1997.

[19] F. Cayré. Méthodes de Volumesfinispour Maillages Non Structutés pour la Simulation

Numérique des Écoulements Incompressibles Monophasiques et Diphasiques. Mémoire
de maîtrise, Université Laval, 1999.

[20] J.P. Chabard. Projet N3S de mécanique des fluides, manuel théorique -version 2.0. Rap-
port technique HE-41/88.09, Électricité de France, Direction des Études et Recherche,
1988.

[21] S. Champier et T. Gallouet. Convergence d'un schéma décentré sur un maillage triangu-
laire pour un problème hyperbolique linéaire. R.A.I.R.O. Mathematical Modelling and
Numerical Analysis, 26 :835-853,1992.
 163

[22] K.C. Chang, C.S. Chen, et C.I. Uang. A hybrid k-e turbulence model of recirculating
flow. International Journal for Numerical Methods in Fluids, 12 :369-382,1991.
[23] K.H. Chen et R.H. Pletcher. Primitive variable, strongly implicit calculation procedure
for viscous flows at all speed. AIAA, 29 :1241-1249,1993.
[24] A.J. Chorin. A numerical method for solving incompressible viscous flow problems.
Journal of Computational Physics, 2 :12-26,1967.
[25] A.J. Chorin. Numerical solution of the Navier-Stokes equations. Mathematics of Com-
putation, 22 :745-762, 1968.
[26] W.L. Chow et CM. Rhie. Numerical study of the turbulent flow past an airfoil with
trailing edge separation. AIAA Journal, 21 :1525-1532,1983.
[27] P.G. Ciarlet. Introduction à l'analyse numérique matricielle et à l'optimisation.
DUNOD, 1998.
[28] Travail collectif. Special Course on Unstructured Grid Methods for Advection Domi-
nated Flow. AGARD, Neuilly Sur Seine, 1992.
[29] M. Crouzeix et P.A. Raviart. Conforming and non-conforming finite elements methods
for solving the stationary Stokes equations. RAIRO, Série Rouge Analyse Numérique,
3 :33-76,1973.
[30] G. De Vahl Davis. An evaluation of upwind and central difference approximations by a
study of recirculation flow. Computers and Fluids, 4 :24—43,1976.
[31] G. De Vahl Davis. Natural convection of air in a square cavity : a benchmark solution.
International Journal for Numerical Methods in Fluids, 3 :249-264,1983.
[32] G. De Vahl Davis. Natural convection of air in a square cavity : a comparaison exercise.
International Journal for Numerical Methods in Fluids, 3 :227-248, 1983.
[33] J. Doormal et G. Raithby. Enhancements of the SIMPLE method for predicting incom-
pressible fluid flows. Numerical Heat Transfer, Part B, 1:147-163,1984.
[34] C. Cuvelier et al. Éléments d'équations aux dérivées partielles pour ingénieurs. Presses
polytechniques romandes, 1988.
[35] A. Brooks et TJ.R. Hughes. A multi-dimensional upwind scheme with no crosswind
diffusion. Proceedings of a Symposium on Finite Element Methods for Convection Dom-
inated Flows, ASME Winter Annual Meeting, pages 19-35, 1979.
 164

[36] R. Eymard, T. Gallouet, et R. Herbin. Finite Volume Methods. Hanbook of Numerical

Analysis, 7,2001.
[37] [Link] et J. Welsh. Numerical calculation of time dependant viscous incompressible
flow with free surface. Physics of Fluids, 8 :2182-2189,1965.
[38] A. Fortin. Analyse numérique pour ingénieurs. Éditions de l'École Polytechnique de
Montréal, 1996.
[39] T. Gallouet, R. Herbin, et M.H. Vigual. Error estimates on the approximate finite volume
solution of convection diffusion equations with Dirichlet, Neumann or Robin boundary
conditions. SIAM, à être publié.
[40] U. Ghia, K. N. Ghia, et C. T. Shin. High-Re solutions for incompressible flow using the
Navier-Stokes equations and a multigrid method. Journal of Computational Physics,
48:387-411,1982.
[41] A. Greenbaum. Iterative Methods for Solving Linear Systems. SIAM, 1997.
[42] P. Hood et C. Taylor. A numerical solution of the Navier-Stokes equations using finite
element technique. Computers and Fluids, 1 :73-100,1973.
[43] K.H. Huebner et E.A. Thornton. The Finite Element Method for Engineers. Wiley and
Sons, 1982.
[44] S. Hussain, C. Ketleborough, et C. Prakash. Solution of fluid flow problems with the
vorticity-stream function formulation and control-volume-based finite-element method.
Numerical Heat Transfer, PartB, 16 :31-58,1989.
[45] F.P. Incropera et D.P. DeWitt. Introduction to Heat Transfer. Wiley and Sons, 1990.
[46] [Link] et P. Moin. Application of a fractional-step method to incompressible Navier-
Stokes equations. Journal of Computational Physics, 59 :308-323,1985.
[47] P. Joyal et A. Mercier. Analyse vectorielle. 1998.
[48] R. Kessler, M. Perie, et G. Scheuer. Comparaison of finite-volume numerical methods
with staggered and colocated grids. Computers and Fluids, 16 :389-403,1988.
[49] L. Landau et E. Lifchitz. Mécanique des fluides. Wiley and Sons, 1994.
[50] B.E. Launder et B.I. Sharma. Application of the energy dissipation model of turbulence
for the calculation of flow near a spinning disc. Letters in Heat and Mass Transfer,
2:131-138,1974.
 165

[51] B. Van Leer. Towards the ultimate conservative difference scheme, ii : Monotocinity and
conservation combined in a second-order scheme. Journal of Computational Physics,
32 .-101-136,1979.

[52] B.P. Leonard. A stable and accurate convective modelling procedure based on quadratic
upstream interpolation. Computational Methods in Applied Mechanical Engineering,
19 :59-98,1979.

[53] B.P. Leonard. A survey of finite differences of opinion on numerical muddling of the
incompressible defective confusion equation. Finite Element Methods for convective
Dominated Flows, AMD-vol.34, 1979.

[54] R.J. LeVeque et J. Oliger. Numerical Methods based on additive splittings for hyper-
bolic partial differential equations. Rapport technique NA-81-16, Computer Science
Département, Stanford University, Stanford, California 94305,1981.

[55] J.M. Lopez et J. Shen. Numerical simulation of incompressible flows in cylindrical

geometries using a spectral projection method. International Journal ofApplied Science
Computing, 5 :25^tO, 1998.

[56] B. Makiola. Experimented Unter-suchungen Zur Strmungber Die Schrge Stufe.

Thèse de doctorat, Institut Fuer Hydromechanik, Universitaet Karlsruhe Kaiserstr, 1992.
http ://[Link].

[57] B. Mohammadi et O. Pironneau. Analysis of the K-Epsilon Turbulence Model. Wiley

and Sons, 1994.

[58] S. Patankar et D. Spalding. A calculation procedure for heat, mass and momentum
transfer in three-dimensional parabolic flows. International Journal of Heat and Mass
Transfer, 15 :1787-1806,1972.

[59] S. V. Patankar. Numerical Heat transfer and Fluid Flow. McGraw-Hill, 1980.

[60] C. Pozrikidis. Introduction to Theoretical and Computational Fluid Dynamics. Oxford

university Press, 1997.
[61] L. Quartarpelle. Numerical Solution of the Incompressible Navier-Stokes Equations.
Birkhauser, 1993.
 166

[62] G. D. Raithby. A critical evaluation of upstream differencing applied to problems in-

volving fluid flow. Computational Methods in Applied Mechanical Engineering, 9 :75-
103,1976.
[63] G. D. Raithby. Skew upstream differencing schemes for problems involving fluid flow.
Computational Methods in Applied Mechanical Engineering, 9 :153-164,1976.
[64] Y. Saad et M. H. Schultz. A generalized minimal residual algorithm for solving non-
symmetric linear systems, research report yaleu/dcs/rr-254. SIAM, 7:856-869,1986.
[65] Y. Saad et H. A. Van Der Vorst. Iterative solution of linear systems in the 20th cen-
tury. Rapport technique umsi-99-152, Minnesota Supercomputer Institute, University
of Minneapolis, 1999.
[66] M. Schafer et S. Turek. Benchmark computations of laminar flow around cylinder. In
Proceedings DFG Priority Research Program, Flow Simulation on High Performance
Computers. Vieweg, 1992.
[67] J. Shen. On error estimates on projection methods for Navier-Stokes equations : first
order schemes. SIAM Journal Numerical, 29 :55-77, 1992.
[68] J. Shen. A remark on the projection-3 method. International Journal for Numerical
Methods in Fluids, 16 :249-253,1993.
[69] J. Shen. On error estimates on projection methods for the Navier-Stokes equations :
second order schemes. Mathematics of computation, 65 :1039-1065,1996.
[70] D. B. Spalding. A novel finite difference formulation for differenctial expressions in-
volving both first and second derivatives,. International Journal for Numerical Methods
in Engineering, 4 :551-559,1962.
[71] C.G. Speziale. On nonlinear K-l and K-e models of turbulence. Journal of Fluid Me-
chanics, 178 :459-475,1987.
[72] P. Wesseling. Principles of Computational Fluid Dynamics. Springer-Verlag, 2000.
[73] F. M. White. Viscous Fluid Flow. McGraw-Hill, 1991.
[74] D.C. Wilcox. Turbulence Modeling for CFD. DCW Industries, 1993.
[75] U.C. Yee. Upwind and symmetric schock capturing methods. Rapport technique TM-
89464, NASA, 1987.