Chapitre 2 : Méthode des volumes finis

Méthode des volumes finis 1 / 38

De la modélisation à la simulation numérique
Les différentes étapes pour modéliser un système complexe :
Problème physique continu est décrit par un modèle mathématique
continue (mise en équations) ;
Elaboration d’un maillage : discrétisation de domaine géométrique en petits
éléments appelés volume de contrôle.
Modèle mathématique continu est discrétisé en s’appuyant sur une(des)
méthode(s) numérique(s) ;
Equations discrétisées (souvent système linéaire) sont approximées à l’aide
des schémas numériques appropriés, l’algorithme de résolution est établie ;

Transcription informatique et programmation des relations discrètes ;

Simulation numérique et exploitation des résultats.
L’ingénieur peut être amené à intervenir sur l’une ou plusieurs de ces différentes
étapes.
Méthode des volumes finis 2 / 38
Discrétisation des EDP
Pour passer d’une problème exact continu régit par une EDP au problème approché
discret, il existe trois grandes familles de méthodes :
Différences finis
La méthode des différences finies est la première méthode utilisée pour résoudre des
équations différentielles. La méthode utilise le développement en série de Taylor pour
la discrétisation des dérivées associées aux variables.

Eléments finis
La méthode des éléments finis consiste à remplacer un problème continu initial par
une formulation variationnelle. La formulation variationnelle est obtenue en intégrant
le produit scalaire de l’équation par une fonction test.

Volumes finis
Dans la méthode des volumes finis le domaine est divisé en un nombre fini de volumes
de contrôle qui recouvrent le domaine de calcul. Sur chaque volume de contrôle les
lois de conservation sont appliquées pour déterminer les différentes variables associées
au problème en des points (ou nœuds) qui ne correspondent pas forcément avec les
points du maillage ou de la grille.
Méthode des volumes finis 3 / 38
Pourquoi utiliser l’approche "Volumes Finis" ?

Différences finis
Bien connue
Mise en œuvre simple pour une géométrie simple
Mise en œuvre difficile pour une géométrie complexe
Pas toujours consérvative
Utilisation dans des codes de "recherche"

Eléments finis
Basée sur des fondements mathématiques rigoureux
S’adapte à une géométrie quelconque
Difficultés pour résoudre les termes non-linéaires
Plus coûteuse
Très utilisée dans le domaine de mécanique des structures.

Méthode des volumes finis 4 / 38

Pourquoi utiliser l’approche "Volumes Finis" ?

Volumes finis
Approche très "physique" : bilan des flux
S’adapte à une géometrie quelconque
Plusieurs schémas pour la résolution des termes non-lineaires hyperboliques
Conservative (par sa formulation)
La base de tout les codes généralistes en Mécanique des Fluides : Fluent
et CFX (ANSYS), StarCCM+, OpenFoam (Libre)...

Méthode des volumes finis 5 / 38

Méthodologie de résolution par la MVF
Les différentes étapes de la mise en œuvre sont :
Le domaine de calcul est discrétisé en un nombre fini de points (les nœuds
du maillage), autour desquels on définit des volumes élémentaires (appe-
lés volumes de contrôle) contigus, juxtaposés et sans discontinuités aux
interfaces ;

Les équations modèles, sous forme conservative, sont intégrées sur chaque
volume de contrôle (VC) ;

Les intégrales sur un volume de contrôle en un nœud donné sont évaluées

en approchant la variation de Φ par des profils ou des lois d’interpolation
entre les nœuds voisins du point considéré ;

Écriture des équations algébriques en fonction des valeurs de Φ aux nœuds

du maillage ;

Résolution du système algébrique linéaire obtenu.

Méthode des volumes finis 6 / 38
Equations de conservation

Pour un écoulement incompressible, on a

∂ui


 =0


 ∂xi
  
∂ui ∂ui ∂p


ρ + uj =− + ρgi + µ∆ui

 ∂t ∂xj ∂xi

∂2T
 
∂T ∂T



 ρcV
 + ui =k +ϕ
∂t ∂xi ∂xi ∂xi

En plus de ces équations, on ajoute les équations d’états des gaz, ρ = ρ (p, T )
et p = p (ρ, T )

Méthode des volumes finis 7 / 38

Forme générale d’une équation de conservation

Les équations de conservation peuvent être réécrites sous la forme générale

suivante :
∂(ρΦ)
+ div (ρΦ⃗u ) = div (Γgrad Φ) + SΦ
∂t
C’est l’équation de convection-diffusion

