pubAchetez de l'or en Suisse en ligne avec Bullion Vault


Home | Publier un mémoire | Une page au hasard

Résolution des équations de Navier Stokes bidimensionnelle par méthode des différences finies et méthode spectrale


par Rachid Benrazouk
Université Hassan 1er - master 2009
Dans la categorie: Sciences
   
Télécharger le fichier original

Disponible en mode multipage

 

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)

dV

001

000000001

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

au Ox 7

Ov

O= 0 (1.16) y

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 (

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.