29 septembre 2006

Rapport de stage ingénieur

5 Juin - 29 Septembre

Création d’une base de données de

Travaux pratiques en calcul numérique

Estelle Etienne

Encadrement : Marc Montagnac.

CERFACS - CFD - 42, avenue Gaspard Coriolis - 31057 Toulouse Cedex 1
Remerciements

Je tiens en premier lieu à remercier Jean-François Boussugue, notre chef d’équipe et Marc
Montagnac, mon maître de stage pour m’avoir accueillie au CERFACS.
Merci aussi à Yann Colin qui a toujours su répondre à mes questions avec compétence et
patience, merci à Guillaume Puigt pour sa gentillesse et sa disponibilité, merci à Antoine Devesa,
mon collègue de bureau, qui m’a suivie dans mon apprentissage de Linux... enfin merci à toute
l’équipe pour la bonne ambiance qui y règne.
Je tiens aussi à remercier mes collocataires pour ces quatre mois agréables passés à Toulouse.
Bonne continuation à vous tous !

2
Table des matières

Introduction 7

1 Moyens mis en oeuvre 8

1.1 Site WIKI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.1 support WIKI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.1.2 Architecture du site . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.2 Code de calcul elsA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.2.1 Script Python . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Déroulement d’un calcul type 12

2.1 TP squelette . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Analyse du cas physique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Pré-traitement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.1 Conditions initiales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.2 Conditions aux limites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.3 Adimensionnement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.5 Simulation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6 post-traitement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3 Points d’étude 20
3.1 Modélisation de la turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.1.1 Modèle de Spalart-Allmaras . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.2 Modèle k-ω . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.3 Couche limite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.1.4 Phénomènes physiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Influence des schémas numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2.1 Schéma centré de Jameson . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.2 Schéma décentré de Roe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3 Techniques d’accélération de la convergence . . . . . . . . . . . . . . . . . . . . . . 34
3.3.1 Méthode explicite / implicite . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3
 3.3.2 Méthode multi-grille . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Conclusion 38

4
Table des figures

1.1 page d’accueil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1 TP squelette . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Conditions initiales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3 Géométrie du profil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 configuration du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.5 Répartition des points sur le profil . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6 repère propre au maillage. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.7 Maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1 Plaque plane (modèle de Spalart) - Loi de paroi . . . . . . . . . . . . . . . . . . . . 23

3.2 Plaque plane - couche limite laminaire et turbulente. Effet du choix du modèle de
turbulence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.3 RAE2822 (cas 9) - coefficient de pression. . . . . . . . . . . . . . . . . . . . . . . . 25
3.4 RAE2822 (cas 9) - coefficient de frottement. . . . . . . . . . . . . . . . . . . . . . . 25
3.5 RAE2822 (cas 10) - vecteurs vitesse suivant une ligne du maillage proche du pied
du choc. Visualisation du retournement des vecteurs vitesse (décollement) . . . . . 26
3.6 RAE2822 (cas 10) - Coefficient de pression. . . . . . . . . . . . . . . . . . . . . . . 27
3.7 RAE2822 (cas 10) - Coefficient de frottement. . . . . . . . . . . . . . . . . . . . . . 27
3.8 NACA - position du choc pour α = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.9 NACA - position du choc pour α = 4 . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.10 Polaire obtenue par calcul Euler : effets de la viscosité numérique d’elsA. . . . . . . 29
3.11 détails du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.12 schéma - détail du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.13 schémas numériques à t=0s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.14 schémas numériques à t=2s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.15 schémas numériques à t=5s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.16 schémas numériques à t=7s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.17 2 niveaux de grilles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.18 comparaison des résidus avec et sans méthode multi-grille . . . . . . . . . . . . . . 37

5
Liste des notations

CFD : Computational Fluid Dynamics

CFL : Courant-Friedrichs-Levy
elsA : Ensemble Logiciel de Simulation Aérodynamique
Rij : tenseur de Reynolds en kg.m.s−2
µt : viscosité turbulente en kg.m−1 .s−1
k : energie cinétique turbulente en m2 .s−2
ω : fréquence caractéristique de la turbulence en s−1
τ : tenseur des contraintes
κ : la constante de Von Karman (sans dimension) qui vaut 0.41
W : le vecteur des variables conservatives
Y : une variable conservative quelconque
α : la vitesse d’advection
f : le flux de Y

6
Introduction

Le stage ingénieur, dont le travail est décrit dans ce rapport, a été effectué au sein de la
société CERFACS (Centre de Recherche et de Formation Avancée en Calcul Scientifique) dans
le département Computational Fluid Dynamics (CFD) au sein de l’équipe AAM (Advanced Ae-
rodynamics and Multiphysics). L’entreprise, ses activité principales, son organisation et celle de
l’équipe AAM sont présentées dans l’annexe A.
Le sujet du stage, "création d’une base de données de Travaux Pratiques en calcul numérique"
s’inscrit dans un projet interne au Cerfacs. Une série de trois Travaux Pratiques avait déjà été
mise au point mais les sujets ne suivaient pas une architecture commune.
Cette base de données a pour but de proposer un ensemble de sujets construits sur le même
modèle. La première étape a été de définir la base de tous les sujets : un TP squelette, qui servira
de canneva pour la construction de tous les TP. Le choix d’un plan chronologique a été fait. A
chaque étape, figurent les points théoriques qui interviennent. Ce sont ces différents points que les
travaux pratiques s’attachent à illustrer. Ils apparaissent ou non dans l’énoncé selon les objectifs
du sujet.
Ces Travaux Pratiques sont à vocation pédagogique. Ils sont destinés à initier les élèves de
troisième année de l’ENSEEIHT à la simulation numérique avec le code de calcul elsA utilisé chez
Airbus. Dans un premier temps, le but est de leur faire intégrer le déroulement d’un calcul CFD et
ensuite d’aborder certaines notions importantes de la modélisation et des méthodes numériques.
Cette base de données a été mise en ligne sur un site internet de type WIKI (www.cerfacs.fr/wiki/cfdwiki),
ce qui offre la possibilité de l’enrichir à tout moment. Vous pouvez aller consulter ce site pour mieux
vous rendre compte du travail effectué durant le stage.

7
Chapitre 1

Moyens mis en oeuvre

1.1 Site WIKI

1.1.1 support WIKI
Principe

PmWiki est un système de type wiki pour la création et l’entretien collectif de sites Internet.
Un WikiWikiWeb est un "système ouvert d’édition" où l’accent est mis sur l’écriture et la
collaboration à l’écriture de documents plutôt que la lecture rapide ou leur visualisation simple. Le
mot "wiki" provient du terme hawaïen "wiki wiki" qui signifie "informel" ou "rapide". Le concept
de base d’un WikiWikiWeb (ou "wiki") est que (presque) tout le monde peut éditer n’importe
quelle page. L’objectif de ce système est de rendre la création et l’édition de contenus en ligne
aussi facile que possible
Les pages PmWiki ont le même aspect et fonctionnent comme des pages Internet ordinaires,
sauf qu’elles possèdent un lien pour "éditer" ce qui permet de modifier ou d’ajouter facilement
des pages à un site, en utilisant les règles d’édition de base. Il n’est pas nécessaire de connaître
le langage HTML. Il existe notamment un émulateur Latex pour faciliter l’écriture des formules
mathématiques. L’édition des pages peut être laissée ouverte à tout public ou restreinte à un
groupe d’utilisateurs.

Pourquoi un site Wiki ?

Le choix d’un support Wiki a été fait déjà pour permettre un travail collaboratif durant le
stage et pour donner l’opportunité à ce site d’évoluer et de s’enrichir. Par la suite, ce site sera
ouvert à d’autres utilisateurs tels que l’ONERA. D’autres sujets de TP pourront être proposés et
rédigés en se servant des travaux déjà effectués. Ils bénéficieront de tout l’environnement du site.

1.1.2 Architecture du site

Le site a été pensé de manière à privilégier la fonctionnalité. Il est organisé sous forme
d’arborescense vers différentes pages auxquelles on accède grâce à des liens. On trouve sur la page
d’accueil du site plusieurs liens :
– la liste des TP à laquelle on accède en ouvrant le lien « TP calcul numérique »
– deux autres parties qui sont détaillées ci-dessous et qui sont disponibles depuis la barre
d’outils : des rappels théoriques et une « boîte à outil ».Ce sont des parties qui encombre-

