Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

∂T
où Ae = AP = ( ∆y / 2) ⋅ 1 et An = ( ∆x / 2) ⋅ 1 . En tenant compte que q = −λ P et
∂x P
∂T
= 0 , on obtient :
∂y P
TE − TP T − TP
λ e Ae + AP q + λ n An N = 0. (4.134)
∆x ∆y

En regroupant les termes, on obtient la forme générale de l’équation discrétisée pour le

noeud de coin 1 :

a P TP = a E TE + a N T N + b , (4.135)

λ e Ae λ n An
où aE = aN = aP = aE + aN b = AP q .
∆x ∆y

Les valeurs numériques des coefficients sont les suivantes :

0.1 0.1
1000 ⋅ ⋅ 0.01 1000 ⋅ ⋅ 0.01
λ e Ae 2 λ n An 2
aE = = =5 aN = = =5 ;
∆x 0.1 ∆y 0.1

0. 1
aP = aE + aN = 5 + 5 = 10 b = AP q = ⋅ 0.01 ⋅ 500 ⋅ 103 = 250 .
2

L’équation à résoudre est la suivante :

10T1 = 5T7 + 5T2 + 250 . (4.136)

Pour obtenir l’équation discrétisée pour le noeud 25 (noeud de coin “E – S”) on intègre
l’équation (4.116) sur le quart de volume de contrôle hachuré à la figure 4.27.

Fig. 4.27 Quart de volume de contrôle (le coin E-S).

∂ ∂T ∂ ∂T
λ dx dy + λ dx dy = 0 . (137)
1 / 4VC
∂x ∂x 1 / 4VC
∂y ∂y

83
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Méthode des volumes finis

Après l’intégration on a :

∂T ∂T ∂T ∂T
λ P AP − λ w Aw + λ n An − λ P AP = 0, (4.138)
∂x P ∂x w ∂y n
∂y P
∂T
où Aw = (∆y / 2) ⋅ 1 et An = (∆x / 2) ⋅ 1 . En tenant compte que =0 et
∂x P
∂T
= 0 , on obtient :
∂y P

T P − TW T − TP
− λ w Aw + λ n An N = 0. (4.139)
∆x ∆y

Finalement on obtient l’équation discrétisée pour le noeud 25 (noeud de coin) :

a P T P = aW TW + a N T N , (4.140)

λ w Aw λ n An
où aW = aN = a P = aW + a N .
∆x ∆y

Les valeurs numériques des coefficients sont les suivantes :

0. 1 0.1
1000 ⋅ ⋅ 0.01 1000 ⋅ ⋅ 0.01
λ A 2 λ A 2
aW = w w = =5 aN = n n = =5 ;
∆x 0. 1 ∆y 0.1

aP = aW + a N = 5 + 5 = 10 .

L’équation à résoudre est la suivante :

10T25 = 5T19 + 5T26 (4.141)

Finalement, le système d’équation à résoudre est le suivant :

10T1 = 5T7 + 5T2 + 250

20T2 = 10T8 + 5T1 + 5T3 + 500
20T3 = 10T9 + 5T2 + 5T4 + 500
20T4 = 10T10 + 5T3 + 5T5 + 500
20T5 = 10T11 + 5T4 + 5T6 + 500
20T7 = 5T1 + 5T13 + 10T8
40T8 = 10T2 + 10T14 + 10T7 + 10T9
40T9 = 10T3 + 10T15 + 10T8 + 10T10
40T10 = 10T4 + 10T16 + 10T9 + 10T11

84
 Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

40T11 = 10T5 + 10T17 + 10T10 + 10T12

20T13 = 5T7 + 5T19 + 10T14
40T14 = 10T8 + 10T20 + 10T13 + 10T15 (4.142)
40T15 = 10T9 + 10T21 + 10T14 + 10T16
40T16 = 10T10 + 10T22 + 10T15 + 10T17
40T17 = 10T11 + 10T23 + 10T16 + 10T18
20T19 = 5T13 + 5T25 + 10T20
40T20 = 10T14 + 10T26 + 10T19 + 10T21
40T21 = 10T15 + 10T27 + 10T20 + 10T22
40T22 = 10T16 + 10T28 + 10T21 + 10T23
40T23 = 10T17 + 10T29 + 10T22 + 10T24
10T25 = 5T19 + 5T26
20T26 = 10T20 + 5T25 + 5T27
20T27 = 10T21 + 5T26 + 5T28
20T28 = 10T22 + 5T27 + 5T29
20T29 = 10T23 + 5T28 + 5T30