Φ est une variable caractéristique généralisée du système (vitesse, tempéra-

ture...).
Γ est le coefficient de diffusion généralisé (propre au phénomène traité). Il
exprime la propriété de transport comme la viscosité ou la diffusivité, conjoin-
tement avec le gradient de la variable Φ.
S est le terme source spécifique à chaque variable. Il exprime le processus de
production interne comme la génération de la chaleur dans le fluide.

Méthode des volumes finis 8 / 38

 Variables
Équations Φ Γ S
Continuité 1 0 0
−∂p
Conser. Quant. Mvt. x u µ + ρgx
∂x
−∂p
Conser. Quant. Mvt. y v µ + ρgy
∂y
−∂p
Conser. Quant. Mvt. z w µ + ρgz
∂z
k
Conservation de l’énergie T Si
Cp
En intégrant sur le volume de contrôle :
Z Z Z Z
∂
ρΦdV + div (ρΦ⃗
u )dV = div (ΓgradΦ)dV + SΦ dV
∂t VC VC VC VC

Méthode des volumes finis 9 / 38

Résolution des problèmes de diffusion stationnaire
1D

L’équation de la diffusion stationnaire s’écrit :

div (Γgrad Φ) + SΦ = 0

− Cas 1D :  
d dΦ
Γ + SΦ = 0
dx dx
Cas typique de la conduction de chaleur dans une barre métallique d’extrémité
A et B.
▶ Étape 1 : Génération du maillage
Le domaine d’étude est subdivisé en un nombre fini de nœuds. Chaque volume
fini entoure un nœud principal "P". Les nœuds voisins sont "E" côté Est et
"W" côté West.

Méthode des volumes finis 10 / 38

Les conditions aux limites sont supposées connues aux extrémités de domaine.
▶ Étape 2 : Intégration de l’équation de diffusion sur le VC
L’équation précédente est intégrée sur le volume de controle limité par les faces
(w,e).
Pour un problème 1D, le volume d’un volume de contrôle est ∆x × 1 × 1.
Z e   Z e
d dΦ
⇔ Γ dx + SΦ dx = 0
w dx dx w

dΦ e
  Z e    
dΦ dΦ
Γ + SΦ dx = 0 ⇒ Γ − Γ + S̄Φ ∆x = 0
dx w w dx e dx w
Méthode des volumes finis 11 / 38
S̄Φ est la valeur moyenne du terme source sur le VC.

C’est la forme discrétisée de l’équation de diffusion. Elle peut être écrite sur
des nœuds P du domaine de calcul.
Signification physique :
Le bilan de ϕ sur le VC ⇒ le flux diffusif de ϕ quittant la face EST du VC - le
flux diffusif de ϕ entrant de la face WEST du VC = génération de ϕ dans le
VC.
   
dΦ dΦ
Γ − Γ + S̄Φ ∆x = 0
dx e dx w
Calcul des gradients (flux diffusif) sur les faces du VC :
Le profil linéaire de température, entre les nœuds du maillage, définie les déri-
vées aux interfaces par :
dΦ ΦE − ΦP ΦE − ΦP
  
 Γ

= Γe = Γe
 dx e xE − x P

δxPE
d Φ Φ P − Φ W Φ P − ΦW
 Γ = Γw = Γw


dx w x P − xW δxPW
Méthode des volumes finis 12 / 38
et on utilise une approximation linéaire pour calculer le terme Γ sur les faces
du VC :
Γ + ΓW

 Γw = P

2
 Γe =
 ΓP + ΓE
2
Linéarisation du terme source :
Le terme source peut être constant ce qui ne pose aucun problème particulier.
Mais souvent il dépend explicitement de la variable Φp . 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 :

S̄Φ = Su + Sp Φp

où Su est la partie constante du terme source et Sp et le coefficient de Φp .

   
ΦE − ΦP ΦP − ΦW
Γe − Γw + Su ∆x + Sp ΦP ∆x = 0
δxPE δxPW

Méthode des volumes finis 13 / 38

      
Γe Γw Γw Γe
+ − ∆xSP ΦP = ΦW + ΦE + Su ∆x
δx δxPW δxPW δxPE
| PE {z } | {z } | {z }
aP aW aE

Ainsi l’équation algébrique finale s’écrit sous la forme :

aP ΦP = aW ΦW + aE ΦE + Su ∆x
▶ Étape 3 : Résolution du système algébrique

La valeur de Φp au centre du VC est exprimée en fonction des valeurs de la