8
 raient l’énoncé et qui ne sont pas forcément nécessaires aux étudiants. Elles permettent
d’obtenir rapidement une aide ponctuelle durant la séance.
– un lexique qui recense les mots difficiles et en donne une brève définition.

Rappels théoriques

On y trouve pour l’instant deux parties, les équations de base et une partie théorique sur
la couche limite, mais au fur et à mesure de l’évolution du site, de nouveaux points théoriques
y seront traités. La partie « Equation de base » rappelle les équations d’Euler, Navier-Stokes, et
Navier-Stokes moyennées (RANS). Ce sont ces équations que résout le code de calcul elsA utilisé
dans les TP. Ce code de calcul est présenté dans la partie suivante.

Outils

Il s’agit d’un inventaire des logiciels que les étudiants vont être amenés à utiliser, avec pour
chaque logiciel, un lien vers un pseudo guide d’utilisation orienté spécifiquement pour ce qui est
demandé en TP. Y figurent les logiciels de visualisation Tecplot et Quickview, l’outil EDM (Editor
of Model) qui permet de générer les fichiers nécessaires à un calcul elsA et de convertir les résultats
du calcul au format souhaité, le code de calcul elsA où l’on trouve la démarche à suivre pour lancer
un calcul et enfin un guide Unix qui répertorie les commandes de base et qui présente l’éditeur de
texte « vi ».

1.2 Code de calcul elsA

ElsA (Ensemble Logiciel de Simulation Aérodynamique) est un logiciel de simulation des
écoulements compressibles, laminaires ou turbulents, en 3D, sur des maillages structurés en blocs.
Ce code de calcul a été développé par l’ONERA depuis l’année 1997 et en collaboration avec le
Cerfacs depuis 2000. ElsA résout les équations de Navier-Stokes moyennées (RANS pour Reynolds
Average Navier-Stokes), dans leur forme intégrale, dans un référentiel fixe ou mobile. La turbulence
peut être modélisée par des modèles algébriques ou avec équations de transport. Des simulations
stationnaires et instationnaires peuvent être réalisées. Les utilisateurs interagissent avec elsA grâce
à un script construit en Python qui sert d’interface.

1.2.1 Script Python

C’est dans ce script que les utilisateurs d’elsA définissent leurs paramètres de calcul. Une
connaissance très limité du language Python suffit pour mettre au point un script. Python est
un langage moderne, orienté objet. Sa syntaxe est facile à lire. En général, on récupère un script
Python que l’on modifie selon les paramètres numériques que l’on souhaite changer ou rajouter.
Un script Python basique s’articule autour de plusieurs grandes parties :
– Création du problème auquel on donne un nom et dont on spécifie la configuration (pro-
blème 2D, 3D).
– Création du modèle physique : cette partie précise tous les paramètres qui définissent le
système d’équations résolu par elsA.On y entre le système d’équation (Euler, Navier-Stokes
laminaire ou turbulent), les caractéristiques du fluide (γ, Cv), les valeurs adimensionnées
µ∞ et T∞ nécessaires au calcul de la viscosité selon la loi de Sutherland, et le modèle de
turbulence souhaité.
– Création de l’algorithme numérique : cette partie précise les paramètres qui définissent
comment le système d’équations va être résolu. On y entre la méthode de discrétisation
en espace (schéma numérique) et d’intégration en temps (nombre d’itérations, méthode
d’accélération de la convergence).

9
Fig. 1.1 – page d’accueil

10
 – Etat amont : On défini le champ infini. On donne les valeurs adimensionnées des variables
conservatives et des quantités turbulentes suivant le modèle de turbulence choisi. On ini-
tialise généralement l’état initial avec les valeurs de l’état infini.
– Création et initialisation des blocs : c’est là que l’on fournit le maillage.
– Conditions aux limites : celles-ci doivent être précisées pour chaque frontière.
– Extraction des résultats : on précise dans cette partie les variables à extraire et le nom des
fichiers dans lesquels on souhaite les mettre.
Il ne s’agit bien évidemment pas là d’une liste exhaustive des possibilités d’un script Python pour
ce qui est de la définition des conditions d’un calcul, mais, sont cités les principaux paramètres
que l’on a été amené à modifier pendant ce stage.

11
Chapitre 2

Déroulement d’un calcul type

2.1 TP squelette
Le TP squelette est, comme il a été dit dans l’introduction, le modèle sur lequel sont et
seront construits tous les autres sujets de TP. Son plan doit pouvoir convenir pour n’importe quel
cas test. C’est pourquoi on a choisi un plan chronologique qui permet aussi de bien mettre en
évidence le déroulement d’un calcul en insistant sur chaque étape. Le TP squelette s’apparente à
un cours succint de CFD organisé de manière chronologique. En effet, à chaque étape, figurent les
points théoriques qui s’y rattachent, sous forme de lien. S’il font partie intrinsèque du déroulement
d’un calcul, ils apparaissent systématiquement dans l’énoncé du TP, sinon, s’il s’agit de points
se rapportant au numérique, ils apparaissent selon l’objectif du sujet. Le plan du TP squelette
(cf Figure 2.1) reprend donc les étapes d’un calcul, à savoir : l’analyse du cas physique, le pré-
traitement, le maillage, le calcul à proprement parler, et le post-traitement.

12
 Fig. 2.1 – TP squelette

Sur ce modèle là, ont déjà été traités trois cas tests : « Plaque plane », « Calcul Euler autour
d’un profil NACA 0012 » et « Ecoulement autour d’un profil RAE2822 ». Nous allons voir comment
chaque étape a été illustrée au sein des TP.

13
2.2 Analyse du cas physique
Cette partie permet de poser le problème. Elle sera plus ou moins importante selon la com-
plexité du cas. On y fixe d’une certaine manière l’objectif du TP : on référence dans cette partie
sous forme de mots clefs les points importants du TP afin d’orienter l’analyse des étudiants.

2.3 Pré-traitement
L’étape de pré-traitement sert à créer le script Python qui sera ensuite « lu » par elsA. On
calcule là les valeurs adimensionnées du champ initial et on fixe les conditions aux limites.

2.3.1 Conditions initiales

Sont données dans l’énoncé les paramètres nécessaires pour démarrer un calcul elsA. On
donne en général le nombre de Mach, le nombre de Reynolds et l’angle d’incidence. les grandeurs
turbulentes sont déjà initialisées dans le script Python. Le lien « Conditions initiales », précise
toutes les conditions initiales que l’on doit fournir à elsA (voir Figure 2.2).

Fig. 2.2 – Conditions initiales

14
2.3.2 Conditions aux limites
L’attribution des conditions aux limites est une étape primordiale, et qui peut ne pas être
évidente. Plusieurs possibilités peuvent convenir sans pour autant qu’elles donnent les même ré-
sultats. C’est pour initier les étudiants à ce choix qu’est donné une liste des conditions aux limites
usuelles. La configuration du maillage leur est fournie avec, en couleur, les frontières des blocs du
maillage pour lesquelles ils ont à choisir les conditions aux limites. Ils n’auront pas cependant à
entrer dans le script Python.

2.3.3 Adimensionnement
Pour résoudre les équations de Navier-Stokes, il est conseillé d’utiliser des quantités adimen-
sionnées. Il n’existe pas de système d’unité dans elsA. C’est donc à l’utilisateur de rentrer des
données utilisant un système d’unité cohérent.
La normalisation est basée sur le choix de 4 quantités indépendantes que l’on prend égales
à l’unité. La longueur de référence est toujours prise égale à 1. Il existe quatre possibilités
qui mènent aux normalisations SRPT (Static-Rho-Pressure-Temperature), SRVT (Static-Rho-
Velocity-Temperature), SRCT (Static-Rho-Celerity-Temperature) et TRCT
(Total-Rho-Celerity-Temperature).
Selon la normalisation choisie, on fixe certaines quantités et les autres sont calculées à l’aide
des relations isentropiques. Le fluide est considéré comme un gaz parfait. Est détaillée ici la nor-
malisation SRPT. On prend donc les trois variables suivantes égales à 1 :

ρ∗∞ = 1 , p∗∞ = 1 , ∗
T∞ = 1.

La loi des gaz parfaits donne la relation :

p∞ = ρrgaz T∞ ,

avec
rgaz = Cv(γ − 1) .
On en déduit :
p∞ = ρ∞ (γ − 1)CvT∞
et il vient
1
Cv = .
γ−1

