Polycope FVM Vfinale Francais V 01

UNIVERSITE KASDI MERBAH OUARGLA
Faculte des Sciences Appliquées

Département Génie Mécanique
Polycopié de cours
Méthode des volumes finis
Réalisé par :
Dr: BELAHYA Hocine
Langue Version
Arabe
English
Francis 01
1
PRÉFACE
L’objectif de ce document pédagogique est de permettre aux étudiants de

première année Master en génie mécanique Énergétique d’acquérir certaines
notions fondamentales d’ingénierie numérique.
Pour une meilleure maitrise de la méthode des volumes finis le polycopié est
répartie en cinq chapitres :
Le premier chapitre illustre les notions des phénomènes de

transport comporte le bilan global et le bilan différentiel, des généralités sur
l’équation de transport convecto-diffusif avec un rappel physique présente
l'analogie entre les trois lois phénoménologiques de la diffusion, tel que le
transport de chaleur, de matière, et de quantité de mouvement.
Le deuxième chapitre présente les notions et l'approche de la
méthode des volumes finis, ainsi que un aperçu sur le calcul des coefficients de
diffusion aux interfaces et la prise en compte des conditions aux limites ont été
développées. La dernière partie met l’accent sur le choix de maillage et leur
indépendance.
Le troisième chapitre consacré à la méthode des volumes finis pour
les problèmes de diffusion. Le contenu du chapitre est suffisant pour
comprendre la pluparts des cas de transport diffusif tel que le transport de
chaleur.
Le quatrième chapitre montre la méthode des volumes finis pour les
problèmes de convection-diffusion. Ce chapitre à pour objectif de connaître les
fondements de l’élaboration des différents schémas de convection tel que le
schéma centré, le schéma avant et le schéma d'ordre supérieur. Puis la prise en
considération du gradient de pression a été introduite pour résoudre les équations
de Navier-Stokes par les algorithmes les plus utilisés tel que l'algorithme
SIMPLER et l'algorithme SIMPLEC.
Le dernier chapitre est conservé pour la résolution des équations

algébriques discrétisées. Ce chapitre contient des généralités sur le facteur de
relaxation, Les méthodes par blocs et les méthodes de résolution itérative telle que
la méthode de Jacobi et la méthode de Gauss-Seidel, ces méthodes contiennent des
définitions ainsi que des informations sur les conditions de convergence.
Je tiens à remercier les collègues qui ont bien voulu juger le manuscrit et
m’aider à l’améliorer. Il est possible que cette première version comporte
quelques imperfections, je serais reconnaissant à tous ceux qui me ferait part de
leurs remarques et suggestions.
2
SOMMAIRE
Chapitre 1 : Notions des phénomènes de transport
1.1 Bilan global et bilan différentiel 05

1.1 Le transport diffusif 06
1.3 Le transport convectif 08
1.4 L’équation de transport convecto-diffusif 08
1.4.1 Cas particuliers de l'équation générale 11
1.4.1.1 Transport de matière 11
1.4.1.2 Transport de quantité de mouvement 12
1.4.1.3 Transport de chaleur 13
1.4.1.4 Transport de turbulence 13
1.5 Conditions aux limites 14
1.5.1 Condition de Neumann 14
1.5.2 Condition de Fourier 15
1.5.3 Condition de Dirichlet 16
Chapitre 2 :Notion de la méthode des volumes finis
2.1 Méthode de conservation du volume fini 17

2.2 L'approche de la méthode des volumes finis 18
2.5 Calcul des coefficients de diffusion aux interfaces 21
2.6 Prise en compte des conditions aux limites 23
2.7 Le Maillage 25
2.7.1 Le choix de maillage 26
2.7.2 Indépendance du maillage 26
Chapitre 3 : Méthode des volumes finis pour les problèmes de diffusion.
3.1 Problème conduction pure 27

3.2 Loi de variation linéaire 28
Chapitre 4: Méthode des volumes finis pour les problèmes de convection-
diffusion.
4.1 Le schema centré (Centered Differencing Scheme, CDS) : 32

4.2 Schéma avant (Upwind Differencig Scheme, UDS) 34
4.3 Le schéma hybride 35
4.4 Schéma la loi de puissance 35
4.5 Schéma d'ordre supérieur 36
4.5.1 Schéma QUICK 38
4.5.2 Schéma TVD (Total Variation Diminishing) 39
4.6 Prise en considération du gradient de pression 41
4.6.1 Le problème du damier 43
4.6.2 Algorithme de Patankar 43
4.6.3 L'algorithme SIMPLER (SIMPLE-Revised) 45
4.6.4 L'algorithme SIMPLEC (SIMPLE-Consistent) 46
3
Chapitre 5: Solution des équations algébriques discrétisées.
5.1 Système matriciel 50

5.2 Les méthodes de résolution directe 51
5.3 Les méthodes de résolution itérative 54
5.3.1 méthode de Jacobi 55
5.3.2 La méthode de Gauss-Seidel 56
5.3.3 Convergence et divergence 57
5.4 Notion de relaxation 59
5.4.1 Facteur de relaxation 60
5.4.2 La relaxation par faux pas de temps 61
5.5 Les méthodes par blocs 61
Références 63
4
Chapitre 1 : Notions des phénomènes

de transport
1.1 Bilan global et bilan différentiel
Les lois de conservation de la matière, de l'énergie, et de la quantité de

mouvement correspondent à trois grandes lois de la physique classique; la loi de
Lavoisier, le premier principe de la thermodynamique, la loi fondamentale de la
dynamique de Newton. La notion de bilan d'une extensité φ(masse ou nombre
de moles d'une espèce chimique, enthalpie, composante selon une direction x
quelconque de la quantité de mouvement) est tout à fait intuitive ; un bilan
s'écrira toujours:
Entre deux temps t1 et t2 dans un volume de référence V de la manière suivante:
entrée de Φ + création de Φ = sortie de Φ +accumulation de Φ

c'est-à-dire:
quantité de Φ qui rentre dans V par une surface entre t1 et t2
+
quantité de Φ produite à l'intérieur de V entre t1 et t2
=
quantité de Φ qui sort de V par une surface entre t1 et t2
+
quantité de Φ qui s'accumule à l'intérieur de V entre t1 et t2
Les premier et troisième termes correspondent à des transports d'extensité à

travers des surfaces: Les second et quatrième termes sont généralement de
nature volumique.
La distinction entre bilan global et bilan différentiel est directement liée à la

taille du volume dans lequel est effectué le bilan. Si ce volume est de dimensions
finies, on parlera de bilan global. Si l'une au moins des dimensions du volume est
infinitésimale, on écrira un bilan différentiel. Dans la pratique, on choisit le plus
souvent un volume qui puisse être considéré comme homogène (notion de
réacteur parfaitement agité) du point de vue de l'extensité qu'on étudie. De ce
fait, l'écriture de bilans globaux semble offrir peu d'intérêt pour une étude fine
des phénomènes de transport. Par contre, les bilans différentiels, qui conduisent
5
à des équations aux dérivées partielles, sont abondamment utilisés. Nous

verrons toutefois que la méthode des volumes finis, qui est exposée dans ce
cours, fait largement appel à la notion de bilan global.
1.2 Le transport diffusif :
Les lois phénoménologiques de la diffusion peuvent se mettre sous la forme:
φ = −λ grad V Eq. 1-1
Où : la densité de flux diffusifs de l'extensité Φ
: le coefficient de conductivité de ce transport
V : le potentiel dont le gradient provoque la diffusion de l'extensité Φ

En fait, Equation 1-1 est déjà un cas particulier. Tout d'abord, elle ne tient aucun
compte du transport diffusif d'une extensité (par exemple, la matière) sous
l'action de gradients de potentiels divers (par exemple, la température : effet
Soret). Une écriture générale nécessiterait l'emploi de l'ensemble des coefficients
phénoménologiques d'Onsager. Ensuite, le transport diffusif de quantité de
mouvement ne peut s'exprimer, sous une forme aussi simple que Equation 1-1
que dans le cas d'un écoulement monodirectionnel ou de celui d'un fluide
newtonien, incompressible, et à viscosité constante. Enfin, la conductivité λ est,
sous sa forme générale, un tenseur a 9 composantes. Seul le cas de matériaux
isotropes est traité ici. Lorsqu'on a affaire à des matériaux orthotropes fibreux
comme le bois ou le graphite, par exemple, les transports dans différentes
directions ont des conductivités très différentes.
L'extensité spécifique (ici, par unité de masse), notée par la suite Φ, et le

potentiel V peuvent le plus souvent être reliés par une relation de type :
dΦ = C dV Eq. 1-2
qui définit C = dΦ / dV, la capacité de l'extensité. On fait généralement apparaître

le coefficient Γ= λ/C ; la densité de flux diffusif devient alors :
Eq. 1-3
ϕ = − Γ grad φ
6
loi de Fick pour la diffusion chimique, loi de Fourier pour la diffusion thermique,
loi de Newton pour le frottement visqueux, auxquelles on peut ajouter la loi
d'Ohm pour le courant électrique et la loi de Darcy pour l'écoulement dans un
milieu poreux, Le tableau 1-1 présente l'analogie entre les trois lois
phénoménologiques de la diffusion pour les transports de chaleur, de matière, et
de quantité de mouvement (ainsi que pour le transport d’électricité). Comme on
l'a déjà précisé, les lois utilisées sont des lois simplifiées, qui supposent entre
autres que le milieu est incompressible. C'est l'analyse physique du problème à
traiter qui permet de remplacer éventuellement ces lois par des lois de même
type, qui peuvent être alors un peu plus complexes.
Tableau 1-1: Lois de la diffusion (analogies entre les transports)
Transport de Transport Transport de Transport

matière quantité de d'électricité
de chaleur mouvement
(constituant i) (pour un fluide
newtonien)
Extensité mi masse de i H enthalpie Px =mvx quantité de Q quantité

Φ mouvement dans la d'électricité (Cb)
(kg) (J) direction x ( kg m s-1)
Densité de φi densité de flux φ densité de flux τ contrainte de j densité de

flux diffusif φ massique chaleur frottement courant
(kg m-2 s-1) (W m-2) ( N m-2) (A m-2)
Potentiel V ρi concentration T température vx vitesse dans la U potentiel

massique (kg m-3) absolue ( K ) direction ( m s-1) électrique ( V )
Conductivité Di coefficient de λ conductivité μ viscosité dynamique σ conductivité

λ diffusion (m2 s-1) thermique ( kg m-1 s-1) électrique
( W m-1 K-1) ( Ώ-1 m-1)
Capacité 1/ρ inverse de la Cp chaleur C capacité

masse volumique spécifique électrique
C totale (m3 kg-1) massique 1
( J kg-1 K-1) (F)
7
1.3 Le transport convectif :
Lorsqu'en un point du système que l'on étudie, la vitesse barycentrique est

différente de zéro, il existe un mouvement d'ensemble du milieu qui transporte
la matière et les extensité associées, chaleur et quantité de mouvement. Ce
monde de transport, appelé transport par convection, est caractérisé par une
densité de flux :
ϕ=ρ φ v Eq. 2-4
Où ρ : est la masse volumique du système
Φ : l'extensité spécifique (par unité de masse)
v : le vecteur vitesse instantané
Ainsi, la densité de flux massique convectif d'un constituant dans un mélange
est égale à ρ v ωi , où ωi : est le titre massique de ce constituant, la densité de
flux thermique convectif à ρ v h , où h est l'enthalpie massique, et la densité
de flux convectif de quantité de mouvement selon x, à ρ v vx . On voit ainsi

que la densité de flux convectif du vecteur quantité de mouvement est un
tenseur, comme l'est d'ailleurs la densité de flux diffusif (tenseur des
contraintes).
1.4 L’équation de transport convecto-diffusif
On a vu que l'écriture d'un bilan différentiel d'une extensité quelconque conduit

au résultat suivant:
accumulation = - ( sortie-entrée) + création
∂ ( ρφ ) r Eq. 1-5
= − div ϕ + S
∂t
Où (ρ ϕ) est l'extensité volumique, et S la vitesse de création de Φ par unité de

volume (quantité de Φ créée par unité de volume et de temps). Ce terme de
création, qui est évidemment négatif lorsqu' il y a une consommation de Φ, peut
correspondre:
8
-pour le transport de matière, à la production (ou la disparition) de l'espèce

chimique i par réaction,
-pour le transport de quantité de mouvement, à une force volumique ou une

force de pression,
- pour le transport de chaleur, a un dégagement de chaleur par effet Joule ou par

frottement visqueux, ou du fait du travail des forces de pression, ou encore à
cause de l'exo- ou endo-thermicité de réactions chimiques
Dans l'équation 1-5, la densité de flux φ correspond à la somme des transports

convectif et diffusif, soit, en remplaçant par leurs expressions l'équation 1-3 et
l'équation 1-4 :
( ) Eq. 1-6
= − div(ρΦ v − Γ grad Φ) + S
Après un léger remaniement, l'équation générale de transport prend sa forme

