0

MODELISATION ET SIMULATION NUMERIQUE DES SYSTEMES

ENERGETIQUES PAR ELEMENTS ET VOLUMES FINIS

Partie : volumes finis

TABLE DES MATIERES

Chapitre : 1
I-INTRODUCTION ..................................................................................................................................... 2
II-FORMULATION DU PROBLEME A RESOUDRE ..................................................................................... 3
II-1. Type de problèmes ....................................................................................................................... 3
II-2. Mise en équations ........................................................................................................................ 3
III-FORMULATION GENERALE DES EQUATIONS DIFFERENTIELLLES DE CONSERVATION ........................ 4
III-1. Remarque .................................................................................................................................... 4
III-2. Forme générale d’une équation de conservation ....................................................................... 4
IV- DISCRETISATION DE L’EQUATION GENERALE DE CONSERVATION .................................................... 6
IV-1. DOMAINE DE DISCRETISATION.................................................................................................. 6
IV-2. DISCRETISATION .......................................................................................................................... 7
2.1- Discrétisation de l’équation de conservation générale. .............................................................. 7
2.2 -Calcul des termes de type valable pour les faces (e et n) ................................ 9
2.3-Calcul des termes de type valable pour les faces (s et w) ............................... 10
2.4 - Calcul exacte du flux généralisé (solution exacte) ................................................................ 11
2.5- Linéarisation du terme source ............................................................................................... 12
2.6- Discrétisation de l’équation à deux dimensions ........................................................................ 13

Chapitre 2
1- DISCRETISATION DES EQUATIONS _________________________________________________ 17
1.1- Choix du maillage ___________________________________________________________ 17
1.2- Discrétisation des équations de conservation sur un volume fini ______________________ 19
2- TRAITEMENT DES CONDITIONS AUX LIMITES ________________________________________ 20
2.1- Conditions aux limites à la paroi _______________________________________________ 20
 1

2.1.1-Vitesse : _______________________________________________________________ 20
2.1.2-Température ___________________________________________________________ 21
2.2- Sortie ____________________________________________________________________ 22
2.3- Entrée ____________________________________________________________________ 22
3- METHODE NUMERIQUE DE RESOLUTION ___________________________________________ 23
4-COUPLAGE PRESSION-VITESSE, ALGORITHME DE RESOLUTION ___________________________ 24
4.1. COUPLAGE PRESSION-VITESSE ________________________________________________ 24
4.2. ALGORITHME SIMPLER _______________________________________________________ 24
4.2.1 Algorithme SIMPLE ______________________________________________________ 24
4.2.2. Algorithme SIMPLER _____________________________________________________ 26
4.3. CONDITIONS AUX LIMITES POUR LA CORRECTION DE PRESSION ET POUR LA PRESSION ____ 27
4.4 SOUS-RELAXATION –CONVERGENCE_____________________________________________ 29
4.5 TEST DE CONVERGENCE ______________________________________________________ 30
4.6 CONCLUSIONS ______________________________________________________________ 30

Chapitre 3 : Applications
Résolution d’un problème de convection libre en 2.D...................................................................... 32
I- Formulation du problème : ........................................................................................................... 32
II- Solution numérique .................................................................................................................. 34
1. Discrétisation des équations gouvernantes par la méthode du volume de contrôle (ou volume
finis) ................................................................................................................................................... 34
a- Equation de mouvement suivant x ......................................................................................... 36
b- Equation de mouvement suivant y ........................................................................................... 36
c- Discrétisation de l’équation de l’énergie................................................................................... 37
d- Discrétisation des conditions aux frontières : ........................................................................... 37
2-Résolution de équations discrétisées............................................................................................. 37
Algorithme de Résolution (SIMPLER →SIMPLEC) .......................................................................... 37
 2

Chapitre : 1

I-INTRODUCTION

Les trois grandes méthodes numériques utilisées dans les codes de calcules numériques sont
les volumes finis, les différences finis et les éléments finis.
La méthode des différences finis consiste à discrétiser les équations de continues aux nœuds
d’un maillage prédéfini en calculant chaque dérivée partielle à l’aide de séries de Taylor
tronquées pour obtenir des équations linéaires reliant la valeur des inconnues en un nœud aux
valeurs de ces mêmes inconnues aux nœuds voisins.
La technique des éléments finis discrétise l’espace à l’aide d’éléments géométriques simples
(triangles ou quadrilatères en général). Comme elle permet de modéliser des géométries très
complexes, elle est parfois préférée à la méthode des différences finis.
La méthode des volumes finis famille de méthodes numériques reposant sur une partition du
domaine en volumes de contrôle, largement décrite par Patankar (1980), consiste à discrétiser
le domaine de l’écoulement en une multitude de volumes de contrôle puis d’effectuer des
bilans (masse, de quantité de mouvement…..) sur ces petits volumes, pour lesquels seule une
valeur moyenne des inconnues est calculée. Pour des lois de conservation, les échanges entre
volumes de contrôle se font par l'intermédiaire de flux, de manière automatiquement
conservative.
L’avantage déterminant des volumes finis par rapport aux autres méthodes est qu’ils sont
conservatifs, tout ce qui sort d’un volume de contrôle entre dans un autre.
Les méthodes de volumes finis sont utilisées depuis longtemps pour la simulation numérique
en mécanique des fluides. Elles permettent d'approcher des solutions presque nécessairement
discontinues, tout en conservant de bonnes propriétés (conservativité , précision, monotonie,
etc...). Ces méthodes ont trouvé une seconde jeunesse avec leur application à
l'électromagnétisme, notamment pour les cas hétérogènes, et les applications de la mécanique
des fluides où les domaines sont déformables (interactions fluide-structure, écoulements
moteur, etc...)
La méthode des volumes finis a été alors proposée pour traiter ces problèmes. Son principe est
simple : on considère comme espace d'approximation les fonctions constantes par morceaux,
ces morceaux ou cellules ou volumes finis étant choisis par l'utilisateur et issus par exemple
d'un maillage de type éléments finis - notamment non structuré- ce qui permet de prendre en
compte des géométries complexes. Cette vision des volumes finis est néanmoins réductrice, et
fortement influencée par l'approche éléments finis (les volumes finis peuvent être vus comme
des éléments finis de type P0)
 3

Cependant, les méthodes en volumes finis prennent un autre sens quand on s'intéresse à une
loi ou à un système de lois de conservation. Les valeurs numériques dans chaque cellule
peuvent être vues comme des approximations de valeurs moyennes sur la cellule, dont les
variations ne dépendent que des flux aux bords de la cellule. Ainsi, il suffit de construire des
fonctions de flux numériques, donnant de bonnes approximations de ce qui passe d'une cellule
à ses voisines. Automatiquement, la méthode produite est conservative (ce qui sort d'une
cellule rentre exactement dans sa voisine).
Ces méthodes peuvent être appliquées dans la grande majorité des domaines d'application
abordés par le projet, des écoulements multi-espèces ou multiphasiques
Enfin, il est important de rappeler ici que les méthodes de volumes finis s'adaptent très
simplement pour la simulation numérique de lois de conservation en formulation
Arbitrairement Lagrangienne Eulérienne (ALE), en maillages dynamiques à topologie
constante, passage presque nécessaire pour la simulation de la plupart des interactions fluide-
structure (où le domaine fluide, complémentaire de la structure, varie avec le temps).

II-FORMULATION DU PROBLEME A RESOUDRE

II-1. Type de problèmes
Il s’agit d’étudier un problème d’écoulements laminaires et du transfert de chaleur convectif.
Exemple :

Entrée Sortie
U0, V0 , T0

Plaque
x

- Domaine fluide limité par deux plans parallèles