La solution numérique du système (4.142) est présentée dans le tableau 4.8. On constate
une bonne correspondance, pour presque le même numéro de noeuds, entre les résultats
obtenus avec la méthode des volumes finis et celles obtenus avec la méthode des élements
finis en utilisant le logiciel QuickField (Student’s QuickField TM 3.40a, Copyright  1993-97
Tera Analysis). Pour apprécier la précision, il faut comparer les résultats numériques avec la
solution analytique. À la figure 4.28 sont présentées les courbes isothermes obtenues à l’aide
du logiciel QuickField.

Tableau 4.8 La solution numérique du système (4.142) (la méthode des volumes finis – MVF,
30 noeuds et 546 noeuds, et la méthode des élements finis – MEF, 31 noeuds et 485 noeuds)

1 2 3 4 5 6 7 8 9 10 11 12 13

T1 T2 T3 T4 T5 T7 T8 T9 T10 T11 T13 T14 T15

319.8 312.9 291.8 254.7 196.3 276.7 270.0 249.8 215.3 165.3 246.9 240.7 222.1

320.6 311.9 289.4 250.8 190.3 277.3 270.2 249.9 216.0 167.4 246.8 240.6 221.6

321.2 314.5 293.7 257.6 201.5 278.0 271.4 251.3 216.9 166.5 247.9 241.6 223.0

321.2 314.5 293.8 257.6 201.6 278.0 271.4 251.2 216.8 166.5 247.9 241.6 223.0

85
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Méthode des volumes finis

14 15 16 17 18 19 20 21 22 23 24 25
Noeud

T16 T17 T19 T20 T21 T22 T23 T25 T26 T27 T28 T29
MVF
191.4 149.7 229.5 223.7 206.4 178.6 141.9 223.8 218.2 201.4 174.6 139.6 30
noeuds

MEF
191.3 148.7 228.9 222.8 205.3 177.4 141.7 223.4 217.6 200.7 173.9 139.2 31
noeuds
MVF
192.0 149.9 230.3 224.4 206.9 178.9 142.0 224.5 218.7 201.8 174.7 139.6 546
noeuds
MEF
192.1 149.9 230.2 224.4 206.9 178.8 142.0 224.4 218.7 201.8 174.7 139.6 485
noeuds

En analysant les résultats du tableau 4.8 on constate une très bonne correspondance
entre les résultats obtenus avec la méthode des volumes finis (546 noeuds) et celles obtenus
avec la méthode des éléments finis (485 noeuds). Comme la matrice des coefficients du
système d’équations algébriques est une matrice de type bande, on peut utiliser l’algorithme
de Thomas adapté au cas bidimensionnel (voir le code source à l’Annexe G – THER2Ds1).

4.3.3 Application de l’algorithme de Thomas aux problèmes 2D (ou ligne par ligne)

L’algorithme de Thomas peut être appliqué itérativement pour résoudre le système

d’équations algébriques dans le cas des problèmes 2D. On combine la méthode directe de
l’algorithme de Thomas dans une direction et la méthode itérative de Gauss-Seidel dans
l’autre direction. Soit l’exemple du maillage de la figure 4.29

Fig. 4.28 Les isothermes.

86
 Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

Fig. 4.29 Application du TDMA dans le cas 2D.

Les équations discrétisées sur la ligne j = 3 , par exemple, sont les suivantes :

* + a T* + b
a32T32 = a31T31 + a33T33 + a22T22 42 42 32
* + a T* + b
a33T33 = a32T32 + a34T34 + a23T23 (4.143)
43 43 33
* + a T* + b .
a34T34 = a33T33 + a35T35 + a24T24 44 44 34

Si 1 ≤ j ≤ m , on résout m problèmes 1D dans la direction “ x ” par le TDMA, en balayant

tous les indices j ∈ [1, m] . Les températures Tij* sont les dernières valeurs de la température
calculées (l’itération précedente) sur les lignes voisines de la ligne où l’on effectue le TDMA.

Remarques
♦ La convergence de la méthode ligne par ligne est plus rapide que celle de Gauss-Seidel car
l’information contenue sur les frontières de la ligne où s’applique le TDMA est transmise
instantanément à l’intérieur du domaine et cela quel que soit le nombre de noeuds sur la ligne.
♦ La vitesse de transmission dans l’autre direction est la même que pour la méthode de
Gauss-Seidel.
♦ En alternant les directions d’application du TDMA, on peut accélérer la propagation de
l’information contenue sur toutes les frontières vers l’intérieur du domaine.
♦ La géométrie et le maillage peuvent donner des coefficients a vs beaucoup plus grands dans
une direction que dans l’autre. À ce moment-là, il faut prendre pour direction du TDMA la