D’après la relation
γ − 1 2
γ−1
−γ
p∞ = pt∞ 1 + M∞ ,
2
on déduit :
γ − 1 2
γ−1
γ
p∗t∞ = 1 + M∞ .
2
De même, pour Tt∞ , d’après la relation :
γ − 1 2
−1
T∞ = Tt∞ 1 + M∞ ,
2
on déduit :
γ−1 2
Tt∗∞ = 1 + M∞ .
2
D’après la relation : q
a∞ = (γrgaz T∞ ) ,

15
soit : p
a∞ = γ(γ − 1)CvT∞ ,
il vient :
√
a∗∞ = γ.
Il en découle :
∗ √
U∞ = M∞ a∗∞ = M∞ γ .
D’après la relation :
1 2
et∞ = Cv T∞ + U∞ ,
2
on déduit :
1 1 2
e∗t∞ = + M∞ γ.
γ−1 2
D’après la relation :
ht∞ = γCvTt∞ ,
il résulte :
γ γ−1 2
h∗t∞ = (1 + M∞ ) .
γ −1 2
Comme on connaît le nombre de Reynolds, et que l’on a la relation :
ρ∞ U ∞
µ∞ =
Re∞
il vient : √
M∞ γ
µ∗∞ = .
Re∞
On accède aux valeurs des variables conservatives ρu, ρv et ρw grâce aux relations suivantes :

ρu = ρ∞ U∞ cos(α)cos(β)

ρv = ρ∞ U∞ sin(β)
ρw = ρ∞ U∞ sin(α)cos(β)

Cet exercice est proposé au début de chaque TP. Les étudiants ont un script Python à remplir.
Dans un souci de facilité, les scripts ont été modifiés, on a créé de nouvelles variables qui renvoient
aux valeurs adimensionnées à entrer et que l’on a placées en début de script.

2.4 Maillage
Le maillage sera, dans la majorité des cas, fourni. C’est une étape longue qui alourdirait
considérablement les TP si elle devait être répétée à chaque séance. Il est cependant souhaitable
que les étudiants en aient l’expérience pour apprendre à ne pas négliger cette étape et se rendre
compte de l’impact d’un maillage plus ou moins bien fait sur le calcul.
Le sujet « calcul Euler autour d’un profil NACA 0012 » est un TP qui s’échelonnera sur
plusieurs séances dont une consacrée au maillage. Le mailleur utilisé s’appelle Megacads. Les
étudiants ne connaissant pas le logiciel, toutes les étapes sont détaillées. Voilà, en général, comment
s’effectue la réalisation d’un maillage :
• On choisit dans quel repère on va travailler. On importe les données géométriques du profil,
un ensemble de points, par lequel on fait passer une courbe. On a ainsi créé la géométrie du
profil (cf. Figure 2.3). On définit ensuite le domaine de calcul en créant les frontières et on crée à
l’intérieur les blocs du maillage qui seront discrétisés (cf. Figure 2.4).

16
 • On répartit les points du maillage aux frontières. On choisit une distribution de Poisson, de
manière à avoir une discrétisation plus fine près du profil. On peut choisir la longueur des premier
et dernier espacements. Pour les blocs qui ne sont pas rectangulaires, on joue sur cette valeur
pour obtenir le maillage le plus orthogonal possible, ce qui est le gage d’un maillage de qualité (cf.
Figure 2.5).
• Le maillage à créer est un maillage structuré : un sommet du maillage est repéré par un
triplet (i,j,k), ce qui impose une conservation du nombre de maille sur les faces « parallèles » des
blocs. Pour construire le maillage, on fonctionne comme dans un repère, en reliant les points des
frontières qui se font face. (cf. Figures 2.6 et 2.7).
Ce maillage sera ensuite utilisé pour le calcul, ce qui permettra de voir (un peu trop tard !) si
il était ou non bien fait. Pour ne pas fausser les résultats attendus dans la suite du TP, un maillage
dont on est sûr est fourni. Ce qui simplifiera aussi la démarche de calcul. En effet le passage d’un
maillage Mégacads à un fichier d’entrée elsA est long et présente de nombreuses occasions de se
tromper. C’est pourquoi, on a mis au point deux scripts qui reprennent les commandes à effectuer.

Fig. 2.3 – Géométrie du profil

Fig. 2.4 – configuration du maillage

17
 Fig. 2.5 – Répartition des points sur le profil

Fig. 2.6 – repère propre au maillage.

Fig. 2.7 – Maillage

2.5 Simulation numérique

La simulation numérique consiste à lancer le calcul elsA. On a créé un lien « Comment lancer
un calcul elsA ? » vers une page qui donne les instructions pour démarrer le calcul. On y fournit
le fichier batch ci-dessous, qui permet de lancer le calcul et dans lequel les étudiants auront à

18
préciser le script Python qu’ils exécutent. Ce fichier batch est un script de commande qui permet
de se « logger » sur le calculateur et d’avoir accès aux processeurs de calcul. Les processeurs
sont attribués selon un système de queue. Nos calculs sont assez courts et s’effectent sur un seul
processeur, on se met donc sur la queue « dev » pour développement. L’avantage de cette queue
est la rapidité de passage, mais le temps alloué par calcul est faible, aux alentours de 20 minutes.

#!/bin/tcsh
# BSUB -n 1
# BSUB -o sortie%J.txt
# BSUB -e sortie%J.txt
# BSUB -q dev
# BSUB -J naca
echo ’------------------------------------------------------------’
setenv MPI_USE_LIBELAN_SUB 0
cd $LS_SUBCWD
source ~etienne/ELSA/setupelsa.csh I3203g dec_r8_so ~etienne/ELSA

elsa -f naca01-1deg_V.epy
echo ’------------------------------------------------------------’

2.6 post-traitement
Le post-traitement consiste à visualiser les résultats. On utilise pour cela le logiciel Tecplot.
Pour les visualisations un peu compliquées, une aide Tecplot est accessible depuis la barre d’outils.
Ce lien répertorie suivant les TP les visualisations à faire et guide l’étudiant. On cherche à valider
les résultats numériques par comparaison avec des données expérimentales, ou à observer des
phénomènes physiques.

19
Chapitre 3

Points d’étude

L’objectif de ce chapitre est de présenter succintement l’ensemble des points théoriques abor-
dés et de montrer comment ils ont été illustrés.

3.1 Modélisation de la turbulence

elsA permet de faire des calculs non-visqueux (équations d’Euler) et visqueux. La turbulence ;
qui se rencontre au quotidient (panaches de fumée de cheminée d’usine,...), est un phénomène
physique très complexe. Pour les écoulements visqueux turbulents, elsA dispose d’un grand choix
de modèles de turbulence conçus grâce à une approche RANS (Reynolds Averaged Navier-Stokes).
Le principe de cette approche est de séparer les champs conservatifs en une partie moyenne et une
partie fluctuante. La partie moyenne est obtenue par application d’un filtre :
– le filtre de Reynolds pour la pression p et la masse volumique ρ,
– le filtre de Favre écrit à partir de la moyenne de Reynolds pour toute les autres grandeurs
physiques.
Avec cette approche, de nouveaux termes issus du passage à la moyenne apparaissent et doivent
etre modélisés. C’est cette modélisation qui fait apparaître la viscosité turbulente µt dont les
effets viennent s’ajouter à ceux de la viscosité laminaire (ou moléculaire). Les modèles à viscosité
turbulente (modèles du premier ordre) sont basés sur l’hypothèse de Boussinesq qui consiste à
exprimer les tensions de Reynolds en fonction de la viscosité turbulente µt que l’on cherche à
déterminer : !
3
∂ui ∂uj 2 X ∂uk 2
Rij = µt + − δij ( ) − ρkδij
∂xj ∂xi 3 ∂xk 3
k=1

avec Rij , le tenseur de Reynolds :

 
u0 u0 u0 v 0 u0 w 0
Rij =  u0 v 0 v0 v0 v 0 w0 
u0 w 0 v 0 w0 w0 w0

et k énergie cinétique turbulente.

De plus, un coefficient de diffusivité thermique turbulent λt permet de prendre en compte les
effets de la turbulence sur la thermique de l’écoulement. L’évolution du coefficient de diffusivité
thermique de la turbulence est liée aux variations de µt . Pour déterminer µt , plusieurs modèles
sont disponibles dont plusieurs sont implémentés dans elsA :
– les modèles algébriques qui expriment directement µt en fonction des variables conserva-
tives,