- Fluide Newtonien incompressible.
Le problème : champs de vitesses, de températures et flux de chaleur échangés.
II-2. Mise en équations
Les équations du mouvement sont celles de l’écoulement plan
- Conservation de la masse

- Conservation de la quantité de mouvement dans le cas bidimensionnelle

( )
 4

( )

- Conservation d’énergie sans dissipation

T T  2T  2T
u v  a(  )
x y x 2 y 2

*Conditions aux limites hydrodynamiques :

- Section d’entrée généralement u,v fixées
- Surface solide u=v=0,
- Section de sortie condition d’écoulement établi

*Conditions aux limites thermiques

- Section d’entrée généralement T fixée
- Section de sortie condition d’écoulement établi
- Les parois solides T ou sont fixées

III-FORMULATION GENERALE DES EQUATIONS

DIFFERENTIELLLES DE CONSERVATION
III-1. Remarque
Plusieurs situations pratiques impliquent en même temps des écoulements de fluide, de
transfert de chaleur, la turbulence…etc. Ces phénomènes sont régis par les principes de conservation
de la masse, quantité de mouvement et d’énergie. Ces principes sont exprimés sous forme
d’équations différentielles, tous possédant une forme commune. L’intérêt de cette forme est de
rendre plus commode la présentation des méthodes numériques… . La présentation mathématique
basée sur une seule équation différentielle générale facilite la construction de la méthode numérique
et rend possible aussi le développement d’un code de calcul général.

III-2. Forme générale d’une équation de conservation

Une équation différentielle de conservation peut être présentée en notation tensorielle cartésienne,
sous la forme suivante :

( ) ( ) ( ) ( ) (1)

Dans cette équation, φ est la variable dépendante généralisée par unité de volume, est coefficient
de diffusion généralisé et ( ) le terme source. Les quantités et ( ) sont
spécifiques à chaque variable φ. L’équation différentielle totale est constituée de quatre termes,
respectivement le terme d’accumulation, le terme de transport convectif, le terme de transport
diffusionnelle et le terme source.
Variable généralisée :
φ représente la variable généralisée, elle peut représenter les différentes quantités physiques. Si
l’équation générale est supposée représenter la conservation de la quantité de mouvement, alors φ
 5

est interprétée comme la composante , de même φ peut représenter la température, la masse ou

autres variables. Chaque φ est accompagnée de l’expression de , de ( ) et de condition
aux limites spécifiques.
Coefficient de diffusion :
est le coefficient de diffusion dans l’équation différentielle générale pour φ, il exprime la
propriété de transport comme la viscosité ou la diffusivité, conjointement avec le gradient de la
variable appropriée φ.
Terme source :

Le terme source ( ) exprime le processus de production (destruction) interne comme la

génération de chaleur dans le fluide. Il contient deux termes de même nature mais d’origines
différentes, le terme ( ) est du aux interactions des couches fluides entre elles, alors que le terme
( ) est dû aux interactions des couches fluides avec des obstacles solides présents dans
l’écoulement.
Le terme source dans l’équation de quantité de mouvement inclut le gradient de pression, les forces
volumiques et la contribution des termes de viscosité qu’on ne peut pas introduire dans le terme de
diffusion. De cette manière, ( ) sert à englober toutes les quantités qu’on ne peut pas
introduire dans les autres termes.
TABLEAU RECAPITULATIF
DES EQUATIONS DE CONSERVATION
Pour un écoulement laminaire de fluide visqueux sans présence d’obstacles ni génération de chaleur

conservation φ
Masse 1 0 0 0

- 0
Vitesse suivant 1
La direction + [ ( )]

Energie dans le T 0 0
fluide

Ce modèle est applicable au calcul des écoulements laminaires et du transfert de chaleur convectif.
Le modèle est constitué de 5 équations aux dérivées partielles non-linéaires et couplées. Leur
résolution nécessite une méthode numérique.
 6

IV- DISCRETISATION DE L’EQUATION GENERALE DE

CONSERVATION

Dans le but d’alléger la présentation de la résolution du système d’équations de conservation nous

nous proposons de résoudre l’équation générale de conservation sur un domaine physique. Pour cela
nous présentons la méthode de discrétisation, la discrétisation et la linéarisation du terme source et
nous proposons une méthode de résolution mathématique.

IV-1. DOMAINE DE DISCRETISATION

La méthode numérique qu’on se propose de décrire ici est basée sur la formulation volume de
contrôle. Le domaine physique est divisé en un nombre fini d’éléments de volume, où les variables
telles que pression, température, vitesse sont calculées. La subdivision du domaine en volumes de
contrôle et l’association d’un réseau (maillage) de points sont schématisés sur la figure 1.
Les lignes en gras indiquent les faces de volumes de contrôle, deux cas sont possibles
 Maillage à volume complet à la frontière, les (o) sont les centres des volumes de contrôle
 Maillage à demi volume à la frontière, les ( ) sont au centre des volumes de contrôle.

Figure 1 : Maillage et volume de contrôle

La plupart des propriétés physiques du fluide sont constantes sur le volume de contrôle. Cependant il
est souhaitable de placer les faces des volumes de contrôle là où il y a discontinuité des propriétés du
fluide ou du terme source et de préférence là où les conditions aux limites sont définies. Après que la
 7

position des faces des volumes de contrôles ait été décidée les nœuds sont placés aux centres des
volumes de contrôle.
Considérons le volume de contrôle type sur la figure ci-dessus placé autour du point P. Il
communique avec les quatre volumes voisins à travers les quatre faces. Les centres des volumes de
contrôle voisins sont notés E, W, N et S pour indiquer Est, Ouest, Nord et Sud. Il n’est souvent
désirable d’utiliser un maillage non uniforme, car pour une utilisation optimale d’un même nombre
de points on peut mettre plus de points là où la variable φ varie rapidement et moins de points où φ
varie peu.

IV-2. DISCRETISATION
La discrétisation est réalisée en intégrant l’équation différentielle sur un volume élémentaire appelé
volume de contrôle (volume fini).

2.1- Discrétisation de l’équation de conservation générale.

La base de la méthode numérique est la conversion de l’équation différentielle (1), en
équation algébrique reliant la valeur de φ au point P aux valeurs de la variable aux points
voisins. Ceci est réalisé par l’intégration de l’équation (1) sur le volume de contrôle centré au
point P et l’approximation des différents termes dans l’intégration par des expressions fonctions
de la variable φ.

Pour l’intégration de cette équation sur le volume de contrôle considéré, il convient

d’introduire la notion de flux généralisé (convection et diffusion).
On appel flux généralisé suivant la direction i la quantité Ji définie par :

(2)

L’équation différentielle (1) peut alors s’écrire :

( ) ( ) ( ) (3)

Le volume de contrôle sur lequel on va intégrer notre équation est le volume de contrôle schématisé
sur la figure 2, centré au point P, il est limité par les quatre faces désignées par e, w, n, s et son
volume est .
Par intégration de l’équation (3) sur le volume de contrôle considéré, on a obtenu l’équation
algébrique (4) en fonction des termes sur les faces. La valeur de ( ) est constante sur les
faces Est et Ouest du volume de contrôle, la valeur de ( )est constante sur les faces Nord et
Sud, et la valeur de (̅̅̅̅ ̅̅̅̅ ) est constant sur l’ensemble du volume de contrôle. Ceci est justifié
par le fait que les valeurs utilisées sont des valeurs moyennes sur la surface ou dans le volume.

(̅̅̅̅ ) (4)
 8

Figure 2 : volume de contrôle type

La discrétisation de l’équation de continuité sur le même élément de volume donne :

( ) ( ) (5)
Avec à l’instant t
: Valeur de à l’instant t +

( ) : représentent les aires des différentes faces du volume de contrôle.

En multipliant (5) par et en la soustrayant de (4) on obtient :