grandeur aux nœuds voisins à l’aide de l’équation discrète.
aP ΦP = aW ΦW + aE ΦE + Su ∆x
Cette équation s’écrit autant de fois qu’il y a de nœuds où la variable Φ est
inconnue.
Dans le cas d’un maillage uniforme, ∆x = (δx)PW = (δx)PE .
Γw Γe
alors on aura aW = , et aE =
∆x ∆x
Méthode des volumes finis 14 / 38
L’écriture de l’équation pour tous les nœuds du domaine conduit à un système
algébrique qui s’écrit sous la forme matricielle suivante :

[M]{ϕ} = {b}
Méthodes de résolution du système algébrique :
◦ Gauss-Seidel,
◦ Tri-Diagonal Matrix Algorithm (TDMA)

◦ Méthode de Gauss-Seidel
Une méthode itérative de résolution d’un système linéaire (de dimension finie)
de la forme [M]{ϕ} = {b}. On calcule les composantes Φk+1 i de Φk+1 de
manière séquentielle pour i = 1, . . . , n par :
 
i−1 n
1  X X
Φk+1
i = bi − mij Φk+1
j − mij Φkj 
mii
j=1 j=i+1

En donnant aux inconnues ϕki des valeurs arbitraires initiales ϕ0i (pour k = 0)
Méthode des volumes finis 15 / 38
Le processus sera arrêté si : |ϕk+1
i − ϕki |< ε, avec ε est la précision absolue
imposée par l’utilisateur.
◦ Tri-diagonal matrix algorithm (TDMA)
Cet algorithme permet de calculer la solution d’un système linéaire lorsque la
matrice est tri-diagonale. C’est notre cas, car les équations discrétisées dans le
système linéaire s’écrivent :

ai Φi = bi Φi+1 + ci Φi−1 + di 2≤i ≤N

Où P correspond au nœud i, E correspond au nœud i + 1, W correspond au

nœud i − 1.

Soit un système d’équation ayant la forme tri-diagonale :

Méthode des volumes finis 16 / 38

 Φ1 = C1
−c2 Φ1 + a2 Φ2 − b2 Φ3 = C2
−c3 Φ2 + a3 Φ3 − b3 Φ4 = C3
−c4 Φ3 + a4 Φ4 − b4 Φ5 = C4
··· ··· ··· ··· =·
−cN ΦN−1 + aN ΦN − bN ΦN+1 = CN
ΦN+1 = CN+1
Φ1 et ΦN+1 sont donnés par les conditions aux limites.
− Élimination de Gauss
L’équation correspondant à i=2 relie ϕ1 ,ϕ2 et ϕ3 et puisque ϕ1 est donnée, alors
cette équation relie ϕ2 et ϕ3 seulement. En continuant ainsi jusqu’à l’équation
correspondant à i = N reliant ϕN et ϕN+1 , avec ϕN+1 donnée.
D’après ce qu’on vient de faire, la nouvelle équation i est de la forme :
Φi = Pi ϕi+1 + Qi
Méthode des volumes finis 17 / 38
Où Pi et Qi à déterminés.
Cette équation s’obtient après avoir déterminé la relation liant ϕi−1 à ϕi :

Φi−1 = Pi−1 ϕi + Qi−1

Où P1 , Q1 , P2 , Q2 ,..., Pi−1 et Qi−1 sont déjà déterminés.

Déterminons maintenant Pi , Qi
Pour ce faire, portons cette dernière équation dans l’équation avant la
dernière, nous obtenons :
bi di + ci Qi−1
ai Φi = bi ϕi+1 +ci (Pi−1 ϕi + Qi−1 )+di ⇒ ϕi = ϕi+1 +
ai − ci Pi−1 ai − ci Pi−1

Alors, nous obtenons :

bi di + ci Qi−1
Pi = et Qi = , i = 2, 3, ..., N
ai − ci Pi−1 ai − ci Pi−1

Méthode des volumes finis 18 / 38

Application 1 : Problème de conduction de la
chaleur en une dimension et stationnaire

Considérons le problème de la conduction thermique dans une barre métallique

de longueur L=0,5 m. En utilisant la méthode des volumes finis, déterminer les
valeurs approchées de T (0, 05), T (0, 15), T (0, 25), T (0, 35) et T (0, 45), où
T (x) est la solution exacte du problème. On donne la conductivité thermique
k = 1000 W/m.K.
 2
d (kT )
 dx 2 = 0 0 ⩽ x ⩽ 0, 5




 T (0) = 100