20
 – les modèles à une seule équation de transport, comme le modèle de Spalart-Allmaras,
permettent de calculer directement l’évolution de µt en fonction des champs conservatifs,
– les modèles à 2 équations de transport, comme le modèle k − ω de Wilcox calculent 2
échelles de la turbulence et ces 2 échelles permettent d’obtenir µt par relation algébrique.
Dans le cadre des TP, seuls les modèles de Spalart-Allmaras et k − ω de Wilcox ont été présentés.

3.1.1 Modèle de Spalart-Allmaras

Le modèle de Spalart-Allmaras est un modèle à 1 équation. Les modèles à 1 équation pro-
posent soit une équation de transport liée à µt soit à une autre grandeur de la turbulence direc-
tement liée à la viscosité turbulente.

3.1.2 Modèle k-ω

Le modèle k −ω est un modèle à deux équations de transport. Pour ces modèles là, la viscosité
turbulente s’exprime en fonction de deux grandeurs de la turbulence. On écrit une équation de
transport pour chacune de ces deux grandeurs, ici :
– pour k, l’énergie cinétique turbulente,
– pour ω, fréquence caractéristique de la turbulence. L’inverse de ω représente l’échelle de
temps caractéristique de la dissipation de l’énergie cinétique k. La viscosité turbulente
s’exprime en fonction de k et ω :
k
µt = ρ
ω

3.1.3 Couche limite

La couche limite est une zone située au voisinage d’un corps immergé dans un fluide en
mouvement, en dehors de laquelle, les effets de la viscosité sont négligeables. Elle est le siège de
forts gradients des champs conservatifs. C’est la capacité à plus ou moins bien prédire cette zone
de l’écoulement qui permet de décider si un modèle de turbulence st précis ou non. Sur le plan
théorique, on distingue deux régions dans la couche limite :
– la première, très près de la paroi, est appelée sous couche laminaire et elle est le siège de
la loi linéaire
– l’autre, loin de la paroi, est appelée zone d’équilibre et elle est le siège de la loi logarithmique.
Les deux régions sont contrôlées par les valeurs respectives de la viscosité laminaire et de la viscosité
turbulente, que l’on retrouve dans le tenseur des contraintes :
2
τ = (µ + µt )(∇u + ∇ut − ∇.uI)
3
On suppose que l’on étudie une couche limite de plaque plane turbulente avec un très faible
gradient de pression longitudinal. Cela permet de simplifier fortement les équations de Navier-
Stokes dans la couche limite. En particulier, on montre que la composante normale du cisaillement
est conservée dans ce cas :
∂u
τ = (µ + µt ) = τw .
∂y
On définit la vitesse de frottement :
τw
r
uτ = .
ρw
On introduit une vitesse et une distance adimensionnelles appelées les « variables de paroi » :
u yuτ
U+ = et y + = .
uτ ν

21
Sous-couche linéaire

Au voisinage immédiat de la paroi, les contraintes turbulentes tendent vers 0 et la viscosité

laminaire domine largement la viscosité turbulente, d’où :
∂u
τ =µ
∂y
Par définition de la vitesse de frottement et en approchant la dérivée de la vitesse par son approxi-
mation linéaire, on peut écrire :
∂u u
τw = ρw u2τ = µ =µ
∂y y
Soit,
u ρw u τ
=y
uτ µ
On écrit une loi de la vitesse en fonction de la distance à la paroi. Si on réécrit cette expression
en fonction des variables de paroi, on obtient la loi linéaire :
u+ = y +
La loi linéaire est vérifiée pour y + < 3.

Zone logarithmique

Dans cette zone, les tensions turbulentes dominent la composante laminaire du tenseur des
contrainte car la viscosité turbulente domine la viscosité laminaire. La viscosité turbulente est
obtenue par produit des échelles caractéristiques de la turbulence. Des essais en souflerie montre
qu’on peut la modéliser comme suit :
µt = ρκuτ y .
La nouvelle constante est appelée constante de Von Karman : κ = 0.41.
On a par ailleurs,
∂u
µt = ρw u2τ ,
∂y
et d’autre part, la contrainte de Reynolds reste en première approximation constante dans la
direction normale à la paroi :
τ = τw = ρw u2τ

On obtient donc,
∂u ∂u
µt = ρu2τ = ρκuτ y
∂y ∂y
où y représente la distance à la paroi. Soit,
uτ ∂u
=
ky ∂y

On peut intégrer cette équation par rapport à y et l’écrire en fonction des variables de paroi,
soit :
1
u+ = ln(y + ) + C
k
où C est une constante d’intégration.
Les très nombreux essais en soufflerie subsonique et transonique montrent que c’est une
constante universelle proche de 5. La loi logarithmique est valable pour des valeurs de y + supé-
rieures à environ 30. La région comprise entre y + = 3 et y + = 30 est appelée région tampon.

22
Comparaison entre calculs et lois théoriques

Les comportements théoriques sont obtenus sur une plaque plane. L’objectif d’un TP est
de les retourver par simulation numérique. Pour cela, on simule avec elsA un écoulement sur
une plaque plane adiabatique à M∞ = 0.85. Pour retrouver les lois linéaire et logarithmique (cf
Figure 3.1), on construit les variables de paroi u+ et y + et on trace la fonction u+ = f (y + ). On se
place en un point assez éloigné du bord d’attaque de la plaque afin de disposer d’une couche limite
suffisamment établie. On récupère la valeur du frottement pariétal en ce point. On fait ensuite
une coupe en ce même point suivant la direction normale à la paroi et on relève les valeurs de la
viscosité et de la masse volumique à la paroi.
Frame 001  21 Sep 2006  Laminar computation -- Flat Plate

20
uplus, uplus_th, uplus_log

15

u+=f(y+)
u+=y+
u+=5.6 log(y+)+5
10

0 0 1 2 3 4 5
10 10 10 10 10 10
yplus

Fig. 3.1 – Plaque plane (modèle de Spalart) - Loi de paroi

Comparaison laminaire / turbulent

Les couches limites pour les régimes d’écoulement laminaire et turbulent sont différentes(cf
Figure 3.2). La couche limite turbulente est plus épaisse que la laminaire. En effet, la loi de
variation de la vitesse dépend de la viscosité du fluide qui induit un frottement entre les couches
voisines. Par rapport à un écoulement laminaire, il existe en plus de la viscosité turbulente qui
renforce le terme de diffusion visqueuse. Comme l’écoulement est plus visqueux, il y a besoin de
plus de place pour que la vitesse locale atteigne la valeur de la vitesse amont.

Comparaison des modèles de turbulence

Les écoulements sur plaque plane sont les plus simples à reproduire par le calcul. On doit
obtenir la même couche limite pour tous les modèles de turbulence (cf. Figure 3.2) dans la mesure
où les modèles ont été mis au point pour coller aux résultats expérimentaux sur la plaque plane.
Sur des configurations plus complexes, comme un profil d’aile d’avion par exemple, les résul-
tats peuvent fortement différer. Dans le cadre du TP « Ecoulement autour d’un profil RAE2822 »

23
 Frame 001  21 Sep 2006  Laminar computation -- Flat Plate

0.7

0.6
laminaire
Spalart
K-Omega
0.5
u

0.4

0.3

0 0.02 0.04
y

Fig. 3.2 – Plaque plane - couche limite laminaire et turbulente. Effet du choix du modèle de
turbulence.

on traite deux cas bien connus d’écoulements transoniques, les cas 9 et 10, pour lesquels on pos-
sède des données expérimentales. Ces cas servent notamment à valider ou comparer des modèles
numériques.
Le cas 9 est un cas considéré sans décollement ou avec un décollement très faible. Pour ce cas,
on obtient des résultats qui collent bien aux données expérimentales. On visualise le coefficient de
pression et le coefficient de frottement (cf Figures 3.3 et 3.4). On constate que les deux modèles
de turbulence donnent des résultats quasiment identiques.

24
 Frame 001  27 Sep 2006  RAE2822 - Mach 0.73

Spalart
1 K-Omega
exp

0.5
Cp, V2

-0.5

0.2 0.4 0.6 0.8 1

x, V1

Fig. 3.3 – RAE2822 (cas 9) - coefficient de pression.

Frame 001  27 Sep 2006  friction on wall_1