( ) ( ) ( ) ( ) (̅̅̅̅ ) (6)

Nous avons présenté la discrétisation dans le cas bidimensionnel mais l’idée peut facilement
s’appliquer à d’autres cas (mono-ou tridimensionnel).
 9

2.2 -Calcul des termes de type valable pour les faces (e et n)

Considérons l’élément linéaire suivant P-E :
𝛿𝑥𝑒 ⬚

P e E x

𝛿𝑥𝑒 𝛿𝑥𝑒

On pose flux généralisé réduit, on utilisant la définition de , on peut écrire

(7)

: Nombre de Peclet.
La discrétisation de sur l’élément de module P-E donne :

( ) (8)

Soit encore
( ) (9)

Où est la valeur de au point e. ce point n’est pas un nœud du maillage de , la valeur de celle-ci
est donc inconnue et il faut le déduire par interpolation entre les valeurs aux nœuds P et E.
Considérons pour simplifier une variation linéaire de entre , la valeur de au point e est
donnée par :
( ) (10)
f : facteur géométrique de pondération

L’équation (9) devient :

( ( ) ) ( ) (11)

avec { (12)
( )
Ces deux termes ont les propriétés suivantes :
- Dépendance :
- Symétrie : ( ) ( ) (c’est-à-dire si on inverse le sens des axes, A jouera le
rôle de B et vice versa)
D’après ces deux propriétés, il suffit de connaitre une des deux fonctions (A ou B) seulement et pour
des nombre de Peclet positifs ( ), l’autre s’en déduit aisément.
 10

En effet si Peclet est négatif ( ), on a d’après la première propriété :

( ) ( ) (13)
Et d’après la seconde :
( ) ( ) (14)
En remplaçant B dans (13) par son expression (14), on obtient :
( ) ( ) (15)
Et comme la valeur absolue d’une valeur négative est égale à moins la valeur elle-même :
( ) (| |) (16)
D’une manière générale on peut mettre la relation (16) sous la forme valable pour tous les nombres
de Peclet Pe (<0 ou >0) :
( ) (| |) ( ) (17)
Revenons au calcul du flux réduit relation (11) en utilisant les propriétés de dépendance et
symétrie plus haut :
( ) ( ) ( )( ) (18)
Soit
( )( ) (19)
Et le flux de quantité de mouvement se déduit facilement par :

( )( ) (20)

Et ( )( ) (21)

Avec Ae section de passage du flux

( ) ( )( ) (22)

On pose ( ) ( )

Avec

( ) (| |) ( ) (23)

Où : conductance diffusionnelle.

Et : conductance convectionnelle.
On constate que les conductances restent encore localisées au point e qui n’est pas un
nœud du maillage et donc on a besoin d’interpoler entre les points P et E pour les calculer.
2.3-Calcul des termes de type valable pour les faces (s et w)
On suivra le même résonnement que précédemment en utilisant un élément linéaire S-P, ce qui
revient à changer le sens de l’axe et à changer la fonction d’interpolation A par B.
( ) (| |) ( ) (24)
 11

Où seul le signe du nombre de Peclet à l’intérieur de la fonction Max a changé. Les conductances
globales sont données par :

( ) (| |) ( ) (25)

La loi d’interpolation est introduite au niveau de (| | ), nous avons effectué le calcul pour une
interpolation linéaire par souci de simplification, alors que cette variation peut prendre n’importe
quelle forme. Pour illustrer cela on se propose de faire un calcul exact de ce flux généralisé.

2.4 - Calcul exacte du flux généralisé (solution exacte)

Les flux généralisés ( ) par conséquent les valeurs de φ et de sa dérivée première sont
calculés en des points (faces des volumes finis) qui ne sont pas des nœuds du maillage de . Il faut
donc exprimer les valeurs de φ au droit des faces en fonctions des valeurs de φ aux centres des
volumes de contrôle. Pour ce faire, on considère le cas simple suivant :
On se propose de résoudre un problème monodimensionnel de convection diffusion en régime
permanent, à coefficient de diffusion constant et sans terme source. L’équation (3) s’écrit dans ce
cas :

( ) (26)

𝛿𝑥𝑒 ⬚

P e E x

𝛿𝑥𝑒 𝛿𝑥𝑒

L’intégration de l’équation (26) entre les points P et E donne :

(27)

Avec suivant la direction i (28)

L’équation (26) s’écrit encore :

(29)

La résolution analytique de l’équation (29) entre les points P et E figure ci-dessus (par intégration
deux fois entre le point P et E), permet de donner la solution exacte de φ(x) entre les points P et E.
Avec les conditions aux limites suivantes :
x=0 (30a)
x=L (30b)
 12

La solution est donnée par :

( )
( )
(31)

En tenant compte de la solution (31) qu’on substitue dans de l’équation(28), on obtient pour i=e, le
flux généralisé J entre les points P et E celui-ci s’écrit :

( ( )
) (32)

Où ; ; (33)

Ce qui donne quantité adimensionnelle appelée nombre de Peclet (traduisant l’importance

de la convection par rapport à la diffusion) de maille.
Les distances et sont indiquées sur la figure 2, le coefficient de diffusion Γ est considéré
uniforme sur chaque volume de contrôle.

( ) ( )
On définie * + (34)

La raison de définir de la sorte est plutôt physique que numérique.

Finalement par substitution de la relation (32) et une expression similaire pour J w dans l’équation (27)
on obtient :

( ( )
) ( ( )
) (35)

Qui peut être écrite sous forme standard :

(36)

Avec ( )

( )
( )
( ) (37)
De façon identique, on peut déterminer les différents flux totaux entre le nœud central P et les
autres nœuds voisins.

2.5- Linéarisation du terme source

Ce terme peut être constant ce qui ne pose aucun problème particulier. Mais souvent il dépend
explicitement de la variable . Cela peut causer des problèmes au niveau du traitement
numérique (instabilités, divergence de la méthode numérique). Pour éviter ce genre de problème il
faut linéariser ce terme en le décomposant en deux parties :
 13

̅̅̅̅ (38)

Où : partie constante du terme source

: Coefficient de
L’opération de linéarisation du terme source est souvent cruciale, car la solution de beaucoup de
problèmes complexes dépend de cette décomposition. Il faut veiller à ce que le coefficient soit
négatif, un coefficient positif peut faire diverger le code de calcul.

La décomposition du terme source peut être quelconque il suffit que le coefficient soit négatif,
une façon de le linéariser est de prendre comme étant la dérivée du terme source au point
( : valeur de φ calculée dans l’itération précédente).

2.6- Discrétisation de l’équation à deux dimensions

A deux dimension l’équation générale de conservation en introduisant la notion de flux

généralisé suivant la direction i la quantité Ji définie par :

L’équation différentielle (1) peut alors s’écrire :

( ) ( ) ( ) ( ) (39)

Où et sont les flux généralisé (convection +diffusion) tel que :

u et v sont les composantes de vitesse suivant x et y.

L’intégration de l’équation (39) sur l’élément de volume schématisé figure 3, on obtient :

( )
( ) (40)

On peut utiliser la solution exacte ce qui donne :

( )
[ ( )] [ ( )]
( ) ( )

+* ( ( )
)+ * ( ( )
)+ = ( )

Sont les valeurs de au point P à l’instant donnée et l’instant précédent.

De la même manière on intègre l’équation de continuité sur l’élément de volume, on obtient :
 14

(41)
sont les flux massiques traversant les interface de volume de contrôle.
Avec

Jn

JW Je
∆y

Js

∆x

Figure 3 : volume de contrôle à deux dimensions

En multipliant équation (41) par et on la soustrayant de (40), on obtient :

( )
( ) ( ) ( ) ( ) (
) (42)
De la même manière que (22) on écrit :
( ) ( )

Ainsi que les autres termes( ), ( ) et ( ) où

(| |) [ ]
(| |) [ ] (43)
(| |) [ ]
 15

(| |) [ ]
L’expression Max [a,b] désigne le maximum de a et b . La fonction d’interpolation entre le point P et
les nœuds voisins, calculée à partir de la solution exacte d’un cas monodimensionnel, présente un
point singulier quand le terme de convection devient nul et nécessite beaucoup de temps dans le calcul
numérique.
Et l’équation finale discrétisée est donnée par :

(44)

Les lois d’interpolation

C’est au niveau de (| |) que s’introduit la loi d’interpolation. Nous avons effectué les calculs
pour une variation linéaire au § 2.2 par souci de simplification, et nous avons présenté au § 2.4 la
solution exacte, donc cette variation peut prendre n’importe quelle forme. Dans la littérature, on
trouve plusieurs, à savoir tableau ci-dessous :

Schéma (| |)
Différence centrale 1-0.5| |
Upwind (en amont) 1
Hybride Max(0, 1-0.5| |)
Puissance Max(0, ( | |) )
Exponentielle | |
(| |)

Le choix d’une loi fait intervenir plusieurs facteurs dont les plus importants sont la précision et le
temps de calcul.

L’équation (41) peut s’écrire sous la forme condensée suivante :

∑ (45)

La notation nb indique les nœuds voisins au nœud P, et la sommation est étendue à tous les nœuds
voisins au nœud P.
 16

Chapitre 2

1- DISCRETISATION DES EQUATIONS _________________________________________________ 17

1.1- Choix du maillage ___________________________________________________________ 17
1.2- Discrétisation des équations de conservation sur un volume fini ______________________ 19
2- TRAITEMENT DES CONDITIONS AUX LIMITES ________________________________________ 20
2.1- Conditions aux limites à la paroi _______________________________________________ 20
2.1.1-Vitesse : _______________________________________________________________ 20
2.1.2-Température ___________________________________________________________ 21
2.2- Sortie ____________________________________________________________________ 22
2.3- Entrée ____________________________________________________________________ 22
3- METHODE NUMERIQUE DE RESOLUTION ___________________________________________ 23
4-COUPLAGE PRESSION-VITESSE, ALGORITHME DE RESOLUTION ___________________________ 24
4.1. COUPLAGE PRESSION-VITESSE ________________________________________________ 24
4.2. ALGORITHME SIMPLER _______________________________________________________ 24
4.2.1 Algorithme SIMPLE ______________________________________________________ 24
4.2.2. Algorithme SIMPLER _____________________________________________________ 26
4.3. CONDITIONS AUX LIMITES POUR LA CORRECTION DE PRESSION ET POUR LA PRESSION ____ 27
4.4 SOUS-RELAXATION –CONVERGENCE_____________________________________________ 29
4.5 TEST DE CONVERGENCE ______________________________________________________ 30
4.6 CONCLUSIONS ______________________________________________________________ 30
 17

Chapitre 2

DISCRETATION DES EQUATIONS DE CONSERVATGION

TRAITEMENT DES CONDITIONS AUX LIMITES
METHODE NUMERIQUE DE RESOLUTION DES EQUATIONS ALGEBRIQUES
COUPLAGE PRESSION-VITESSE, ALGORITHME DE RESOLUTION

Les équations de bases sont discrétisées spatialement en utilisant l’approche des volumes de
contrôle (ou volume finis). Cette technique, à base physique, s’appuie sur le principe de
conservation ; elle consiste à faire des bilans, dans chaque volume fini, des différentes quantités
physiques telles que : la masse, la quantité de mouvement, l’énergie …., etc. Mathématiquement,
cette technique se traduit par l’intégration des équations, écrites sous forme conservative, sur
chaque volume.

1- DISCRETISATION DES EQUATIONS

1.1- Choix du maillage
De nombreux schémas aux volumes finis ont été utilisés auparavant, la principale différence entre
ceux-ci est les stockages des variables dépendantes (composantes de la vitesse, pression,
température) ainsi que les différents volumes pour le calcul de ces variables.
Le choix du maillage, doit être basé sur des considérations physiques, à savoir : les propriétés de
conservation des variables physiques doivent être maintenues quand on passe d’un volume de
contrôles au volume voisin, la discontinuité des propriétés physiques doit être respectée, et
l’implémentation des conditions aux limites doit être aisée. Dans notre cas, deux problèmes doivent
être considérés :
- la présence de la pression sous forme de gradient dans le terme source.
- découplage entre les nœuds pairs et impairs dans la discrétisation de l’équation de
continuité.
1°/ la présence de la pression sous forme de gradient dans le terme source des équations
de conservation de la quantité de mouvement: L’intégration du gradient de pression sur le volume
de contrôle suivant la direction x, donne :
∫ ( ) (1)

: est la différence de pression exercée sur le volume de contrôle.

Or sont inconnues et doivent être calculées à partir des valeurs nodales par interpolation.

(2)

L’équation de conservation de la quantité de mouvement contient donc la différence de pression entre

deux points alternés et non pas entre deux points adjacents. Un champ de pression uniforme ne peut
donc être différencié d’un champ de pression en « zig-zag » oscillant entre un niveau constant pour
 18

tous les nœuds pairs et un autre niveau pour tous les nœuds impairs. Cela permettrait au programme
de calcul de converger vers une solution numériquement exacte mais qui ne présente aucune réalité
physique.
2°/ La discrétisation de l’équation de continuité dans le cas monodimensionnel suivant la
direction x permet d’écrire :
(3)
Sont des valeurs inter-nodales et doivent être déterminées à partir des vitesses nodales
par interpolation (linéaire pour simplifier). On aboutit à :
( ) ( ) ( )
(4)
Une fois de plus, il y a découplage entre les nœuds pairs et impairs, ce qui peut provoquer des
oscillations non physique du champ de vitesse. Cela remet en cause la propriété de conservation de
la masse en passant d’un volume de contrôle au volume voisin, et en conséquence celle des autres
variables.
Les deux problèmes évoqués ci-dessus, liés directement à la discrétisation, ont conduit Harlow et
Welch (1965) à repenser autrement le problème de la localisation des variables et utiliser des
maillages décalés pour les vitesses. On considère donc une grille pour la pression (et toutes les
variables scalaires), et on met les nœuds pour les composantes de la vitesse sur les faces des
volumes de contrôle figure 1. Cela permet d’éviter les deux problèmes cités ci-dessus provocant
l’oscillation des champs de pression et de vitesse.

Figure 1 : Localisation des différentes variables sur les volumes de contrôle

 19

1.2- Discrétisation des équations de conservation sur un volume fini

L’intégration des équations de conservation sur un volume fini (volume de contrôle) dans un système
physique bidimensionnel s’obtient facilement en appliquant la méthode de discrétisation de l’équation
générale présentée dans le chapitre 1. La pression P, la température T, sont définies au centre du
volume de contrôle centré au point P( ) figure 2. Les composantes de la vitesse (U,V) sont
définies respectivement sur les faces Est et Ouest localisées à (i± ) et Nord et Sud localisées à
(i,j 1/2). La figure 1 illustre les différents éléments dans le cas d’un maillage cartésien
bidimensionnel.
Les coefficients ainsi que le terme source sont spécifiques à chaque variable. Nous développons dans
ce qui suit la discrétisation et l’introduction des conditions aux limites et aux interfaces relatives à
chaque variable.

N
()
𝛿𝑦(𝑗)

W P(xi ,yj)
E

dy(j)

𝛿𝑥(𝑖)

y dx(i)

S
x

Figure 2: Position et dimensions du volume de contrôle

On a
() , ( )

() ( )

() , ( ) (5)

() ( )
 20

Il y a deux possibilités de mailler un domaine physique figure (1 chap. 1).

1- Volume complet à la frontière
2- Demi-volume à la frontière
Le choix de l’une ou de l’autre façon de mailler un domaine ne peut être déterminé que par le bon
conditionnement des conditions aux limites. En effet les conditions aux limites sont souvent
génératrices d’erreurs dans le calcul numérique et peuvent causer des d’instabilités. De plus la
distribution du champ de vitesse a une influence considérable sur les autres variables qui en
dépendent, à savoir la température et le transfert de chaleur. Un maillage à volume complet à la
frontière pour les variables scalaires ( P, T, …) s’adapte bien à beaucoup de problèmes et intègre bien
les conditions aux limites dans le code numérique, là où elles sont définies.

2- TRAITEMENT DES CONDITIONS AUX LIMITES

La qualité de toute simulation numérique est fortement liée à la façon avec laquelle ces conditions
aux frontières sont calculées ou imposées.
L’équation de conservation discrétisée (6 chap. 1) fait intervenir la valeur du flux généralisé sur la
face du volume de contrôle confondue avec la paroi. Ayant choisi le maillage à volume complet à la
frontière pour les variables scalaires telles que P, et T comment peut-on implémenter les conditions
aux limites pour chaque variable dans l’équation discrétisée correspondante ? On va essayer de
répondre à cette question dans ce qui suit.

2.1- Conditions aux limites à la paroi

Comme nous l’avons introduit dans le chapitre 1,
*Conditions aux limites hydrodynamiques :
- Section d’entrée généralement u,v fixées
- Surface solide u=v=0,
- Section de sortie condition d’écoulement établi

*Conditions aux limites thermiques

- Section d’entrée généralement T fixée
- Section de sortie condition d’écoulement établi
- Les parois solides T ou sont fixées

Près d’une paroi solide la sous-couche visqueuse est très mince et caractérisée par un gradient de
vitesse tangentielle considérable. En général elle ne peut être décrite dans le maillage utilisé. Un
maillage fin est nécessaire pour pouvoir tenir compte des effets de cette sous-couche visqueuse. Cela
peut devenir coûteux en temps de calcule et en mémoire.

2.1.1-Vitesse :
a- cas d’une valeur imposée à la frontière
La discrétisation de l’équation de conservation de la quantité de mouvement suivant la direction x,
sur un volume de contrôle dont la face Nord est une paroi, s’écrit :
( ) ( ) ( ) ( ) (6)
 21

Où on a remplacé le terme ( ) par ( ) sur la frontière Nord c’est le seule terme

à calculer, le calcul des autres termes, reste inchangé, comme pour un volume central.

( ) (7)

( )

Par différence finie :

( )

Donc :

( )( )

En multipliant par et en utilisant les conductances de diffusion D et de convection F, on a :

( ) ( )( ) ( )( )

( ) ( )

Comme pour la conductance globale intérieure, on met la conductance globale à la frontière sous la
forme :

[ ] (8)

Ce qui revient à prendre : (| |)

Cette conclusion est tout à fait logique puisque la valeur de sur la frontière qui coïncide avec la
paroi est donnée et on n’a pas besoin d’interpoler pour l’obtenir.

b- cas d’un gradient imposé à la frontière

La discrétisation de l’équation permet d’écrire :

=0 ou (9)

On effectue alors le calcul normalement et on adapte à chaque fois les valeurs à la frontière en
fonction des valeurs calculées. Ceci revient à imposer une conductance nulle à la frontière .

2.1.2-Température
a- Cas du gradient imposé à la frontière
L’équation à la frontière reste la même, il suffit de recalculer le terme

( ) ( ) (10)

L’expression – ( ) exprime le flux de chaleur à la paroi.

 22

On constate qu’en plus du flux imposé à la frontière, on a besoin de la valeur de sur cette même
frontière. La manière la plus générale pour résoudre le problème est de développer une équation pour
chaque à la frontière, tenant compte du flux imposé, que l’on résout avec le système global pour
obtenir .

b- Flux imposé constant

Ce cas présente seulement pour l’équation d’énergie et dans le cas où la frontière est une paroi
solide. On a alors la vitesse à la frontière est nulle et l’équation à la frontière se réduit à :
cte (11)
Qui s’introduit dans le terme source de l’équation discrétisée.
Pour tous autres cas, avec un peu de réflexion on établira l’expression aux limites.

c- Paroi adiabatique
Si la paroi est adiabatique, il suffit d’assurer un coefficient nul à la frontière. Si par exemple c’est la
face Est du volume de contrôle qui est confondue avec la paroi on pose alors :
. (12)

2.2- Sortie
Le fait d’imposer un gradient nul à la sortie, se traduit numériquement par un coefficient nul sur la
face aval du volume de contrôle à la sortie.
(13)

Cela peut être justifié par le fait que, généralement à la sortie le phénomène de convection devient
prépondérant par rapport à la diffusion, par conséquent le nombre de Peclet devient important. Cela
donne un coefficient nul à l’Est si cette face est confondue avec la frontière à la sortie.

2.3- Entrée
Les valeurs des différentes variables sont imposées, leur introduction dans les calculs ne pose aucun
problème.
 23

3- METHODE NUMERIQUE DE RESOLUTION

L’équation discrétisée (44 chap.1) exprime la valeur de φ au point P en fonction des valeurs de φ aux
nœuds voisins. Elle s’écrit autant de fois qu’il y a de nœuds où la variable φ est considéré comme
inconnue. Les coefficients etc, dépendent généralement de la variable φ.
En effet sur l’ensemble du domaine physique nous obtenons un système d’équations algébriques
non-linéaires de la variable φ. La solution de ce système d’équations non-linéaires, peut être obtenue
par l’intermédiaire de solution d’équations algébriques linéaires, en utilisant une méthode de
résolution itérative. Cela consiste à calculer les coefficients en fonction des précédentes estimations
des valeurs de φ, l’équation discrétisée (44 chap.1) est ainsi traitée comme une équation linéaire,
permet de calculer une nouvelle estimation de la variable recherchée φ. Les coefficients pour
l’itération suivante sont calculés en fonction des nouvelles estimations.
En pratique une méthode de résolution doit être peu coûteuse à la fois en place mémoire, et en
temps de calcul (avoir une convergence rapide). Cela n’est possible que par la combinaison de
plusieurs techniques complémentaires qui sont par exemple :
Méthode de balayage rangée par rangée
Algorithme de Thomas (méthode de résolution directe)
Procédure de correction par bloc
La combinaison de ces trois techniques constitue l’algorithme de résolution.
Le principe de base est le suivant :
L’idée de base de la méthode de balayage rangée par rangée, valable dans le cas d’un domaine de
dimension supérieure ou égale à deux, est transformer l’équation (44 chap.1), en une pseudo-
équation monodimensionnelle par la substitution des valeurs de φ dans les autres directions par des
approximations. L’équation obtenue pour une ligne est alors résolue par l’Algorithme de Thomas. On
parcourt le système ligne par ligne jusqu’à ce que toutes les lignes soient traitées de cette manière
dans une direction. On parcourt ensuite de la même façon l’autre direction. Pour illustrer cela,
considérons un système bidimensionnel, l’équation (44 chap.1) s’écrit dans ce cas :
(14)

Si la méthode de balayage rangée par rangée est appliquée suivant la direction x, on doit introduire des
approximations des valeurs de ϕ aux nœuds voisins suivant la direction y, soit le résultat est :
( ) (15)

Le signe (*) indique que la valeur est approximative, et le terme entre parenthèse dans l’équation (15)
est considéré comme connu. En regroupant les termes inconnus dans le même membre, cette équation
peut s’écrire :
(16)
Equation dans laquelle i indique la position du nœud suivant la direction x, et les coefficients
sont donné par :

(17)
 24

On obtient donc pour chaque ligne un système d’équations de structure tridiagonale qui peut être
résolu par l’Algorithme de thomas.
La méthode de résolution rangée par rangée consiste donc à écrire une équation semblable à
l’équation (15) pour chaque nœud du maillage (volume fini) suivant la direction x, on obtient un
système d’équations auquel on applique une méthode de résolution directe (par exemple
l’Algorithme de Thomas), après on explore la ligne suivante. Le domaine est parcouru de la même
façon suivant les autres directions suivant que le domaine est bi ou tridimensionnel. Le processus est
répété jusqu’à la convergence de la solution. On trouve dans l’annexe A bref rappel de l’algorithme
de Thomas et de la procédure de correction bloc par bloc.

4-COUPLAGE PRESSION-VITESSE, ALGORITHME DE

RESOLUTION

Dans cette partie on se propose de décrire succinctement la procédure que nous avons adoptée pour
la résolution du système des équations de conservation. Après avoir décrit dans un chapitre
précédent l’algorithme de résolution d’une équation de conservation, nous présenterons
l’algorithme que nous avons choisi d’utiliser pour tenir compte du couplage les différentes variables.

4.1. COUPLAGE PRESSION-VITESSE

Par couplage pression-vitesse, on entend la liaison qui existe entre le champ de pression et le champ
de vitesse. En effet l’équation de conservation de la quantité de mouvement, contient le gradient de
pression.
Cependant cette équation ne peut être résolue que si le terme de pression est spécifié, c’est-à-dire
on fixe un champ de pression. Alors pour différents champs de pression, qu’on aura fixée, on aura
différents champs de vitesse. Si le champ de pression est exact, le champ de vitesse résultant de la
résolution de l’équation de conservation de la quantité de mouvement doit satisfaire également
l’équation de continuité. La contrainte qui permet la détermination du champ de pression est donc
l’équation de continuité.
La manière dont il faut résoudre le couplage vitesse-pression a fait l’objet de nombreuses études. Il
en a résulté un certain nombre d’algorithmes (SIMPLE, SIMPLER, SIMPLEC, PISO …), qui ne différent
que par la technique utilisée pour le traitement de ce couplage.

4.2. ALGORITHME SIMPLER

L’algorithme SIMPLER, est une amélioration de l’algorithme SIMPLE proposé par PATANKAR. Nous
commençons donc par expliciter ce dernier. Nous introduisons les différentes améliorations
apportées par les différents auteurs pour construire l’algorithme SIMPLER.

4.2.1 Algorithme SIMPLE

Les équations algébriques de conservation de la quantité de mouvement permettent de déterminer
le champ de vitesse. Ces équations sont fonctions du gradient de la pression et d’un terme source qui
dépend de la vitesse.
∑ ( ) (18)

∑ ( ) (19)
 25

Où les termes désignent respectivement le terme source correspondant à chaque

composante de la vitesse, excepté le gradient de pression.
Il est possible de résoudre ces équations en se donnant un champ de pression P*. Il en résulte un
champ de vitesse approximatif U*, donné par :
∑ ( ) (20)
∑ ( ) (21)
Le champ de pression étant approximatif, les solutions obtenues après une itération ne satisfont pas
à l’équation de conservation de la masse. Il faut corriger les vitesses et les pressions pour que les
trois champs satisfassent simultanément les trois équations de conservations.
On écrit alors :
(22)

En substituant les quantités définies ci-dessus dans l’équation (18,19) et en soustrayant l’équation
(20,21) des résultats, on obtient :
∑ ( ) (23)
∑ ( ) (24)
La caractéristique principale de l’algorithme SIMPLE est de considérer comme négligeables les
termes ∑ ∑ qui constituent la contribution des corrections de vitesse des nœuds
voisins.
Cet approximation est justifiée par le fait que, quand la solution aura convergé vers la solution finale
les corrections U’ et V’ deviennent nulles et par conséquent les contributions des nœuds voisins sont
nulles. Ce qui permet d’écrire les corrections de vitesse comme suit :

( ) (25)

( ) (26)

On pose et

On obtient donc les vitesses corrigés par :

( ) (27)
( ) (28)
Equation de correction de pression
En introduisant ces expressions dans l’équation de continuité écrite pour le volume associé au nœud
P, on obtient l’équation suivante pour les corrections de pression.
(29)
Avec (30)
(31)
(32)
(33)
 26

(34)
(35)
Nous avons ainsi obtenu l’équation de correction de la pression. Elle présente la même forme que
l’équation générale discrétisée. Elle peut donc facilement être intégrée dans le code de calcul
numérique. Sa résolution permet de calculer le champ de correction de pression P’, ainsi que ceux
des corrections des vitesses par application des relations (25,26).
Remarquons que n’est rien d’autre que le bilan de masse correspondant au champ de vitesse
approximatif. Lorsque le terme source devient nul il n’y aura plus de génération du champ de
correction de pression, U* et V* seront les solutions exactes U et V.
Le champ de correction de pression obtenu à l’aide de l’équation (29), est généralement grossier : il a
le rôle de corriger les vitesses, mai entraîne aussi des corrections exagérées pour le champ de
pression. Cela peut faite diverger la procédure de résolution. Il est donc nécessaire d’apporter
certaines modifications à savoir :
 Introduction d’un facteur de sous-relaxation pour la pression
(36)
 Introduction d’un facteur de sous-relaxation pour les vitesses
Cette sous relaxation pour les vitesses se fait à la fin de la résolution de toutes les équations du
système, à l’aide de la valeur de l’itération en cours de la variable à relaxer et de sa valeur lors de
l’itération précédente :
( )
= ( ) (39)
* : indique la valeur de φ à l’itération précédente.
Les facteurs de sous-relaxation ont une grande influence sur la convergence et sont parfois difficiles à
ajuster. Ils varient d’un problème à un autre. Pour remédier à ce problème une première
amélioration introduire par Patankar, consiste à trouver une équation pour la pression. Celle-ci- sera
calculée directement ce qui permet d’éliminer le coefficient de sous relaxation . C’est cette
amélioration qui est à l’origine de l’algorithme SIMPLER.

4.2.2. Algorithme SIMPLER

Considérons l’équation de conservation de la quantité de mouvement que nous pouvons écrire sous
∑ ( )
la forme suivante : (40)

∑ ( )
(41)

On définit les composantes de vitesse ̂ ̂ , appelées pseudo-vitesses, comme étant les vitesses
que l’on aurait aux faces du volume fini, en l’absence de gradient de pression.
Elles sont calculées juste en fonction des vitesses aux nœuds voisins.
∑
̂ (42)

∑
̂ (43)
 27

Equation de la pression
En substituant les vitesses données par les relations (40,41) dans l’équation de continuité discrétisée,
on obtient une équation analogue à celle obtenue pour la correction de pression P’(29), mais cette
fois-ci en P.
(44)
Avec les coefficients donnés par les relations (30-34. Quant au terme source
il est défini par :
̂ ̂ ̂ ̂ (45)
La seule différence entre l’équation de correction de pression et l’équation de pression est le terme
source b, calculé à partir de vitesses U* dans un cas et des vitesses ̂ dans l’autre.
L’équation obtenue nous permet de calculer le champ de pression exact qui satisfait l’équation de
conservation de la masse correspondant aux vitesses estimées. Par conséquent on a plus besoin du
facteur de sous-relaxation .
Nous avons ainsi décrit la modification introduire sur l’algorithme ‘Simple ». Il en résulte un nouvel
algorithme, « SIMPLER », donc voici les étapes de calcul.
Etapes de calcul « Algorithme SIMPLER »
1- Estimation des vitesses (valeurs initiale,)
2- Calcul des coefficients des équations de quantité de mouvement et calcul des
pseudo-vitesses ̂ ̂ par les équations (42 ,43)
3- Evaluation du terme source par la relation (45) et résolution de l’équation
(44) ce qui donne le champ de pression P= P*
4- Résoudre des équations de la quantité de mouvement, on obtient U* et V*
5- Evaluation du terme source à partir de la relation (35) et résolution de
l’équation (29) en P’
6- Correction du champ de vitesse (U,V) relations (22) , en calculant les corrections
par les relations (27,28)
7- Résoudre des autres équations, énergie … si c’est couplé
8- Test de convergence :
* oui passer à l’étape suivant
* Non revenir à l’étape n°2 avec les champs de vitesse (U,V) comme
nouvelles estimations
9- Résolution des équations relatives aux autres variables non couplées ni avec la
pression ni avec les vitesses.

4.3. CONDITIONS AUX LIMITES POUR LA CORRECTION DE PRESSION ET POUR

LA PRESSION

Le maillage que nous avons adopté (maillage à volume complet à la frontière), permet d’avoir une
face du volume de contrôle pour la pression qui coïncide avec la frontière, (figure 1).
La vitesse sur la paroi est, spécifié d’après les conditions aux limites, c’est-à-dire que est fixe. Par
conséquent la correction de vitesse est nulle D’après la relation (25) on a :
( ) (46)
On en déduit que :
(47)
 28

Cela se traduit donc au niveau de la procédure numérique par l’introduction d’un coefficient nul à la
frontière Est ( ) dans l’équation de correction de pression.
Pour la pression, comme la vitesse est connue sur la frontière, donc on n’a pas besoin de valeur de
pression. Mais la condition aux limites reste posée, car la pression intervient sous forme de gradient
et on a besoin de sa valeur sur la frontière. Pour détourner cette situation on adopte la même
condition que pour la correction de pression c'est-à-dire .

N
Volume fini V
𝛿𝑦𝑛 N’
n
𝛿𝑦(𝑗)
Volume fini U

W P(xi ,yj) E
w e
dy(j)
𝛿𝑥𝑒

𝛿𝑥(𝑖)
s
y S’
dx(i)

S
x

Figure 3 : Volume de contrôle frontière

Cette façon de faire justifie par le fait que près de la paroi, l’écoulement est dominé par les forts
gradients de vitesse tangentielle. L’erreur est négligeable en adoptant l’approximation suivante :
(48)
La pression n’intervient dans les calculs que sous forme de gradients. La pression n’est définie qu’à
une constante près. Il faut donc fixer le niveau de pression en un point au moins de l’écoulement.
Cela s’intègre facilement dans le code numérique, en posant :
Soit I le point où la pression et imposée.
(49)
(50)
(51)
A : constante quelconque non nulle
: Niveau de pression imposé
 29

Il est souhaitable, pour de raisons numériques, de prendre une valeur de A proche des coefficients
aux nœuds voisins qui entourent le point I.

4.4 SOUS-RELAXATION –CONVERGENCE

La résolution des équations algébriques de conservation est réalisée à l’aide d’une méthode itérative
décrite par l’algorithme SIMPLER pour traiter le couplage et la non-linéarité de ces équations. Le
processus peut converger vers la solution, comme il peut osciller autour sans jamais l’atteindre ou
diverger complètement. L’introduction d’un coefficient de sous-relaxation ( ) permet de
contrôler la variation de la variable φ d’une itération à une autre (non-linaire), et par suite assurer la
convergence du schéma numérique.
Plusieurs techniques sont utilisées dans la littérature, nous présentons celle que nous avons utilisée
due à Patankar. Elle consiste à opérer la sous-relaxation pour les vitesses (ou une autre variable) au
sein de la boucle non-linéaire au moment du calcul des coefficients. Cela a l’avantage d’accélérer la
convergence et de diminuer le temps de calcul.
Soit l’équation de conservation pour la variables φ
∑ (52)

La variable φ est donnée par la relation suivante

∑
(53)

Si on ajoute et on retranche au deuxième membre de l’équation ci-dessus, la valeur de φ à l’itération

précédente, notée φ*, on obtient
∑
( ) (54)

Les quantités sont calculée en fonction de *.

L’expression entre parenthèse de l’équation (54) représente la variation de φ par rapport à l’itération
précédente. Cette variation peut être modifiée par l’introduction du facteur de sous-relaxation ,
ce qui donne :
∑
( ) (55)

Ce qui s’écrit encore

∑ (56)

Cette sous-relaxation s’introduit par une modification du coefficient et du terme source b qui
deviennent
∑( )
(57)

(58)

Quand la solution aura convergé, la valeur de φ devient égale à la valeur de φ*, et la valeur de φ
satisfait à l’équation (52).
 30

Le choix du coefficient de sous-relaxation dépend de la nature du problème et de la variable φ. Il

peut être variable d’une itération à une autre pour un même problème, comme il peut avoir des
valeurs différentes d’une maille à l’autre.

4.5 TEST DE CONVERGENCE

La résolution du système d’équations différentielle se fait d’une manière itérative, il faut donc avoir
un critère pour juger de la convergence de la solution, c'est-à-dire un critère pour arrêter les calculs.
Le premier test qui vient à l’esprit est de mesurer la variation des différentes variables entre deux
itérations successive es en chaque nœud. C'est-à-dire que l’erreur maximale détectée en un nœud du
maillage ne soit pas supérieure à une valeur E .
Soit
( ) (59)
E . : Nombre très petit choisi suivant la précision désirée. Tant que la relation (59) n’est pas
satisfaite, le calcul se poursuit. Mais cela est insufflant, car une faible vitesse de convergence peut
donner l’illusion que la variable considérée converge vers la solution exacte.
Un test plus réalise physiquement pour contrôler la convergence est de calculer à la fin de chaque
itération, les résidus en chaque nœud du maillage et pour chaque variable φ. est calculé
d’après par la relation :
(∑ ) (60)

Le test d’arrêt et alors que les associés à chaque variable φ tendent vers zéro :

| | 61)

Nous réalisons ces tests de convergence dans notre code numérique comme suit :
- La résolution du système d’équations correspondant à une variable φ est considérée comme
convergente si la variable φ entre deux itérations successives satisfait à la relation (59), en
chaque nœud du maillage.
- La solution convergente est obtenue par l’algorithme SIMPLER si l’ensemble des variables
satisfont à la relation (61) séparément et en chaque nœud du maillage.

4.6 CONCLUSIONS
Dans ce chapitre nous avons décrit l’algorithme SIMPLER, pour la résolution du système des
équations de conservation. Cet algorithme traite le couplage des variables, permettant ainsi l’écriture
d’un système d’équations algébriques indépendantes. La solution de chaque équation peut être
obtenue alors par simple application de l’algorithme présenté.
 31

Chapitre 3 : Applications

Résolution d’un problème de convection libre en 2.D...................................................................... 32

I- Formulation du problème : ........................................................................................................... 32
II- Solution numérique .................................................................................................................. 34
1. Discrétisation des équations gouvernantes par la méthode du volume de contrôle (ou volume
finis) ................................................................................................................................................... 34
b- Equation de mouvement suivant x ......................................................................................... 36
b- Equation de mouvement suivant y ........................................................................................... 36
c- Discrétisation de l’équation de l’énergie................................................................................... 37
d- Discrétisation des conditions aux frontières : ........................................................................... 37
2-Résolution de équations discrétisées............................................................................................. 37
Algorithme de Résolution (SIMPLER →SIMPLEC) .......................................................................... 37
 32

Chapitre 3

Résolution d’un problème de convection libre en 2.D

I- Formulation du problème :

Il s’agit d’étudier un problème de convection naturelle bidimensionnelle dans une cavité

rectangulaire de hauteur H et de longueur L, inclinée d’un angle ϕ, par rapport à l’horizontale,
deux parois sont adiabatiques les autres parois sont symétriquement chauffées et refroidie à
l’aide d’un flux de chaleur constant qui les maintient à des températures constantes et
différentes.

Adiabatique
Tf=cte

Tc=cte
Adiabatique
ϕ

Figure 2 : système étudié

En supposant que les propriétés thermiques du fluide sont constantes à l’exception de la densité
dont le terme de poussé qui est supposé obéir à l’approximation de Boussenesq.
Les équations gouvernantes pour un phénomène permanent, en utilisant les variables primitives
T, p, u, et v sont définies par :
- Conservation de la masse
(62)

- Conservation de la quantité de mouvement dans le cas bidimensionnelle

( ) (62)

( ) (64)

- L’équation d’énergie :
T T  2T  2T
u v  a(  ) (65)
x y x 2 y 2
 33

 Approximation e Boussenesq :
( ( ))

 Adimensionnalisation des équations :

Pour la réduction adimensionnelle de ces équations on introduit les échelles caractéristiques

tel que :

L= longueur de référence.
a
U0  : Vitesse de référence, , [ ]
[ ]
[ ][ ]
L
P0  U 02 : pression caractéristique

T  Tc  T f : Écart de température de référence

On définit, alors, les variables adimensionnelles suivantes :

u v x y
U , V , X , Y
U0 U0 L L

T  Tc

Tc  T f

Avec P  [ p  0 g ( y sin   x cos  )] / U 02 (66)

Par l’introduction de ces variables réduites dans les équations (62), (63) et (64) ci-dessus, il en résulte
les équations adimensionnelles gouvernantes suivantes :

- Conservation de la masse
(67)

- Conservation de la quantité de mouvement dans le cas bidimensionnelle

-
( ) (68)

( ) (69)

- L’équation d’énergie :
   2  2
U V   ) (70)
X Y X 2 Y 2

 Avec : nombre de Rayleigh