Fig. 4.30 Choix de la direction du TDMA.

87
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Méthode des volumes finis

direction des grands coefficients car cela accélère la convergence. Soit l’exemple suivant de la
figure 4.30. Les coefficients de l'
équation discrétisée sont les suivants :

λ w ∆y λ e ∆y λ s ∆x λ n ∆x
aW = aE = aS = aN = ,
δx w δx e δy s δy n

où δx w(e) >> δy s ( n ) ∆x >> ∆y .

Ceci implique que aS ( N ) >> aW ( E ) et donc on prend pour direction du TDMA la direction
“ y ”. L’équation discrétisée est la suivante :

a P TP = a S TS + a N T N + aW TW* + a E T E* + b . (4.144)

♦ Dans certains cas, la direction de balayage pour les itérations de la méthode de Gauss-
Seidel est imposée par les conditions aux limites ou par la direction de l’écoulement (de
l’amont vers l’aval).

Exemple 2
On considère une plaque métallique rectangulaire ( 0.05 x 0.04 m ). La conductivité
thermique du matériau de la plaque est λ = 4 W/m/K . Toutes les frontières de la plaque sont
maintenues à la température de 0 °C et le terme source est S = 40 ⋅ 10 6 W/m 3 .
Calculer la distribution stationnaire de la température dans la plaque en utilisant le
maillage de la figure 4.31 ( ∆x = ∆y = 0.01 m ).

Solution
L’équation à résoudre est la suivante :

∂ ∂T ∂ ∂T
λ + λ +S =0. (4.145)
∂x ∂x ∂y ∂y

L’équation discrétisée pour un noeud intérieur (le noeud 16 par exemple) est la suivante :

a P T P = aW TW + a E T E + a S TS + a N T N + b , (4.146)

λAw λAe λAs λAn

où aW = aE = aS = aN = ,
∆x ∆x ∆y ∆y

a P = aW + a E + a S + a n b = S ∆x∆y .

Le terme source, S étant constant, ce n’est pas nécessaire d’être linéarisé. Les valeurs
des coefficients voisins, pour le maillage choisi, sont :

4 ⋅ 0.01 ⋅ 1
aW = aE = aS = a N = =4 ;
0.01
a P = aW + a E + a S + a n = 4 + 4 + 4 + 4 = 16 ;

88
 Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

b = S ∆x∆y = 40 ⋅ 10 6 ⋅ 0.01 ⋅ 0.01 = 4000 .

Fig. 4.31 Maillage et conditions aux limites pour l’exemple 2

Pour les noeuds de frontière ( noeuds 1, 2, 3, 4, 5, 6, 7, 12, 13, 18, 19, 24, 25, 26, 27, 28, 29, 30 )
la température étant connue ce n’est pas nécessaire d’écrire les équations discrétisées. Alors le
nombre d’équations à résoudre est égal à 12, les équations discrétisées sont les suivantes :

4T8 = T14 + T9 + 1000

4T9 = T15 + T8 + T10 + 1000
4T10 = T16 + T9 + T11 + 1000
4T11 = T17 + T10 + 1000
4T14 = T8 + T20 + T15 + 1000
4T15 = T9 + T21 + T14 + T16 + 1000 (4.147)
4T16 = T10 + T22 + T15 + T17 + 1000
4T17 = T11 + T23 + T16 + 1000
4T20 = T14 + T21 + 1000
4T21 = T15 + T20 + T22 + 1000
4T22 = T16 + T21 + T23 + 1000
4T23 = T17 + T22 + 1000

En regroupant les inconnues, le système d’équations à résoudre est le suivant :

89
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Méthode des volumes finis