0.01
Spalart
Spalart
K-Omega
K-Omega
exp
Cf, frictionmodulus, V2

0.005

-0.005

0 0.2 0.4 0.6 0.8

x, V1

Fig. 3.4 – RAE2822 (cas 9) - coefficient de frottement.

Le cas 10 est un cas avec décollement en pied de choc (cf Figure 3.5). Là, la physique est
beaucoup plus complexe et les modèles ne sont pas initiallement conçus pour la reproduire par-
faitement. On touche là aux limites d’une modélisation RANS. Les calculs effectués dans le cadre
du TP montrent que la position du choc et la taille du décollement sont dépendants du choix du

25
modèle de turbulence. Aucun ne donne la même chose et aucun non plus n’est capable de prédire
correctement la position du choc. On remarque que la simulation prédit le choc trop tard par
rapport à la solution expérimentale (cf Figure 3.6) et que l’évolution du frottement (variation de
u) change fortement entre les 2 modèles (cf. Figure 3.7). On retrouve la présence d’un décollement
quand le coefficient de frottement change de signe le long du profil, car il est directement lié aux
variations de la vitesse par la relation Cf = µ ∂u
∂y .

Frame 001  26 Sep 2006  RAE2822 - Mach 0.73

0.0535

0.053
z

0.0525

0.052

0.659 0.66
x

Fig. 3.5 – RAE2822 (cas 10) - vecteurs vitesse suivant une ligne du maillage proche du pied du
choc. Visualisation du retournement des vecteurs vitesse (décollement)

26
Frame 001  21 Sep 2006  RAE2822 - Mach 0.73

K-Omega
1 Spalart
exp
Cp, V2

0.5

-0.5

0.2 0.4 0.6 0.8 1

x, V1

Fig. 3.6 – RAE2822 (cas 10) - Coefficient de pression.

Frame 001  21 Sep 2006 

0.008 Spalart
Spalart
K-Omega
0.006 K-Omega
exp
Cf, frictionmodulus, V2

0.004

0.002

-0.002

-0.004

0 0.2 0.4 0.6 0.8

x, V1

Fig. 3.7 – RAE2822 (cas 10) - Coefficient de frottement.

27
3.1.4 Phénomènes physiques
Déplacement du choc avec l’incidence

L’écoulement autour d’un profil est fortement dépendant des conditions de vol. Pour le vérifier
à moindre coût numérique, un TP propose d’effectuer 2 calculs Euler autour du profil NACA0012
pour 2 angles d’incidence différents. On constate alors que la position du pied de choc sur le profil
change fortement (décalage de 5% vers le bord de fuite). Ce résultat est montré sur les figures 3.8
et 3.9.
Frame 002  20 Sep 2006  NACA0012 - Mach 0.85

0.2 Frame 003  20 Sep 2006  NACA0012 - Mach 0.85

0.3

0.1
0.2
z

0.1

0
0

-0.1

0.8 0.9 1 1.1

-0.1 x

0.8 1
x

Fig. 3.8 – NACA - position du choc pour α = 2

Frame 003  20 Sep 2006  NACA0012 - Mach 0.85

0.2

Frame 002  20 Sep 2006  NACA0012 - Mach 0.85

0.1
0.2
z

0.1
z

0 0

-0.1

0.8 1
x

-0.1

0.8 0.9 1 1.1

x

Fig. 3.9 – NACA - position du choc pour α = 4

28
Phénomène de décrochage

Le phénomène de décrochage est lié à une baisse de la portance Cz et peut être obtenu
par le calcul. Pour preuve, une polaire est tracée en scéance de TP pour le profil NACA 0012
en Euler. Normalement, le décrochage n’existe pas avec les écoulements obtenus par résolution
des équations d’Euler car le décrochage est lié à la viscosité du fluide dans la couche limite (en
particulier, présence de décollement). Les équations d’Euler ne permettent pas de modéliser une
couche limite et c’est la viscosité numérique du solveur elsA qui permet ici de construire une telle
polaire.

Polaire : Cz = f (α)
1.4

1.2

0.8
Cz

0.6

0.4

0.2

0
0 2 4 6 8 10 12 14 16
angle d’attaque en degrés
Fig. 3.10 – Polaire obtenue par calcul Euler : effets de la viscosité numérique d’elsA.

3.2 Influence des schémas numériques

Un schéma numérique sert à discrétiser les équations de Navier-Stokes et se ramène au calcul
du flux à l’interface. On définit donc un « pas » d’espace qui représente la taille des mailles du
maillage.
Soit W le vecteur des variables conservatives :
 
ρ
 ρu 
 
W =  ρv  .

 ρw 
ρet
Les schémas sont très généralement construits par une analyse 1D puis étendus aux autres dimen-
sions. Ainsi, les explications seront donc faites en 1D.

29
 Fig. 3.11 – détails du maillage

On se propose de résoudre l’équation modèle d’advection suivante :

∂Y ∂Y
+α = 0,
∂t ∂x
où Y représente une variable conservative et α une constante appelée vitesse d’advection.
L’équation d’advection est résolue par une méthode volumes finis qui permet de prendre en
compte des discontinuités dans la solution. Ainsi, on intègre l’équation précédente sur une cellule
indicée i du maillage :
∂Y ∂Y
Z Z
+α = 0.
Ci ∂t Ci ∂x

On peut réécrire :
∂Y
Z Z
= Y.n.ds ,
Ci ∂x ∂Ci
avec
– ∂Ci bord de la cellule Ci ,
– n vecteur normal dirigé vers l’extérieur de Ci .
Ainsi, en 1D, :
∂Y
Z
α = fi+1/2 − fi−1/2
Ci ∂x

où fk représente le flux à l’interface k et f peut s’écrire sous la forme f = αY . Pour les deux
schémas expliqués ci-dessous, on ne s’interessera plus qu’au flux numérique permettant le calcul
des flux advectifs aux interfaces.

3.2.1 Schéma centré de Jameson

On cherche donc à déterminer la valeur du flux f aux interfaces i − 1/2 et i + 1/2.
Pour le schéma de Jameson, on utilise un développement en série de Taylor. Si on développe
f(x) en séries de Taylor au point d’abcisse z :

f (x) = f (z) + f 0 (z)(x − z) + o(x − z)2 .

Si on développe à l’interface z = i + 1/2 :

0
fi = fi+1/2 + fi+1/2 (xi − xi−1/2 ) + o(xi − xi−1/2 )2
0
fi+1 = fi+1/2 + fi+1/2 (xi+1 − xi−1/2 ) + o(xi+1 − xi−1/2 )2
∆x ∆x
0
=⇒ fi + fi+1 = 2fi+1/2 + fi+1/2 (− + ) + o(∆x)2
2 2
fi + fi+1
⇐⇒ = fi+1/2 + o(∆x)2 .
2

30
 Il s’agit donc d’un schéma centré d’ordre 2. Mais un schéma purement centré est instable,
c’est pourquoi on utilise le schéma centré de Jameson qui lui rajoute de la dissipation numérique
ou dissipation artificielle. L’équation d’advection s’écrit alors :

∂u ∂u ∂4u
+a = 4
∂t ∂x ∂x
et le flux à l’interface :
1
fi+1/2 = (fi + fi+1 ) − k 4 (ui+2 − 3ui+1 + 3ui−1 + ui−2 ) .
2

3.2.2 Schéma décentré de Roe

Comme son nom l’indique, le schéma décentré ne fait pas une moyenne des deux valeurs des
cellules voisines de l’interface, mais prend la valeur de la cellule qui précède selon le signe de α
(pour un schéma décentré d’ordre 1).

si α > 0

fi
fi+ 12 =
fi+1 si α < 0

a>0 a<0

i i+1

i−1/2 i+1/2 i+3/2

a>0

a<0

Fig. 3.12 – schéma - détail du maillage

Y se propage à la vitesse α, suivant l’équation d’advection :

∂u ∂u
+a =0
∂t ∂x
On pose :
α+ = max(α, 0)

31
 α− = min(α, 0)
|α| = α+ − α−
α = α+ + α−

On réécrit l’équation d’advection :

∂u ∂u
+ (α+ + α− ) =0
∂t ∂x
∂u ∂u ∂u
+ α+ + α− =0
∂t ∂x ∂x
   
∂u ui − ui−1 ui+1 − ui
+ α+ + α− =0
∂t ∆x ∆x
   