définitive:
( ) Eq. 1-7
+ div(ρΦ v ) = div (Γ grad Φ) + S
L'universalité de l'équation 1-7 générale de transport convecto-diffusif a pour

conséquence directe que la même méthode numérique peut être utilisée pour
traiter l'ensemble des phénomènes de transport. En fait, toute équation aux
dérivées partielles qui peut se mettre formellement sous la forme de l'équation
1-7, ou de l'équation 1-14 que nous verrons plus loin, pourra être résolue par la
méthode numérique présentée dans ce cours, même si aucune extensité n'est
réellement transportée. Il suffira d'établir une correspondance terme à terme
pour définir une “ extensité ”, une “ conductivité ”...
L'équation 1-7 est composée de quatre termes qui correspondent aux

phénomènes physiques bien connus : accumulation, convection, diffusion
et création. Selon le problème considéré, chacun de ces termes peut disparaître
de l'équation. Par exemple, en régime permanent, le terme ∂(ρϕ)/ ∂t est nul
pour un fluide idéal (viscosité nulle), ou un constituant ne diffusant
pratiquement pas, le terme div (Γ grad ϕ) disparaît. . . Toute combinaison de ces
phénomènes est bien sûr possible. En cas de doute, on pourra toujours refaire la
9
démarche complète en remplaçant les équations générales 1-2 à 1-4 par les
équations correspondant au mieux au phénomène étudié.
Après un bref examen, il apparaît clairement que les dérivées du second ordre
(de type ∂²ϕ/∂ x²) seront toujours dues à la présence des termes diffusifs. À
l'inverse, les termes convectifs sont à l’ origine des dérivées du premier ordre
(type ∂ϕ/ ∂x) par rapport aux variables d'espace. Si, pour une variable d'espace
donnée, le terme diffusif peut être négligé (ce qui signifie en pratique que le
transport convectif dans la direction de cette variable est très important),
l'équation sera dite parabolique par rapport à cette variable. Dans le cas
contraire, où la dérivée de second ordre est présente, l'équation sera dite
elliptique par rapport à la variable d'espace. Un cas typique d'équation de
transport parabolique dans une direction et elliptique dans une autre se
rencontre lorsqu'on veut traiter un problème d'écoulement monodirectionnel
avec transport de chaleur, comme le représente sur la figure 1-2:
T2
x=0 x =L
Figure 1-2: Equation elliptique ou parabolique
Entre deux plaques planes portées à des températures différentes, un fluide
s'écoule parallèlement aux plaques. En régime permanent et sans terme de
création, l'équation de transport de chaleur au sein de ce fluide sera, si les
propriétés physiques (p,Cp,λ) sont indépendantes de la température (voir plus
loin la discussion de l'équation 1-22 :
² ² Eq. 1-8
ρCp v x = λ div (grad T) =λ ( + )
² ²
Si La vitesse v x est suffisamment importante, ce qui peut être caractérisé par la
valeur du nombre adimensionnel de Péclet ρ Cp v x L / λ, le terme diffusif

(dérivée du second ordre) dans la direction x, λ ∂² T/ ∂x² , est négligeable
10
par rapport au terme convectif (dérivée du premier ordre), ρ Cp v x ∂² T/ ∂x².

L'équation de transport de chaleur devient alors:
² Eq. 1-9
ρ Cp v x =
²
Selon nos définitions, il s'agit bien d'une équation parabolique suivant x et

elliptique suivant y.
L'intérêt d'une telle distinction sera particulièrement mis en évidence quand

nous en viendront à l'étude des conditions aux limites. On peut immédiatement
remarquer que l'équation de transport 1-7, en régime instationnaire, est toujours
parabolique par rapport au temps.
Par ailleurs, tous les bilans que nous écrirons dans le cadre de ce cours seront
différentiels du point de que du temps, c'est-à-dire que les instants t1 et t2 sont
infiniment voisins et séparés par un intervalle de temps infinitésimal dt.
1.4.1 Cas particuliers de l'équation générale
Nous allons à présent étudier quelques cas particuliers de l'équation générale.

correspondant aux transports de matière, de quantité de mouvement, de chaleur
, ou d'autres extensité.
1.4.1.1 Transport de matière
Lorsque l'extensité considérée est la masse d’un constituant i, l'équation

générale de transport devient, en remplaçant terme à terme :
∂ ( ρω i )
∂t
+ div (ρω v ) =
i div ( D i grad (ω i )) + S i
Eq. 1-10
où ωi est le titre massique en constituant i et si le taux volumique de production

de i par réaction chimique.
Si l'extensité considérée Φ est la masse totale, l'extensité spécif que ϕ est

alors égale à 1, et le flux de diffusion est nul. Par ailleurs, si l'on exclus le cas des
11
réactions nucléaires, le terme de création est également nul. L'équation ,générale

de transport I-7 se réduit alors tout simplement a l'équation de continuité:
∂(ρ )
∂t
+ div (ρ v ) = 0
Eq. 1-11
Il est très courant d'utiliser à la place de l'équation générale 1-7 une équation
obtenue en retranchant de l'équation 1-7 le produit de l'équation de continuité
par la variable ϕ.
∂ (φ ) Eq. 1-12
ρ + ρ div (φ ) = div ( Γ grad (φ )) + S
∂t
Alors que l'équation 1-7 provient directement d'un bilan différentiel et exprime
la conservation de l'extensité ϕ, l'équation 1-12 est le résultat d'une
manipulation mathématique et n'a pas de signification physique absolue. Lors
de la résolution numérique par ta méthode des volumes finis, on essaiera le
plus souvent de discrétiser la forme conservative de l'équation 1-7 générale de
transport.
1.4.1.2 Transport de quantité de mouvement
Le transport de la quantité de mouvement selon une direction x s'écrit, à partir

de l'équation générale :
∂ ( ρv x )
∂t
( )
+ div ρφ v = div ( µ grad ( v x )) + S i
Eq. 1-13
Dans le terme S interviennent alors principalement les forces volumiques (et en

particulier la pesanteur ρ gx) et les forces de pression - "P / "x. En fait, pour un
fluide non-newtonien en particulier, des termes de frottement visqueux
supplémentaires peuvent apparaître, que nous négligerons dans la suite de ce
cours. L'équation 1-13 devient alors la “ célèbre ”équation de Navier-Stokes:
Des équations similaires régissent évidemment le transport des autres
composantes de la quantité de mouvement.
∂ ( ρv x )
∂t
+ div (ρ vv ) = − ∂∂px
x div ( µ grad v x ) + ρg x
Eq. 1-14
12
1.4.1.3 Transport de chaleur

L'extensité transférée est ici l'enthalpie, ou quantité de chaleur (si on néglige le
travail des forces de pression). Et h est l'enthalpie massique, on écrit directement
l'équation de transport :
∂(ρh)
∂t
( )
+ div ρhv = div(
λ
Cp
grad h) + q
Eq. 1-15
Si, dans le domaine de températures considéré, on peut supposer que la

chaleur spécifique Cp est constate, on peut alors poser:
Eq. 1-16
h = CpT+ h0
L'équation 1-16 devient alors

Cp
∂ ( ρT )
∂t
+ Cp div (ρT v ) = div ( λ grad T ) + q
Eq. 1-17
ou. sous la forme non conservative (c'est-a-dire en retranchant le produit de T

par l'équation de continuité 1- 13) :
ρ Cp
∂ (T )
∂t
( )
+ ρ Cp div T v = div (λ grad T ) + q
Eq. 1-18
En fait, cette équation reste correcte même si la chaleur spécifique

varie avec la température. La démonstration devient simplement beaucoup plus
difficile. Par contre, elle n’est valable que s’ il existe clairement une chaleur
spécifique, c’est-à-dire en particulier que le système ne subit aucun changement
de phases.
Si la conductivité thermique λ est également constante, l'équation 1-18 prend

une forme condensée:
∂T
∂t
( )
+ div T v =
λ
ρCp
div( grad T ) +
q
ρCp
Eq. 1-19
qui est peut-être la forme la plus connue de l'équation de la chaleur
1.4.1.4 Transportde turbulence

A ce stade, notons simplement que la résolution numérique des problèmes de
transport dans un fluide en écoulement turbulent impose le plus souvent la
résolution d'équations convecto-diffusives supplémentaires, telles que celle qui
13
exprime le transport de l'énergie cinétique turbulente par unité de masse k. Cette

équation s'écrit :
∂ ( ρk )
∂t
( )
+ ρ div φ vk = div (Γk grad k ) + S k
Eq. 1-20
Les termes Γket Sk ne serons pas explicités pour l'instant. Il nous suffit de savoir
que sur le plan numérique, la nature laminaire ou turbulence de l'écoulement ne
change rien a la méthode fondamentale de résolution.
1.5 Conditions aux limites
résolution par une méthode analytique ou ,numérique, de l'équation générale

de transport (convecto-diffusif ) nécessite la détermination de ses conditions
initiales et spatiales (aux limites). Il existe fondamentalement trois types de
conditions aux limites (ou en abrégé C.L) qui portent sur la détermination de la
densité de flux d'extensité à la frontière du système étudié. Par ailleurs, la nature
de l'équation, elliptique ou parabolique, par rapport aux différentes variables
d'espace et de temps, impose la nécessité de ces conditions. Nous allons tout
d'abord passer en revue les différentes catégories des conditions aux limites,
puis nous examinerons les exigences mathématiques de l'équation aux dérivées
partielles en ce qui concerne ces conditions aux limites.
Les trois sortes de conditions aux limites à la frontière du système étudié, la

densité de flux d'extensité normale φ#. % que nous noterons φ est donnée par
l'une des trois conditions aux limites suivantes.
1.5.1 Condition de Neumann
La condition de Neumann exprime la connaissance de la valeur de φ en un

point de la frontière et s'écrit par conséquent en ce point : ϕ=cte
Eq. 1-21
ϕ = cte
La condition de Neumann est fréquemment utilisée :
-lorsque des mesures expérimentales de flux ont été réalisées
-sous la forme φ = 0, pour représenter les conditions de symétrie,
14
-également sous la forme φ = 0, pour simuler :
• l'adiabatisme (transport de chaleur),
• l'imperméabilité (transport de matière)
• le frottement nul (ou glissement parfait) à la surface libre d'un

écoulement (transport de quantité du mouvement)
1.5.2 Condition de Fourier
La condition de Fourier exprime l'existence d'une relation entre densité de flux

normale et extensité spécifique en un point de la surface, et s'écrit de la manière
générale:
Eq. 1-22
φ = f (φ)
Dans la majorité des cas, cette relations est linéaire. Un exemple typique est de
transfert de chaleur à la surface d'un solide pur convection de l'atmosphère
ambiante. La densité de flux de chaleur sortant du solide peut s'exprimer sous la
forme :
Eq. 1-23
φ = h(T − Text)
Où Text est la température de l'air ambiant et h un coefficient de transfert, ou

conductance. qui exprime la résistance thermique correspondant à
l'établissement d'une couche limite d’air au voisinage de la surface. Ce coefficient
de transfert thermique ne doit bien sûr pas être confondu avec l'enthalpie
massique, également notée h. L'expression de l'équation 1-23 a l'aide de la
température plutôt que de l'enthalpie spécifique est très courante l'équation
1-18.
Une autre condition aux limites classique, pour le transport de chaleur, est le
rayonnement d'une surface:
Eq. 1-24
φ = σ ε (T 4 −Text4 )
15
ou σ est la constante de Stefan-Boltzmann et ε l'émissivité de surface (dans

l'hypothèse du corps gris).
1.5.3 Condition de Dirichlet
La condition aux limites de Dirichlet peut être considérée comme une forme
asymptotique de la C. L. de Fourier linéarisée, avec une valeur infinie du
coefficient de transfert. Elle exprime en fait la connaissance de la valeur de
l'extensité en un point de la frontière :
Eq. 1-25
φ = cte
Cette condition aux limites se rencontre classiquement:
-en transport de chaleur, lorsque des mesures de température sont disponibles,

ou pour exprimer un contact parfait avec un milieu dont la température est
connue,
-en transport de matière, pour exprimer l'équilibre thermodynamique avec une

phase de composition connue,
-en mécanique des fluides, pour simuler une vitesse imposée, ou une vitesse
nulle en un point au contact d'une paroi solide immobile.
16
Chapitre 3: Méthode des volumes finis pour les problèmes de diffusion.
Chapitre 2: Notions de la Méthode des

volumes finis
2.1 Méthode de conservation du volume fini
Si nous utilisons l’approche de la différences finies et les éléments finis pour

l’équation de Navier-Stokes discrétisée, nous devons contrôler manuellement la
conservation de la masse, de la quantité de mouvement et de l’énergie. Mais avec la
méthode des volumes finis, nous pouvons facilement découvrir que, si l’équation de
Navier-Stokes est satisfaite dans chaque volume de contrôle, elle sera automatiquement
pour l’ensemble du domaine. En d'autres termes, si la conservation est satisfaite dans
chaque volume de contrôle, elle sera automatiquement satisfaite dans tout le domaine.
C'est la raison pour laquelle le volume fini est préféré dans la dynamique des fluides
computationnelle.
l’approche de la méthode des volumes finis consiste à diviser le domaine de calcul

en plusieurs petits volumes qui ne se chevauchent pas et dont la somme fait exactement
le volume du domaine de calcul à étudier. C’est très important pour assurer le principe
de conservation et surtout la conservation des flux entre l’entrée et la sortie du domaine
de calcul.
Pour simplifier les équations de Navier-Stokes, nous pouvons les réécrire sous la forme
générale.
∂ (ρ Φ ) ∂  ∂Φ  Eq. 2-1
+  ρ U i Φ − ΓΦ  = qΦ
∂t ∂xi  ∂ x i 
Quand: Φ = 1,U j , T ,nous pouvons respectivement obtenir l’équation de continuité,
l’équation du moment et l’équation de l’énergie.
17
2.2 L'approche de la méthode des volumes finis

Cette méthode basant sur l'intégration de la forme générale de l'équation de Navier-
Stokes sur un volume de contrôle et appliquez la théorie de Gauss
∂ Eq. 2-2
∫ ∂x ΦdV = ∫ Φ ⋅ n dS
V i S
i
Nous pouvons obtenir la forme intégrale de l'équation de Navier-Stoke
∂ (ρ Φ )  ∂Φ  Eq. 2-3
∫V ∂t dV + ∫S  ρU i Φ − Γ ∂xi  ⋅ ni dS = V∫ qΦ dV
Pour approximer l'intégrale du volume, nous pouvons multiplier le volume et la valeur

au centre du volume de contrôle. Par exemple, nous avons un domaine 2D comme sur la
figure 2-1. Pour approximer la masse et la quantité de mouvement du volume de
contrôle P, nous avons:
m = ∫ ρ dV ≈ ρ pV Eq. 2-4
Vi
mu = ∫ ρ i u i dV ≈ ρ P u PV
Vi
Eq. 2-5
Pour approcher l’intégrale de surface, par exemple la force de pression, nous avons
∫ P dS ≈ ∑ Pk S k k = n, s, e, w Eq. 2-6
Si
k
Normalement, nous stockons nos variables au centre du volume de contrôle. Nous
devons donc les interpoler pour obtenir Pk , qui sont situées à la surface du volume de
contrôle.
En règle générale, nous avons plusieurs types d’interpolations (l’interpolation au avant,

l’interpolation centrale…et autre pour plus de détail voir Ch. 4 transport convectif ).
18
UP UE
Ue
Figure 2-1 Domaine du maillage structuré en 2D
Pour rendre cette présentation plus concrète, nous appliquons maintenant la

méthode des volumes finis dans le cas de l'équation de diffusion de la chaleur pour une
géométrie monodirectionnelle cartésienne, en régime permanent. Dans le cadre de ces
hypothèses, l'équation de la chaleur prend la forme suivante :
d  dT  Eq. 2-7
Γ =S
dx  dx 
Comme il a été déjà montré pour un volume de contrôle interne, l'équation 2- 7

est discrétisée sous la forme où Γ est la conductivité thermique, T la température et S
terme source représentant la création ou l'absorption d'énergie par unité de volume.
Figure 2-2: conduction dans une géométrie monodirectionnelle
Nous appliquons la méthode pour le volume de contrôle construit autour du

point P Figure 2-2. Nous appelons 'East' le premier voisin du point P dans la direction
des x croissants et 'West' le premier voisin dans la direction des x décroissante, les deux
19
traits verticaux en pointillé délimitent le volume de contrôle associé à P. Les frontières

du volume de contrôle sont représentées par les lettres minuscules e et w.
e e
d  dT 
∫
w
Γ  dx =
dx  dx  ∫
w
Sdx
Eq. 2-8-1
On suppose que est constant Γ = Γe = Γw
e e w
 dT   dT   dT  Eq. 2-8-2
Γ  − S δx ew =  Γ  − Γ  − S δ x ew = 0
 dx  w  dx   dx 
e
 dT   T − TP  Eq. 2-8-3
Γ  =  Γ E 
 dx   (δx ) e 
w
 dT   T − TW  Eq. 2-8-4
Γ  =  Γ P 
 dx   (δx ) w 
 TE − TP   TP − TW  Eq. 2-8-5
 Γ  −  Γ  + S ∆ x = 0
 (δx ) e   (δx ) w 
et après arrangement :
aPTP = aETE + aWTW + b Eq. 2-9
où
aP = aE + aW Eq. 2-9-1
Γ
aE = Eq. 2-9-2
(δx) e
Γ
aW =
(δx) w Eq. 2-9-3
Eq. 2-9-4
b = S ∆x
20
2.5 Calcul des coefficients de diffusion aux interfaces
Dans l'équation 2-7, si la valeur du coefficient de diffusion est variable en fonction de x.
Γe représente la valeur du coefficient de diffusion à l'interface e du volume de contrôle. Dans

ce cas où Γ dépend de x, sa valeur est connue uniquement aux points W, P et E. Ceci se
produit par exemple avec un domaine composé de différents milieux, ou si Γ dépend de φ
Nous décrivons maintenant une méthode simple pour déterminer Γe
Une technique possible consiste à considérer une variation linéaire de Γ entre les points P
et E, ce qui permet d'exprimer Γe ,
Γe = f e ΓP + (1− f e ) ΓE Eq. 2-10
(δx)e +
avec fe = Eq. 2-11
(δx)e
Figure 2-3 :Distances associées à l'interface
Si, comme sur la Figure 2-3, l'interfacée est située au milieu du segment (P,E), fe prend la
valeur 0,5 et Γe est la moyenne arithmétique de ΓP et ΓE .
Nous allons montrer que cette technique s'avère mauvaise, notamment dans le cas de
fortes variations de Γ . Notre objectif principal n'est pas en fait de déterminer Γe . Mais
d'estimer le flux ϕe de l'extensité φ à l'interface e par l'expression
21
Γe (φP − φE ) Eq. 2-12

ϕe =
(δx)e
La “bonne" valeur Γe est celle qui conduit à une estimation correcte de ϕe .

Considérons que les volumes de contrôle construits autour des points P et E sont
constitués par des milieux dont le coefficient de diffusion Γ est constant. Quand le
terme source est nul, la valeur de ϕe peut être déterminée à partir de la solution exacte
de l’équation (2-7)
(φP − φE )
ϕe =
(δx ) e − (δx ) e +
+ Eq. 2-13
ΓP ΓE
les équations 2-11 à 2-13 permettent d'écrire

−1
 1 − fe 1 − fe 
Γe =  + 
 P Γ ΓE  Eq. 2-14
Quand l'interface e est placée au milieu du segment (P,E), fe = 0,5 et
1  1 1 
= 0.5  + 
Γe  ΓP ΓE  Eq. 2-15
Γe est la moyenne harmonique de ΓP et ΓE .
Les équations 2-9-1 et 2-14 permettent d'exprimer aE sous la forme:

−1
 (δx) (δx) 
aE =  e− + e+ 
 ΓP ΓE  Eq. 2-16
Ainsi si ΓE →0, l'équation (2-14) permet de voir que:

Eq. 2-17
Γe →0
Le flux ϕe déterminé a partir des relations 2-12 et 2-17, est nul à l'interface e, ce
qui est normal puisque le coefficient de diffusion de la cellule centrée sur E est nul.
Dans le cas d'une prise de moyenne arithmétique équation 2-10, ϕe n'aurait pas
été nul. L'équation 2-14 permet un calcul exact de ϕe sa justification a été menée pour
une géométrie simple monodirectionnelle, en régime permanent, sans terme source

22
et avec un coefficient de diffusion constant dans chaque cellule. Même dans des cas
moins triviaux (prise en compte d'un terme source, variation continue du coefficient de
diffusion); la moyenne harmonique donne de meilleurs résultats que la moyenne
arithmétique.
2.6 Prise en compte des conditions aux limites
Considérons le maillage de la figure 2-6. A chacune des frontières correspond un

point du maillage. Les autres points son appelés des points internes. Autour de chaque
point interne est construit un volume de contrôle. Pour ces points, l'équation 2-1 est
discrétisée sous la forme de l'équation 2-8-5. Dans deux de ces équations, intervient la
valeur de φ sur une frontière.
Figure 2-6:volumes de contrôle pour un point interne et pour un point frontière.
Comme le traitement d'un point frontière ne dépend pas de sa localisation, nous

allons nous intéresser au point B. Les trois types de conditions aux limites couramment
utilisées sont:
condition de Dirichet φB = φ0 Où : φ0 est connu Eq. 2-18-1
condition de Neumann  dφ  Eq. 2-18-2

ϕ 0 = − Γ  Flux imposé
 dx  B
condition de Fourier Eq. 2-18-3

h
(ϕext − ϕB ) = − Γ dφ 
Cp  dx  B
23
Si la valeur de φ est imposée au point B, il n'est pas nécessaire d'écrire des

relations supplémentaires. Dans le cas contraire, nous obtenons une nouvelle équation
discrétisée par intégration de l'équation (2-7) sur le volume de contrôle frontière du
point B (figure 2-5).
figure 2-5:volume de contrôle frontière construit autour du point B
Cette intégration s'exprime sous la forme:

 dφ   dφ  Eq. 2-19
 Γe  −  Γ  + (SC + S PφB ) ∆x = 0
 dx e  dx  B
soit encore
Γe
(φE − φB ) −  Γ dφ  + (S + S φ ) ∆x = 0
Eq. 2-20
 
(δx )e  dx  B C P B
La suite de la discrétisation dépend du type de condition à la frontière.
a) condition de Neumann
 dφ  Eq. 2-21
−  Γ  = ϕ0
 dx  B
Dans ce cas prend la forme discrétisée suivante :
a B φB = a E φE + b Eq. 2-21-1
Γe
aE = Eq.2-21-2
(δx)e
b = SC∆x +ϕ0
Eq.2-21-3
aB = aE − SP∆x Eq.2-21-4
24
b) condition de Fourier
 dφ  h