: Nombre de Prandtl
 34

 Conditions aux limites :

* Conditions aux limites hydrodynamiques :
U = V = 0 sur toutes les parois solides (71)
* Conditions aux limites thermiques :
X=0 et X=1 =1 flux imposé
Parois adiabatiques (72)
Adiabatique
Tf=cte

H Y,V

0
X,U

Tc=cte
Adiabatique
L ϕ

Figure 3 : système étudié

Conditions aux
Transfert de chaleur :
Le transfert de chaleur à travers le système exprimé en termes de nombre de Nusselt local est défini
par :

( ) ( )
(73)
Le modèle est constitué de 4 équations aux dérivées partielles non-linéaires et couplées. Leur
résolution nécessite une méthode numérique.

II- Solution numérique

L’intégration numérique des équations gouvernantes consiste à déterminer les valeurs des variables
primitives U,V,P et en un certains nombres de points du domaine appelés nœuds , la méthode
numérique conduit à un certain nombre d’équations algébriques où les inconnues sont U,V, P et .
Il s’agit de diviser le domaine physique en un certain nombre de volume de contrôle entourant
chacun un nœud principale P, en chaque volume, différents intégrales sont effectués.

1. Discrétisation des équations gouvernantes par la méthode du volume de

contrôle (ou volume finis)
La formulation plus générale des équations de conservations est utilisée. Le domaine de calcul est
divisé en un nombre fini de volume de contrôle, chaque volume caractérisé par son nœud « P » et
ses faces qui sont « est «, « ouest », « nord » et « sud » désignés respectivement par « e » , « w »,
« n » et « s » figure 3.
Les équations (7, 8, 9) peuvent s’écrire sous la forme condensée (voir chapitre 1) suivante :