∂u α + |α| ui − ui−1 α − |α| ui+1 − ui
+ + =0
∂t 2 ∆x 2 ∆x
   
∂u ui+1 − ui−1 ui+1 − 2ui + ui−1
+α = |α|∆x
∂t 2∆x 2∆x2
"    #    
∂u ui+1 + ui ui + ui−1 ui+1 − ui ui − ui−1
+α − − |α| + |α| =0
∂t 2∆x 2∆x 2∆x 2∆x
"    # "    #
∂u ui+1 + ui ui+1 − ui ui + ui−1 ui − ui−1
+ α − |α| − α − |α| =0
∂t 2∆x 2∆x 2∆x 2∆x

Par identification avec l’équation d’advection en volumes finis

fi+1/2 − fi−1/2
 
∂u
+ =0
∂t ∆x

On obtient :    
ui+1 + ui ui+1 − ui
fi+1/2 =α − |α| ,
2 2
soit :
1 |α|
fi+1/2 = (fi+1 + fi ) − (ui+1 − ui ) .
2 2
Mettre en évidence l’influence des schémas numériques sur une simulation est assez difficile
car on constatera seulement que le calcul aboutit ou n’aboutit pas. C’est ce que l’on montre
aux étudiants dans un premier temps en leur faisant lancer un calcul avec un schéma purement
centré puis avec un schéma décentré de Roe. Ces calculs introduisent une activité Matlab qui
vient en démonstration et explication du comportement des schémas numériques. Trois schémas
numériques différents ont été codé sous Matlab. Il s’agit du schéma de Lax-Wendroff, que l’on
peut rapprocher d’un schéma centré, d’un schéma « upwind method » qui correspond à un schéma
décentré et d’un schéma « superbee limiter » que l’on peut comparer à un schéma décentré d’ordre
2. On les applique à un maillage, où l’on a initialisé une moitié en forme de porte à la valeur u0 .
On recalcule la valeur de l’écoulement à chaque pas de temps avec les trois schémas. On observe
l’évolution de la solution exacte.
On constate effectivement que le schéma purement centré est instable, que le schéma décentré
est trop dissipatif alors que le schéma décentré d’ordre deux reste assez proche de la solution exacte.

32
Fig. 3.13 – schémas numériques à t=0s

Fig. 3.14 – schémas numériques à t=2s

33
 Fig. 3.16 – schémas numériques à t=7s

3.3 Techniques d’accélération de la convergence

3.3.1 Méthode explicite / implicite
Les équations de Navier-Stokes couplées aux modèles de turbulence contiennent un terme
d’advection, un terme de diffusion et un terme source. Pour faciliter l’exposé, ces 3 termes seront
regroupés sous la dénomination de R (pour « Right-Hand-Side « ).
Le but des phases explicite et implicite est l’intégration en temps de l’équation suivante :

W n+1 − W n
V +R=0
∆t
avec V le volume d’une cellule et W vecteur des variables conservatives.

Phase explicite

On réécrit l’équation :
W n+1 − W n
V + Rn = 0
∆t
avec Rn calculé au temps n, c’est-à-dire avec les champs conservatifs pris au temps n. On parle
alors de schéma explicite car l’état n + 1 se calcule directement à partir de l’état n :
∆t n
W n+1 = W n − R
V

La phase explicite à l’avantage d’être facile à implémenter. Cependant, elle n’est pas très
robuste. Sa condition de stabilité impose un pas de temps très faible et par conséquent un grand
nombre d’itérations pour atteindre la convergence. Le temps et le coût de calcul peuvent ainsi être
très importants.

34
Phase implicite

On réécrit l’équation :
∆W
V + Rn+1 = 0
∆t
On développe Rn+1 en séries de Taylor :
 ∆R n
Rn+1 = Rn + ∆W
∆W
On réinjecte dans l’équation :
−
→  ∆R n
∆W
V + Rn + ∆W = 0
∆t ∆W
  ∆R n 
V
⇐⇒ I+ ∆W + Rn = 0
∆t ∆W
⇐⇒ K∆W + Rn = 0
⇐⇒ ∆W = K −1 Rn

La phase implicite revient à faire une inversion de matrice, ce qui est assez lourd ; cependant,
il n’y a pas de limitation sur le CFL, ce qui implique une réduction du temps de calcul importante.
Dans la pratique, on ne calcule pas exactement la matrice inverse, mais on l’approxime. Il
existe plusieurs méthodes. On peut soit l’approximer par un scalaire, soit essayer d’en calculer une
valeur aprochée par des méthodes numériques d’inversion de matrice (méthode LU par exemple).

3.3.2 Méthode multi-grille

La méthode multi-grille est un moyen d’accélérer la convergence.
Elle se base sur plusieurs niveaux de grille, de la grille la plus fine (maillage initial) à des
grilles de plus en plus grossières (on divise par deux le nombre de points dans chaque direction du
maillage). On cherche d’abord une solution un peu convergée sur la grille fine. On transpose cette
solution sur la grille grossière. Là, le temps de calcul jusqu’à la convergence est beaucoup plus
court car directement lié au nombre de points du maillage. On interpole cette solution "grossière"
sur le maillage fin et le calcul converge alors beaucoup plus rapidement.
Il existe deux cycles de multi-grille, le cycle en V et le cycle en W. Le plus souvent, les calculs
se font sur deux niveaux de grille.
– Le cycle en V décrit un V ! Il part d’une grille fine, passe sur 2 niveaux de grilles grossières
et remonte de la même manière à la grille fine
– Le cycle en W décrit, comme son nom l’indique, un W. Il part d’une grille fine, descend
sur 2 niveaux de grilles, remonte à la grille intermédiaire, puis redescend à la grille la plus
grossière pour enfin remonter à la grille fine
Le cycle en V est le plus souvent utilisé. Le cycle en W est pourtant plus puissant mais
pour atteindre la valeur des résidus souhaités, ils sont équivalents. Alors que le cycle en V stagne
presque passé cette valeur, le cycle en W continue de décroître.

35
Fig. 3.17 – 2 niveaux de grilles

36
 Frame 001  28 Sep 2006  NACA0012 - Mach 0.85 - Residus

0
10

10-2 monogrille
multi-grille
10-4
residual_ro

-6
10

-8
10

-10
10

10-12

10-14

1000 2000 3000 4000 5000

iteration

Fig. 3.18 – comparaison des résidus avec et sans méthode multi-grille

37
Conclusion

Le site propose donc aujourd’hui trois TP cohérents avec le modèle de base que l’on s’est fixé.
Grâce au système Wiki, il va être amené à evoluer, en s’enrichissant de l’apport de chacun : de
nouveaux sujets de TP, des points théoriques supplémentaires ou plus approfondis, un lexique plus
complet, etc. Le fait d’avoir défini un modèle à suivre impose quelques contraintes de rédaction
selon le sujet mais garantit cependant une unité des sujets de TP et une évolution structurée du
site.
On peut se rendre compte avec le recul que les TP ne présentent pas aux étudiants une vision
très objective de ce qu’est le numérique. En effet, les exercices proposés ont été « prémachés »
pour que les calculs soient rapides et aboutissent. En réalité, le travail d’ingénieur en simulation
numérique est une boucle comprennant :
– le maillage,
– le prétraitement du calcul,
– le calcul lui-même,
– le post-traitement.
et plusieurs passages de cet algorithme sont souvent nécéssaires pour obtenir une solution de bonne
qualité.

38
Annexe A

présentation du Cerfacs
Créé en 1987, le CERFACS est une organisation de recherche qui a pour but de développer
des méthodes avancées pour la simulation numérique et la solution algorithmique de problèmes
scientifiques et technologiques complexes intéressant aussi bien la recherche que l’industrie. Cette
mission requiert l’accès aux ordinateurs les plus puissants disponibles.
Le CERFACS a cinq actionnaires :
1. le CNES (Centre National d’Études spatiales)
2. EADS (European Aeronautic and Defence Space company)
3. EDF (Électricité de France)
4. Météo-France
5. SNECMA (Société Nationale d’Étude et de Construction de Moteurs d’Avions)
Le CERFACS abrite des équipes inter-disciplinaires, aussi bien pour la recherche que pour
la formation avancée dans les domaines de la physique, les mathématiques appliqués, l’analyse
numérique et l’informatique.
Environ 100 personnes travaillent au CERFACS, dont 90 chercheurs et ingénieurs venant de
10 pays différents. L’activité est répartie en six groupes de recherche :
1. algorithmique parallèle
2. aérodynamique
3. combustion
4. climat et environnement
5. fusion de données
6. électromagnétisme.