− Γ  = (φ f − φB ) Eq.2-22
 dx  B C p
L'équation 2-20 est alors discrétisée sous la forme suivante :
aB φB = aE φE + b Eq. 2-23-1
Γe Eq. 2-23-2
aE =
(δx)e
h Eq. 2-23-3
b = SC ∆x + φext
CP
h Eq. 2-23-4
a B = a E − S P ∆x +
Cp
2.7 Maillage
La distance entre deux points du maillage n'est pas obligatoirement constante. Si

cette distance est constante, on parle de maillage régulier, si non de maillage irrégulier. Il
existe trois types de maillages: les maillages structurées, les maillages non structurées et
les maillages structurées en blocs. Le plus simple est la maillage structurée figure 2-7-1.
Dans ce type de maillage, tous les nœuds ont le même nombre d'éléments voisins. Nous
pouvons les décrire et les stocker facilement. Mais ce type de maillage ne concerne que
le domaine simple. Si nous avons un domaine complexe, nous pouvons utiliser une
maillage non structurée. Par exemple, la figure 2-7-2 profil aérodynamique. La structure
de la surface portante est très complexe. Les gradients des différents grandeurs telles
que la vitesse, la pression et la température à proximité de l'objet sont très important et
complexe, nous avons besoin d'un maillage raffiné très fin dans cette région. Loin du
profil aérodynamique, le gradient est relativement simple, nous pouvons donc utiliser
une maillage grossière. Généralement, la maillage non structurée convient à toutes les
géométries. Il est très populaire dans les CFD. L'inconvénient est que, la structure des
données étant irrégulière, Il est plus difficile de décrire et de stocker. La maillage de
structure en blocs est un compromis entre la maillage structurée et non structurée.
L'idée est tout d'abord de diviser le domaine en plusieurs blocs, puis d'utiliser
différentes maillages structurées dans les différents blocs.
25
Figure 2-7-1:maillage structurés Figure 2-7-2:maillage non structurés
2.7.1 Le choix de maillage:
Le choix du maillage dépend du problème posé ; dans une zone où Φ varie

fortement, il sera nécessaire d'employer des maillages fines, tandis que des maillages plus
larges pourront être utilisées dans des zones de variations plus faibles. Bien qu'il n'existe
pas de règle stricte, il ne faut pas passer brusquement d'une maillage très fine à une
maillage beaucoup plus large. En pratique, le rapport des dimensions entre deux mailles
voisines doit être compris entre 1/3 et 3.
2.7.2 Indépendance du maillage:
Si l'on ne connait pas à l'avance les zones de fort gradient, une méthode
consiste à faire un premier calcul avec un maillage grossier, puis si nécessaire à
raffiner localement le maillage. il faut par ailleurs toujours s'assurer que la solution
est indépendante du maillage.
Erreur
Résultats Indépendance du maillage
Nombre d'éléments
Figure 2-8: Indépendance du maillage
26
Chapitre III : Méthode des volumes