( ) ( ) ( ) (74)

: est une variable physique générale qui représente U,V et

 35

: Terme source correspondant à la variable .

Donc c’est les méthodes de discrétisation et de résolution proposées et développer dans les
chapitres 1 et 2, qui vont être appliqué dans ce cas. On utilisera un maillage décalé, pour les, les
vitesses ne seront donc pas évaluées aux nœuds principaux des volumes mais plutôt aux interfaces
de ceux-ci.

N
Volume fini V
𝛿𝑦𝑛 N’
n
𝛿𝑦(𝑗)
Volume fini U

W P(xi ,yj) E
w e
dy(j)
𝛿𝑥𝑒

𝛿𝑥(𝑖)
s
y S’
dx(i)

S
x

Figure 3: Position et dimensions du volume de contrôle

On a () , ( )

() ( )

() , ( )

() ( )
 36

a- Equation de mouvement suivant x

L’équation algébrique de conservation de la quantité de mouvement permet de déterminer la

composante du champ de vitesse suivant la direction X. L’équation du mouvement discrétisée dans la
direction (x) prend la forme suivante :
∑ ( ) (75)

Où le terme désigne le terme source correspondant à la composante de la vitesse U, excepté le

gradient de pression.
On peut noter d’après la figure 3 que la composante U de la vitesse est calculée sur la face e qui
représente le centre du volume de contrôle et les faces sont : (P≡w, E≡e, N’≡n, S’≡s)
Le terme source est donné par :
() () (76)
Comme la température sera calculée aux nœuds principaux la valeur aux interfaces sera calculée
de manière suivante :
( ) (77)
f : coefficient de pondération géométrique définie comme suit figure 3 :
𝛿𝑥𝑒