4 −1 0 0 −1 0 0 0 0 0 0 T8
0 1000
−1 4 −1 0 0 −1 0 0 0 0 0 T9
0 1000
0 −1 4 −1 0 0 −1 0 0 0 0 T10
0 1000
0 0 −1 4 0 0 0 −1 0 0 0 T11
0 1000
−1 0 0 0 4 −1 0 0 −1 0 0 T14
0 1000
0 −1 0 0 −1 4 −1 0 0 −1 0 T0 1000
⋅ 15 = (4.148)
0 0 −1 0 0 −1 4 −1 0 0 −1 0 T16 1000
0 0 0 −1 0 0 −1 4 0 0 0 −1 T17 1000
0 0 0 0 −1 0 0 0 4 −1 0 0 T20 1000
0 0 0 0 0 −1 0 0 −1 4 −1 0 T21 1000
0 0 0 0 0 0 −1 0 0 −1 4 −1 T22 1000
0 0 0 0 0 0 0 −1 0 0 −1 4 T23 1000

La solution du système (4.148) est :

T8 746.479
T9 1028.000
T10 1028.000
T11 746.479
T14 957.746
T15 1338.000
= (4.149)
T16 1338.000
T17 957.746
T20 746.479
T21 1028.000
T22 1028.000
T23 746.479

Fig. 4.32 Les courbes isothermes (obtenues à l’aide du logiciel QuickField)

90
 Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

La solution analytique du problème pour un quart du domaine d’analyse (fig. 4.34) est
la suivante, [28] :

πx πy
∞ (−1) n cos (2n + 1) cosh (2n + 1)
S (a 2− x2 ) 16 Sa 2 2a 2a
T ( x, y ) = − (4.150)
2λ π 3λ n=0 (2n + 1)3 cosh (2n + 1) πb
2a

Fig. 4.34 Le domaine de calcul analytique

Fig. 4.34a Les courbes isothermes (solution analytique)

91
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Méthode des volumes finis

À la figure 4.34-a on présente la solution analytique obtenue à l’aide de l’expression

(4.150) pour le domaine de calcul hachuré de la figure 4.34. À la figure 4.34-b on présente la
solution numérique. On constate une très bonne correspondance entre la solution numérique et
la solution analytique (mais il faut tenir compte que même la solution analytique est le résultat
d’un calcul numérique d’évaluation d’une série tronquée). Dans le tableau ci-dessous on
présente une comparaison des résultats obtenus en utilisant la méthode des volumes finis (le
programme THERM2D2), la solution analytique et à l’aide du logiciel QuickField.

Fig.4.34b Les courbes isothermes (solution numérique)

x=0 T[K ] T[K ] T[K ] Erreur [ % ]

y[mm] QuickField Analytique numérique MVF
1 1429.10 1430.000 1428.26 0.12
3 1413.70 1416.000 1414.27 0.12
5 1386.80 1388.000 1385.97 0.14
7 1341.60 1345.000 1342.71 0.17
9 1281.60 1286.000 1283.49 0.19
11 1206.50 1209.000 1206.97 0.16
13 1113.30 1113.000 1111.47 0.14
15 996.47 996.558 994.92 0.16
17 852.53 856.310 854.90 0.16
19 683.60 689.757 688.63 0.16
21 486.86 493.798 493.01 0.16
23 259.76 265.067 264.66 0.15
485 noeuds 546 noeuds

92
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

93
 Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique

4.4 Conduction thermique instationnaire en deux dimensions (2D)

4.4.1 Forme générale de l’équation discrétisée

L’équation de conduction thermique instationnaire pour un problème 2D est la suivante

∂T ∂ ∂T ∂ ∂T
ρc p = λ + λ + S. (4.151)
∂t ∂x ∂x ∂y ∂y

Fig. 4.35 Volume de contrôle en deux dimensions.

L’intégration de l’équation (4.151) sur le volume de contrôle schématisé à la figure

4.35 donne :
t + ∆t t + ∆t t + ∆t t + ∆t
∂T ∂ ∂T ∂ ∂T
ρc p dx dy dt = λ dx dy dt + λ dy dx dt + Sdx dy dt
t VC
∂t t VC
∂x ∂x t VC
∂y ∂y t VC
(4.152)
t + ∆t t + ∆t
∂T ∂T ∂T ∂T
ρc p (TP − TP0 )∆x∆y = λe − λw ∆y dt + λn − λs ∆x dt +
t
∂x e ∂x w t
∂y n
∂y s
t + ∆t
S ∆x∆y dt
t

En remplaçant les gradients de température on obtient :

t + ∆t
TE − TP T P − TW
(
ρc p T P − T P0 ∆x∆y = ) λe
δx e
− λw
δx w
∆y dt +
t
t + ∆t
T N − TP T P − TS
λn − λs ∆x dt + (4.153)
t
δy n δy s
t + ∆t
S ∆x∆y dt
t

93