finis pour les problèmes de diffusion.
3.1 Problème conduction pure
Pour rendre cette présentation plus concrète, nous appliquons maintenant la méthode
des volumes finis dans le cas de l'équation de diffusion de la chaleur pour une géométrie
monodirectionnelle cartésienne, en régime permanent. Dans le cadre de ces hypothèses,
l'équation de lu chaleur prend la forme :
d  dT  Eq. 3-1
Γ =S
dx  dx 
Où Γ est la conductivité thermique, T la température et S terme source
représentant la création ou l'absorption d'énergie par unité de volume. Nous appliquons
la méthode pour le volume de contrôle construit autour du point P (figure3-1)
Figure3-1: Volume de contrôle construit autour du pointP
L'intégration de l'équation 3-1 sur le volume de contrôle donne:
 dT   dT  e Eq. 3.2
Γ  − Γ  + ∫w Sdx = 0
 dx  e  dx  w
27
3.2 Loi de variation linéaire

Nous représentons sue la figure 3-3-a, la loi de la variation dans le cas où l’on suppose
que la température est constante dans chaque volume de contrôle, dT / dx n'est pas
définie aux points frontières du volume de contrôle (w et e), tandis que la loi de la
figure (3-2-b) où l'on suppose une variation linéaire de la température entre deux
points du maillage, permet le calcul de dT/ dx
a) T constante dans chaque volume de b) variation linéaire de T entre deux

contrôle points du maillage.
figure 3-2:Types de loi de variation locale pour T
En adoptant cette loi de variation linéaire, l équation (3-2) devient
 TE − TP   TP − TW  Eq. 3-3
 Γe  −  Γw  + S ∆x = 0
 (δx ) e   (δx ) w 
 dT   TE − TP 

Eq. 3-3-1
  =
 dx  e  (δx ) e 
 dT   T − TW  Eq. 3-3-2
  =  P 
 dx  w  (δx ) w 
28
Où Γe et Γw sont les valeurs de la conductivité thermique sur les faces e et w du volume
de contrôle et S une valeur moyennée de S sur le même volume.
Pour simplifier l’exposé nous allons utiliser un maillage uniforme(∆X e=∆Xw=∆X)et
appliquer un schéma centré d’ordre deux pour remplacer les dérivés premières sur les
facettes du volume de contrôle.
Il est pratique d'écrire l'équation discrétisée sous la forme:
aPTP = aETE + aWTW + b Eq. 3.4
Où
aP = aE + aW Eq. 3-3-1
Γe
aE = Eq. 3-3-2
(δx) e
Γw
aW =
(δx) w Eq. 3-3-3
b = S ∆x
Eq. 3-3-4
1- L'équation 3-4 est la forme standard que nous adopterons pour l'écriture d'une
équation discrétisée. La température Tp au point P apparaît à gauche de l'équation,
tandis que la température des nœuds voisins et la constante b apparaissent à droite.
Comme nous le verrons plus loin (chapitre 4), le nombre des points voisins augmente
avec la complexité de la géométrie du problème (2D-3D). Nous noterons de façon
générale une équation discrétisée sous la forme:
aPTP = ∑ anbTnb + b Eq. 3-5

nb
Où l'indice nb désigne un voisin du point P, la sommation portant sur tous les voisins de
P.
2-Pour discrétiser l'équation (3-4), nous avons utilisé une loi de variation locale linéaire
pour calculer les dérivées aux frontières du volume le contrôle. Il est bien sur possible
d'utiliser une autre forme de loi.
3- Il est possible de choisir une loi de variation locale différente selon la grandeur à
29
exprimer. Par exemple, pour calculer S , il n'est pas nécessaire de prendre une loi
linéaire, une loi constante dans chaque cellule pouvant convenir.
4- La liberté de choix en ce qui concerne les lois de variation locale montre qu'il est
possible d'obtenir un grand nombre de forme d'équations discrétisées pour une même
équation aux dérivées partielles et un même maillage. Il est cependant important de
respecter deux règles.
a- Le choix de la loi de variation locale d’une grandeur doit avoir un

sens physique. Par exemple dans un problème de diffusion de la chaleurs sans
terme source, les extremums de la température se situent sur les frontières du
domaine (figure 4-3).
b- Il faut toujours s'assurer que la conservation globale est vérifiée. Nous

développerons ce point dans le paragraphe suivant.
Pas de sens physique
Solution approchée
ayant un sensphysique
Solution
exacte
Pas de sens physique
Figure 3-3 : Exemples de lois de variation locale avec ou sans sens physique.
5-traitement du terme source. Le terme source qui apparaît dans l'équation (3-4) peut
dépendre de la température. Nous allons supposer que cette dépendance est linéaire.
Dans ce cas, une bonne approximation pour S est :
S = SC + SPTP Eq. 3-6
30
Cette expression est obtenue en considérant une loi de variation locale ou la

température est, constante dans chaque cellule(figure3-2-b).La forme discrétisée de
l'équation(3-3)devient
Où
aP = aE + aW − SP∆x Eq. 3-7-1
Γe
aE = Eq. 3-7-2
(δx) e
Γw
aW =
(δx) w Eq. 3-7-3
b = SC ∆x
Eq. 3-7-4
31
Chapitre 4: Méthode des volumes finis pour les problèmes de convection-diffusion.
Chapitre 4: Méthode des volumes finis

pour les problèmes de convection-
diffusion.
On reprend ici, l’équation stationnaire de transport d'une variable par convection

diffusion sous sa forme générale, cette équation peut être interprétée comme une
équation d'équilibre de flux convecto-diffusif:
∂ (ρuiφ ) ∂  ∂φ  Eq. 4-1
= Γ +S
∂xi ∂xi  ∂xi 
La forme de l'équation (4-1) pour un problème à une seule dimension (1D), donne:
∂ (ρu φ ) ∂  ∂φ  Eq. 4-2
= Γ +S
∂x ∂x  ∂x 
(ρuφ )w Aw −  Γ dφ   dφ 
Aw − (ρuφ )e Ae −  Γ  Ae + S∆V = 0 Eq. 4-3
 dx  w  dx  e
Le terme source sera linéarisé suivant la forme suivante:

S = SU + SPφP Eq. 4-3-1
Pour le terme de la diffusion nous avons :
 dφ  φ −φ Eq. 4-3-2
Γ  = Γ P W
 dx  w ∆x
 dφ  φ −φ Eq. 4-3-3
Γ  = Γ E P
 dx e ∆x
On définit les deux coefficients suivants :

Eq. 4-3-4
F :flux convectif ; F = ρu A
ΓA Eq. 4-3-5
D : flux diffusif ; D =
∆x
32
Pour estimer les valeurs des termes convectifs (ρuφ)w et (ρuφ)e nous avons plusieurs
schémas, le choix est dépond du nombre adimensionnel de Peclet qui représente le
rapport du flux convectif sur le flux diffusif:
ρu F Eq. 4-4
Pe = =
Γ / ∆x D
Si Pe → 0 dans ce cas la vitesse est nulle, le fluide réagi comme un corps solide
Si Pe → ∞ le flux convectif est dominant
4.1 Le schema centré (Centered Differencing Scheme, CDS) :
Pour ce schéma en prenant la valeur moyenne des φw et φe comme suit :
φ P + φW Eq. 4-4-1
φw =
2
φE + φP Eq. 4-4-2
φe =
2
L'équation 4-2 devient

Fe Eq. 4-6
(φ E + φ P ) − Fw (φW + φ P ) + S ∆V = De (φ E − φ P ) − Dw (φ P − φW )
2 2
Et après arrangement:
 Fw   Fe    Fw   Fe  Eq. 4-7
 Dw − 2  +  De + 2  − S p φ P =  Dw + 2 φW +  De − 2 φ E + SU
        
 Fw   Fe   Fw  Fe  Eq. 4-8
 
 Dw + 2  +  De − 2  + (Fe − Fw ) − S p φP =  Dw + 2 φW +  De − 2 φE + SU
        
l'équation de continuité : (Fe − Fw ) = 0
Enfin:
où
aP = aE + aW + (Fe − Fw ) − SP Eq. 4-9-1
 F 
aE =  De + e  Eq. 4-9-2
 2
 F 
aW =  Dw − w 
 2 Eq. 4-9-3
Eq. 4-9-4
b = SU
33
4.2 Schéma avant (Upwind Differencig Scheme, UDS)
Ce schéma suit le sens de la vitesse, si u est supposée positive dans la direction de x:
Eq. 4-10
Fw ≥ 0 et Fe ≥ 0 : φw = φW et φe = φP
L'équation 4-2 devient:

 φP − φw   φ −φ  Eq. 4-11
ρuw AwφW − Γ  Aw − ρue AeφP + Γ E P  Ae + SU Aδx + SPφP Aδx = 0
 dx   dx 
Fwφw − Dw(φP −φw ) − Fwφw − De (φP −φw )e + SU Aδx + SPφP Aδx = 0 Eq. 4-11-1
Eq. 4-11-2
(Fw + Dw + De − SPδx)φP = (Fw + Dw)φw + (De )φe + SU Aδx = 0
et après arrangement:
où
aP = aE + aW − SP Aδx Eq. 4-12-1
Eq. 4-12-2
aE = De
aW = Fw + Dw Eq. 4-12-3
b = SU Eq. 4-12-4
si la vitesse u est supposée négative :
Fw ≤ 0 :φw = φW et Fe ≤ 0 :φe = φE
aP = aE + aW − SP Aδx Eq. 4-14-1
Eq. 4-14-2
aE = De − Fe
aW = Dw Eq. 4-14-3
b = SU Eq. 4-14-4
Une forme de notation générale pour les coefficients qui couvre les deux sens
d’écoulement est donnée ci-dessous:
34
Eq. 4-15-1
aE = De + max [0,−Fe ]
Eq. 4-15-2
aW = Dw + max [Fw,0]
4.3 Le schéma hybride (de Spalding1972)

Il repose sur une combinaison de deux schémas, la différence centré, qui est
précis au second ordre, est utilisé pour les petits nombres de Peclet (Pe< 2) et le schéma
Upwind, qui est précis au premier ordre, est utilisé pour les grands nombres de Peclet
(Pe ≥2).
Comme précédemment, nous développons la discrétisation de la convection
unidimensionnelle dans l'équation 4-2 de diffusion sans terme source il donne.
  F 
aW = max Fw ,  Dw + w ,0 
  2  Eq. 4-16-1
  F 
aE = max− Fe ,  De − e ,0 
  2 
Eq. 4-16-2
4.4 schéma la loi de puissance
Le schéma de la loi de puissance de Patankar (1980) est une approximation plus précise
de la solution exacte unidimensionnelle et donne de meilleurs résultats que le schéma
hybride. Dans ce schéma, la diffusion est fixée à zéro lorsque Pe supérieur à 10.
Si 0 < Pe< 10, le flux est évalué à l'aide d'une expression polynomiale. Par exemple, le
flux net par unité de surface au volume de contrôle ouest est évalué à l’aide des
propriétés du schéma de la loi de puissance sont similaires à celles du schéma hybride.
Le schéma de la loi de puissance est plus précis pour les problèmes monodimensionnels
puisqu'il représente la solution exacte.
[ ]
aW = Dw max 0, (1 − 0.1 Pew ) + max[Fw ,0]
5 Eq. 4-17-1
Eq. 4-17-2
[ ]
aE = De max 0, (1 − 0.1 Pee ) + max[Fe ,0]
5
35
4.5 Schémas d'ordre supérieur
4.5.1 Schéma QUICK (quadratic upstream interpolation for convective kinetics)le

schéma d'interpolation amont quadratique de Leonard (1979) utilise une interpolation
quadratique pondérée en amont à trois points pour les valeurs de face de cellule. La
valeur faciale de φw;e est obtenue à partir d'une fonction quadratique passant par deux
nœuds successifs (de chaque côté de la face) et un nœud du côté amont (Figure 4-1).
Figure 4-1:Profils quadratiques utilisé dans le schéma QUICK
si la vitesse est supposée positive Fw ≥ 0 et Fe ≥ 0 ; les coefficients est donné par la
formule suivante:
6 3 1 Eq. 4-18
a face = φ i −1 + φ i − φ i − 2
8 8 8
6 3 1
aw = φW + φ P − φWW Eq. 4-18-1
8 8 8
6 3 1
ae = φ P + φ E − φW
8 8 8 Eq. 4-18-2
et si nous utilisons les équations 4.18-1 et 4.18-2 pour les termes convectifs et la
différence centrale pour termes de diffusion, la forme discrétisée de l’équation
unidimensionnelle de transport de convection-diffusion Eq. 4-2 peut être écrite ainsi:
 6 3 1  6 3 1  Eq. 4-19
 Fe  8 φP + 8 φ E − 8 φW  − Fw  8 φW + 8 φ P − 8 φWW  = De (φ E − φ P ) − Dw (φP − φW )
    
36
 3 6   6 1   3  1
 D − F + D + F φ = D + F + F φ + D − F φ − FwφWW
