|
|
ion d'ENS2D : la méthode des
diffé
U n i v e r s i t é H a s s a n 1 e
r
Faculté des Sciences et Techniques
Settat
|
Mémoire de Fin d'Eudes présentée par
:
R. BENRAZOUK
pour l'obtention du
master
« Ingénierie et modélisation des
systèmes mécaniques »
Sous le thème :
Résolution Numérique des équations
de Navier Stokes
bidimensionnelle par Méthode des Différences
finies et
Méthode Spectrale
Soutenu le : 24/11 /2009
Les membres du jury :
- Mr J .DABOUNOU : Professeur à la FST de Settat
Examinateur
- Mr M.LOUZAR : Professeur à la FST de Settat
Examinateur
-Mr E. SEMAA : Professeur à la FST de Settat
Examinateur
-Mme G. MANGOUB : Professeur à la FST de Settat
Encadrante
Je dedie ce modeste travail a tous ceux qui ont contribue a son
aboutissement de pres ou de loin.
A mes parents pour leur disponibilite, conseil, aide materiel et
morale.
A mon frere, mes soeurs pour leurs gentillesse.
A mes cheres amies avec qui j 'ai passe des meilleurs moments,
sans oublie bien sur mes formateurs, mes formatrices de la faculte des sciences
et Techniques de Settat qui mon donnee beaucoup de renseignements durant toute
ma formation.
Remerciements
Je tiens à exprimer mes remerciements les plus
vifs à mon encadrante Mme MANGOUB Ghita professeur à la FST de
Settat, qui n'a épargné aucun effort pour m'assister et
m'orienter durant ma période de projet, pour avoir accepté
d'encadrer ce projet de fin d'études, pour la confiance qu'il m'a
témoigné mais aussi pour ses encouragements, sa
disponibilité
Je remercie également Mr SEMMA El Alami
responsable de notre Master et directeur du laboratoire de mécanique
à la FST de Settat.
Mes remerciements vont à tous les membres de ma
famille et mes parents sans qui je ne serais pas là. Merci de m'avoir
toujours soutenu tout au long de ces années d'études
Enfin, que tout ceux qui ont contribué de
prés ou de loin à la réalisation de mon projet, trouvent
ici l'expression de mes sentiments les plus distingués.
Je remercie, enfin, les membres du jury d'avoir
accepté d'examiner ce modeste travail.
Résumé :
Dans ce travail nous présentons une
méthode numérique basée sur une approximation spectrale
collocation-Tchebychev pour résoudre les équations de
Navier-stokes qui régissent l'écoulement d'un fluide visqueux,
incompressible et bidimensionnel dans une cavité carrée
entrainée.
Les équations de Navier-stokes sont
formulées en termes de vitesse pression ou encore en termes de fonction
courant-tourbillon dans un espace à deux dimensions. La
discrétisation temporelle des équations de Navier-Stokes se fait
par un schéma d'intégration de deuxième ordre. Ce dernier
est une combinaison de deux schémas : le premier de Grank-Nicolson
appliqué sur le terme de diffusion et l'autre d'AdamsBaschforth de
second ordre qui est appliqué sur le terme d'advection. Les
résultats numériques sont présentés,
analysés et confrontés à d'autres résultats
numériques trouvés par l'autre méthode numérique
(comparaison avec la méthode des différences finies).
Abstract:
In this work we present a numeric method based on a
spectral approximation Collocation-Chebyshev of internal recirculating flows
encompassing a two dimensional viscous incompressible flow generated inside a
regularized square driven cavity.
The equations of Navier-stokes are formulated in terms
of speed pressure or in terms of function current-whirlwind in a space to two
measurements. The temporal discrétisation of the equations of
Navier-Stokes makes itself by a diagram of integration of second order. This
last is a combination of two diagrams: the first of applied Grank-Nicolson on
the term of diffusion and the other of Adams-Baschforth of second order that is
applied on the term of advection. The numeric results are presented, analyzed
and confronted to other numeric results found by the other numeric method
(comparison with the finished method difference).
SOMMAIRE
NOTATIONS 5
LISTE DES FIGURES 6
INTRODUCTION 7
1. Equations générales 9
1.1
Définition du problème 9
1.2 Equations de base
10
1.2.1 Les équations de Navier-stokes
10
1.2.2 L'équation de continuité
11
1.3 Forme adimensionnelle des équations de
Navier-Stokes . 12
2. Différentes formulations des équations
de Navier-stokes ... 14
2.1 Formulation en variables primitives (u, v, P)
14
2.2 Formulation (fonction de courant-fonction
tourbillon) 15
2.3 Formulation (vitesse-tourbillon) 17
2.4 Résolution numérique des
équations de Navier-stokes en formulation (iji,
w) 18
2.5 Les conditions aux limites ... 18
3. Résolution des équations de
Navier-stokes par la méthode des différences
finies .. 21
3.1 Introduction 21
3.2 Principe générale de la méthode
d'ordre O(h4) 22
3.3 Résolution de l'équation de poisson de
la fonction de courant . 23
3.4 Calcul des composantes de la vitesse 25
3.5 Résolution de l'équation de transport
du Tourbillon 25
4. Résolution des équations de
Navier-stokes par la méthode spectrale........... 28 4.1 Introduction .
28
4.2 Méthodes spectrales pour l'approximation
d'une fonction 29
4.3 Généralité sur la
méthode collocation-Tchebychev 29
4.3.1 Les propriétés principales des
polynômes de Tchebychev 29
4.3.2 Quadrature de Gauss-Lobatto 30
4.3.3 Représentation d'une fonction sur les
polynômes de Tchebychev 30
4.3.4 Matrice de passage pour les points de collocation
de Gauss-Lobatto 31
4.3.5 Dérivation numérique 31
4.4 Exemple 32
4.5 Projection de méthode spectrale
34
4.5.1 L'intégration temporelle 34
4.5.2 Les conditions aux limites 35
4.5.3 Procédures de résolutions des
équations 36
5. Les résultats et discussion 38
5.1 Formulation (ijj,
w): méthode des différences finies d'ordre
(O(H2)-O(H4)) 38
5.2 Formulation (vitesse -pression): méthode
spectrale . 40
5.3 Formulation fonction de courant fonction tourbillon :
méthode spectrale 43
Conclusion générale 44
Références bibliographiques .
45
NOTATIONS
Vo La vitesse de référence
Lo La longueur de référence
Champ de la pression non dimensionnel
~
p
p
=
PVo2
P Champ de la pression dimensionnel
PVoLo
Re =
Nombre de Reynold
11
At Pas de temps
AT Pas de temps fictif
h Pas de l'espace
Ax Pas de l'espace suivant la direction x
Ay Pas de l'espace suivant la direction y
Coordonnée horizontale non
dimensionnelle
~
x
=
x
Lo
Coordonnée verticale non
dimensionnelle
~
y
=
y
Lo
|
~
t
|
Vo
=
Lo
|
t Temps non dimensionnel
|
Composante horizontale non dimensionnel de la
vitesse
~
~
=
~
Vo
Composante verticale non dimensionnel de la
vitesse
~
V
=
V
Vo
xi les points de collocation
fN(x, t) L'approximation polynomiale de la
fonction
f~ic(t) Coefficients spectraux
N Nombre de polynômes ou points de collocation
Tk(x) = cos kx Polynômes de Tchebychev d'ordre k
N la viscosité cinématique
% Fonction tourbillon
p La masse volumique de fluide
(Pic(x) Les fonctions de base
' Fonction courant
Sii Symbole Kronecker
~)* Eléments de tenseur de contrainte
M La viscosité dynamique
A Coefficient de poisson
A L'operateur Laplace
ü L'operateur gradient
LISTE DES FIGURES
Figure (1.1) : cavité à paroi
supérieure entrainée avec une vitesse horizontale u=1 Figure
(4.1) : Procédures de résolutions des équations
Navier-Stokes : méthode spectrale
Figure (5.1) : les lignes iso valeurs de la fonction
de courant (Le calcul est fait avec la méthode des différences
finies d'ordre (O(H2)-O(H4)), formulation fonction de courant
tourbillon)
Figure (5.2) : les lignes iso valeurs de la fonction de
courant (méthode des différences finies d'ordre
(O(H2)-O(H4))
Figure (5.3) : Composante U de la vitesse en x=0.5 en
fonction de Y
Figure (5.4) : Composante U de la vitesse en x=0.5 en
fonction de Y, résultats obtenu Katuhiko Goda et R.
Burggraf.
Figure (5.5) : les lignes iso valeurs de la fonction de
courant (Calcul global spectral, vitesse pression).
Figure (5.6) : les lignes iso valeurs de la fonction de
courant (Calcul de fonction de courant pour chaque itération spectral,
vitesse pression).
Figure (5.7) : la composante horizontale de la vitesse en
x=0.5 en fonction de y en formulation vitesse pression avec la méthode
spectrale de collocation Tchebychev
Figure (5.8) : les lignes iso valeurs de la fonction de
courant (Le calcul est fait avec la méthode de collocation de
Tchebychev, formulation fonction de courant tourbillon).
INTRODUCTION
Dans ce travail, on s'intéressera à la
résolution numérique des équations de Navier-Stokes qui
régissent l'écoulement d'un fluide visqueux, incompressible, et
bidimensionnel.
Les Simulations numériques des équations
de Navier-Stokes pour les écoulements de fluide visqueux, incompressible
et bidimensionnel sont généralement basées sur une
formulation en termes de variables primitives (vitesse et pression) ou sur une
formulation en termes de variables fonction courant- fonction tourbillon. La
difficulté majeure qui survient avec la première des deux
formulations précédentes provient du fait d'associer la pression
avec la vitesse pour satisfaire la condition d'incompressibilité.
L'équation de la continuité contient seulement les composantes de
la vitesse, et il n'y a aucun lien direct avec la pression. Avec la formulation
de fonction courant-tourbillon on évite ce problème. Pour le cas
où on a trois dimensions il est préférable d'utiliser la
formulation en variables primitives.
Plusieurs méthodes ont été
proposées pour vaincre la difficulté qui survient dans la
formulation des variables primitives. Parmi celles-ci, la méthode des
différences finies est sans doute la plus utilisée en raison de
son aspect universel, de sa simplicité relative et de la facilité
de sa mise en oeuvre. Les Inconvénients de cette méthode sont :
limitation à des géométries simples, difficultés de
prise en compte des conditions aux limites de type Neumann...etc.
Dans le cadre de la méthode des
différences finies, le degré d'approximation locale est fixe
(typiquement d'ordre 2 ou 4) et la convergence avec l'augmentation du nombre de
points de collocation, vers la solution exacte est due au fait qu'une
approximation de degré donné est d'autant plus précise que
l'intervalle (la distance entre points de collocation voisins) sur lequel elle
est construite est petit.
C'est sur ce point que diffère la
méthode pseudo-spectrale : au lieu de figer le degré
d'approximation local en chaque point de collocation et d'augmenter le nombre
de ceux-ci, on construit une approximation globale (c.-`a-d. basée sur
l'ensemble de points de collocation). Cette approximation (polynomiale), de
degré d'autant plus élevé que le nombre de point sur
lequel elle est construite est grand, sera ainsi d'autant plus précise
que le nombre de points de collocation employés est grand.
L'objectif principal de ce projet est de
développer une méthode numérique effective pour la
résolution des équations de Navier-stokes qui régissent
l'écoulement d'un fluide visqueux,
incompressible et bidimensionnel dans une
cavité carrée. Pour ce but, la solution numérique de ces
équations est basée sur la méthode spectrale avec un
polynôme Tchebychev (nommée aussi méthode de
pseudo-spectral, collocation-Tchebychev). La motivation pour utiliser ce type
des méthodes est lorsque les méthodes spectrales ont de haute
précision et des erreurs très basses pour la prédiction de
l'écoulement dans un régime instationnaire. L'intégration
temporelle du système des équations est exécutée
par un schéma de second ordre semiimplicite (Adams-Bashforth et
Crank-Nicolson).
Les méthodes spectrales ont été
utilisées dans la combinaison temporelle avec un schéma de haut
ordre. Par exemple, [2] Johnny (2007) a utilisé un schéma
temporel d'ordre trois, pour améliorer la exactitude de son algorithme,
Dans ce projet on présente un algorithme basé sur une
méthode spectrale (méthode collocation-Tchebychev). Cet
algorithme numérique utilise une technique de la diagonalisation complet
(technique du non-itérative) lequel est très efficace et rapide
pour la solution directe des équations résultantes après
la discrétisation spatial et temporel.
Les résultats numériques sont
comparés et évalués avec des résultats
numériques précédemment publiés par d'autres
auteurs par exemple ( [4]-Maciej MATYKA -2004- et [2].Johnny de Jesús
Martinez and Paulo de Tarso T -2007- ).
Le rapport est organisé comme suit. En premier,
les formulations mathématiques sont présentées, suivi par
les différentes formulations des équations Navier-stokes. La
section suivante est consacrée à l'étude numérique,
qui consiste en la discrétisation temporelle et spatiale des
équations résultantes. On commence par la méthode des
différences finies puis après celle spectrale (méthode
avec collocation-Tchebychev). Dans la section suivante les résultats
numériques sont présentés. Nous terminerons notre travail
par une conclusion générale.
Chapitre1 : Equations générales
1.1 Définition du problème:
On considère un écoulement
instationnaire et bidimensionnel d'un fluide visqueux et incompressible dans
une cavité carrée. L'écoulement est engendré sous
l'influence de la viscosité par le mouvement de la paroi
supérieure qui est animée d'une vitesse de translation uniforme,
les autres parois étant immobiles. On a pris comme vitesse de
référence la vitesse de translation de la paroi en mouvement et
comme longueur de référence celle d'un coté.
Fig. (1.1) : cavité carrée
entrainée
Un fluide emplit un carré. En chaque point du
carré, la vitesse du fluide est V = (u, v). Le fluide est mis en
mouvement par le déplacement vers la droite de la paroi
supérieure. Les trois autres parois (gauche, droite et
inférieure) sont fixes. Le fluide adhère aux parois. Un
tourbillon central prend naissance puis se stabilise ou oscille en fonction de
la viscosité.
Pour étudier cet écoulement, on se base sur
les équations complètes de Navier-stokes.
1.2 Les équations de base:
1.2.1 Les équations de Navier-stokes
Les équations de Navier-Stokes s'obtiennent en
appliquant le principe fondamental de la dynamique, le principe de conservation
de la masse au mouvement du fluide. Elles se présentent sous la forme
d'un système d'équations aux dérivées partielles
non-linéaire.[1]
001
dV
P .~ ~ 2
3 456.
000000001
P 7 V-T (1.1)
Dans l'hypothèse d'un fluide newtonien,
hypothèse correspondant au comportement de la plupart des liquides et
gaz usuels nous avons les éléments du tenseur des contraintes qui
sont donnés par:
iii = XEkk6ii 7 4E4 (1.2)
Avec:
iii: Éléments du tenseur de
contrainte
Eii: Éléments du tenseur de
déformation
t et A: Coefficients de Poisson caractérisant les
propriétés visqueuses des fluides. Ils représentent
respectivement la viscosité de cisaillement et la viscosité de
dilatation.
Les éléments du tenseur de
déformation sont donnés par :
|
:)* ~
|
(1.3)
; >?~)
= 7 ?~* @
?~* ?~)
|
En reportant ces expressions dans l'équation
(1.1), on obtient:
P .~ ~ 2
3 456.
~ 7 ~~/ 7 9 7 ~456.000000001
divV (1.4)
it et A Vérifient l'équation suivante
:
3A 7 21.1. = 0 (1.5)
En l'absence de forces extérieures, lorsque le
fluide est incompressible l'équation (1.4) devient :
dV
(1.6)
P dt ~ 3456.
000000001 7
~/01
En explicitant la dérivée particulaire
à gauche, c'est-à-dire :
dV
=
dt
aVi_,
at
+ (V. v). V (1.7)
On obtient alors l'équation de Navier-stokes pour
un fluide incompressible:
p(OV
a
t
Avec :
p: La masse volumique
u : la viscosité dynamique du fluide
V_1(u, v): Le champ de vitesse du fluide P: la pression
local
1.2.2 L'équation de continuité:
La masse m d'un domaine fluide quelconque D, que l'on
suit au cours du temps reste constante :
dm
dt
= 0 (1.9)
En explicitant la dérivée particulaire de
m, on obtient l'équation globale de conservation de la masse
:
|
dm
dt
|
= Ill p dV = Ill(
|
|
at + div(pV)) dV = 0 (1.10)
|
Cette équation devant être
vérifiée pour tout domaine D, le théorème de
l'intégrale nulle permet d'obtenir l'équation locale de
conservation de la masse dite équation de continuité
:
Op
+ div(pV_1) = 0
(1.11)
at
Si le fluide est incompressible : p = pa
div(V_1) = 0 (1.12)
p à
= 0
Si l'écoulement est stationnaire :
at
div(pV) = pdiv(V) + grad(pV) = 0 (1.13)
+ V(0 V)) = --gradP + u0V(1.8)
La conservation de la masse d'un fluide en mouvement, de
masse volumique p et de champ de vitesse V(M ,t) est traduite par
l'équation locale.
0 p (1.14)
at
7 div(pV) = 0
Les deux équations (1.8) et (1.12)
définissent le système à résoudre (1.15 et 1.16)
suivant :
7 u 7 v =3 V2u
{au au au 1 OP t
at Ox Oy p Ox p
Ov Ov Ov 1 OP [t 2
ay =
at
7 u
Ox
7 v
P
(1.15)
ay7 V v P
Nous avons donc trois équations et trois inconnues
u ,v ,p considérées comme fonctions de x, y ,t . Cette forme des
équations est dite forme dimensionnelle.
La détermination physique du problème
c'est-à-dire l'intégration des équations pour la
détermination des inconnues u, v et p nécessite en outre la
connaissance des conditions initiales et des conditions aux
limites.
1.3 Forme adimensionnelle des équations de
Navier-Stokes
L' emploi de variables adimensionnelles dans les
équations permet d'approcher de plus la réalité des
phénomènes physiques. Cette nouvelle forme de mise en
équations est basée sur deux données une vitesse de
référence vo (nous rapporterons toute vitesse ou toute composante
de vitesse à cette valeur vo ) Une longueur de référence
Lo . Les temps seront
rapportés à la quantité to = Lo
qui est homogène à un temps et les pressions motrices à
la
vo
quantité Po = pv02qui est
homogène à une pression. Le changement de variables est
par
suite défini par les relations suivantes :
[1]
1 _ X y
t
vo =
_ v
VY =
-
vo
Lo
x= ,
Lo
U
_
U = ,
-
t
Lo
Vo
Lo Pvo
_ p
t , p = 2
vo
(1.17)
Ces nouvelles variables xp yp up vp p , sont sans
dimension. On les appelle aussi des variables réduites.
Le changement de variables nous conduit donc aux
relations suivantes :
Ou
?~ ~
Oft
Of
d(v011) Vo2
~
Lo
d
(L°
vo ~~[
|
OP
~
Ox
|
(PPv02) a(Rl0)
|
Pvo2
~
Lo
|
OP
R
|
Ou 0(Vail) Vo Oft
Ox ~ (1,0) ~Lo
OR
|
02u Ox2 =
|
02(you)
|
Vo
|
02i1
|
|
|
|
|
0(LoR)2
|
L2
|
OR2
|
|
02u 0y2 ~
|
02(Vo~~~
|
Vo
|
02i1
as,r,2
etc ...
|
|
a 0SO2
|
L2
|
Considérons la première des
équations de Navier-stokes (éq.1.15) elle s'écrit
maintenant :
2
Vo
Oil Vo 2
L0 at Lo
Soit :
a(u2) Vo2 031
7Lo
OP- uv0
(1.18)
7 -
L0 OR p La
1
~
0(UV)
~
ay
Pvo2
|
au
at 7
|
a(u2)
R 7
|
0071i-0
|
OP 1.1. 1
~ 3Ox 7 ~~~
~ ~~~~
|
|
ay
|
Le système d'équations (1.15) et (1.16)
devient donc :
V?
?~~ 7 ?~~~~~
?~~ 7 ?~~~~~~
?~~ ~ 3 ?~\
?~~ 7 = ~~~
T ~~
(1.19)
U ?~~ ?~~~~~~ ?~~~~~ 3 ?~\
at 7
?~~ 7
T 1
ay si 7
~~~
Re
-- a
Oft OV
aR 7 ay = 0 (1.20)
En posant
Re
PvoLo
=
VoLo
~
^
(1.21)
11
Sous les formes (1.19) et (1.20) les équations
indéfinies du mouvement sont appelées équation
réduites, elles sont sans dimensions et ne font intervenir maintenant
qu'un seul coefficient dépendant des données de
l'écoulement. Ce coefficient Re est le nombre de Reynolds de
l'écoulement, ce nombre est sans dimension, sa valeur numérique
ne dépend donc pas du système d'unités
choisies.
Chapitre2:Différentes formulations des
équations de
Navier-stokes
L'équation adimensionnelle qui définit
l'écoulement d'un fluide visqueux incompressible instationnaire
s'écrit sous la forme suivante :
|
Oui
Ot +
|
O(uiu;)
Ox*
|
=
|
OP
, +
axi
|
1
Re
|
O
Oxi
|
O
Oxi
|
ui (2.1)
|
Oui
Oxi
= 0 (2.2)
Avec i, j =1,2
On peu distinguer trois formulations différentes
de ces équation suivant les inconnues qu'on
considère.
2.1 Formulation en variables premières (u, v, P)
Equation dynamique
|
Equation de continuité
|
Ou
Ot + u
Ov
+ u
Ot
|
Ou
+ v
Ox
Ov
+ v
Ox
|
Ou
=
Oy
Ov
=
Oy
|
OP Ox + OP Oy +
|
RV2u (2.3)
e
RV2v (2.4)
e
|
Ou Ox +
Ov
O = 0 (2.5)
y
Où sous la forme conservative:
Ou
Ot +
OP
+
Ox
=
O(uv)
Oy
O(u2)
+
Ox
RV2u (2.6)
e
Ov
Ot +
OP
+
Oy
=
O(u2)
Oy
O(uv)
+
Ox
1
V2v (2.7)
Re
La formulation en variables primitives et la plus
naturelle et elle se prête directement à l'extension au cas
trimensionnel. Dans la plupart des méthodes, on substitue à
l'équation de continuité, une équation de poisson. Cette
équation s'écrit sous la forme convective suivante :
V2P = 2
(Ou Ov ) Ov Ou
(2.8)
Ox Oy + Ox Oy
U = 3
v=
Ox
OtIJ
Oy
OtIJ
(2.12)
(2.13)
Les conditions aux limites
Pour le fluide visqueux en contact avec une paroi, on
impose la condition d'adhérence à tout instant t,V( , y, t) pour
x et y appartenant à la paroi égal à Vp vitesse du point
lié à la paroi .Pour un fluide parfait, on impose la condition de
glissement, qui traduit uniquement que le fluide ne traverse pas la paroi. Il
est naturel de se donner les conditions initiales pour V, comme pour une
équation différentielle du premier ordre en t. Par contre, les
conditions que doit vérifier p sont moins claires .On remarque que la
pression est solution d'une équation de la forme
Grad(p)=A(V). La richesse des conditions aux limites du
domaine de résolution des équations à tout instant, est
plus grande que celle des conditions initiales. Pour le liquide visqueux en
contact avec une paroi, on impose la condition d'adhérence pour la
vitesse à tout instant. Cette condition est type Dirichlet et est
donnée par :
u = v = 0 (2.9)
Pour la pression, on impose une condition de type Newman,
soit :
OP
=
On
1 O2vi
(2.10)
Re Ox2
Dans le cas bidimensionnel, bien qu'il existe
différentes méthode pour résoudre numériquement le
système des équations de Navier Stokes en termes de variables
primitives, c'est-à-dire u, v, p, on peut le réduire à un
système plus simple de deux équations différentielles pour
la fonction de courant et la fonction tourbillon.
2.2 Formulation (fonction de courant-fonction
tourbillon)
Pour un fluide incompressible p = cte g div
V01 =0
|
01
gV
|
= roti; (2.11)
|
|
-,
Le vecteur vitesse V
|
01
dérive d'un potentiel vecteur V)
|
. Pour un écoulement plan
xOy, on a:
|
Les lignes de courant sont définies
par
dx dy
=
U v
|
?'
gOx dx 7
|
NJ
Oy
dy = 0 = 1(x, y, t) = cte
|
* Est appelé fonction de courant
La fonction tourbillon est définie en fonction de
V par la relation to = rot(V), elle se réduit à la seule
composante scalaire t( , y, t).
?~ ?~
% ~ 5"/ = 3(2.13)
?~ ?~
Elle s'exprime au moyen de ij(x, y, t) par la formule
:
|
02*
co = 7
Ox2
|
02*
|
= V2* (2.14)
|
|
0y2
|
On dérive les deux équations (2.6) et (2.7)
respectivement par rapport à x et y
|
a au at ay7
|
Ou Ou
Oy Ox7 u
|
02u
Oy Ox7
|
Ov Ou
Oy Oy7 v
|
02u 0y2 =
|
= N ? ?~ >?~~
?~~ 7 ?~~
?~~@O~~ (2.15)
|
|
a av at ax 7
|
au av
Ox Ox7 v
|
02v
Ox Oy 7 u
|
02v Ox2 7
|
av
Ox
|
av Oy =
|
= N ? ?~ >?~~
?~~ 7 ?~~
?~~@O~~ (2.16)
|
En soustrayant l'une de l'autre on obtient :
?~ _?~
? ?~ 3 ?~
?~` 7 ~ >?~~
?~~ 3 ?~~
?~ ?~@ 7 ~ > ?~~
?~ ?~ 3 ?~~
?~~@ 7 ?~
?~ _?~
?~ 7 ?~ 3 ?~
?~ _?~
?~ 7 ?~
?~` ?~`
= N ? ?~ >?~~
(2.17)
?~~ 7 ?~~
~ ?~ ~@O
~~
D'après les équations (2.5) et (2.13)
l'équation (2.17) devient.
Ow
7 u
at
OW
7
ax v
1702W ?~~ 7 ?~%
?~~ @
~~
Ow
=
Oy
(2.18)
0(to)
7
Ox
0(vc) 1
= V26) (2.19)
Oy Re
OW at 7
C'est l'équation de transport du tourbillon, sous
forme conservative l'équation peut s'écrire comme suit
:
Et d'après les équations (2.12) et (2.13),
on a :
|
? _?'
?% ?~ %`
3 at ax 7
|
o _?' ?~ w)= 1 V2)
(2.20)
Oy Re
|
a2*
(0 = 7
ax2
a2*
(2.21)
0y2
Les conditions aux limites pour la fonction courant et
sa dérivée normale sont fournies dans le cas de frontières
solides par les conditions d'adhérence. Pour l'équation de
transport du rotationnel, il n'existe pas de conditions aux limites physiques.
Il ya des relations basées sur des développements en série
de Taylor. Parmi ces relations, on trouve la relation de Woods donnée
par :
COP =
F F 3 =
~~ H'pqr 3 'pJ 7 ~ _?'
?d` ; 'pqr 7 G ~ (2.22)
p
El la relation de Jensen donnée par :
COP =
3 ;~~ H3'pq~7s'pqr 3 t'pJ 7 ~ _?'
F ?d` 7 G ~ p (2.23)
2.3 formulation (vitesse-tourbillon)
A partir de la définition de la fonction
tourbillon et de l'équation de continuité, on déduit deux
équations de Poisson pour la vitesse :
|
Ow
=
Ox
|
02v Ox2 7
|
02v
0y2
|
(2.24)
|
0(i)= 3>?~~ (2.25)
?~~ 7 ?~~
?~~@
?~
Le système (2.24 ,2.25) avec l'équation
(2.19) constitue les équations de Navier-stokes en formulation
(vitesse-tourbillon) sous forme conservative.
La plupart des méthodes proposées pour
la résolution numérique des équations régissant les
écoulements bidimensionnels d'un fluide visqueux incompressible,
utilisent comme inconnues, soit la fonction vitesse et la pression, soit la
fonction de courant et le tourbillon cette dernière formulation
présente des avantages importants sur les autres formulations. Ces
avantages sont bien connus:
(1) le champ de la vitesse est automatiquement à
divergence
(2) les propriétés mathématiques
des équations permettent de construire des méthodes de la
solution de manière plus simple. Le nombre d'équations est
inférieur à celui avec la formulation précédente
donc la résolution demande moins de temps de calcul. La
difficulté
classique associée à cette formulation est
le manque de conditions aux limites pour la fonction tourbillon.
2.4 Résolution numérique des
équations de Navier-stokes en formulation (4,, w)
|
Equation de transport du tourbillon
aco a (a* co)
ay
3 at ax +
|
a Cat co)
|
=
|
1
|
V2 co (2.26)
|
|
|
|
ay
|
Re
|
Equation Poisson de la fonction de courant
a2*
co = +
ax2
a2*
(2.27)
ay2
Les composantes de la vitesse sont définies par
:
|
u =
|
a*
ay
|
, V =
|
a*
(2.28)
ax
|
2.5 Les conditions aux limites :
- Condition initiales
On suppose le fluide au repos, à l'instant t=0
;
co = 0 , u = 0 , v = 0 (2.29)
A t=0 la paroi supérieure est mise en mouvement
avec une vitesse de module, IIVII = 1
(u = 1, v = 0) les autres parois sont fixes à
chaque instant (conditions d'adhérence). Les conditions aux limites sur
la fonction de courant (V)) sont déduites à partir les conditions
aux limites sur la vitesse.
- Conditions aux limites sur la fonction de
courant
a*
=
ax
a*
a (pour x = 0 ou x = 1) et y E [0,1] (2.30)
y
1 a* = 0
ax
a* _ _
- 1
y
pour x E [0,1] et y = 1 (2.31)
a* a*
=
ax ay
= 0 pour x E [0,1] et y = 0 (2.32)
Les conditions (2.30),(2.31)et(2.32) impliquent que V) =
cte sur le contour de la cavité. On peut sans perdre en
généralités poser V) = 0 sur le contour du
carré.
De même l'équation de poisson (2.27)
écrite sur le contour de la cavité donne immédiatement
:
Résolution d'ENS2D par : méthode des
différences finies et méthode spectrale
V?;'
UT ?~; ~ % ?;' T ?~; ~ G S
pour (x = 0 oux = 1 ety E [OM (2.33)
|
V?;'
UT ?~; ~ G
T?;' ?~; ~ % S
|
pour (y = 0 ou y = 1 et x E [OM (2.34)
|
A ce stade nous avons toutes les conditions aux limites
nécessaires pour la fonction de courant et ses dérivées
premières et secondes.
- Conditions aux limites pour la fonction
tourbillon
Les conditions aux limites pour la fonction tourbillon
sont calculées implicitement à partir d'un développement
limité de la fonction * au voisinage de la paroi. Soit 11Jp
la valeur de la fonction courant au noeud (i, 1) , et 14+1 la valeur de la
fonction courant au noeud (i, 2) .
On aura alors pour la fonction * la relation
:
Ay
ipp+i . ipp
7
_?' 7 ~~~~~
; >?~' 7 ~~~~~
>?~'
?~` ?~~@ ?~~ @ 7 G~~~~~€~
(2.35)
p p p
On sait que sur la paroiinférieuree on a
:
|
?'
TV ?~ ~ G
SUp ?' T ?~ ~ G p
|
g
|
02*
|
?~'
~ G g <~ .É6~5#~;IF...~ ?~~ ~ %p P
P P
|
|
Ox2
|
D'autre part on a :
V?~' ?~ >% 3 ?~'
?~~ @ ~ ?% ?~ 3 ? ?~ >?~'
U
a 02* 02 alp
Oy(Ox2). Ox2
(0y)
partout dans le domaine
T ?~~ ~ ? ?~~ @
S
Or
?;
_?' ?
?~` ~ G g g ?~ N?;'
?~; >?'
?~@ ~ G ?~; O ~ G
p 2
P P
Ce qui donne
a (0
ipp+i 3 ipp = 0 2302 cop 7 06303
lay) 7 0(0(0y)4)
P
On discrétise ZP‡
Pà[p par différence finie, on obtient :
aop 7 cop+i =
(AY)2 (16+1 3 16) 7
0(0302) (2.36)
6
|
Sur la paroi supérieure
:13 = --1
aop 7 cop+i =
|
6 6
~~~~~ H'p}r 3 'pJ 3 7 G ~
(2.37)
~~
|
Les relations (2.36) et (2.37) permettent de calculer
implicitement la fonction detourbillon n sur les deux parois
(supérieure,inférieure)..
On se base sur la formulation (fonction courant-
fonction tourbillon)puisqu'ellee comporte uneéquationn et une inconnue
en moins ce qui implique un temps de calcul moins importantet t
une économie en occupation mémoire. Plus le fait
quel'équationn de da la fonction courantest t munie de
conditions aux limites de types Dirichlet, ce qui permet une
convergencebeaucoup p plus rapide que celle
del'équationn de Poisson de la pression, qui est munie de condition de
type Newman.
Chapitre 3:Resolution des equations de
Navier-stokes
par la Methode des differences finies
3.1 Introduction
Les equations de Navier-Stockes utilisees pour
modeliser le comportement d'un fluide sont des equations aux derivees
partielles dont on ne connaît pas de solution analytique. Il est donc
usuel de recourir à la simulation numerique pour calculer des solutions
approchees de ces equations. Nous etudierons ici comment obtenir une
approximation numerique des equations de Navier-Stockes par la methode des
differences finies.
Les equations de Navier-Stokes contiennent donc bien
des difficultes. Avant d'en aborder la resolution, il convient de bien
maitriser les specificites des equations scalaires presentees ci-dessus et des
methodes numeriques associees.
Les methodes numeriques les plus utilisees en
mecanique des fluides sont les differences finies, et les methodes spectrales.
Il en existe d'autres (volumes finis, les elements finis....) dont nous ne
parlerons pas dans ce projet. Ces methodes transforment le problème
continu en un problème discret.
La methode des differences finies à laquelle
nous nous interessons est une combinaison de deux schemas d'ordre
O(h2) et O(h4) respectivement pour le tourbillon et pour
la fonction de courant.
Le schema à deux dimensions s'obtient en
generalisant le schema à une dimension et en
considerant les deux variables ( , y) d'une
manière separee. Le maillage du domaine se fait dans les deux directions
x et y
La methode des differences finies consiste à
remplacer les derivees partielles aux points du maillage par des developpements
de Taylor :
|
~ }r~ ~
~ qr~ ~ ~
|
~
3
|
OE~
|
7 OE~ ~
~ ~
|
~
~ 3
|
~
OE
|
of( i)
7 OE 7
|
hn a(n)f( i)
é 7 7
|
'~OEè~
7 '~OEè}r~
|
|
ç
ç~ ~
7 é 7
|
êç è
~3=~è OEè
ç~è~~ ~
|
|
ç
|
êç
è
|
Par combinaisons linéaires des
développements de Taylor, on exprime les dérivées
partielles en fonction des valeurs aux points de discrétisation. Ainsi,
en négligeant les erreurs de troncature.
Dérivée première
=
2h (fii 3 fi-1,j) (3.1)
aft; =
oXi
1 ;~ H~)~*}r 3 fi,j-1)
(3.2)
Oki
=
ay;
|
Dérivée seconde
|
a2fij
oXi2
a2fij
aYi2
|
1
? ~~ H~)}r~* 3 ;~)~* 7
~)qr~*J (3.3)
=
? ~~ H~)~*}r 3 ;~)~* 7
~)~*qrJ (3.4)
|
3.2 Principe de la méthode d'ordre O(h4)
C'est un schéma de différences finies
hermitien compact. La précision d'ordre 4 est obtenue avec seulement 3
points de discrétisation, en considérant les
dérivées comme des inconnues supplémentaires. La fermeture
du système est assurée en utilisant des relations additionnelles
obtenues à partir des développements en série de Taylor de
ces dérivées, ainsi que les expressions de l'équation
différentielle en trois points du maillage, au lieu d'un dans chaque
direction d'espace successivement.
É
Si on désigne par h le pas d'espace de la
discrétisation et par ~)~ ~) , fr, les valeurs de la fonction et ses
dérivées premières et secondes au noeud (i), on peut
écrire les relations tri diagonales suivantes :
F (3.5)
~)qr
É 7 ...~) É 7 )}r
É ~ ~ )}r 7 ~)qr~ 7
G~~€~
12
fi" 1 7 10fr 7 )}r
ÉÉ = h2 (fi-Fi 3
2fi 7 fi-1) 7 0(h4) (3.6)
Afin que le système soit bien
déterminé, il est nécessaire d'imposer des conditions aux
limites non pas seulement pour f, mais aussi pour ses drivées
premières et secondes.
Si N est le nombre de noeuds du maillage, on aura donc
en général un système de 3N équations à 3N
inconnues à résoudre, ce qui peut apparaître comme un
inconvénient mais étant donné que la précision
d'ordre O(h4) permet de diminuer le nombre de noeuds du maillage, le gain en
temps de calcul et en occupation mémoire reste encore important. Pour
résoudre l'équation de transport de la fonction du tourbillon, on
utilise la méthode O(h2), par contre pour résoudre
l'équation de poisson de la fonction de courant, on utilise une
méthode d'ordre 4.
3.3 Résolution de l'équation de poisson
de la fonction de courant
Pour résoudre l'équation de Poisson de la
fonction de courant, on utilise une méthode
d'ordre 0(h4) combinée avec un
schéma aux directions alternées implicites (ADI).
Par
conséquent, il a été nécessaire de
transformer l'équation elliptique de la fonction de courant
1p
en une équation parabolique en temps fictif T, en
ajoutant le terme a au second membre de
aT
l'équation de Poisson (éq.2.27)
|
02* Ox2 7
|
02*
|
= co 7
|
011J
|
(3.7)
|
|
0y2
|
OT
|
La solution stationnaire de cette équation
parabolique représente la solution de l'équation de la fonction
de courant. Chaque pas de temps fictif AT est résolu en deux
étapes ou demipas. Dans le premier demi- pas, l'intégration se
fait suivant une direction, et dans le second demi -pas dans l'autre
direction.
En considérant T = kAT et en désignant par
par IIJk+1/2et IIJk+1 les valeurs de IIJ
à chaque demi- pas fictif, on aura à résoudre les
équations suivantes :
|
~}r»
~'
?~~ @
|
~
7
|
~
~'
7 >?
?~~ @
|
~
~
|
%)*
%)*
|
7
7
|
=
|
|
>?)*
~}r»~
~'
?~~@
|
)*
}r
~'
?~~ @
|
ï~-
=
|
|
>?)*
|
>?)*
|
ï~à
|
~[ (3.8)
Z')*
~}r»~ 3 ')*
Z')* }r 3 ')* ~}r»~[ (3.9)
Premier demi- pas :intégrationn suivant
xOn suppose que les valeurs de ZP%o ~
Pà[)* dans l'équation (3.8) sont connues,
et on pose : i
k
) ~ %)* 3 >?~' 3 =
?~~ @ ')* ~
ï~-
)*
On aura alors :
~}r»~
>?~' ~}r»~
~ ™i 7 =
?~~ @ (3.10)
ï~- ')*
)*
A partir de cette relation on peut déduire les
valeurs de ZP%o ~}r»~ ~}r»~
et ZP%o
P-[)qr~*
P-[)}r~*
En portant dans la relation hermitienne (3.6), on
obtient une équation linaire qui permet de
k+1
2
.
déterminer les valeurs nodales iij
|
12 1 \ .c+1(2 24 1 \ 12 1 \++1/2
..rik+11(2
k(0x)2 Aix)1
k(AX)2 Aix) k(Ax)23 Aix)
= Si-i
7 113Si 7 Si+1
|
(3.11)
|
Les conditions aux limites étant de type
Dirichlet, on a à résoudre un système tridiagonal
pour
k+1
k+1
déterminer les valeurs iij
2 sur chaque ligne j = cte.une fois qu'on
connaît ltrij
2 on déduit
k+i
la valeur de (a211') 2 à
partir de la relation (3.10). ax2ij
Deuxième demi- pas : intégration suivant
y
k1
2 et moe
2
k+i
L'intégration suivant x nous a permis de connaître
*if
A partir de l'équation (3.9), on peut
déduire la relation :
}r
>?~' }r
~ * 7 =
?~~@ ')* (3.12)
ï~à )*
Avec
~}r»~ 3 =
* ~ %)*3 >?~'
?~~ @ ')* ~}r»~
ï~à
)*
De la même manière que
précédemment, en utilisant la relation hermitienne correspondant
à la dérivée suivant y, on obtient la relation
:
> =; @ ')*
~ï~~~ 3 = }r 7 > =;
@ ')~*qr
}r 3 > ;...
~ï~~~ 3 = ~ï~~~ 3 = @ ')*}r
}r ~ *qr 7 =G* 7 *}r
ï~à ï~à
ï~à
Pour i= cte, et on a un système tridiagonal
à résoudre pour déterminer *71 }r
La dérivée second Zùú
ù·[oe est calculée à partir de
la relation (3.12).
3.4 Calcul des composantes de la vitesse
Le champ de * étant déterminé, les
composantes u et v de la vitesse sont déduites de la relation
hermitienne (3.5), laquelle nous permet d'écrire les schémas
suivants :
}r
>?~'
?~~@ ')* }r
~ * 7 =
ï~à
)*
a*
Calcul de = V
ax
_?' 7 ... _?' 7 _?' ~ F
?~` ?~` ?~` ï~ H')}r~*
3 ')qr~*
J
)}r* )~* )qr~*
Calcul de alp =
--u
ay
(aa 1P)n7 +
4 (191P)n7 + (19
a1P)n = 3 * 4+1 Vi 4-1)
y i,]+1+i ay . Li . y ti-1
A
y 3.5 Résolution de
l'équation de transport du tourbillonnUn schéma directionnel
d'ordre 2 a été choisi pour résoudre l'équation de
transport du utourbillon..C'est un schéma semi-implicite du type
directions alternées. Ce choix est motivé par la
arecherche d'une simplicité dans l'affichage des conditions
aux limites, des critères de stabilité emoins
restrictifs et d'un gain en temps de calcul. Ce schéma sera
précis d'ordre 2 en espace, et tprécis d'ordre 1 en
temps..Sous forme conservative l'équation de transport du
tourbillon est donnée par ::
|
aw a (ay ,
4(0)
\
at ax7+
|
°Cat
6))
|
=
|
1>?~% (3.13)
?~~ 7 ?~%
?~~ @
~~ j
|
|
ay
|
aw
at : Terme
''évolution
|
a_?' ?~ %` /
7+
ax
|
a Cat
6))
|
: Terme de convection
|
|
ay
|
1 (026)
02(o)
Re aX2
: terme de diffusion
7
+ ay2
On note :
|
S. =
|
a
ox ' 8 Y =
|
0 02 02
~ (- ~ ~ ~ (à ~ ~ ?~~
?~ ?~~
|
Le calcul de la fonction tourbillon se fait en deux
etapes ou demi -pas. Si on designe par cn+1/2 la valeur du tourbillon calculee
au premier demi- pas, et par con' la valeur de w au deuxième
demi- pas (valeur de to au (n 7 1)At ), le schema A.D.I s'ecrit de la
manière suivante :
1er demi-pas :
COn+1/2 3 con
3 (- >%}r»~ _?' Z(-
~H%}r»~J 7 (à ~~%~[
?~` @ 7 (à >% _?' @ ~ =
?~`~~ ~~
;
2ème demi-pas :
1
con+1 3 COn+ 2
(- >%}r ~ _?' @ 7 (à >%}r
_?' @ ~ = >(- ~ _%}r ~%}r@
?~` ~` 7 (à
?~`~~ ~~
;
Au premier demi- pas, les valeurs inconnues
COn+1/2 n'apparaissent
que sous les operateurs
8x et 4 .
Or ces operateurs ne font intervenir que les points
immediatement voisins et situes sur la même horizontale (seul x
varie).
On pourra donc resoudre sur chaque ligne horizontaley =
cte.
De même au deuxième demi -pas les valeurs
inconnues co' n'apparaissent que sous les operateurs
Sy etSy2, on pourra donc resoudre sur chaque
ligne verticale x=cte. La precision d'ordre 2 en espace est obtenue de deux
manières differentes suivant la façon de discretiser les termes
de convection.
Schema directionnel d'ordre 2 differences centrees
pour les termes de convection Soient f et g deux fonctions
S. (resp Sy) le pas de discretisation en
s(resp en y)
fmn (respgm,n) la valeur de f(resp
g) au noeud (m,n)
Alors l'operateur de convection discretise approchant
a ù (resp ù ù· ) prendra la forme
suivante.
|
6x (fg) =
|
fi+Li
· gi+Li 3 fi-Li.
gi-Li
|
|
2.6a
|
|
6y(fg) =
|
~)}r~*I gi+i,j 3
~)qr~*I 4)qr~*
|
|
2Ay
|
Schema directionnel d'ordre 2 utilisant des differences
decentrees dans le sens de la vitesse pour les termes de convection
Cette technique de discretisation des termes de
convection permet d'assurer la preponderance de la diagonale principale des
matrices resultantes. Ce genre de schema conduit à un temps de calcul
moins important que le schema centre, et ces conditions de stabilite sont moins
restrictives. Par consequent, il a ete choisi pour la resolution numerique de
l'equation de transport du tourbillon.
On definit :
~)* ~ ZP%o Pà[)* pi; = 1 si uji < o et 0
ailleurs
~)* ~ ZP%o P-[)*
aii = 1 si vij > 0 et 0 ailleurs
Les operateurs de convection discretises approchant
a ù et ù ù· prendront la forme
suivante : 1er demi -pas :
}r»~[ }r»~3)*%)* }r»~[
(- >%}r ~ _?' §)* Z)*%)*
}r»~3~)qr~*%)qr~* H= 3 §)*J Z~)}r*%)}r*
?~` @ ~ 7 Ax Ax
(à >% _?' H= 3
#172;)*JH~)*%)*
3~)~*qr%)~*qr
J J
#172;)*H~)~*}r%)~*}r
3)*%)*
?~` @ ~ 7 Ay Ay
2ème demi - pas :
(-
(à
>%}r ~
|
>%}r
|
_?'
|
@
@
|
~
~
|
|
?~`
_?'
|
|
?~`
|
}r»~[ }r»~[
}r»~3)*%)*
H= 3 §)*J Z~)*%)r*
}r»~3~)qr~*%)qr~*
§)* Z~)}r~*%)}r~*
7 Ax Ax
}r J }r3)*%)* }rJ
#172;)*H)*%)* }r3~)~*qr%)~*qr H= 3
#172;)*JH~)~*}r%)*}r
7
~~ ~~
Chapitre 4:Résolution des équations de
Navier-stokes
par la Méthode spectrale
4.1 Introduction
Les méthodes spectrales sont un outil
très puissant dans le cadre de la résolution numérique des
équations aux dérivées partielles. Elles consistent
à utiliser des changements d'espaces sur la fonction
étudiée. Nous présentons une méthode spectrale
colocationTchebychev pour la résolution numérique des
équations de Navier-stokes pour le cas de l'écoulement
instationnaire d'un fluide visqueux (de viscosité cinématique
v et densité p) dans une cavité
à paroi supérieure entrainée. L'écoulement est
supposé bidimensionnel. La vitesse sur la paroi supérieure est
bien définie u=1 et v=0. La difficulté majeure quand on approche
numériquement les problèmes avec les méthodes spectrales
est liée à la complexité du domaine de calcul qui varie au
cours du temps. Une autre difficulté qui se pose lors de la
résolution d'une équation non linéaire est liée
à un problème technique dû au phénomène
d'aliasing. Les méthodes spectrales, dans leur formulation d'origine
sont mal adaptées pour ce genre de situations
Les méthodes spectrales font partie des
méthodes à résidus pondérés dans lesquelles
les approximations sont définies en terme d'une série
(développement en série) de telle manière qu'une certaine
quantité appelée erreur ou résidu qui doit être nul
est forcé d'être nul de manière approximative ceci est
obtenu en utilisant le produit scalaire définit par:
#177;
(f,g)0 = f f 4 o dx
²
Ou f(x) et g(x) sont deux fonctions définies sur
l'intervalle [a, f3] et w(x)
est une fonction poids donnée.
4.2 Méthodes spectrales pour l'approximation
d'une fonction
L'approximation d'une fonction donnée f(x) est
représentée sous forme d'un développement en série
tronquée de la manière suivante :
N
fN(xP t) = 1 f~k(t)
(Pk(x)
k=0
Avec
f~k(t): Coefficients spectraux qui sont des
inconnues à déterminer (pk( x) : Les fonctions de base qui sont
données
Le choix de la fonction poids et de la fonction de base
dépend de la nature du problème, en effet, si la solution n'est
pas périodique à tout ordre et si le domaine est borné
(par exemple
-1<,x <1), dans
ce cas on utilise le polynôme de Tchebychev . 4.3
Généralité sur la méthode
collocation-Tchebychev
4.3.1 Les propriétés principales des
polynômes de Tchebychev
|
Les polynômes de Tchebychev (de première
espèce) Tk(x), sont tels que :
|
|
|
Tk(x) = cos kx
|
|
(4.1)
|
|
Les Tk sont des polynômes de degrés k
liés par triplets par récurrence :
|
|
|
To(x) = 1
|
|
(4.2)
|
|
Ti(x) = x
|
|
(4.3)
|
|
Tk+1 = 2xTk(x) -- Tk-1(x)P pour k
|
1
|
(4.4)
|
Ces polynômes de Tchebychev sont orthogonaux sur
l'intervalle [-1; 1] avec la fonction poids qui est donné
par:
co(x) = 1 (4.5)
il -- x2
la propriété d'orthogonalité est
:
0 si k * m
11 Tk(x)Tm(x) dx = TC si k = m =
0
(4.6)
J-1 V1 - x2 TC
2 si k = m * 0
4.3.2 Quadrature de Gauss-Lobatto
En décomposition Tchebychev, la quadrature de
Gauss-Lobatto, exacte pour des polynômes de degré 2N - 1 ou moins,
impose que les points de collocation soient les zéros de (1 -
x2)TN '(x).Ces points ont pour abscisses
:
iTC
xi = - cos (N) pour i E [0, N]
(4.7)
Et les poids de cette quadrature sont alors :
TC
CO0 = CON = 2N (4.8)
COi =
TC
N pour i E [1, N - 1] (4.9)
4.3.3 Représentation d'une fonction sur les
polynômes de Tchebychev
Sur (N + 1) points de collocation repartis sur
l'intervalle [-1; 1], le polynôme d'interpolation fN(x) d'une fonction
f(x) est :
N
fN(x) = f~~(t) Tk(x)
(4.10)
k=0
fN(x) : Fonction donné
f~k(t) : Les coefficients spectraux sont
fixés,
En particulier, les opérations de
dérivation, intégration et interpolation numériques, se
ramène à des manipulations des polynômes de Tchebychev
Tk(x), dont les propriétés et particularités sont bien
connues.
Connaissant les valeurs prises par f(x) aux points de
collocation xi, yi = f(xi), on voit que la relation ci-dessus peut s''ecrire
sous la forme d'un produit matrice vecteur : P C = Y, où C est le
vecteur des coefficients spectraux f~~(t)et P est la matrice de
passage dont les éléments sont P(i, j) = Tj(xi).
En d'autres termes, connaissant les valeurs nodales yi,
on a accès aux coefficients spectraux ~~~~~~ (par multiplication du
vecteur y par l'inverse de P) et réciproquement.
On qualifié la matrice P de matrice de passage
de l'espace spectral à l'espace physique (puisque son application sur
les coefficients spectraux permet d'obtenir les valeurs nodales). Son inverse,
P-1, permettant l'opération inverse (l'obtention des
coefficients spectraux à partir des valeurs nodales) est
qualifiée de matrice de passage de l'espace physique à l'espace
spectral.
4.3.4 Matrice de passage pour les points de collocation de
Gauss-Lobatto
Les éléments de la matrice de passage P et
de son inverse P-1 sont tels que :
Pii = Tj(xi) (4.11)
~)* qr ~co
·
Ti(x
·) (4.12)
(Ti, TON
Ou les xi sont les points de collocation, les wi sont les
poids et (Ti, Ti)N le produit scalaire discret de la quadrature
considérée.
4.3.5 Dérivation numérique via le
développement sur les polynômes de Tchebychev
Si on note par u le vecteur dont les composantes sont les
valeurs de uN aux points de collocation et par U(p) le vecteur des
valeurs de la dérivée Pième de u, on a
:
N N
|
141(xi) = uk
|
r
£~ É~~)~ ~ ~»~
|
Tk(Xi) (4.13)
|
k=0 k=0
Avec
N
p=k+1
(p+k)pair
~»~ r ~ ;
ck P up
(4.14)
ck donné par :
t2 si k = 0
ck =
1 si k > 1
4.4 Exemple :
On considère une fonction f(x) avec x E
[--l,l]
Soit
f(x) = x-il -- x2
Dans ce cas oil on a une fonction qui n'est pas
périodique donc on utilise un développement polynomial et puisque
le problème est posé en domaine borné, on utilise
polynôme de Tchebychev.
Soit
N
fN = / f~k Tk
k=0
Pour la méthode de collocation, le résidu
est exactement nul en certains points (points de collocation) par contre la
méthode de Gallerkin, le résidu est nul au moyenne.
RN(xi) = 0
Avec : xi = cos (sN) j = 0, ....,
N
Ce qui implique :
N
/ f~k Tk(xi) = f(xi)
k=0
On obtient donc un système algébrique de
N+1 équations à N+1 inconnues f~k , l'existence de la
solution pour ce système implique que :
det (Tk(xj)) * 0
C'est une première condition que doivent
satisfaire les points de collocation. Faisons une application numérique
pour le cas oil on prend N=4 :
4
/ f~k Tk(x*) = f(xJ)
k=0
Sous la forme matricielle : El Y = b Avec :
|
£r
|
£~~~~~
|
£~~~~~
|
£€~~~~
|
|
£r~~r~
|
£~r
|
£~r
|
£€~~r~
|
|
£r ~
|
£~~~~~
|
£~~~~~
|
£€~~~~
|
|
£r~
|
£~~~~~
|
£~~~~~
|
£€~~~~
|
|
£r~~€~
|
£~~~€~
|
£~~~€~
|
£€~~€~
|
Æfa)
~~~ É
É Æ Æ£
É
r Å ~ Z!"# À €[
È
È r £r
Å È Å È Å
~~~ È , ~ , et Á ~
Å ~~~~~ È ~ È £~~~~~
Å
Y=
Å ~ Z!"# À ~[
~~~ È Å ~~~~~ È Å
£~~~~~
È ~ Z!"# ~À È
Å
€ [
~~€Ç Ä ~~~€~Ç Å
È Ä £~~~€~
Ä ~~!"# ° Ç
Application numérique :
|
1
0.7071
0
0.7071
--1
|
1
0
--1
0
1
|
1
0.7071
0
0.7071
--1
|
1
1
Å
El = 1
Å 1
Ä 1
1 É
--1.7071
È
0 È
--0.2929 È 3 Ç
Et puis que : Y=M-1b Implique et finalement
:
~
f ~ É
Æ G
r È GIFËF É
~~~ È ~ G
ÅÅÅ
ÈÈÈ
~~~ È 3GIFËF
È Ä
0 Ç
f4
Dans le cas oil on ne veut pas inverser la matrice du
système, on peut utiliser la relation discrète
d'orthogonalité pour les polynômes de Tchebychev pour exprimer
explicitement les valeurs desfk.
Soit
|
Avec
|
~
fk =
|
N
!Ì =
; Ck
J=0
|
f(ci)Tk(xi) k = 0,
·
· , N
|
2 si j = 0
ci = Í1 si 1 j < --1 2 si j =
N
Appliquons maintenant la méthode spectrale
à la résolution des équations de Navier
Stokes.
4.5 Projection de la méthode spectrale
Problème à résoudre ;
0(uw)
7
Ox
0(vc) 1
= V2c (4.15)
ay Re
aco at 7
4.5.1 L'intégration temporelle
|
co = a2ip7
ax2
|
a2*
(4.16)
0y2
|
Le schéma d'intégration temporelle que nous
avons utilisés est de type semi-implicite à pas multiples. Pour
illustrer ceux-ci, considérons l'équation:
Ow
at
= ,C(w) 7 N(w) (4.17)
Ou L et N sont des operateurs respectivement
linéaire et non linéaire. Les
déférentes
discrétisations des trois termes de
l'équation, permettant de construire le système linéaire
à
résoudre pour déterminer c"l(c.-`a-d.: le
champ au temps présent: (n+1).At, en
fonction
des valeurs aux temps précédents
n.At, (n - 1).At, . . .). Adams-Bashforth
Crank-Nicholson (ordre 2 en temps)
Ow
at Ð
1 1
co(n+1) 3 co(n)
(4.18)
At At
1 1
£(w)2 L(6)(n+1)) 7 2 46)(n)) (4.19)
3 1
~~%~ Ð ; ~H%~~J 3 ; ~H%~qr~J
(4.20)
d d d7= d3=
= =
~~ %d7= 3 ~~ %d ~ ^ ; %d7= 7 ^ ; %d 3 F ;
>~ ?% 7 ~ ?% @ 7 ; = N ?% 7 ~ ?% O
?~ ?~ Ox ay
|
_~
|
d d d7= d3=
^~~` %d7= ~ 3 ;
; ^~~ %d 3 %d 7 F ^ >~ ?% 7 ~
?% @ 3 ^ = N ?% 7 ~ ?% O
?~ ?~ ?~ ay
|
Après application du schéma temporel, le
champ c0(n+1) est tel que:
0 3 oco(n+i) = s(l+i)
2
Oil k =vAT est une constante (positive) et
Sn+1 un terme source (connu).
Ce système est d'abord réduit aux points
de collocation intérieurs (procédure d'injection des Conditions
aux limites), puis résolu (via les diagonalisations successives de
l'opérateur), avec reconstruction des valeurs aux
frontières.
4.5.2 Les conditions aux limites
Les conditions aux limites pour chaque variable, suivant
chaque direction peuvent s'écrire sous la forme :
OC
aC + 13 an = y
Les six coefficients sont, dans l'ordre: a_,
a+,13_, 13+, y_ety+ Les indices -- et + font
|
référence aux frontières en --1 et
+1.
-les conditions aux limites sur la vitesse :
|
|
|
|
|
a_ = 1, 13_
|
= y_ =
|
0
|
--> U(x = --1)
|
= 0
|
|
a+ = 1, 13+
|
= y+ =
|
0
|
--> U(x = +1)
|
= 0
|
|
a_ = 1, 13_
|
= y_ =
|
0
|
--> U(y = --1)
|
= 0
|
|
a+ = 1,13+
|
= y+ =
|
1
|
--> U(y = +1)
|
= 1
|
|
a_ = 1, 13_
|
= y_ =
|
0
|
--> V(x = --1)
|
= 0
|
|
a+ = 1,13+
|
= y+ =
|
0
|
--> V(x = +1)
|
= 0
|
|
a_ = 1, 13_
|
= y_ =
|
0
|
--> V(y = --1)
|
= 0
|
|
a+ = 1,13+
|
= y+ =
|
0
|
--> V(y = +1)
|
= 0
|
Les conditions aux limites pour la fonction courant
sont les même que celle de V. Par contre pour la fonction tourbillon ig
est toujours nul et y change à chaque itération puisqu'il est
fonction de * .
-paroi supérieure x E [--1,1], Y =
1
|
o(i, NY) =
|
2*(i, NY) -- *(i, NY --
1)
|
|
Dye
|
|
-paroi inférieure E [-1,1],y = --1
|
|
|
w(i3 O) =
-paroi droite
|
2m(i3O) 3 m(i,1)
|
|
ïb~
|
|
w(NX,j) =
|
2m(NX,j) 3 m(ÙÚ 3 1,])
|
|
-paroi gauche
|
LX 2
|
|
w(O,j) =
|
2m(O,j) 3 m= Û
|
|
LX 2
|
4.5.3 Procédures de résolutions des
équations de Navier-Stokes
Condition initiales pour la
La fonction courant et
tourbillon
Résolution de l'équation
de
poisson
Obtention les composantes
u et v
L'équation de transport
Visualisation des résultats
Fig. (4.2) : Procédures de résolutions
des équations de Navier-stokes
Les résultats que nous allons présenter
dans le chapitre qui va suivre concernent :
1) La résolution des équations de
Navier Stokes bidimensionnelles dans une cavité carrée à
paroi supérieure entrainée en variables fonction de courant
fonction tourbillon par la méthode des différences finies d'ordre
(O(H2)-O(H4)).
2) La résolution des équations de
Navier Stokes bidimensionnelles dans une cavité carrée à
paroi supérieure en variables vitesse pression par méthode
spectrale de collocation Tchebychev.
3) La résolution des équations de
Navier Stokes bidimensionnelles dans une cavité carrée à
paroi supérieure en variables fonction de courant fonction tourbillon
par méthode spectrale de collocation Tchebychev.
Chapitre 5 : Résultats et discussion
5.1.Formulation fonction de courant fonction tourbillon
: méthode des différences finies d'ordre (O(H2)-O(H4)).
Dans cette section, on présente les lignes iso
valeurs de la fonction de courant pour le cas d'une cavité à
paroi supérieure entrainée avec une vitesse horizontale u=1. Le
calcul est fait avec la méthode des différences finies
(O112-O114) avec un schéma ADI. Ces résultats sont
comparés à ceux donnés pour le cas d'une formulation
vitesse pression.
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
1
0 0.2 0.4 0.6 0.8 1
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
Fig.(5 .1) : Re=100, LX=LY=1, NX=NY=29,
DT=0.05, t=25
Champs de vitesse avec la formulation vitesse pression pour
Re=100
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig.(5.2) : Re=100, LX=LY=1, NX=NY=29, DT=0.05,
t=25
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1
1
0.9
0.8
0.7
0.6
Y
0.5
0.4
0.3
0.2
0.1
0
Avec la formulation Fct de courant- Tourbillon Avec la
formulation Vitesse- Pression
Composante U de la vitesse en X=0.5
Fig.(5.3) :Re=100, LX=LY=1, NX=NY=29, DT=0.05,
t=25
Composante U de la vitesse en x=0.5 en fonction de Y
Fig.( 5.4) : Re=100, LX=LY=1, NX=NY=29, DT=0.05,
t=25
Sur la figure (5.3) on compare le résultat sur
la composante horizontale de la vitesse en x=0.5 en fonction de y. Les calculs
sont faits en formulation vitesse pression et en formulation fonction de
courant tourbillon avec la méthode des différences finies. Ces
résultats sont comparables à ceux obtenus par d'autres auteurs
(figure 5.4) [12] , [13]
5.2.Formulation vitesse pression: méthode
spectrale de collocation Tchebychev
Dans cette section, on présente les lignes iso
valeurs de la fonction de courant pour le cas d'une cavité à
paroi supérieure entrainée avec une vitesse horizontale u=1. Le
calcul
est fait avec la méthode de collocation de
Tchebychev pour le cas où les équations sont exprimées en
variables vitesse pression. Une fois qu'on a le champ de vitesse à
l'instant voulu et pour le nombre de Reynolds voulu, on applique le rotationnel
à ce champs de vitesse pour en déduire le tourbillon, puis on en
déduit la fonction de courant par résolution de l'équation
:
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0.1
0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
Fig.(5.5) : Re=100, LX=LY=1, NX=NY=20, DT=0.005,
t=25
Calcul global spectral, vitesse pression
Résolution d'ENS2D par : méthode des
différences finies et méthode spectrale
0
0 0.1 0.2
0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
Fig.(5.6) : Re=100, LX=LY=1, NX=NY=20, DT=0.005,
t=25
Calcul de fonction de courant pour chaque
itération spectral, vitesse pression
0
-1 -0.8 -0.6 -0.4
-0.2 0 0.2 0.4 0.6 0.8 1
Y
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
1
programme avec psi à chaque pas programme avec psi
global
Composante U de la vitesse
Fig.(5.7) : Re=100, LX=LY=1, NX=NY=20, DT=0.005,
t=25
Composante U de la vitesse en x=0.5 en fonction de Y
Sur les figures
(5.5) et (5.6), on compare le résultat sur la fonction de courant. Les
calculs
sont faits en formulation vitesse pression avec la méthode
spectrale de collocation
Tchebychev respectivement pour un calcul de fonction de
courant à chaque pas de temps puis seulement pour la dernière
itération.
Sur la figure (5.7), on donne la composante
horizontale de la vitesse en x=0.5 en fonction de y en formulation vitesse
pression avec la méthode spectrale de collocation Tchebychev
respectivement pour un calcul de fonction de courant à chaque pas de
temps puis seulement pour la dernière itération.
5.3.Formulation fonction de courant tourbillon:
méthode spectrale de collocation Tchebychev
Dans cette section, on présente les lignes iso
valeurs de la fonction de courant pour le cas d'une cavité à
paroi supérieure entrainée avec une vitesse horizontale u=1. Le
calcul est fait avec la méthode de collocation de Tchebychev pour le cas
où les équations sont exprimées en variables fonction de
courant tourbillon.
|
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
|
|
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig.(5.8) : Re=100, LX=LY=1, NX=NY=20, DT=0.005,
t=25
On peut remarquer que nos résultats pour les
deux premières sections sont relativement comparables par contre les
résultats obtenus avec la méthode spectrale avec la formulation
fonction de courant tourbillon présentent un problème que nous
n'avons pas pu lever faute de temps mais que nous essayerons de résoudre
par la suite.
Conclusion générale
Dans ce travail nous avons présenté la
résolution des équations de Navier stokes dans différentes
formulations et avec différentes méthodes numériques
:
Pour les formulations, il s'agit des formulations vitesse
pression et fonction de courant tourbillon.
Pour les méthodes numériques il s'agit
des méthodes de différences finis d'ordre O(H2) - O(H4) et de la
méthode spectrale avec une collocation Tchebychev. Le problème
test qu'on a considéré est celui de la cavité à
paroi supérieur entrainée et les résultats concernent les
isovaleurs de la fonction de courant pour les cas suivantes :
- Différences finies O(H2) - O(H4) en (fonction de
courant -tourbillon)
w = klt0000001 a01
- Spectrale en (vitesses-pression) et
iji calculé par et à chaque
itération.
m ~ 3w
w = klt0000001 a01
- spectrale en vitesse pression iji
calculé par et seulement après
m ~ 3w
la fin des calcules sur la vitesse
- spectrale en (fonction de courant-
tourbillon)
Les résultats pour ce dernier cas ne sont pas
satisfaisants et doivent être revus.
Références bibliographiques
[1] R.COMOLET -Mécaniques Expérimentales
des Fluides Tome II -- Dynamique des Fluides Réels, Turbomachines 1994
,4éme édition.
[2] Johnny de Jesús Martinez and Paulo de Tarso
T. Esperança, A Chebyshev Collocation Spectral Method for Numerical
Simulation of Incompressible Flow
Problems, July-September 2007.
|
|
[3] John P. BOYD, Chebyshev and Fourier spectral method,
Dover Publications, 2001, Second Edition.
[4] Maciej MATYKA, Solution to two-dimensional
Incompressible Navier-Stokes Equations with SIMPLE, SIMPLER and
Vorticity-Stream Function Approaches Driven-Lid Cavity Problem: Solution and
Visualization. 30 czerwca 2004 roku.
|
|
[5] Eric GONCALVES, Résolution numérique,
discrétisation des EDP et EDO, septembre 2005
[6] William RAPHA, Simulation numérique des
équations de la mécanique des fluides à l'aide de la
méthode des différences finies. Ecole Centrale Paris
|
|
[7] David LOUREIRO, Solving the Laplacian equation with
Boundary conditions through the pseudospectral method , 29 mai 2007
[8] Olivier Botella , Résolution des
équations de Navier-Stokes par des schémas de projection
Tchebychev ,Rapport de recherche n°3018 ,Octobre 1996
|
|
[9] Frédéric DABBENE et Henri PAILLERE,
Initiation à la simulation numérique en mécanique des
fluides : Eléments d'analyse numérique. Cours ENSTA MF307, 6 juin
2003
[10] Claude Delannoy. Programmer en Fortran 90. Guide
complet. Eyrolles, 2002.
|
|
[11] Berrada Mohamed, Schéma numérique de
différences finis, Septembre 13, 2006
[12] Katuhiko Goda,`Amultistep technique with implicit
difference schemes for calculating two or three dimensional cavity flow', J. of
Comput. Phys, Vol. 30. P. 76-95, 1979
|
|
[13] R. Burggraf, `Analytical and numerical studies of
the structure of steady separated flow', J. of Fluid Mechanics, Vol. 24, P.
113-151, 1966.
| 