T (0, 5) = 500


Méthode des volumes finis 19 / 38

Application 2 : Problème de conduction de la
chaleur en une dimension et stationnaire en
présence de terme source

En utilisant la méthode des volumes finis, déterminer les valeurs approchées

de T (0, 002), T (0, 006), T (0, 01), T (0, 014) et T (0, 018), où T (x) est la
solution exacte du problème.
 2
d T 6
 dx 2 + 2 × 10 = 0 0 ⩽ x ⩽ 0, 02




 T (0) = 100


T (0, 02) = 200


Méthode des volumes finis 20 / 38

VF pour les problèmes de diffusion stationnaire 2D
L’équation de la diffusion stationnaire dans le cas 2D est :
   
∂ ∂ϕ ∂ ∂ϕ
Γ + Γ + Sϕ = 0
∂x ∂x ∂y ∂y

Le discétisation de domaine 2D

Méthode des volumes finis 21 / 38

L’intégrale de l’équation de la diffusion stationnaire sur le volume de contrôle
donne
Z   Z   Z
∂ ∂Φ ∂ ∂Φ
Γ dx · dy + Γ dx · dy + SΦ dV = 0
VC ∂x ∂x VC ∂y ∂y VC

Cela conduit à
     
∂Φ ∂Φ
Γe ∆y − Γw ∆y
∂x e ∂x w
    
∂Φ ∂Φ
+ Γn ∆x − Γs ∆x + S̄Φ ∆V = 0
∂y n ∂y s

- Calcul des gradients (flux diffusif) sur les faces du VC

Flux à travers la face West

∂ϕ (ΦP − ΦW )
Γw = Γw
∂x w δxWP

Méthode des volumes finis 22 / 38

 Flux à travers la face Est
∂Φ (ΦE − ΦP )
Γe = Γe
∂x e δxPE
Flux à travers la face Sud
∂Φ (ΦP − ΦS )
Γs = Γs
∂y s δySP
Flux à travers la face Nord
∂Φ (ΦN − ΦP )
Γn = Γn
∂y n δyPN
- Système algébrique
En remplaçons les flux dans l’équation :
     
∂Φ ∂Φ
Γe ∆y − Γw ∆y
∂x e ∂x w
    
∂Φ ∂Φ
+ Γn ∆x − Γs ∆x + S̄Φ ∆V = 0
∂y n ∂y s
Méthode des volumes finis 23 / 38
Et en linéarisant le terme source, on obtient le système algébrique suivant :

aP Φp = aW ΦW + aE ΦE + aN ΦN + aS ΦS + Su ∆x∆y

avec :
Γw ∆y Γe ∆y Γs ∆x Γn ∆x
aW = ; aE = ; aS = ; aN =
δxWP δxPE δySP δyPN
ap = aW + aE + aS + aN − Sp ∆x∆y

Méthode des volumes finis 24 / 38

VF pour les problèmes de convection-diffusion
Dans les problèmes où le mouvement du fluide est mis en évidence et joue un
rôle significatif, il faut tenir en compte l’effet de convection. Ceci est pris en
considération dans l’équation de convection-diffusion (stationnaire) :
div(ρΦ−
→
v ) = div(Γ grad Φ) +SΦ
| {z } | {z }
terme convectif terme diffusif

Pour la diffusion, la propagation de la propriété Φ s’effectue équitablement

dans toutes les directions,
Pour la convection, la propagation de la propriété Φ s’effectue seulement
suivant le sens de l’écoulement du fluide ; elle dépend fortement du sens
et vitesse de l’écoulement.
□ Discrétisation de l’équation de convection-diffusion
Le terme diffusif : déjà discrétisé et montre une grande stabilité.
Le terme convectif : des difficultés apparaissent lors de l’approximation
de la propriété Φ sur les faces de chaque VC ainsi que son flux convectif.
Méthode des volumes finis 25 / 38
C’est pourquoi, un grand nombre de schémas, se différenciant par leur précision
sont utilisées pour traiter ce problème. Parmi ces schémas, on peut citer :

Schéma Centré
Schéma décentré (Upwind)
Schéma exponentiel
Schéma hybride

▶ Application : convection-diffusion sans terme source 1-D :

 
d d dΦ
(ρΦu) = Γ
dx dx dx

Le fluide doit aussi satisfaire l’équation de continuité :

d
(ρu) = 0
dx

Méthode des volumes finis 26 / 38

En intégrant sur le VC, on obtient :
    