8   8   8 
w w e e P w w e W e e E
 8  8  8 Eq. 4-19-1
La form compact :
aPφP = aWφW + aEφE + aWWφWW Eq. 4-19-2
aW aE aWW aP
 6 1   3  1 aW + aE + aWW
- Fw
 Dw + 8 Fw + 8 Fe  De − 8 Fe  8
   
Pour le flux Fw<0 et Fe <0, entre les limites sont donnés par les expressions
6 3 1
aw = φ P + φW − φ E
8 8 8
Eq. 4-20-1
6 3 1
ae = φ E + φ P − φ EE
8 8 8
Eq. 4-20-2
aPφP = aWφW + aEφE + aWWφWW Eq. 4-21
aW aE aEE aP
 3   6 1  1 aW + aE + aEE
Fw
 Dw + 8 Fw   De − 8 Fe − 8 Fw  8
   
L'expression générale, valables pour les directions de flux positives et négatives,
peuvent êtreobtenu en combinant les deux ensembles de coefficients ci-dessus
aPφP = aW φW + aEφE + aWWφWW + aEEφEE Eq. 4-22
et après arrangement
aW aWW aE aEE aP
6 1 1 3 6 1
Dw + α w Fw + α e Fe − α w Fw De − α e Fe − (1 − α e )F (1 − α e )Fe aW + aWW + aE + aEE
8 8 8 8 8 8
3 1
+ (1 − α w )Fw − (1 − α w )Fw
8 8
37
Où
Eq. 4-23-1
α w = 1 pour Fw > 0 et α e = 1 pour Fe > 0
Eq. 4-23-2
α w = 0 pour Fw < 0 et α e = 0 pour Fe < 0
-Condition aux limites pour le schéma QUICK
On peut facilement montrer que la valeur

extrapolée linéairement au nœud en question (A)
est donné par φ0 = 2φA −φP
En remplaçant la valeur φ0
6 3 1 Eq. 4-24-1
φe = φP + φE − (2φA − φ p )
8 8 8
7 3 2
= φP + φE − φ A Eq. 4-24-2
8 8 8
On peut montrer que le flux diffusif à travers la limite ouest (nœud A) est donnée par
∂φ DA* Eq. 4-24-3
Γ = (9φP − 8φA − φE )
∂x A 3
Où
Γ 2Γ
D*A = et DA = 2D =
∂x ∂x
Le tableau récapitulatif des coefficients
Nœud aWW aW aE SP Su
1 0 0 1 3 8 2  8 2 
De + D A* − Fe −  DA* + Fe + FA   D*A + Fe + FA φ A
3 8
3 8  3 8 
3 1 1
De − Fe Fw − Fwφ A
2 ..n-1 0 7 1 8 4 4
Dw + Fw + Fe
8 8
1 * 6 8  8 * 
Dw + D B + Fw −  DB* − FB   DB − FB φB
3 8 3  3 
n 1 0
− Fw
8
38
4.5.2 Schéma TVD (total variation diminishing)
TVD est un schéma compacte en tenant compte des propriétés de base des schémas
standard et de leurs carences. Il donne des résultats pour les problèmes de dynamique
des gaz. le schéma upwind est le plus stable et limité sans conditions, mais il introduit un
niveau élevé d'erreur de diffusion due à son faible ordre de précision (premier ordre).
Pour les schémas d'ordre supérieur tels le schéma QUICK peuvent donner des
oscillations lorsque le nombre de Peclet est élevé.
φ −φ  Eq. 4-25
ψ = ψ (r ) et r =  P W 
 φE − φP 
Nous supposons que le flux est positif x-direction, donc u> 0, en développant le concept
de TVD en tant que généralisation de expressions de la quantité transportée φe d'un
volume de contrôle unidimensionnel:
1 Eq. 4-26
φ e = φ P + ψ (r )(φ E − φ p )
2
En remplaçant φ0 par sa valeur pour chaque schéma
schémas ψ (r ) φe =ψ (r )
schéma Upwind(UD) ψ (r ) = 0 φe = φP
1
Schema Centre (CD) ψ (r ) = 1 φe = φ P + (φ E − φ p )
2

1 φ −φ 
Schema Upwind linear(LUD) ψ (r ) = r φe = φP +  P W (φE − φ p )
2  φE − φ p 

1
schéma (QUICK) ψ (r ) = (3 + r ) / 4 φe = φ P + (3φ E − 2φ p − φW )
8
- Implementation des schemes TVD

l'équation 4-1 pour un problème à une seule dimension (1D), donne:
∂ (ρu φ ) ∂  ∂φ  Eq. 4-27
= Γ +S
∂x ∂x  ∂x 
39
Pour un écoulement à une seule dimension (1D), si u positive dans la direction x (u>
0) et les valeurs de φe et φ w à l'aide d'un schéma TVD peut être écrit :
1 Eq. 4-28
φ e = φ P + ψ (re )(φ E − φ p )
2
1 Eq. 4-28-1
φ e = φ P + ψ (rw )(φ P − φW )
2
 φ −φ   φ − φWW 
Où re =  P W  et rw =  W 
φ −
 E pφ  φ P − φW 
En remplaçant dans l'équation 4-27

1 1 Eq. 4-29
Fe φP + ψ (re )(φE − φ p ) − Fw φW + ψ (rw )(φP − φW ) = De (φE − φP ) − Dw (φP − φW )
   
 2   2 
et après arrangement
Eq. 4-30
[De + Fe + Dw ]φ p = [Dw + Fw ]φW + DeφE − Fe  1ψ (re )(φE − φ p ) + Fw  1ψ (rw )(φP − φW )
2  2 
aPφP = aWφW + aEφE + SuDC Eq. 4-31
La forme compacte
aW = Dw + Fw Eq. 4-31-1

Où aE = De Eq. 4-31-2
a = a + a + (F − F ) Eq. 4-31-3
 P W E e w
1 1
SuDC = −Fe  ψ (re )(φE − φ p ) + Fw  ψ (rw )(φP − φW )
 
Eq. 4-31-4
2  2 
Pour noter la direction du flux, nous utilisons l' indice ‘+’ qui représente le sens positif de
l'écoulement comme rw+ ; re+ et l'indice ‘−’ pour le sens négatif.
En remplaçant le terme source par cette notation:
1 1 Eq. 4-32
( ) ( )
SuDC = −Fe  ψ re+ (φE − φ p ) + Fw  ψ rw+ (φP − φW )
 
2  2 
Le tableau récapitulatif des coefficients de TVD
aW aE SuDC
Dw + max(Fw,0) De + max(− Fe ,0) 1

[ ( ) ( )]
Fe (1 − α e )ψ re− − α e .ψ re+ (φ E − φ p )
2
1
[ ( ) ( )]
+ Fw α w .ψ rw+ − (1 − α w )ψ rw− (φ P − φW )
2
40
Où
Eq. 4-34-1
α w = 1 pour Fw > 0 et α e = 1 pour Fe > 0
α w = 0 pour Fw < 0 et α e = 0 pour Fe < 0

Eq. 4-34-2
-Condition aux limites pour les schémas TVD
En supposant une valeur donnée en point A àl' entrée φ = φA et une quantité convective
du flux par unité de surface: F = F A . L'équation discrétisée des TVD est
1 Eq. 4-35-1
Fe φP + ψ (re )(φE − φ p ) − FAφA = De (φE − φP ) − D*A (φP − φA )
 
 2 
1 Eq. 5-34-2
(D e ) ( )
+ Fe + D A* φ P = Deφ E + D A* + FA φ A − Fe ψ (re )(φ E − φ p )
2
Γ r =  φP − φW 
 et r =  φW − φWW 
Où D A* = ; e   w  φ −φ 
δx  φE − φ p   P W 
pour le terme de correction l'extrapolation des nœuds miroirs φ = φW
 φ −φ  2(φ − φ ) Eq. 4-35

φ0 = 2φA −φP et re =  P 0  = P A
 φE − φ p  φE − φ p
5.6 Prise en considération du gradient de pression

Normalement le terme source SU dans les équations du mouvement contient un gradient
de pression, l'équation 4-2 s'écrit donc:
d (ρu φ ) dp d  dφ  Eq. 4-36

=− + Γ +S
dx dx dx  dx 
Malheureusement la pression n'apparaît pas dans l'équation de continuité. Pire encore
pour un fluide incompressible où la masse volumique est constante, l'équation de
continuité qui traduit le principe de conservation de la masse se trouve complètement
découplée des équations du mouvement. Pour φ = 1 ; Γ = 0 et S = 0 dans l'équation 4-36,
on obtient l'équation de continuité.
dui Eq. 4-37
=0
dxi
et pour un problème à une seule dimension:
41
du Eq. 4-37-1
=0
dx
ce qui donne:
ue − uw = 0 Eq. 4-37-2
4.6.1 Le problème du damier

Supposons qu'on veut intégrer l'équation 4-36 par rapport au volume de contrôle de
centre P. Le gradient de la pression sera discrétisé comme suit:
dp  p − pe Eq. 4-38
−  = w
dx  p ∆x
et en utilisant une interpolation linéaire:
dp  1  pW − pP pP − pE  pW − pE Eq. 4-38-1
−  =  − =
dx  p ∆x  2 2  2∆x
Où
Figure 4-2 :champ de pression non uniforme
Ce qui fait que l'information de la pression au point P a simplement disparue de notre

équation. En conséquence, un champ de pression non uniforme de type (50, 100, 50,
100, 50, …). Figure 4-2 comme étant un champ uniforme. L'addition d'un tel champ à la
solution exacte sera aussi solution des équations discrétisées. C'est le problème très
connue sous le nom du problème du damier (Maillage décalé)
Figure 4-3 : Maillage décalé
42
Pour remédier au problème du damier, il convient d’utiliser plusieurs grilles de calcul
décalées. Une première grille passe par les centres des volumes de contrôle. Dans cette
grille sera stocké la pression et toutes les autres quantités scalaires (température,
concentration, k, e…). Ensuite d’autres grilles seront construites au niveau des facettes
des volumes de contrôle où seront stockées les composantes de la vitesse.
4.6.2 Algorithm de Patankar
L'intégration de l'équation 4.36 sur le volume de contrôle de centre e et de limites P et

E donne:
aeue =∑anbunb +b +Ae (pP−pE ) Eq. 4.39
Soit un champ de pression initial p* . La solution provisoire de l'équation précédente

sera noté u* (notons que u* ne vérifie pas l'équation de continuité).
aeuee* =P∑ anbu* + b + Ae ( pP* − p*E )

E nb
Eq. 4.40
A ce stade, aucune des deux variables n'est correcte. Toutes les deux nécessitant une
correction.
u = u *+ u ' Eq. 4.40.1
p = p * + p' Eq. 4.40.2
Où u 'et p ' sont les corrections qu'il faut estimer.

L'introduction des équations 4-40-1 et 4-40-2 dans l'équation 4-36 et en tenant compte
de l'équations 4-40, il s'en suit:
(
ue = ue* + de pP' − pE' ) Eq. 4.41
Où
Ae Eq. 4-41-1
de =
ae
43
Notons ici que pour linéariser l’équation, le terme ∑a '

u
nb nb a été tout simplement
négligé.
Normalement, ce terme doit s'annuler lors de la convergence de la procédure. C'est à
dire que cette omission n'influe par sur le résultat final, mais elle fausse un peu le
résultat temporaire. C'est d'ailleurs la seule simplification faite dans l'algorithme
SIMPLE. Elle a été corrigée dans les variantes plus évoluées (SIMPLER et SIMPLEC).
L'introduction de l'expression corrigée de l'équations 4-41 dans l'équation de

continuité 4-40, donne l'équation de correction de la pression, qu'on écrira sous la
forme suivante:
Eq. 4-42
' ' '
aP p = aE p + aW p + b
p E W
Où
Eq. 4-42-1
aE = (d A )e
Eq. 4-42-2
aW = (d A )w
Eq. 4-42-3
aP = aE + aW
Eq. 4-42-4
b= u A( *
) + (u A )
w
*
e
D’après l'équation 4-42, le terme b représente le terme source de masse présent à cause
du champ de pression aléatoire initial. Normalement, l'algorithme de résolution doit
annuler ce terme.
Enfin, l'algorithme SIMPLE sera résumé comme suit:
1. Choisir un champ de pression initial p* .
2. Résoudre les équations de transport 4-39 pour déduire un champ de vitesse u*.
3. Calculer le terme source de la masse b de l'équation 4-42 et résoudre l'équation 4-
40 de correction de la pression.
4. Corriger les champs de pression et de vitesse via les équations 4-40-1 et 4-40-2.
5. Résoudre les autres équations de transports d'autres scalaires du problème, tel
que la température ou les quantités turbulentes
44
6. Remplacer l'ancien champ de pression par le nouveau et revenir à l'étape 2.
Répéter les calculs jusqu'à convergence de toutes les variables.
Comme il a été mentionné plus haut la simplification du terme ∑a '

u n'affecte en rien
nb nb
la solution finale, puisque si la convergence est atteinte ce terme devra s'annuler.
Toute fois le taux de convergence est modifié par cette simplification. Il se trouve que la
correction p' est surestimée par SIMPLE et le calcul a tendance à diverger. Le remède
pour stabiliser les calculs est d'utiliser un coefficient de sous relaxation aussi bien pour
les équations du moment que celle de la pression. L'équation 4-40-2 devient:
p = p* + α p p' Eq. 4-43
4.6.3 L'algorithme SIMPLER (SIMPLE-Revised)
Dans la version révisée de SIMPLE, seule la correction de la pression p' est utilisée pour
corriger la vitesse. Une autre équation est utilisée pour estimer la nouvelle pression.
L'équation 4-41 est reécrite:
^
ue = ue + de ( pP − pE ) Eq. 4-44
Où uˆe est une pseudo vitesse définit par:
^
u= ∑a u nb + b
nb Eq. 4-45
ae
La ressemblance marquante entre l'équation 4-41 et l'équation 4-42-44, nous permet

d'injecter cette dernière dans l'équation 4-37-2 comme on l'a fait pour l'équation 4-41.
Le résultat sera une équation similaire à l'équation 4-42, mais en p et non en p'.
aP pP=aE pE +aW pW +b Eq. 4-46
Où les coefficients aE , aW et aP sont donnée par l'équation 4-42-1, l'équation 4-42-2) et