Moyens informatiques

En interne, le CERFACS possède :

– un super-calculateur Cray XD1 muni de 120 processeurs optérons ;
– un Compaq AlphaServer SC45 de 40 processeurs ;
– une SGI Origin 2000 avec 32 processeurs.
En externe, nous avons la possibilité d’accéder au super-calculateur Fujistsu de Météo-France,
aux moyens du CCRT (Centre de Calcul Recherche et Technologie) du Commissariat à l’Énergie
Atomique (CEA) et au CINES (Centre Informatique National de l’Enseignement Supérieur). En
mars 2005, l’équipe combustion a même eu la possibilité de porter son code de recherche sur l’IBM
BlueGene/L du Thomas Watson Research Center qui était à ce moment la 8ème machine la plus
puissante au monde selon le Top500. SUr un maillage de 40 millions de cellules, ils ont observé
un speed-up de 2037 sur une configuration de 2048 processeurs et un speed-up de 4078 sur une

39
configuration de 4096 processeurs. Démontrant ainsi l’excellente stalabilité de leur code puisque
proche d’un speed-up linéaire.

40
Annexe B

Nous donnons ici un exemple de script python pour elsA.

!/usr/bin/env EpMuse.py
# ==========================================================================
# Project: elsA - DSNA/ELSA - Copyright (c) 1998-2000 by ONERA
# Type : <17470693 7474> Test file
# File : script/naca01-EXTRACT/naca01-EXTRACT.epy
# Vers : $Revision: 1.13 $
# Chrono : No DD/MM/YYYY Author V Comments
# 1.1 15/05/2000 MCLP 2.0 User Python, tex text, nb it, join.
# 1.0 30/08/1999 AJ/SB 1.0 Creation
# ==========================================================================
# Name naca01-EXTRACT
# Memo transonic inviscid flow around a 2-D profile - 1 domain
# See
# void
# Physics
# 1- Configuration
# > 2-D profile
# 2- Computational domain
# > The computational domain is a C-shape domain around the profile.
# 3- Equations
# > Euler equations with perfect gas EOS
# 4- Wall conditions
# > Slip condition on the walls
# 5- External conditions
# > Non-reflective conditions
# 6- Initial conditions
# > Uniform flowfield
# Numerics
# 1- Computational domain
# > See physics for the definition of the domain
# > Mesh : 1 domain with 257x33x2 mesh points
# 2- Space discretization
# > Jameson centered scheme with artificial dissipation
# 3- Time discretization
# > Explicit Runge-Kutta (RK4) with local time stepping
# 4- Resolution method
# > Explicit
# 5- Wall conditions
# > boundary type "wallslip" with extrap 0
# 6- External conditions
# > boundary type "nref"
# 7- Initial conditions
# > At t=t0, the far-field state is defined in each cell.
# Computer
# 1- Type : SGI
# 2- Processors: 1
# 3- Memory: ?
# 4- CPU: ?
#
# Files
# 1- Data
# > not used
# 2- Mesh
# > ./naca_eu.mai
# 3- t0 fields
# > not used
# 4- t1 fields
# > not used
# Design
# void

41
# References
# void
#
# ===========================================================================
#
# --- HighLigths
# This par is used into the Validation document
# Each line is a LaTeX item, thus you can insert pure LaTeX text
#
#@ NACA0012 2D wing profile.
#@ Infinite flow : ${M_\infty}~=~0.85$, $\alpha~=~1\deg$.
#@ Steady transonic perfect gas flow.
#@ Mesh:~ $257 \times 33 \times 2$ .
#@ Initial field : far-field state.
#@ Jameson centered fluxes with artificial dissipation ‘‘dissca’’ : $[~0.5~,~0.032~,~1~]$.
#@ Time integration : Runge-Kutta 4 steps with freezing.
#@ Implicit residual smoothing.
#@ Local time step.
#@ 1500 iterations, CFL=15.
#@ Flux extraction.
#
import time
#
# ---------------------------------------------------------------------------
# PRELIMINARY COMPUTATION
#
# Reference used by normalization: RoInf, Vinf
# Calculation of the far-field state corresponding to the values
# of the Mach number and the angle of attack
# (flow is 2-D in coordinates (x,z))
#
from math import *
#
MachInf=0.85
Gamma=1.4
Alpha = 1.0 * pi /180.
Pinf = 1. / (Gamma*MachInf*MachInf)
#
RoInf = 1.
RouInf = cos(Alpha)
RovInf = 0.
RowInf = sin(Alpha)
RoeInf = Pinf/(Gamma-1.) + 0.5
#
IBF1=33
IBF2=225
#
IMIN=1
IMAX=257
#
JMIN=1
JMAX=33
#
KMIN=1
KMAX=2
#
#-----------------------------------------------------------
# PROBLEM CREATION AND NAMING
#----------------------------
#
naca = cfdpb(name=’naca’)
naca.set(’config’,’2d’)
#
#-----------------------------------------------------------
# MESH CREATION, NAMING, DESCRIPTION
#-----------------------------------
#
mshnaca = mesh(name=’mshnaca’)
#
mshnaca.set(’file’,’naca_eu.mai’)
mshnaca.set(’format’,’fmt_tp’)
#
mshnaca.submit()
#
#-----------------------------------------------------------
# BLOCK CREATION, NAMING, DESCRIPTION
#------------------------------------
#
blknaca = block(name=’blknaca’)
#

42
blknaca.set(’mesh’,’mshnaca’)
#
blknaca.show()
blknaca.submit()
#
#-----------------------------------------------------------
# MODEL CREATION, NAMING, DESCRIPTION
#------------------------------------
#
mod = model(name=’mod’)
#
mod.set(’fluid’,’pg’)
mod.set(’phymod’,’euler’)
mod.set(’gamma’,1.4)
#
mod.show()
mod.submit()
#
#-----------------------------------------------------------
# NUMERICS CREATION, NAMING, DESCRIPTION
#---------------------------------------
#
num = numerics(name=’num’)
num.set(’limiter’,’none’)
#num.set(’multigrid’,’none’)
num.set(’restrict_type’,’synchronous’)
#
num.set(’time_algo’,’steady’)
num.set(’ode’,’rk4’)
num.set(’implicit’,’irs’)
num.set(’flux’,’jameson’)
num.set(’artviscosity’,’dissca’)
num.set(’av_base’,[0.5, 0.032, 1.0])
num.set(’freezing’,1)
num.set(’multigrid’,’w_cycle’)
num.set(’nbcoarsegrid’,2)
num.set(’nitercoarse’,2)
num.set(’niter’,700)
num.set(’inititer’,1)
#
#
num.set(’avcoef_k2’,0.5)
num.set(’avcoef_k4’,0.032)
num.set(’avcoef_sigma’,1.)
#
num.set(’freqcomptimestep’,1)
num.set(’cfl’,11.)
#
#
num.set(’freqcompres’,1)
num.set(’convergence_level’,1e-10)
#
num.show()
#
#-----------------------------------------------------------
# BOUNDARY CREATION, NAMING, DESCRIPTION
#---------------------------------------
# No boundary on K-faces due to 2-D calculation
#
#--- jmin - join
#
w1 = window(’blknaca’,name=’w1’)
#
w1.set(’iw1’,IMIN)
w1.set(’iw2’,IBF1)
w1.set(’jw1’,JMIN)
w1.set(’jw2’,JMIN)
w1.set(’kw1’,KMIN)
w1.set(’kw2’,KMAX)
#
b1 = boundary(’blknaca’,’w1’,’join’,name=’b1’)
#
b1.set(’jtype’,’match’)
b1.set(’blkrac’,’blknaca’)
b1.set(’wndrac’,’w3’)
#
#--- jmin - wallslip
#
w2 = window(’blknaca’,name=’w2’)
#

