Chapitre 4 - A
Chapitre 4 - A
Chapitre 4 - A
d dT
+ S = 0. (4.1)
dx dx
L’équation discrétisée sur le volume de contrôle présenté à la figure 4.1 est la suivante :
a P TP = a E TE + aW TW + b ; (4.2)
e w
aE = aW = a P = a E + aW - S P x b = Sc x .
( x )e ( x )w
x e+ xe
e = fe P + (1 fe ) E avec fe = et 1 fe = . (4.3)
xe x
44
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
Si l’interface " e " (fig. 4.2) est située à mi-distance de P et E , alors f e = 1 / 2 et on obtient :
P + E
e = . (4.4)
2
Si l’on considère le flux à l’interface " e " (voir fig. 4.3) on peut écrire :
TP TE TP TE
qe = = . (4.6)
xe x e+ xe
+
P E e
45
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
1
e = . (4.8)
1 fe fe
+
P E
2 P E
e = . (4.9)
P + E
qui est la moyenne harmonique des conductivités thermiques des noeuds voisins du maillage.
Remarque
La formule (4.9) est plus recommandée que la formule (4.4) :
• Si E 0 (dans la formule (4.4)), ce qui signifie qu’on a un matériau isolant en x e+ , ce
qui implique qe 0 ), alors E f e P 0 , et la densité du flux thermique est donnée, à
la limite, par la relation suivante :
fe P (TP TE )
qe 0, (4.10)
xe
E (TP TE ) E (TP TE )
qe = , (4.11)
f e xe x e+
b = b(T ) a P = a P (T ) a vs = a vs (T ) , (4.12)
46
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
( ) ( )
S (T ) = S c T P* + S P T P* T P , (4.13)
où TP* est la température obtenue à l’itération précédente. En même temps il faut respecter la
règle No 3, c’est-à-dire S P (TP* ) 0 .
Lors de la linéarisation de S (T ) en S c + S P T P on peut avoir :
• S P optimum donné par la pente de la courbe S = f (T ) en TP* (si cette pente est négative) ;
• S P > 0 , on a la divergence du procesus itératif (la règle No 3 n’est pas vérifiée) ;
• 0 > S P > S P optimum , il y a le risque de divergence ;
• SP S P optimum < 0 , convergence plus lente mais assurée.
En général, il y a plusieurs possibilités de linéarisation du terme source, illustrées par les
exemples suivants, [50] :
Exemple 1
Soit S = 5 4T donné. Alors plusieurs linéarisations sont possibles :
1. S c = 5 et S P = 4 , c’est la forme optimum recommandée ;
2. S c = 5 4T P* et S P = 0 , utile lorsque l’expression de S est plus compliquée, ce qui
n’est pas le cas ici ;
3. Sc = 5 + 7TP* et S P = 11 , permet d’avoir une pente négative plus grande et de
ralentir la convergence. Ceci est nécessaire s’il y a trop de non-linéarités dans le
problème, car assure la convergence.
Exemple 2
Soit S = 3 + 7T . Dans ce cas, les linéarisations suivantes sont possibles :
1. Sc = 3 et S P = 7 , on a la divergence car la règle No 3 n’est pas assurée ;
2. S c = 3 + 7T P* et S P = 0 , recommandé lorsqu’on a malgré tout S P 0 ;
3. S c = 3 + 9T P* et S P = 2 , la création artificielle de S P < 0 , qui ralentit la
convergence.
Exemple 3
Soit S = 4 5T 3 . Les linéarisations suivantes sont possibles :
1. S c = 4 5T P*3 et S P = 0 , dans ce cas, on ne profite pas de la relation connue entre
S et T ;
2. S c = 4 et S P = 5TP*2 , linéarisation correcte mais la pente n’est pas assez grande ;
3. La méthode recommandée. On cherche la tangente à la courbe S = f (T ) en TP*
47
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
*
dS S S*
= (4.14)
dT TP TP*
En utilisant la relation (4.15) pour le terme source, considéré comme exemple, on obtient :
S =4 5TP*3 + ( 15TP*2 TP )( )
TP* = 4 + 10T P*3 15TP*2TP , (4.16)
S c + S P TP 0, (4.17)
Sc
TM = . (4.18)
SP
Si la valeur de la pente S P est grande, la solution TM tende vers TP* , tandis que pour
une valeur plus petite de la pente, celle-ci implique un grand changement de T , de la valeur
TP* à la valeur TM . Dans certaines situations, le cas limite du terme source dominant est
utilisé pour linéariser un terme source fortement non-linéaire et pour lequel on ne connaît pas
la courbe analytique. Si l’on suppose que pour la valeur courante TP* on désire qu’à l’itération
suivante la variable inconnue ne dépasse pas la valeur limite TM , la linéarisation se fait
comme suit :
- on exprime la dérivée au point P , ainsi
*
dS S P*
= ; (4.19)
dT TM T P*
48
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
L’intégration de l’équation (4.1) sur un volume de contrôle (VC) autour d’un noeud
intérieur "i " donne :
dT dT
+ S i (Ti ) x = 0 . (4.21)
dx i +1 / 2 dx i 1/ 2
Ti + 1 Ti Ti Ti 1
i +1 / 2 i 1/ 2 + Si (Ti ) x = 0 i = 2, K , N 1, (4.22)
x x
où S i est la valeur moyenne de S i sur le volume de contrôle autour du point " i " .
En regroupant les termes l’équation (4.22) peut être écrite ainsi :
Ti 1 2Ti + Ti + 1 = Gi i = 2, K , N 1, (4.24)
où
S i ( x) 2
Gi = .
Pour les points 1 et N situés aux frontières du domaine de calcul on intègre l’équation
(4.1) sur un demi-volume de contrôle.
Les trois cas typiques de conditions aux limites rencontrées pour les problèmes de
conduction thermique sont :
1. Température imposée (connue) à la frontière (condition de type Dirichlet) ;
2. Densité du flux thermique imposée, donc dT / dx connue (condition de type Neumann) ;
3. Densité du flux thermique spécifiée par un coefficient d’échange ( h ) et une température
du fluide environnant ( T f ) ou par un flux radiatif (condition mixte ou de type Fourier)
dT h Tf( TN )
QN =
dx N
=
(T f4 T N4 ), (4.24a)
49
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
Dans ce cas, si la température T1 ou/et T N est connue, il n’est pas nécessaire d’écrire
une équation discrétisée supplémentaire au noeud 1 ou/et N . Quand on a une seule condition
de type Dirichlet (au noeud 1 ou N ) le nombre d’équations à résoudre est N 1 . Si l’on a
deux conditions aux limites de type Dirichlet (dans les points de frontière 1 et N ) alors le
nombre d’équations à résoudre est N 2 . Pour i = 2 de l’équation (4.23) on obtient :
S 2 ( x) 2
2T2 + T3 = T1 . (4.27)
( N 1) 1 / 2 T N 2 ( ( N 1) + 1 / 2 + ( N 1) 1 / 2 )T N 1 + ( N 1) + 1 / 2 T N = SN 1( x) 2 . (4.28)
SN 1( x)2
TN 2 2TN 1 = TN . (4.30)
En intégrant l’équation (4.1) sur le demi-volume de contrôle (VC) illustré à la figure 4.5
on obtient :
Fig. 4.5 Traitement d’une condition à la limite de type flux imposé (Neumann)
50
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
d dT
dV + SdV = 0 , (4.31)
VC
dx dx VC
N N
d dT
dx + Sdx = 0 , (4.32)
n
dx dx n
dT dT x
+ SN = 0. (4.33)
dx N dx n 2
dT
QN = , (4.34)
dx N
l’équation (4.33), en exprimant la dérivée au point " n " avec différences centrales, devient :
TN TN 1 x
QN n + SN = 0. (4.35)
x 2
En regroupant les termes dans l’équation (4.35), on obtient l’équation discrétisée, valable pour
le point de frontière N , pour la condition à la limite de type Neumann (flux imposé) :
S N ( x)2 QN x
TN 1 TN = + . (4.36)
2 N 1/ 2 N 1/ 2
Cet algorithme permet de calculer la solution d’un système linéaire lorsque la matrice
est tridiagonale. C’est notre cas, car les équations discrétisées dans le système linéaire
s’écrivent :
a i Ti = bi Ti + 1 + c i Ti 1 + di 1 i N, (4.37)
où le maillage est présenté à la figure 4.6.
51
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
c1 = 0 et bN = 0 . (4.38)
Dans cette étape on cherche les relations de type Ti = f (Ti + 1 ) sous la forme,
Ti = Pi Ti +1 + Qi , (4.40)
Ti 1 = Pi 1Ti + Qi 1 (4.41)
En regroupant les termes dans l’équation (4.42) sous la forme générale (4.41) on obtient
les coefficients Pi et Qi en fonction des coefficients Pi 1 et Qi 1 :
bi d i + ci Qi 1
Pi = Qi = . (4.43)
ai ci Pi 1 ai ci Pi 1
b1 d1
P1 = et Q1 = . (4.44)
a1 a1
Il est bien de préciser que les relations (4.44) sont obtenues si l’on remplace c1 = 0
dans les relations (4.43).
À la fin du processus de récurrence on constate que b N = 0 et donc PN = 0 et de
l’équation (4.40) on obtient :
TN = Q N . (4.45)
Résumé de l’algorithme
1. Calculer P1 et Q1 en utilisant les relations (4.44) ;
2. Calculer Pi et Qi , pour i = 2 ÷ N , avec les relations de récurrence (4.43) ;
3. Poser T N = Q N ;
4. Utiliser l’équation Ti = Pi Ti +1 + Qi de i = N 1 à 1 pour obtenir TN 1 , TN 2, ..., T1 .
4.1.7 Exemples
Exemple 1
On considère une barre cylindrique, sans source de chaleur, ayant l’aire transversale
A = 10 2 m 2 et la longueur L = 0.5 m . Les extrémités, A et B de la barre sont maintenues
aux températures constantes de 100 °C et de 500 °C respectivement.
Calculer la distribution de la température le long de la barre. On connaît la conductivité
thermique = 1000 W/mK .
Solution
La distribution de la température est gouvernée par l’équation :
d dT
= 0. (4.46)
dx dx
53
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
où aP = aE + aW aW = aE = .
x
Pour les noeuds 2 et 5 on utilise la même équation que pour un noeud intérieur mais on
tient compte que pour les noeuds voisins 1 et 6 les températures sont connues, TA et TB
respectivement.
Pour le noeud 2 l’équation discrétisée est la suivante :
a P TP = a E T E + aW T A . (4.48)
a P T P = aW TW + a E T B . (4.49)
Les équations à résoudre sont les suivantes (le nombre d’équations est égal à 4) :
1000
= = 10000 a P = aW + a E = 10000 + 10000 = 20000 ,
x 0 .1
2T2 = T3 + TA
2T3 = T2 + T4 (4.51)
2T4 = T3 + T5
2T5 = T4 + TB
54
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
'T2 $ '180 $
%T " %260"
% 3" = % " (4.52)
%T4 " %340"
% " % "
&T5 # &420#
La solution analytique de l’équation (4.46) est une distribution linéaire entre les valeurs
de la température des points A et B , c’est-à-dire :
TB TA
T ( x) = x + T A = 800 x + T A . (4.53)
L
À la figure 4.9 sont représentées les solutions analytique et numérique qui correspondent.
55
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
Exemple 2
On considère une plaque très longue d’épaisseur L = 20 mm , ayant la conductivité
thermique constante = 0.5 W/m/K et une source de chaleur uniforme, S = 1000 kW/m 3 .
Les faces de la plaque se trouvent à la température constante de 100 °C et 200 °C
respectivement.
En supposant que les dimensions de la plaque dans les directions “y” et “z” soient très
grandes et donc le gradient de la température est significatif dans la direction “x” seulement,
calculer la distribution de la température et comparer les résultats numériques avec la solution
analytique.
Solution
L’équation différentielle qui gouverne la distribution de la température est la suivante :
d dT
+ S = 0. (4.54)
dx dx
Le domaine d’analyse est divisé en six noeuds, comme dans l’exemple précédent, avec
x = 0.004 m . L’aire A = 1 est considérée dans le plan y z .
d dT
dV + SdV = 0 . (4.55)
VC
dx dx VC
' dT dT $
% A A "+S V =0 ; (4.56)
& dx e dx w#
' TE TP TP TW $
% eA wA " + SA x = 0 . (4.57)
& x x #
56
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
a P TP = aW TW + a E TE + b , (4.58)
w e
où : a P = aW + a E aW = aE = b = S x.
x x
Pour le noeud 2 on utilise la même équation discrétisée que pour un noeud intérieur (les
noeuds 3 et 4) mais on tient compte que le noeud voisin “W” corresponde au noeud “1” où la
température est connue TW ( T1 ( TA et passe comme un terme source supplémentaire.
L’équation discrétisée est donc :
a P TP = a E T E + b + aW T A . (4.59)
Pour le noeud 5 on utilise la même équation discrétisée que pour un noeud intérieur (les
noeuds 3 et 4) mais on tient compte que le noeud voisin “E” correspond au noeud “6” où la
température est connue TE ( T6 ( TB et passe comme un terme source supplémentaire.
L’équation discrétisée est :
a P T P = aW TW + b + a E T B . (4.60)
0 .5
= = 125 a P = aW + a E = 125 + 125 = 250 ,
x 0.004
57
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
'T2 $ '184 $
%T " %236"
% 3" = % " (4.62)
%T4 " %256"
% " % "
&T5 # &244#
'T TA S $
T ( x) = % B + (L x )" x + T A . (4.63)
& L 2 #
La comparaison des résultats numériques, obtenus avec la méthode des volumes finis,
avec la solution analytique est présentée au tableau 4.1 et à la fig. 4.11.
Numéro du noeud 2 3 4 5
x [m] 0.004 0.008 0.012 0.016
Solution numérique 184 236 256 244
Solution analytique 184 236 256 244
Erreur % 0 0 0 0
On constate que malgré un maillage très grossier les solutions numérique et analytique
correspondent. Le code source en Fortran est présenté à l’Annexe G.
58
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
Exemple 3
On considère une barre cylindrique (fig. 4.12) de l’aire A avec une extrémité maintenue
à la température constante de 100 °C ( TB ) et l’autre extrémité est isolée (le flux de chaleur est
nul). Sur le long de la barre il y a un échange de chaleur par convection dépendante de la
température. La température du milieu extérieur est de 20 °C .
Calculer la distribution de la température et comparer les résultats avec la solution
analytique.
On connaît : L = 1 m , hP /( A) = 25 m -2 .
Solution
L’équation diférentielle qui gouverne le transfert thermique dans ce cas est :
d dT
A hP(T T) ) = 0 , (4.64)
dx dx
cosh[n( L x )]
T ( x) = (TB T) ) + T) , (4.65)
cosh (nL )
59
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
d dT
n 2 (T T) ) = 0 où n 2 = hP / A . (4.66)
dx dx
d dT
dV n 2 (T T) )dV = 0 . (4.67)
VC
dx dx VC
La première intégrale de l’équation (4.67) sera traitée comme dans les exemples 1 et 2.
La deuxième intégrale, à cause du terme source, est évaluée en supposant que la quantité à
intégrer est localement constante sur chaque volume de contrôle et donc :
A
dT
dx
A
dT
dx
[n 2 (TP ]
T) ) A x = 0 . (4.68)
e w
TE TP TP TW
+ n 2 (T P T) ) x = 0 . (4.69)
x x
1 1 1 1
+ TP = TW + T E + n 2 xT) n 2 xT P . (4.70)
x x x x
a P T P = aW TW + a E T E + b , (4.71)
1
où aW = a E = a P = aW + a E SP b = S c = n 2 xT)
x
SP = n2 x .
a P TP = a E T E + b + aW T B , (4.72)
1
où aW = a E = a P = aW + a E SP b = S c = n 2 xT) .
x
60
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
d dT
dV n 2 (T T) )dV = 0 ; (4.73)
1 / 2VC
dx dx 1 / 2VC
dT dT x
A A n 2 (TP T) ) A = 0.
dx P dx w 2
dT
Parce que le flux dans le point P est nul ( q P = A = 0 ), on obtient l’équation
dx P
discrétisée suivante :
TP TW x
0 A n 2 (T P T) ) A = 0. (4.74)
x 2
a P TP = aW TW + b , (4.75)
1 x x
où aW = a P = aW SP b = Sc = n2 T) SP = n2 .
x 2 2
1 1
= =5 n 2 xT) = 25 0.2 20 = 100 n 2 x = 25 0.2 = 5 ,
x 0 .2
61
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
'T2 $ '50.569$
%T " %31.707"
% 3" % "
%T4 " = % 24.553" (4.77)
% " % "
%T5 " % 21.951"
%&T6 "# %& 21.301"#
Dans le tableau 4.2 on compare les résultats numériques avec la solution analytique.
Tableau 4.2 Comparaison : solution numérique – solution analytique
Noeud x [m] Numérique Analytique Erreur %
2 0.2 50.569 49.439 -2.28
3 0.4 31.707 30.853 -2.76
4 0.6 24.553 24.056 -2.06
5 0.8 21.951 21.663 -1.33
6 1.0 21.301 21.078 -1.05
La précision de la solution numérique peut être augmentée en utilisant un maillage plus fin.
Dans le tableau 4.3 et à la figure 4.14 sont présentés les résultats numériques et les erreurs
pour un maillage de 21 noeuds ( x = 0.05 m, quatre fois plus petit ). Le code source en Fortran
est présenté à l’Annexe G.
Tableau 4.3 Comparaison : solution numérique – solution analytique
Noeud x [m] Numérique Analytique Erreur %
2 0.2 49.515 49.439 -0.150
3 0.4 30.910 30.853 -0.180
4 0.6 24.088 24.056 -0.130
5 0.8 21.682 21.663 -0.087
6 1.0 21.092 21.078 -0.066
62
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
,T , ,T
-c p = +S, (4.78)
,t ,x ,x
t+ t t+ t t+ t
,T , ,T
- cp dV dt = dV dt + S dV dt . (4.79)
t VC
,t t VC
,x ,x t VC
e 't + t t+ t t+ t
,T $ ' ,T ,T $
% -c p dt " dV = % A A " dt + S V dt . (4.80)
% t
w&
,t "# t & ,x e ,x w# t
't + t ,T $
% -c p
,t
dt " dV = -c p T P ( T P0 ) V, (4.81)
%
VC & t "#
63
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
t+ t t+ t
' TW $
(
-c p T P T P0 ) x = % e
TE
xe
TP
w
TP
xw #
" dt + S x dt , (4.82)
t & t
t+ t
TP0 t
TP dt = TP t , (4.83)
t
( fTP + (1 f )TP0 ) t
' $
%-c p
x
t
+ f
xe
e
+ w
xw
"T P =
e
xe
[
fT E + (1 f )T E0 + w fTW + (1 f )TW0
xw
] [ ]
%& "#
(4.85)
' x $
+ %-c p (1 f ) e (1 f ) w "T P0 + S x
& t xe xw #
a P TP = aW fTW + (1 [ ]
f )TW0 + a E fT E + (1 [ f )T E0 ], (4.86)
[
+ a P0 (1 f )aW (1 ]
f )a E T P0 + b
x
où a P = f (aW + a E ) + a P0 a P0 = -c p ;
t
w e
aW = aE = b=S x .
xw xe
64
Chapitre 4 : Méthode des volumes finis appliquée aux problèmes de conduction thermique
Dans le cas du schéma explicite le terme source est linéarisé par l’expression
b = S c + S P TP0 . En remplaçant f = 0 dans l’équation (4.86) on obtient la discrétisation
explicite de l’équation de conduction thermique 1D instationnaire :
a P T P = aW TW0 + a E T E0 + a P0 [ (aW + aE ]
S P ) T P0 + S c , (4.87)
x w e
où a P = a P0 a P0 = - c p aW = aE = .
t xw xe
La règle No 2 n’est pas toujours satisfaite. Le coefficient de TP0 peut être regardé comme le
coefficient du “voisin” de TP dans la “direction” temps ou un coefficient voisin qui fait la
liaison entre les valeurs de T à l’instant t et celles à l’instant t + t . Pour que le coefficient
de TP0 soit positif il faut aP0 aW aE 0 . Dans le cas général la condition devient :
x e w
- cp 0. (4.88)
t xe xw
x 2 -c p ( x ) 2
- cp ou t (4.89)
t x 2
Si l’on note . = / - c p , on obtient le nombre de Fourier qui doit être inférieur ou égal à 1 / 2 :
. t 1
(4.90)
( x) 2 2
Remarques
• le critère de convergence du schéma utilisé, pour l’intégration dans le temps, résulte d’une
considération physique (la règle No 2).
• si l’on réduit x pour améliorer la précision spatiale, il faut diminuer beaucoup t
( / 1 /( x) 2 ).
65
MODÉLISATION NUMÉRIQUE DU TRANSFERT THERMIQUE • Métode des volumes finis
1 1 x
où aP = (aW + a E ) + a P0 SP a P0 = -c p .
2 2 t
w e 1
aW = aE = b = Sc + S P T P0
xw xe 2
À l’instant t + t plusieurs inconnues sont présentes dans l’équation (4.91), le schéma est
donc implicite et les équations doivent être résolues simultanément pour tous les noeuds à
chaque pas dans le temps.
Mathématiquement le schéma Crank-Nicolson est inconditionnellement stable, mais
numériquement la convergence vers une solution physiquement acceptable n’est pas assurée
(par exemple, des solutions oscillantes d’amplitude constante ou décroissante).
La règle No 2 est satisfaite uniquement lorsque
a E + aW x e w
a P0 = -c p 0. (4.92)
2 t 2 xe 2 xw
-c p ( x ) 2 . t
t ou 1. (4.93)
( x) 2
Remarque
La relation (4.93) est moins restrictive que la relation (4.90) associée au schéma
explicite. La précision du schéma Crank-Nicolson est de second ordre dans le temps donc
pour le même pas dans le temps la précision des résultats est plus grande que dans le cas du
schéma explicite.
a P T P = aW TW + a E T E + a P0 T P0 + S c , (4.94)
w e
où a P = a P0 + aW + a E SP aW = aE = .
xw xe
La règle No 2 est toujours vérifiée, donc le schéma totalement implicite (TI) est
inconditionnellement stable. La précision du schéma TI est de premier ordre dans le temps,
donc un petit pas dans le temps est nécessaire pour augmenter la précision des résultats.
66