l'équation 4-42-3
45
 ^   ^  Eq. 4-47
b =  uA −  uA 
  w  e
Ici encore b est considéré comme un terme source de la masse qui doit s'annuler.
Enfin, l'algorithme SIMPLER sera résumé comme suit:
1. Choisir un champ de vitesse initial.
2. Calculer les coefficients des équations de quantité de mouvement et déduire uˆ à
partir de l'équation 4-45.
3. Evaluer le terme source de la masse b de l'équation 5-47 et résoudre l'équation de la
pression 4-46.
4. Utiliser le champ de pression pour résoudre les équations de quantité de mouvement
uˆ que l'équation 4-40 pour déduire un champ de vitesse u* .
5. Calculer le terme source de la masse b de l'équation 4-42-4 et résoudre l'équation
4-42 de correction de la pression.
6. Corriger le champ de la vitesse via l’équation 4-41 (ne pas corriger la pression).
7. Résoudre les autres équations de transports d'autres scalaires du problème, tel que la
température ou les quantités turbulentes
7. Revenir à l'étape 2. Répéter les calculs jusqu'à convergence de toutes les
variables.
La supériorité de l'algorithme SIMPLER par rapport à SIMPLE réside dans le fait que la
déduction de l'équation de la pression 4-46 ne fait intervenir aucune simplification.
Dans SIMPLE, la déduction de l'équation de correction de la pression l'équation 4-42
passe par l'annulation du terme ∑ *
a nb u nb .Par conséquent le champ de pression dans
SIMPLER est plus proche de la réalité que celui de SIMPLE, puisqu'en général
l'estimation d'un champ de vitesse initial est plus facile que celle d'un champ de
pression. Notons, ici que l'algorithme SIMPLER ne nécessite pas de champ de pression
initial. La pression est directement générée à partir de l'initialisation de la vitesse. Par
conséquent des coefficients de sous relaxation plus consistants peuvent être utilisée
pour les vitesses. Mieux encore, aucune sous relaxation n'est nécessaire pour la
pression. Il est vrai qu'une itération suivant l'algorithme SIMPLER nécessite environ
30% de temps plus que celle de SIMPLE, mais cet effort est largement compensé par la
réduction consistante en nombre d'itérations nécessaires pour la convergence.
46
4.6.4 L'algorithme SIMPLEC (SIMPLE-Consistent)
Cette amélioration a été apportée par Van Doormal et Raithby, (1984) pour rectifier la
négligence du terme ∑ *
a nb u nb ' dans SIMPLE.
La correction de la vitesse obtenue précédemment par les équations 4-41 et 4.41-1 sera
écrite:
(
ue = ue* + de pP' − pE' ) Eq. 4-48
Où
Ae Eq. 4-49
de =
ae − ∑ anb
Au lieu de faire comme dans SIMPLE et négliger complètement le terme ∑ *

a nb u nb ' , on
préfère garder la partie connue et négliger seulement ce qui est inconnu. Ce qui donne la
formulation de l’équation 4-44. Les étapes de SIMPLEC restent les mêmes que ceux de
SIMPLE.
L'algorithme PISO (Pressure Implicit with Splitting of Operators of Issa, 1986)

Cet algorithme a été développé initialement pour les calculs non itératifs des
écoulements non stationnaires et compressibles. Il a été ensuite adapté avec succès pour
les calculs itératifs des problèmes stationnaires. L’algorithme est similaire à SIMPLE
avec une amélioration qui consiste à faire deux corrections successives au lieu d’une
seule.
Etape de prédiction:
En utilisant un champ de pression initial p* qu’on utilise pour résoudre l'équation de
transport 4-40. Un champ de vitesse u* est déduit.
Première étape de correction:

Le champ de vitesse u* obtenu ne vérifie pas l'équation de continuité, sauf si le champ de
pression p* est correct. Les mêmes étapes de SIMPLE sont suivit pour obtenir une
47
première correction du champ de la vitesse. Puisque dans PISO, on doit faire deux
corrections successives, le champ de vitesse obtenue sera noté u** .
Eq. 4-50
** * '
u =u +u
Eq. 4-51
p** = p* + p'
Eq. 4-52
ue** = ue* + de ( pP' − pE' )
Deuxième étape de correction:
L'équation 4-40 pour u **
aeue** = ∑anbunb
*
+ b + Ae pP** − pE**( ) Eq. 4-53
Un deuxième champ corrigé de vitesse sera noté u*** et calculé par
aeue*** = ∑anbunb
**
(
+ b + Ae pP*** − pE*** ) Eq. 4-54
La sommation dans l'équation précédente est faite avec les vitesses issues da la
correction précédente. La soustraction de l'équation 4-53 de 4-54 donne:
***
u =u **
+
∑a nb
**
(unb *
− unb )
(
+ de pP'' − pE'' ) Eq. 4-55
e e
ae
où p'' est une deuxième correction de la pression:
p*** = p** + p'' Eq. 4-56
L’injection de u*** dans l'équation de continuité donne une deuxième équation de

correction de la pression:
aP p'p' = aE pE'' + aW pW'' + b'' Eq. 4-57
48
où:
les coefficients aE , aW et aP sont donnée par les équations 4-42-1, 4-42-2 et 4-42-3.
Aw A
b'' =
aW
∑ **
anb(unb *
− unb ) − e ∑anb(unb
aE
** *
− unb ) Eq. 4-58
Enfin, l'équation 4-56 sera écrite:
p*** = p** + p'' = p* + p' + p'' Eq. 4-59
D'après plusieurs études comparatives, il semblerait que SIMPLER est le meilleur

algorithme de couplage, puisqu'il ne nécessite pas de pression initiale. Il paraît aussi
qu'il est de 30 à 50% plus rapide que SIMPLE (Anderson et al, 1984). Par contre la
comparaison entre SIMPLEC et PISO dépend des conditions de l'écoulement considéré,
du degré de couplage entre les équations du moment et celles des autres scalaires
(combustion par exemple). Une étude comparative entre PISO, SIMPLER et SIMPLEC
pour des problèmes stationnaires non couplés à des scalaires a été faite par Jang et al
(1986). La conclusion est que l'algorithme PISO était le plus stable et le plus rapide. Pour
les problèmes où le scalaire est très lié aux équations de la vitesse l'avantage de PISO
n'est pas significatif.
49
Chapitre 5: Solution des équations

algébriques discrétisées.
5.1 Système matriciel: Nous avons vu dans le chapitre précédent que la discrétisation
d'une équation de transport permet d'établir, pour tout point P du maillage, une
équation algébrique. Le forme de cette équation est semblable, mais légèrement
différente selon que le transport de l'extensité ϕ ait lieu en régime stationnaire ou
transitoire. Rappelons les différentes formes possibles de l'équation algébrique, qui
restent valables 'aussi bien si le volume de contrôle associé au nœud P est une cellule
interne que s'il s'agit d'une cellule frontière:
-régime transitoire, schéma explicite (à chaque pas de temps) :
aPφP = ∑ anb
0 0
φnb + b Eq. 5-1
nb
-régime permanent ou régime transitoire, schéma implicite (de chaque pas de temps)
aPφP = ∑ anbφnb + b Eq. 5-2

nb
Nous allons immédiatement écarter de la discussion le cas du schéma explicite, qui ne pose,
du point de vue de la résolution numérique, aucun problème particulier, à partir du moment
où le critère de stabilité aP ≥ ∑a nb est vérifié en tout point. En effet; l'équation 5-l ne
comporte à droite aucune inconnue et le calcul est direct (explicite, en somme). Dans les deux
autres cas (régime permanent et schéma implicite), la présence d'inconnues φnb dans les
termes de droite impose une résolution matricielle. Celle-ci doit être effectuée à chaque pas
de temps en régime transitoire .
Il peut paraître étrange que, du point de vue numérique, la résolution d'une équation de
transport en régime permanent pose plus de difficultés qu'en régime transitoire
(schéma explicite).En fait, la nature parabolique de l'équation aux dérivée partielles par
rapport au temps rend nécessaire une condition initiale de Dirichlet, c'est-à-dire la donnée
du champ d'extensité ϕ en tout point à 1' instant t =0. Cette connaissance exacte de l'état du
50
système à un moment donné permet de traiter avec beaucoup de précision et de rapidité

les problèmes en régime transitoire. Ceci est vrai à un tel point qu'une résolution en
régime permanent fréquemment à peine moins coûteuse en temps de calcule qu'une
résolution complète en régime transitoire jusqu'à obtention de l'état quasi-stationnaire.
Comme nous le verrons plus loin. la résolution par itérations successives de l'équation de
transport en régime permanent a d'ailleurs beaucoup de points communs, numériquement
parlant, avec un calcul direct effectué à des “ pas de temps ” fictifs, qui sont en fait des “ pas
de calcul ”
Hormis le cas du schéma explicite (et des contraintes qui y sont associées pour assurer la
stabilité du calcul), il est donc nécessaire de résoudre une équation matricielle de la forme :
(φ ) = [A] (φ ) + (B) Eq. 5-3
[ A ] est une matrice carrée, de dimension N ou N est le nombre de points du maillage, qui
contient les termes anb a p . Un problème tridimensionnel discrétisé a l'aide d'un maillage
relativement grossier (20 x 20 x 20 points) conduit donc à une matrice de 6000 x 8000, soit
64.000.000 termes! Ces quelques chiffres montrent bien l'importance de la méthode de
résolution choisie en termes d'encombrement mémoire et de temps de calcul. Remarquons
toutefois que la matrice [ A ] est très creuse (pleine de "0"), ce qui facilite grandement son
stockage.
La résolution de l'équation matricielle 5-3 est basée, soit sur une méthode directe, soit sur
une méthode itérative.
5.2 Les méthodes de résolution directe

L'équation matricielle 5-3) peut se réécrire sous la forme :
([I ] − [A]) + (φ ) = (B ) Eq. 5-4
où[I]est la matrice identité de dimension N², N étant le nombre de points du maillage.

Résoudre le système 5-4 revient pratiquement: à inverser la matrice([I]-[A]) et calculer
exactement le champ des inconnues (ф) par:
(φ ) = ([I ] − [A])−1(B) Eq. 5-5
En fait, ceci ne " marche " rigoureusement que si les coefficients a p et anb qui forment la
matrice [A] , sont indépendants de φ c'est-à-dire si le système est réellement linéaire.

51
Dans le cas contraire, quelques itérations sont théoriquement nécessaires, en actualisant
les coefficients a partir des valeurs nouvellement calculées des φ p .
L'inversion de la matrice à partir du calcul de la matrice des cofacteurs est bien entendu
une méthode a prohiber absolument ! On a pu évaluer qu'elle n'était " rentable ”, en
terme de temps de calcul, que pour des valeurs de N inférieures ou égales à 5. On est loin
de la valeur 5000 précédemment évoquée ! les seules méthodes directes classiques
susceptibles de permettre la résolution de l'équation 5-4 sans impliquer de temps de
calcul prohibitif sont les techniques d'élimination de Gauss ou " méthodes du pivot ”. Il
en existe un certain nombre de variantes (pivot maximum, méthode de Gauss-Jordan,
réduction de Crout, méthode de Cholesky,...). Il n'est pas du ressort de ce cours (ni
d'ailleurs de la compétence des auteurs !) de disserter plus avant sur les différentes
méthodes d'élimination de Gauss, qui sont traitées par ailleurs dans la plupart des
ouvrages de référence d'analyse numérique. En effet, bien que la matrice ([ I ] - ( A ]) soit
très creuse, ses dimensions, dès que le maillage n'est pas particulièrement grossier,
rendent ces méthodes très lourdes d'utilisation. Pour fixer les idées, notons qu'une
méthode itérative est généralement préférable lorsque le nombre de nœuds du maillage
dépasse environ 200, ce qui est extrêmement faible si on a affaire à un problème bi- ou
tri-dimensionne1. En outre, la nécessité d'actualiser les coefficients dépendants de φ ,
c'est-a-dire les non-linéarités, (nécessité dont beaucoup de numériciens se moquent
avec une tranquillité d'esprit surprenante !) joue également contre l'emploi des
méthodes directes. Si, de toute manière, il faut itérer à un moment, pourquoi ne pas le
faire systématiquement ? On n'utilisera dans la pratique la résolution directe que dans
deux cas:
a) Transfert 1-D: le T.DM.A.
Comte tenu de la géométrie, la matrice du système linéaire à inverser est tridiagonale