P e x
E
𝛿𝑥𝑒 𝛿𝑥𝑒
b- Equation de mouvement suivant y
L’intégration de l’équation de mouvement dans la direction Y est faite dans le volume de contrôle
décalé vers le haut le nœud centrale est n.
∑ ( ) (78)

Où le terme désigne le terme source correspondant à la composante de la vitesse V, excepté le

gradient de pression.
On peut noter d’après la figure 3 que la composante V de la vitesse est calculée sur la face n qui
représente le centre du volume de contrôle.
Le terme source est donné par :
() () (79)
Comme la température sera calculée aux nœuds principaux la valeur aux interfaces sera calculée
de manière suivante :
( ) (80)
f : coefficient définie comme suit figure 3 :
 37

c- Discrétisation de l’équation de l’énergie

L’intégration est faite dans les volumes de contrôle correspondants aux nœuds principaux de
l’équation (70), avec un terme source nulle on aura donc :

∑ (81)

Les coefficients sont définis au chapitre 1 par les équations (43).

d- Discrétisation des conditions aux frontières :

La qualité de toute simulation numérique est fortement liée à la façon avec laquelle ces conditions
aux frontières sont calculées ou imposées (voir chapitre 2, §2).

2-Résolution de équations discrétisées

La méthode numérique de résolution est présentée au chapitre 2, §3. Avec une modification au
niveau des étapes de l’algorithme qui font inclure le couplage vitesse-température, soit :

Algorithme de Résolution (SIMPLER →SIMPLEC)

1°/ Estimer les champs de vitesse et de Température
2°/ Calculer les coefficients de mouvement à partir du champ obtenu en (1) et les vitesses ̂ ̂ en
utilisant l’équation correspondantes.
3°/ Intégrer l’équation discrétisée de pression pour obtenir du champ de pression P*
4°/ Intégrer les équations de mouvements pour obtenir le champ de vitesse U* et V*
5°/ Intégrer l’équation discrétisée de la correction de pression pour obtenir le champ de pression P’
6°/ Corriger le champ de vitesse U et V.
7°/ Intégrer l’équation discrétisée de l’énergie pour obtenir le champ de température (T)
8°/ Test de convergence :
* oui passer à l’étape suivant
* Non revenir à l’étape n°2 avec les champs de vitesse (U,V) comme nouvelles estimations
9°/ fin des calculs
****************************************