dΦ dΦ
(ρuΦ)e − (ρuΦ)w = Γ − Γ

dx e dx w
(ρu)e − (ρu)w = 0


Γ
On pose F = ρu et D =
δx
donc sur les faces du VC on a :

 Fw = (ρu)w
 ; Fe = (ρu)e
Γw Γe
 Dw =
 ; De =
δxWP δxPE

Nous pouvons réécrire l’équation précédente comme suit :


Fe Φe − Fw Φw = De (ΦE − ΦP ) − Dw (Φp − ΦW )
Fe − Fw = 0

Méthode des volumes finis 27 / 38

1. Schéma Centré
Ce schéma a été utilisé pour représenter le terme convectif. On exprime la
variable Φ sur la face du volume de contrôle par une approximation linéaire
comme suit :
ΦE + ΦP ΦW + Φ P
Φe = et Φw =
2 2
En remplaçant dans l’équation :

Fe Φe − Fw Φw = De (ΦE − ΦP ) − Dw (ΦP − ΦW )

On obtient
        
Fw Fe Fw Fe
Dw + + De − + (Fe − Fw ) ΦP = Dw + Φw + De − Φ
2 2 2 2

=⇒ aP ΦP = aW ΦW + aE ΦE

Méthode des volumes finis 28 / 38

2. Schéma décentré (Upwind Differencing Scheme)
Le schéma dit Upwind prend en considération la direction du flux lors de la
détermination de la propriété Φ sur la face du VC.
Lorsque l’écoulement est dans le sens positive, on a :

uw > 0, ue > 0 (Fw > 0, Fe > 0)

Alors l’équation précédente :

Fe Φe − Fw Φw = De (ΦE − ΦP ) − Dw (ΦP − ΦW )

devient :

Fe ΦP − Fw ΦW = De (ΦE − ΦP ) − Dw (ΦP − ΦW )
Méthode des volumes finis 29 / 38
 (Fe + De + Dw )ΦP = (Fw + Dw )ΦW + De ΦE

aP ΦP = aW ΦW + aE ΦE
Résumé du schéma Upwind :

aW aE
Fw > 0, Fe > 0 Dw + Fw De
Fw < 0, Fe < 0 Dw De − F e

Ou bien sous la forme compacte suivante en tenant compte les deux directions
de l’écoulement,

aW aE aP
Dw + max(Fw , 0) De + max(0, −Fe ) Dw + max(−Fw , 0) + De + max(Fe , 0) = aW + aE + Fe − Fw

Méthode des volumes finis 30 / 38

Application 3 : Problème de Convection et
Diffusion en régime stationnaire
Soit le problème aux limites suivant :
  
dϕ

 d Γ
 d (ρuϕ) = dx


0⩽x ⩽1

dx dx


 ϕ (0) = 1


ϕ (1) = 0


En utilisant la méthode des volumes finis et le schéma centré, déterminer les

valeurs approchées de ϕ (0, 1), ϕ (0, 3), ϕ (0, 5), ϕ (0, 7), ϕ (0, 9) pour les trois
cas ci-dessous , où ϕ (x) est la solution exacte du problème. On prend Γ =
0, 1Kg /m.s et ρ = 1Kg /m3 pour les applications numériques.
− Cas 1 : u=0,1 m/s
− Cas 2 : u=2,5 m/s
− Cas 3 : u=2,5 m/s, ∆x = 0, 05
Méthode des volumes finis 31 / 38
Corrigé 3 : résolution par le schéma centré
− Cas 1 : u=0,1 m/s

Méthode des volumes finis 32 / 38

− Cas 2 : u=2,5 m/s

Méthode des volumes finis 33 / 38

− Cas 2 : u=2,5 m/s, ∆x = 0, 05

La comparaison des résultats des cas 2 et 3 montre qu’un maillage plus fin
à déterminé la réduction du nombre de Peclet de 5 à 1,25. Le schéma de
différences centrales fournis de bons résultats quand le rapport F/D est petit.

Méthode des volumes finis 34 / 38

Résoudre le problème considéré à l’exemple 3 en utilisant le schéma de discré-
tisation "upwind" pour les cas :
− Cas 1 : u=0,1 m/s ;
− Cas 2 : u=2,5 m/s,

Méthode des volumes finis 35 / 38

Corrigé 3 : résolution par le schéma upwind
Cas 1 : u=0,1 m/s

Méthode des volumes finis 36 / 38

Cas 1 : u=2,5 m/s

Méthode des volumes finis 37 / 38