(seules les 3 diagonales principales sont non nulles)
De fait, la méthode du pivot prend une forme simplifiée, appelée algorithme de Thomas,
ou double balayage de Cholesky, ou plus couramment TDMA (Tri Diagonal Matrix
Algorithm).
52
Pour une présentation plus facile du TDMA, nous utilisons ici une notation différente :
les points sont numérotés de 1 à N dans le sens des x croissants. La forme générale d'une
équation discrétisée prend alors la forme :
&' (' = )' ('*+ + ,' ('-+ + .' pour i allant de 1 à N Eq. 3-18
b
Pour les points frontières, nous avonsnécessairementc1= 0et N= 0Lesystèmelinéaire est
résolu par la méthode suivante:
1º étape : initialisation de la récurrence pour l'étape de descente
)1 .1
/1 = 11 =
&1 &1
2º étape : descente, calcul des Pi et Qi par l es relations de récurrence
)2
/2 =
&2 − ,2 /2−1
.2 + ,2 12−1 i allant de 2 à N
12 =
&2 − ,2 /2−1
3º étape : initialisation on de la récurrence pour l'étape de remontée
(3 = 13
4º étape : remontée, les i sont déterminés par les relations de récurrence
(2 = /2 (2+1 + 12 i allant de N-1 à 1
Cette méthode de résolution requiert:
• 3 (N- 1 ) additions (on ne calcule qu'une fois &' − ,' /'-+ ),

• 3 (N- 1 ) multiplications,
• 2N divisions
soit (8N - 6) opérations au total, ce qui constitue une réduction considérable du nombre
d'opérations par rapport à une méthode du pivot classique où le nombre d'opérations
est de l'ordre de 2N3/3.
Notonsquecetalgorithme,extrêmementperformantentermedetempsdecalcul,permetde
ne jamais stocker en tant que telle la matrice ([I ] - [ A ]). Celle-ci est remplacée par trois
vecteur semi-colonnes qui correspondent aux trois diagonales non nulles. Remarquons
une fois de plus que le TDMA ne permet de résoudre que des systèmes linéaires. Si les
coefficients a4 , b4 , c4 , d4 dépendent de ф, il faut alors utiliser une procédure itérative de la
forme:
53
1) estimer la valeur de Φ en tout point du maillage
2) calculer les coefficients ai, bi, ci , di avec les dernières valeurs estimées
3) appliquer le T.D.M.A. pour résoudre le système
4) avec les nouvelles valeurs de Φ, retourner à l’étape 2 jusqu'à convergence.
b)Transport 2-D ou 3-D en régime transitoire: le schéma A.D.I.
Le schéma A.D.I. (Alternating Direction Implicit), ou schéma aux directions alternées, est
surtout utilisé dans le cas d'une résolution en régime transitoire, pour un problème 2-D
(ou 3-D). La méthode est basée sur la division du pas de temps en 2 ( ou3). Dans un
premier temps dt / 2 (ou dt / 3), le schéma est totalement implicite dans une direction
privilégiée et explicite dans l'autre. A chaque 1 / 2 (ou 1 / 3) pas de temps, la “ direction
implicite " est changée. Cette alternance rend le schéma A.D.I., comme le schéma
implicite, inconditionnellement stable pour un système d'équations totalement linéaire.
Les matrices (B) et ([ I ] - ( A )) sont réordonnées à chaque alternance de manière que
([ I ] - [ A ]) reste toujours tri-diagonale. Le T.D.M.A. est alors appliqué 2 (ou 3) fois par
pas de temps.
5.3 Les méthodes de résolution itérative
Dans le cadre de la méthode des volumes finis, qui permet de traiter des
phénomènes de transport souvent comptés et fortement non-linéaires, notre préférence
va nettement aux méthodes itératives. Nous appliquerons systématiquement une
méthode itérative
-dans tous les problèmes 3-D,
-dans tous les problèmes fortement non-1inéaires (transport convectif de quantité de

mouvement en particulier),
-dans les problèmes 2-D en régime permanent (a moins d'un maillage très grossier
méthode d'élimination de Gauss) ou en régime transitoire (la concurrence avec le
schéma A.D.I. est alors très forte).
54
5.3.1 méthode de Jacobi
La méthode de Jacobi, également appelée méthode de Liebmann ou méthode de

déplacements simultanés, est la plus “ évidente ” des méthodes itératives. A partir de
l'équation matricielle 5-3, qui est en fait la représentation condensée d'un système
d'équations algébriques de forme générale :
anb b Eq. 5.6
φP = ∑ φnb +
nb ap ap
on procède par itérations successives, en commençant par une valeur
"supposée" du résultat. Que nous noterons (φ )0 . Dans le cas d'un problème en régime
permanent, seules des considérations physiques permettent de choisir un champs
d'extensité (φ )0 pas trop éloigné de la solution réelle, afin de réduire le nombre total
d'itérations nécessaires. Lorsque les considérations physiques sont inexistantes, ou que
le numéricien est très paresseux, un champ (φ)0 aussi artificiel que :
Pour tout P
φP,0 = 0 Eq. 5-7
est ou vent utilisé. Dans le cas d'un problème en régime transitoire, le choix le plus
évident, et le meilleur, consiste à utiliser, comme point de départ des itérations au pas de
temps N, la valeur calculée de (φ ) au pas de temps N- 1
(φ )0N = (φ )0N −1 Eq. 6-8
La suite matricielle (φ)0 , (φ)1 = [A](φ)0 + (B)......(φ)i = [A](φ)i−1 + (B)....est calculée jusqu'à ce
qu'un critère de convergence (voir plus loin) ait été atteint. La valeur courante du champ
(φ )i à cette itération est alors le résultat de l'équation matricielle.

Cette méthode implique de garder à tout moment en mémoire deux vecteurs uni
colonnes pour (φ ) , pour l'itération courante et l'itération précédente. Bien entendu, a
chaque itération, les matrices [A]et [B], c'est-à-dire les coefficients ap , a nb et b de

l'équation (5-6), sont actualisées dans le cas de problèmes non-linéaire.
Nous allons a présent examiner d'un peu plus près le problème de transports couplée
(chaleur et masse, pour fixer les idées). Si (φ ) et (φ ) sont les extensité spécifiques
'
transférées, nous avons affaire a deux équations matricielles :
55
(φ ) = [A] (φ ) + (B) Eq. 5-9
(φ ) = [A ] (φ ) + (B )
' ' ' ' Eq. 5-10
Le couplage peut s'exprime par le fait que [ A] dépende de (φ ), par exemple,

'
ou que
(B') défende de (φ ) , ou par toute relation de ce genre.
La méthode de Jacobi s'écrira séquentiellement de la manière suivante:
(φ)1 = [A]0 (φ)0 + (B)0

(φ ) = [A ] (φ ) + (B )
'
1
'
0
'
0
'
0
.
. Eq. 5-11
.
(φ)i = [A]i−1(φ)i−1 + (B)i−1
(φ ) = [A ] (φ ) + (B )
'
i
'
i−1
'
i−1
'
i−1
jusqu'à convergence simultanée des: deux équations. L'ordre de traitement des

équations n'est pas indifférent, et peut jouer un rôle important en ce qui concerne la
vitesse de convergence du processus global.
La méthode de Jacobi est une méthode particulièrement facile à comprendre. Son intérêt
est en fait surtout pédagogique car sa vitesse de convergence est généralement très
faible. Elle a été totalement supplantée dans la pratique par la méthode voisine de
Gauss-Seidel.
5.3.2 La méthode de Gauss-Seidel
Cette méthode, aussi appelée méthode des déplacements successifs ou méthode point-
par-point, se présente comme une variante de la précédente. De la même manière, on
part d'une valeur "supposée" pour l'ensemble du champ d'extensité (φ )0 , puis , sans
garder deux vecteurs en mémoire, on calcule, en “ balayant ” le maillage tout entier, les
valeurs de (φ )1 par
anb * b Eq. 5-12
φP ,1 = ∑ φnb +
nb ap ap
Si le point voisin nb a déjà été traité lors du balayage, on a φnb
*
= φnb,1 . Dans le cas
contraire, φnb
*
= φnb,0 . De toute manière, on utilise la valeur de φ nb disponible (il n'y en a
56
qu'une) au moment du calcul. Une itération consiste en un balayage complet du maillage.

Comme pour la méthode de Jacobi, on répète les itérations jusqu'à convergence. La
résolution d'équations couplées s'effectue de la même manière que pour la méthode de
Jacobi.
La méthode de Gauss-Seidel sévère en général bien plus rapide que celle de

Jacobi. En outre, elle demande moins de stockage-mémoire. De plus, la confusion
mentale entre itérations et pas de temps (Ces deux notions sont tout de même
différentes. Une méthode itérative peut parfaitement utilisée en régime permanent) est
moins susceptible de se produire que lors de l'utilisation de la méthode de Jacobi.
Cette méthode réunit à première vue toutes les qualités. Toutefois, elle reste
encore souvent trop lente, et de plus sujette à divergence (voir le paragraphe suivant)
lorsque les non- linéarités sont trop fortes. C'est pourquoi on l'emploie peu sous sa
forme initiale. La méthode de Gauss-Seidel est généralement utilisée en association avec
des techniques de relaxation ou de traitement par blocs (voir plus loin).
5.3.3 Convergence et divergence
La notion de convergence d'un schéma itératif est intuitive et correspond

mathématiquement à la propriété selon laquelle toute suite de Cauchy converge dans un
espace de Banach. De manière plus pratique, cela signifie que, si la différence entre les
champs d'extensité (φ )i−1 et (φ )i aux itérations i-1 et i est inférieure à une certaine
limite, on peut considérer qu'ils sont égaux ct qu'ils vérifient par conséquent l'équation
matricielle 5-3). La définition du critère de convergence adéquat représente un
problème important. La première idée est de calculer une grandeur telle que
Max p φ p ,i − φ p ;i −1 et de poser que, pour que la convergence soit atteinte, cette
grandeur doit être inférieure à un nombre suffisamment petit ε :
Max p φ p ,i − φ p ;i −1 ≤ ε
Eq. 5-13
L'utilisation d'un tel critère de convergence possède un grave inconvénient : si,

en un point, la convergence est lente, alors qu'elle est très rapide pour presque tout le
domaine, le calcul itératif général continuera à s'effectuer alors même qu'on pourrait
souhaiter s'arrêter. Ce cas est extrêmement courant. L'utilisation de critères tels que :
57
∑ p
φ p ,i − φ p;i −1 ≤ ε Eq. 5-14-1
Eq. 5-14-2
∑V
p
p φ p ,i − φ p ;i−1 ≤ ε
Où VP est le volume de contrôle associé au nœud P, est par conséquent bien plus
satisfaisante.
Toutefois, e n'est pas un nombre adimensionnel. L'ordre de grandeur des valeurs
de φ va par conséquent jouer un rôle anormal : A la limite, si toutes les valeurs des φp
sont très petites, parce qu'expirées dans des unités inadéquates, le critère sera satisfait
dès la première itération, on en vient donc à l'utilisation de critères normés tels que
∑ p
φ p ,i − φ p ;i −1
Eq. 5-15-1
≤ε
∑
p
φ p ,i
∑Vp
P φ p ,i − φ p;i −1
Eq. 5-15-2
≤ε
∑V
p
P φ p,i
De tels critères sont très couraient utilisés. Une technique alternative, que nous
conseillerons malgré sa moins grande facilité de programmation, n’utilise pas la
notion de suite de Cauchy. Elle consiste a calculer, pour chaque point et chaque
itération (et pour chaque variable dans le cas transports couplés, bien en tendu) le
résidu de l'équation RP par:

 a b 
RP,i = φ p,i −  ∑ nb φnb +  Eq. 5-15
 nb a a p 
 p
Le vecteur (R)i représente en fait "l'écart à l'exactitude (φ ) "du champ (φ )i .Des

critères analogues à 5-13, 5-14 et 5-15 peuvent alors être définis simplement en
remplaçant la quantité φ p ,i − φ p ;i −1 par RP,i . Le choix de ε reste bien entendu dans tous les cas
a la discrétion du numéricien.
A ce stade, il est bon de se poser une question capitale un schéma itératif est-il
automatiquement convergent ? Evidemment, la réponse est non ! Nous n'aurons pas
cette chance !
Dans le cas d'un système linéaire (les coefficients a et b sont totalement

indépendants des valeurs de φ ), on peut démontrer mathématiquement, et nous l'admettrons
58
sans le discuter, qu'une condition suffisante pour qu'un schéma itératif soit convergent est,
avec les notations de l'équation 5-6:
anb
∑
nb ap
≤ 1pour tout point P
Eq. 5-17-1
anb
∑
nb ap
〈 1 pour au moins un point P
Eq. 5-17-2
Cette condition, appelée critère de Scarborough, signifie que la matrice ([ I ] - [ A ])