43
w2.set(’iw1’,IBF1)
w2.set(’iw2’,IBF2)
w2.set(’jw1’,JMIN)
w2.set(’jw2’,JMIN)
w2.set(’kw1’,KMIN)
w2.set(’kw2’,KMAX)
#
b2 = boundary(’blknaca’,’w2’,’wallslip’,name=’b2’)
#
b2.set(’extrap’,0)
#
#--- jmin - join
#
w3 = window(’blknaca’,name=’w3’)
#
w3.set(’iw1’,IBF2)
w3.set(’iw2’,IMAX)
w3.set(’jw1’,JMIN)
w3.set(’jw2’,JMIN)
w3.set(’kw1’,KMIN)
w3.set(’kw2’,KMAX)
#
b3 = boundary(’blknaca’,’w3’,’join’,name=’b3’)
#
b3.set(’jtype’,’match’)
b3.set(’blkrac’,’blknaca’)
b3.set(’wndrac’,’w1’)
#
#--- jmax - nref
#
w4 = window(’blknaca’,name=’w4’)
#
w4.set(’iw1’,IMIN)
w4.set(’iw2’,IMAX)
w4.set(’jw1’,JMAX)
w4.set(’jw2’,JMAX)
w4.set(’kw1’,KMIN)
w4.set(’kw2’,KMAX)
#
b4 = boundary(’blknaca’,’w4’,’nref’,name=’b4’)
#
b4.set(’state’,’stateInf’)
#
#--- imin - nref
#
w5 = window(’blknaca’,name=’w5’)
#
w5.set(’iw1’,IMIN)
w5.set(’iw2’,IMIN)
w5.set(’jw1’,JMIN)
w5.set(’jw2’,JMAX)
w5.set(’kw1’,KMIN)
w5.set(’kw2’,KMAX)
#
b5 = boundary(’blknaca’,’w5’,’nref’,name=’b5’)
#
b5.set(’state’,’stateInf’)
#
#--- imax - nref
#
w6 = window(’blknaca’,name=’w6’)
#
w6.set(’iw1’,IMAX)
w6.set(’iw2’,IMAX)
w6.set(’jw1’,JMIN)
w6.set(’jw2’,JMAX)
w6.set(’kw1’,KMIN)
w6.set(’kw2’,KMAX)
#
b6 = boundary(’blknaca’,’w6’,’nref’,name=’b6’)
#
b6.set(’state’,’stateInf’)
#
#-----------------------------------------------------------
# INITIAL FIELD CREATION, NAMING, DESCRIPTION
#--------------------------------------------
#
stateInf = state(name=’stateInf’)
#
stateInf.set(’ro’, RoInf)

44
stateInf.set(’rou’,RouInf)
stateInf.set(’rov’,RovInf)
stateInf.set(’row’,RowInf)
stateInf.set(’roe’,RoeInf)
#
winit = window(’blknaca’,name=’winit’)
#
winit.set(’all’,’yes’)
#
init1 = init(’winit’,name=’init1’)
#
init1.set(’state’,’stateInf’)
#
#-----------------------------------------------------------
# EXTRACTION CREATION, NAMING, DESCRIPTION
#-----------------------------------------
#
#--- Cell extraction
#
wext1 = window(’blknaca’,name=’wext1’)
#
wext1.set(’all’,’yes’)
#
ext1 = extract(’wext1’,name=’ext1’)
#
ext1.set(’file’,’naca01_1deg_W_cell.tp’)
ext1.set(’var’,’conservative ’)
ext1.set(’title’,’NACA0012 - Mach 0.85 ’)
ext1.set(’loc’,’cell’)
#
#--- Node extraction
#
wext2 = window(’blknaca’,name=’wext2’)
#
wext2.set(’iw1’,1)
wext2.set(’iw2’,IMAX)
wext2.set(’jw1’,1)
wext2.set(’jw2’,JMAX)
wext2.set(’kw1’,1)
wext2.set(’kw2’,1)
#
ext2 = extract(’wext2’,name=’ext2’)
#
ext2.set(’file’,’naca01-1deg_W_node.tp’)
ext2.set(’var’,’xyz ro rou rov row roe mach’)
ext2.set(’title’,’NACA0012 - Mach 0.85 ’)
ext2.set(’loc’,’node’)
#
#--- Cell + boundary extension extraction
#
ext3 = extract(’wext1’,name=’ext3’)
#
ext3.set(’file’,’naca01-1deg_W_cellfict.tp’)
ext3.set(’var’,’xyz ro mach’)
ext3.set(’title’,’NACA0012 - Mach 0.85 ’)
ext3.set(’loc’,’cellfict’)
#
#--- Residus
#
extRes = extract_group(name=’extRes’)
#
extRes.set(’title’,’NACA0012 - Mach 0.85 - Residus’)
extRes.set(’windows’,’winit’)
extRes.set(’var’,’residual_ro residual_rou residual_row residual_roe’)
extRes.set(’file’,’./residual_1deg_W.tp’)
extRes.set(’norm’,1)
#
#--- Lift (standard output)
#
extLift1 = extract_group(name=’extLift1’)
#
extLift1.set(’windows’,’+w2’)
extLift1.set(’var’,’convflux_row’)
extLift1.set(’fluxcoeff’,-2.)
#
#--- Lift (file)
#
extLift2 = extract_group(name=’extLift2’)
#
extLift2.set(’title’,’NACA0012 - Mach 0.85 - Lift (fluxcoeff -2)’)

45
extLift2.set(’windows’,’+w2’)
extLift2.set(’var’,’convflux_row’)
extLift2.set(’file’,’./lift_1deg_W.tp’)
extLift2.set(’fluxcoeff’,-2.)
#
#--- Test : mode = 2
#
extMode = extract_group(name=’extMode’)
#
extMode.set(’title’,’Test : mode=2’)
extMode.set(’windows’,’-w2’)
extMode.set(’var’,’convflux_row’)
extMode.set(’file’,’./allfluxes_1deg_W.tp’)
extMode.set(’fluxcoeff’,-2.)
#
#--- Test : multiple observors period/mode
#
extObs = extract_group(name=’extObs’)
#
extObs.set(’windows’,’-w2’)
extObs.set(’var’,’convflux_row’)
extObs.set(’period’,2)
#
#--- Test : lift on two windows
#
w7 = window(’blknaca’,name=’w7’)
#
w7.set(’iw1’,IBF1)
w7.set(’iw2’,100)
w7.set(’jw1’,JMIN)
w7.set(’jw2’,JMIN)
w7.set(’kw1’,KMIN)
w7.set(’kw2’,KMAX)
#
w8 = window(’blknaca’,name=’w8’)
#
w8.set(’iw1’,100)
w8.set(’iw2’,IBF2)
w8.set(’jw1’,JMIN)
w8.set(’jw2’,JMIN)
w8.set(’kw1’,KMIN)
w8.set(’kw2’,KMAX)
#
extWnds = extract_group(name=’extWnds’)
#
extWnds.set(’title’,’Test : lift on two winodws’)
extWnds.set(’windows’,’-w7-w8’)
extWnds.set(’var’,’convflux_row’)
extWnds.set(’file’,’./liftWindows_1deg_W.tp’)
#
#-----------------------------------------------------------
# CONNECTING NUMERICS AND MODEL WITH PROBLEM
#-------------------------------------------
#
naca.set(’global_numerics’,’num’)
naca.set(’global_model’,’mod’)
#
naca.submit()
naca.show()
#
#-----------------------------------------------------------
# FLOW COMPUTATION
#-----------------
#
t1=time.clock()
#
naca.compute()
#
t2=time.clock()
#
#-----------------------------------------------------------
# RESULT EXTRACTION
#------------------
#
naca.extract()
#
#-----------------------------------------------------------
# PRINT CPU Time
#---------------
#

46
print ""
print "============================="
print "CPU Time"
print "-----------------------------"
print "total for compute : ", t2-t1
print "============================="
print ""
#
#-----------------------------------------------------------
# PRINT Memory
#-------------
#
naca.statistics()
#

47
Bibliographie

[1] « Manuel théorique elsA »

[2] « beginner guide 1 » : documentation elsA
[3] « elsA User’s reference manual »
[4] Yann Colin « Overview of preconditionning techniques », septembre 2005
[5] Yann Colin et Jean-Philippe Boin, « Driving state definition and normalisation within
elsA-Python script », fevrier 2006
[6] Marcel Lesieur « La turbulence », Collection Grenoble Sciences,
[7] Yann Colin « TP3 » et « TP4 », novembre 2005
[8] « Reynolds stress transport modelling of transonic flow around the RAE2822 airfoil »,
AIAA-paper, 94-0309, 32nd Aerospace Science Meeting, janvier 1994
[9] J.L. Estivalezes « Notes de cours sur la turbulence », octobre 2005

48