de l'équation 5-4 est à diagonale dominante. Il est particulièrement intéressant de'
remarquer que, si la discrétisation par volumes finis a été effectuée correctement, le critère
de Scarborough est automatiquement vérifié ! Une fois de plus, l'équivalence entre principe
physique et condition mathématique se vérifie, comme dans le cas de l'analyse de la stabilité
conditionnelle du schéma explicite.
Le critère de Scarborough n'est toutefois pas suffisant dans le cas d'équations non-
linéaires ou couplées. Dans de telles circonstances, il n'existe pas de critère simple
permettant d'assurer la convergence d'un schéma itératif. Si le schéma manifeste une:
tendance à la divergence (la suite (φ )i ne se comporte pas comme une suite de Cauchy), il
est nécessaire de sous-relaxer l'équation matricielle.
5.4 Notion de relaxation

La relaxation est une technique fréquemment employée pour accélérer la
convergence (c'est-à-dire réduire le nombre d'itérations nécessaire) d'un schéma itératif
ayant naturellement tendance à converger, aussi bien que pour empêcher la divergence d'un
système “rebelle”. Les méthodes les plus classique sont l'utilisation d'un facteur multiplicatif
de relaxation, et la relaxation par faux pas de temps.
5.4.1 Facteur de relaxation
La technique consiste simplement, à chaque itératif i, à remplacer la valeur qu'on

vient de calculer φ p,i par une quantité :
φ ' p ,i = α φ p ,i + (1 − α ) φ p ;i −1 Eq. 5-18
Et d'utiliser φ p,i à la place de φ p,i à partir de cette itération (méthode de Jacobi) ou

'
immédiatement a partir de ce calcul (méthode de Gauss-Seidel).
59
Dans l'équation 5-18, α est le facteur de relaxation. Plusieurs commentaires

peuvent être formulés à ce stade
• α =1 correspond à une relaxation sans effet : φ ' p ,i = α φ p ,i
=0 signifie que le processus itératif est Stoppé: φ = φ p,i−1 .Si un tel schéma
'
• α p ,i
ne risque pas de diverger, il ne va pas non plus permettre d'avancer beaucoup

dans la résolution.
• 0 〈α 〈1 correspondàlanotiondesous-relaxation.Leprocessusitératifestfreinéafin de limiter
les risques de divergence. Des valeurs aussi faibles que α = 0,01 doivent parfois être
utilisées dans le cas d'équations non-linéaires fortement couplées (turbulence, par
exemple).
• l < α <2 correspond à la notion de sur-relaxation. Le schéma itératif est cette fois "accéléré”
afin d'obtenir une convergence plus rapide. Très efficace dans le cas d'équations linéaires ou
presque linéaires, la sur-relaxation ne doit être employée qu'avec précaution du fait du
risque de divergence du processus.
Par ailleurs, une valeur supérieure à 2 conduit inéluctablement a la divergence des
calculs. Remarquons que, bien entendu, quelle que soit la valeur du facteur de
relaxation, les vecteurs successifs (φ ), si

'
convergent, convergent bien vers la solution
de l'équation matricielle 5-3.
Il n'existe pas de règle absolue permettant de choisir un facteur de relaxation

optimal Certaines équations doivent impérativement être violemment sous-relaxées.
D'autres convergent plus vite si une légère sur-relaxation est introduite. Dans le cas de
transports couplée, chaque équation peut évidemment être traitée avec un facteur de
relaxation différent. De même, il n'est pas nécessaire de garder le même facteur au
cours des itérations. Il est souvent intéressant de visualiser l'évolution d'un paramètre
de convergence tel que Rp,i l'équation 5-16 et d'ajuster le facteur de relaxation au mieux
lors du processus itératif. Cela implique que le programme "rende périodiquement la

main" à l'utilisateur.
L'utilisation d'un facteur de relaxation couplée à la méthode itérative de Gauss-

Seidel est probablement une des plus courantes des techniques de résolution numérique
des équations convecto-diffusives discrétisées par la méthode des volumes finis. Ce
60
schéma a reçu le nom, une fois de plus anglo-saxon, de méthode S.O.R. (Successive Over-
Relaxation). Nous le conseillons fortement dans le cadre de ce cours. Il est en effet
efficace et relativement facile à programmer.
5.4.2 La relaxation par faux pas de temps

L'idée de base est ici, non pas de relaxer l'équation 5-3 de manière multiplicative
en calculant un barycentre entre les valeurs φ p,i et φ p ,i −1 comme dans l'équation 5-8,
mais directement de manière additive au niveau de l'équation discrétisée, par exemple
5-12:
(a + I )φ p ,i =
Eq. 5-19
p ∑
nb
a nb φ * nb + b + φ p ;i −1
où I est une inertie. La sur-relaxation correspond à des valeurs négatives de I et la sous-

relaxation à des valeurs positives.
La ressemblance entre les équations 5-1, 5-2 et 5-19 montre bien, une fois de
plus, l'analogie formelle entre pas de temps et itération. La relaxation par inertie revient
à résoudre, éventuellement à l'intérieur même d'un “ vrai ” pas de temps (régime
transitoire - schéma implicite), un problème en “ faux régime transitoire ”, avec un faux
0
pas de temps ∆ t * " tel que, par analogie entre I et ap
ρ p VP Eq. 5-20
∆t * =
I
On voit que l'inertie peut aussi dépendre du point considéré, si l'on préfère se
fixer à l'avance le faux pas de temps pour la relaxation.
5.5 Les méthodes par blocs

Les méthodes par blocs sont relativement récentes (20 ans environ). L'idée de
base est de partager le domaine en différents sous-domaines (blocs) qui contiennent.
Chacun un certain nombre de cellules de contrôle. On résout le système d'équations
dans Un bloc en supposant connu le résultat dans les blocs voisins, puis on change de
bloc de calcul, en itérant sur l'ensemble du domaine jusqu'à convergence totale.
La technique la plus courante est la résolution ligne-par-ligne. Dans cette méthode, le

bloc de calcul est une ligne de nœuds du maillage dans une direction x, par exemple,
pour un problème bidimensionnel en coordonnées cartésiennes. L'intérêt principal de
la méthode est que le T.D.M.A peut être appliqué pour résoudre l'équation de transport
61
et calculer φ en chaque point de la line de calcul. Un balayage selon les y croissants

permet de calculer le champ d'extensité sur tout le domaine. Un avantage de la méthode
ligne-par-ligne sur la méthode point-par-point (Gauss - Seidel) est que l'il formation en
provenance des conditions aux limites en x est transmise beaucoup plus rapidement sur
l'ensemble du domaine du fait de la résolution directe par leT.D.M.A. Comme dans la
méthode point-par-point d'ailleurs, le sens de balayage selon y peut être changé a
chaque itération afin d'accélérer encore la convergence.
La méthode ligne-par-ligne présente conceptuellement beaucoup de points

communs avec le schéma A.D.I. :
• utilisation répétée du T.D.M.A. bien que le problème soit multi dimensionnel,

• choix d'une direction privilégiée pour l'application du T.D.M.A.
Elle en diffère toutefois fondamentalement en ce se ns que le processus de base
reste itératif, et non direct. Quoi qu'il en soit, cette ressemblance formelle est la aussi
due à l'analogie numérique entre l'idée d'itération et l'idée de pas de temps.
L'ensemble formé par les chapitres 1 à 5 de ce cours doit permettre au lecteur

de discrétiser par la méthode des volumes finis, puis de résoudre numériquement, des
équations de transport diffusif, éventuellement couplées, a I exception (notable) des
équations de transport de quantité de mouvement, généralement appelées équations de
Navier-Stokes. Bien que quelques notions restent à acquérir, les idées des plus
importantes ont été introduites.
62
Références
Références
- An Introduction to CFD Finite volume method , H K Versteeg and W

Malalasekera - Second edition published, 2007
- Abbès AZZI, "Méthodes numériques, la méthode des volumes finis"
- Launder B.E. et Splanding D.b. 1974, the numerical computation of

turbulent flows, comuter Methods in Aplied Mechanics and Engineering,
vol.3,pp.269-289.
- Numerical Heat Transfer and Fluid Flow, Suhas Patankar - CRC Press,
1980
- Simulation numérique des phénomènes de transport, A. JARDY et H.

COMBEAU, 2001
63
64

Polycope FVM Vfinale Francais V 01

Transféré par

Droits d'auteur :

Formats disponibles

Polycope FVM Vfinale Francais V 01

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Polycope FVM Vfinale Francais V 01

Transféré par

Droits d'auteur :

Formats disponibles

UNIVERSITE KASDI MERBAH OUARGLA

Faculte des Sciences Appliquées

L’objectif de ce document pédagogique est de permettre aux étudiants de

Le premier chapitre illustre les notions des phénomènes de

Le dernier chapitre est conservé pour la résolution des équations

1.1 Bilan global et bilan différentiel 05

2.1 Méthode de conservation du volume fini 17

3.1 Problème conduction pure 27

4.1 Le schema centré (Centered Differencing Scheme, CDS) : 32

5.1 Système matriciel 50

Chapitre 1 : Notions des phénomènes

Les lois de conservation de la matière, de l'énergie, et de la quantité de

entrée de Φ + création de Φ = sortie de Φ +accumulation de Φ

Les premier et troisième termes correspondent à des transports d'extensité à

La distinction entre bilan global et bilan différentiel est directement liée à la

à des équations aux dérivées partielles, sont abondamment utilisés. Nous

Les lois phénoménologiques de la diffusion peuvent se mettre sous la forme:

φ = −λ grad V Eq. 1-1

Où : la densité de flux diffusifs de l'extensité Φ

: le coefficient de conductivité de ce transport

V : le potentiel dont le gradient provoque la diffusion de l'extensité Φ

L'extensité spécifique (ici, par unité de masse), notée par la suite Φ, et le

qui définit C = dΦ / dV, la capacité de l'extensité. On fait généralement apparaître

Tableau 1-1: Lois de la diffusion (analogies entre les transports)

Transport de Transport Transport de Transport

Extensité mi masse de i H enthalpie Px =mvx quantité de Q quantité

Densité de φi densité de flux φ densité de flux τ contrainte de j densité de

(kg m-2 s-1) (W m-2) ( N m-2) (A m-2)

Potentiel V ρi concentration T température vx vitesse dans la U potentiel

Conductivité Di coefficient de λ conductivité μ viscosité dynamique σ conductivité

Capacité 1/ρ inverse de la Cp chaleur C capacité

1.3 Le transport convectif :

Lorsqu'en un point du système que l'on étudie, la vitesse barycentrique est

ϕ=ρ φ v Eq. 2-4

Où ρ : est la masse volumique du système

Φ : l'extensité spécifique (par unité de masse)

v : le vecteur vitesse instantané

Ainsi, la densité de flux massique convectif d'un constituant dans un mélange

est égale à ρ v ωi , où ωi : est le titre massique de ce constituant, la densité de

flux thermique convectif à ρ v h , où h est l'enthalpie massique, et la densité

de flux convectif de quantité de mouvement selon x, à ρ v vx . On voit ainsi

1.4 L’équation de transport convecto-diffusif

On a vu que l'écriture d'un bilan différentiel d'une extensité quelconque conduit

Où (ρ ϕ) est l'extensité volumique, et S la vitesse de création de Φ par unité de

-pour le transport de matière, à la production (ou la disparition) de l'espèce

-pour le transport de quantité de mouvement, à une force volumique ou une

- pour le transport de chaleur, a un dégagement de chaleur par effet Joule ou par

Dans l'équation 1-5, la densité de flux φ correspond à la somme des transports

Après un léger remaniement, l'équation générale de transport prend sa forme

L'universalité de l'équation 1-7 générale de transport convecto-diffusif a pour

L'équation 1-7 est composée de quatre termes qui correspondent aux

Si La vitesse v x est suffisamment importante, ce qui peut être caractérisé par la

valeur du nombre adimensionnel de Péclet ρ Cp v x L / λ, le terme diffusif

par rapport au terme convectif (dérivée du premier ordre), ρ Cp v x ∂² T/ ∂x².

Selon nos définitions, il s'agit bien d'une équation parabolique suivant x et

L'intérêt d'une telle distinction sera particulièrement mis en évidence quand

1.4.1 Cas particuliers de l'équation générale

Nous allons à présent étudier quelques cas particuliers de l'équation générale.

1.4.1.1 Transport de matière