REPUBLIQUE ALGERIENNE DEMOCRATIQUE ET POPULAIRE

Ministère de l’Enseignement Supérieur et de la Recherche Scientifique

Centre Universitaire Abdel Hafid Boussouf – Mila

Institut des Sciences et de la Technologie

Département des Sciences et Techniques

Filière: Génie Mécanique

Spécialité : Energétique

POLYCOPIE
Méthodes numériques I
Méthodes des différences finies
Ce polycopié est destiné aux étudiants de la 1ére année Master,
Mécanique Energétique

Auteur : Dr. BOUCHOUCHA Abd El malik

(Maître de Conférences B)
Sommaire Page

CHAPITRE I : Formulation mathématique de l‟équation de la conduction de chaleur

unidimensionnelle (1-D) et bidimensionnelle (2-D)

I.1INTRODUCTION…………………………………………………………………........ 1
I- Equation de la conduction de chaleur unidimensionnelle (1-D)………………………... 1
I- 1.Equation de la conduction de chaleur dans une paroi plane large…………………….. 1
I.2 Equation de la conduction de chaleur dans un cylindre long………………………….. 4
I.3.Equation de la conduction de chaleur dans une sphère…………………………........... 8
I.4.Equation générale de la conduction de chaleur…………………………………… 9
I.5.Conditions aux limites…………………………………………………………………. 11
I.5.1 Conditions de Direchlet „Température imposée‟…………………………………….. 12
I.5.2 Flux imposé „condition de Neuman‟ ……………………………………................... 12
I.5.3 Condition aux limites du 3eme type.‟ Condition de Cauchy…………….................... 12
CHAPITRE I : Classification des équations aux dérivées partielles.

II.1 Définition…………………………………………………………………………….. 15
II.2. Classification mathématique des EDP linéaires du second ordre (cas de deux

variables indépendantes) :………………………………………………………………… 15

II.3. Classification mathématique dans le cas général (n variables indépendantes)………. 16
II.4. Classification physique des EDP……………………………………………………... 17
II.4.1.‐ Equation de la chaleur……………………………………………………………... 17
II.5. Résolution des systèmes d‟équations linières………………………………………… 19
II.5.2. Méthode de Gauss-seidel…………………………………………………………… 21
II.5.3 Amélioration de la convergence en utilisant la relaxation…………………….. …... 24
 II.5.4 Méthode de Jacobi…………………………………………………………………... 24
II.5.5. Méthode de résolution directe ……………………………………………………... 25
II.5.5.1 Méthode d‟élimination de Gauss………………………………………………….. 25
II.5.5.2 Algorithme de Thomas (Méthode direct)…………………………………………. 27
Chapitre III : Méthode des différences finies

III.1 Approximation des opérateurs différentiels………………………………………...... 31

III.2 Maillage………………………………………………………………………………. 31
III.4 Résume……………………………………………………………………………….. 33
III.5 Etude numérique de la conduction unidimensionnel stationnaire……………………. 33
Chapitre IV: Discrétisation de l‟équation de la conduction de chaleur (2-D)

IV.2.Discrétisation de l‟équation de la conduction de chaleur (2-D) en Régime

permanent par la méthode des différences finies………………………………………….. 39

Chapitre V. Discrétisation de l‟équation de la conduction de chaleur en régime
transitoire par la méthode des différences finies
V.1 Discrétisation de l‟équation de la conduction de chaleur en régime transitoire par la

méthode des différences finies……………………………………………………… 47

V.2.Conduction de chaleur transitoire dans une paroi plane……………………………… 48
V.3. Formulation explicite des différences finies…………………………………………. 49
V.3. 1 Critères de la stabilité pour la méthode explicite…………………………………... 50
a)la Méthode explicite……………………………………………………………………... 55
b) la Méthode implicite……………………………………………………………………. 55
V.4.Conduction thermique transitoire 2-D………………………………………………… 54
I.2.1 I.5.CONCLUSION............................................................................................................... 47
 Chapitre 1
FORMULATION MATHEMATIQUE DE L’EQUATION DE LA
CONDUCTION DE CHALEUR UNIDIMENSIONNELLE (1-D) ET
BIDIMENSIONNELLE (2-D)
I.1 Introduction

Le processus de transfert de chaleur par la conduction s‟appuie sur un milieu matériel

sans mouvement de matière est dû à des phénomènes physiques microscopiques
(agitation des atomes ou des molécules, flux d‟électrons libres…). Il peut être vu
comme un transfert d‟énergie des particules les plus énergétiques (les particules
chaudes qui ont une énergie de vibration élevée) vers les particules les moins
énergétiques (les particules froides d‟énergie de vibration moins élevée), dû aux
collisions entre particules. Dans les solides, le transfert d‟énergie peut également se
produire sous l‟effet du déplacement d‟électrons libres dans le réseau cristallin (par
exemple pour les métaux). Ainsi les bons conducteurs d‟électricité sont en général
également de bons conducteurs de la chaleur.
I.2Equation de la conduction de chaleur unidimensionnelle (1-D)
I.2.1Equation de la conduction de chaleur dans une paroi plane large :
.
Egen 𝐥′ é𝐥𝐞𝐦𝐞𝐧𝐭 𝐝𝐞 𝐯𝐨𝐥𝐮𝐦𝐞

Q.x Q.x − Q.x+dx

x dx

FigI.1 conduction de chaleur unidimensionnelle à travers un élément de volume

dans une paroi plane large.
Nous considérons un élément mince d‟épaisseur dx dans une paroi plane et
large (Fig I-1) supposé que la densité ρ, de la paroi, la chaleur spécifique est C, et la
surface de la paroi normale à la direction du transfert de chaleur est A. En l'absence de

1
génération de chaleur, un bilan énergétique sur cet élément mince d'épaisseur ∆x
pendant un petit intervalle de temps ∆t peut être exprimé comme
Taux de la conduction Taux de la conduction
de − de
chaleur à x chaleur à x + dx
Taux de generation
de Taux de changement
′ de
+ chaleur à l intérieur = ′
de l énergie contenue dans
l′ élement l′ élement de volume
de voulume
∆E élément
Q.x − Q.x+dx +Egen
.
. element = …………………… I. 1
∆t
∆Eélément = Et+∆t − Et = mc Tt+∆t − Tt = ρ. dx. A. C. Tt+∆t −Tt I. 2
∆Eélément = e.gen élément . Vélément = e.gen élément . A. dx … … … … … . I. 3
En remplaçant I. 2 et I. 3 dans l‟equation I. 1
T t+∆t −T t
Q.x − Q.x+dx + e.gen élément . A. dx= ρ. dx. A. C. … … … … . . I. 4
dt
e.gen = e.gen élément

.
∆Eélément = 𝑒𝑔𝑒𝑛 . 𝐴. 𝑑𝑥 … … … … … … … … … … … … … … … … … … I. 5

. T t+∆t −T t
𝑄𝑥. − 𝑄𝑥+𝑑𝑥 .
+ 𝑒𝑔𝑒𝑛 . 𝐴. 𝑑𝑥 = 𝜌. 𝑑𝑥. 𝐴. 𝐶. … … … … … … … ..(1.1)
dt

. T t+∆t −T t
−𝑄𝑥. + 𝑄𝑥+𝑑𝑥 .
+ 𝑒𝑔𝑒𝑛 . 𝐴. 𝑑𝑥 = 𝜌. 𝑑𝑥. 𝐴. 𝐶. … … … … … … . …(1.1)
dt
On divise l‟équation (I.4) par A.dx
.
1 𝑄𝑥+𝑑𝑥 −Qx . T t+∆t −T t
− + 𝑒𝑔𝑒𝑛 = 𝜌. 𝐶. … … … … … … … … … … … … … . I. 6
𝐴 𝑑𝑥 dt
Qx+ dx −Qx 𝜕𝑄 𝜕 𝜕𝑇
= = −𝐾. 𝐴 … … … … … … … … … … … … … … … . I. 7
𝑑𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥
La loi de fourrier de la conduction de chaleur
Limite dx→ 0
En prenant, dx→ 0, dt→ 0 l‟équation I. 6 devient.
𝑇𝑡+∆𝑡 − 𝑇𝑡 𝜕𝑇
lim =
∆→0 𝑑𝑡 𝜕𝑡
1 ∂ ∂T ∂T
. K. A = e.gen = ρ. C. … … … … … … … … … … … … … … … . . I. 8
A ∂x ∂x ∂t

2
En notant que la surface de la paroi plane est constante. L‟équation de la conduction
de chaleur transitoire unidimensionnelle (1-D) dans une paroi plane devient

*Conductivité thermique variable

∂ ∂T ∂T
K. A = e.gen = ρ. C. … … … … … … … … … … … … … … … … . … I. 9
∂x ∂x ∂t

*Conductivité thermique constante

∂2T e .gen ρ.C ∂T ∂2T e .gen 1 ∂T
+ = . = + = . … … … … … … … … … … … … . . I. 10
∂2t K k ∂t ∂2t K α ∂t
*Régime permanent :
.
𝜕 2 𝑇 𝑒𝑔𝑒𝑛
+ =0
𝜕2 𝑡 𝐾 … … … … … … … … … … … … … … … … … … … … … . I. 11
𝜕2 𝑇
=0
𝜕2 𝑡
*Régime transitoire pas de génération de chaleur
∂2T 1 ∂T
= … … … … … … … … … … … … … … … … … … … … … … … … I. 12
∂2t 2 ∂t
.
Génération de chaleur (𝑒𝑔𝑒𝑛 =0)
*Régime permanent pas de génération de chaleur
∂2T
= 0… … … … … … … … … … … … … … … … … … … … … … … … … … I. 13
∂2t
∂T
( = 0, e.gen =0)
∂t
Exemple(I.1)

Une paroi plane large est soumise à une température spécifiée sur la surface gauche et
à une convection sur la surface droite. La formulation mathématique, la variation de la
température et le taux de transfert de chaleur doivent être déterminés pour un transfert
de chaleur unidimensionnel régulier.

Hypothèses k
T1=80°C
T=15°C
1- La conduction thermique est stable et unidimensionnelle. .
A=20 m2m2
A=20 2
L=0.4 m h=24 W/m2 .°C
h=24 W/m .°C
2 -La conductivité thermique est constante.
A=20 m2
3- Il n'y a pas de génération de chaleur.
x
Propriétés La conductivité thermique est donnée comme étant
k = 2,3 W / m. ° C. 2
h=24 W/m .°C

3
En prenant la direction normale à la surface du mur comme étant la direction x avec x
= 0 à la surface gauche, la formulation mathématique de ce problème peut être
exprimée comme suit:
d 2T dT ( L)
 0 Et T (0)  T1  80C  k  h[T ( L)  T ]
dx 2 dx

(b) Intégration de l'équation différentielle deux fois par rapport aux rendements x
dT
 C1 T ( x)  C1x  C2
dx

Où et C2 sont des constantes arbitraires. L'application des conditions aux limites

donne

x = 0: T (0)  C1  0  C2  C2  T1

h(C2  T ) h(T  T )
x = L:  kC1  h[(C1 L  C2 )  T ]  C1    C1   1 
k  hL k  hL

En se substituant à la solution générale, la variation de température est

déterminée comme :
h(T1  T )
T ( x)   x  T1
k  hL
(24 W / m2  C)(80  15) C
 x  80 C
(2.3 W / m C)  (24 W / m2  C)(0.4 m)
 80  1311
.x

(c) Le taux de conduction thermique à travers le mur est

dT h(T  T )
Q wall  kA  kAC1  kA 1
dx k  hL
(24 W/m 2  C)(80  15)C
 (2.3 W/m  C)(20 m 2 )  6030 W
(2.3 W/m  C)  (24 W/m 2  C)(0.4 m)

Notez que dans des conditions stables, le taux de conduction thermique à travers un
mur est constant.

4
 I.2.2Equation de la conduction de chaleur dans un cylindre long

.
Egen . élément

𝑞𝑟.
l′ élément de volume
.
𝑞𝑟+𝑑𝑟 𝑟

Fig I.2 conduction de chaleur unidimensionnelle à travers un élément de volume

dans un cylindre long

On considère un élément de coque cylindrique mince d'épaisseur dr dans un cylindre

long. La densité du cylindre est ρ, la chaleur spécifique est C, et la longueur est L. la
surface de cylindre à la direction normale de transfert de chaleur à n'importe quel
endroit est où r est la valeur du rayon à cet en droit est A= π2rL. Notez que la surface
de transfert de chaleur A dépend de r dans ce cas, et donc elle varie selon
l'emplacement. Un bilan énergétique sur cet élément de coque cylindrique mince
d'épaisseur dr pendant un petit intervalle de temps ∆t peut être exprimé comme:

En remplaçant l‟équation (I, 15) et (I, 16) dans l‟Eq (I,14)

Taux de generation
de
Taux de la conduction Taux de la conduction
− chaleur à l′ intérieur
de de +
de
chaleur à x chaleur à x + dx ′
l élément
de volume
Taux de changement
de
= ′
l énergie contenue dans
l′ élément de volume
. ∆𝐸é𝑙é𝑚𝑒𝑛𝑡
𝑄𝑟. − 𝑄𝑟+𝑑𝑟 .
+ 𝐸𝑔𝑒𝑛 . é𝑙é𝑚𝑒𝑛𝑡 = … … … … … … … … … I. 14
∆𝑡

∆Eélément = Et+∆t − Et = mc Tt+∆t − Tt = ρ. dr. A. C. Tt+∆t − Tt … . I. 15

. .
∆𝐸é𝑙é𝑚𝑒𝑛𝑡 = 𝑒𝑔𝑒𝑛 é𝑙é𝑚𝑒𝑛𝑡 . 𝑉é𝑙é𝑚𝑒𝑛𝑡 = 𝑒𝑔𝑒𝑛 é𝑙é𝑚𝑒𝑛𝑡 . 𝐴. 𝑑𝑟 … … … … I. 16

5
 T t+∆t −T t
Q.r − Q.r+dr + Egen
.
. A. dr = ρ. C. A. dr. … … … … … … … . . I. 17
dt

On divise l‟Eq (I,17) par A.dr

1 𝑄 . 𝑟+𝑑𝑟 −𝑄 . 𝑟 . 𝑇𝑡+∆𝑡 −𝑇𝑡

− + 𝑒𝑔𝑒𝑛 . = 𝜌. 𝐶. … … … … … … … … … … … … I. 18
𝐴 𝑑𝑟 𝑑𝑡

Dans la limite, dt→ 0 , dr→0, l‟équation (I, 18) devient

Q .r +dr −Q .r 𝜕𝑄 . T t+∆t −T t 𝜕𝑇
= , =
𝑑𝑟 𝑑𝑟 dt 𝜕𝑡

Et en utilisant la loi de Fourier

1 𝜕 𝜕𝑇 . 𝜕𝑇
. 𝐾. 𝐴 + 𝑒𝑔𝑒𝑛 = 𝜌. 𝐶. … … … … … … … … … … … I. 19 − I. 15
𝐴 𝜕𝑟 𝜕𝑟 𝜕𝑡
La loi de Fourier
Q .r +dr −𝑄 . r 𝑄 . r 𝜕 𝜕𝑇
lim𝑛→∞ = = −𝐾. 𝐴 … … … … … … … … … … … … I. 20
𝑑𝑟 𝑑𝑟 𝜕𝑟 𝜕𝑟
L‟équation de la conduction de chaleur transitoire unidimensionnelle dans un
cylindre devient

*conductivité thermique variable

1 ∂ ∂T ∂T
K. r + e.gen = ρ. C. … … … … … … … … … … … … … … … … . I. 21
r ∂r ∂r ∂t

*Conductivité thermique constante

1 ∂ ∂T e .gen 1 ∂T
r + = . … … … … … … … … … … … … … … … … … … . I. 22
r ∂r ∂r K α ∂t

𝒌
α=
𝜌. 𝐶
*Régime parementant
1 𝜕 𝜕𝑇 .
𝐾. 𝑟 + 𝑒𝑔𝑒𝑛 =0…………… … … … … … … … … … … … … … . I. 23
𝑟 𝜕𝑟 𝜕𝑟
𝜕𝑇
=0
𝜕𝑡
*Régime transitoire pas de génération de chaleur (𝒆.𝒈𝒆𝒏 =0)
1 𝜕 𝜕𝑇 1 𝜕𝑇
𝑟 = . … … … … … … … … … … … … … … … … … … … … … . I. 24
𝑟 𝜕𝑟 𝜕𝑟 𝛼 𝜕𝑡

*Régime parementant pas de génération de chaleur

𝜕 𝜕𝑇
𝐾. 𝑟 = 0 … … … … … … … … … … … … … … … … … … … … … I. 25
𝜕𝑟 𝜕𝑟
𝜕𝑇 .
=, 𝑒𝑔𝑒𝑛 =0
𝜕𝑡

6
Exemple (I.2)
Les surfaces supérieure et inférieure d'une tige cylindrique solide sont maintenues à
des températures constantes de 20°C et 95°C tandis que la surface latérale est
Insulated

D = 0.05 m T2=95°C
Parfaitement isolée. Le taux de transfert de chaleur T1=25°C

à travers la tige doit être déterminé pour les cas de tige de cuivre, L=0.15 m

D‟acier et de granit.

Hypothèses

1- La conduction thermique est stable et unidimensionnelle. .

2 -La conductivité thermique est constante.

3- Il n'y a pas de génération de chaleur.

Les propriétés des conductivités thermiques sont données à k = 380 W/m .°C pour le
cuivre, k = 18 W/m.°C pour l'acier et k = 1,2 W/m.°C pour le granit.

Notant que la surface de transfert de chaleur (la surface du cylindre à la direction

normale du transfert de chaleur) est constante, le taux de transfert de chaleur le long
de la tige est déterminé à partir de :

T T
Q  kA 1 2
L

Où L = 0,15 m et la surface de transfert de chaleur A est

A  D2 / 4   (0.05 m) 2 / 4  1964

.  103 m2

Ensuite, le taux de transfert de chaleur pour chaque cas est déterminé comme suit:
T T (95  20) C
(a) Le cuivre: Q  kA 1 2  (380 W / m C)(1.964  103 m2 )  373.1 W
L 0.15 m

T T (95  20) C
(b) l'acier: Q  kA 1 2  (18 W / m C)(1.964  103 m2 )  17.7 W
L 0.15 m
T T (95  20) C
(c) granit: . W / m C)(1.964  103 m2 )
Q  kA 1 2  (12  1.2 W
L 0.15 m

7
Discussion:
Le taux constant de conduction thermique peut différer de plusieurs ordres de
grandeur, selon la conductivité thermique du matériau.

I.2.3 Equation de la conduction de chaleur dans une sphère

Elément de volume
Fig I.3 conduction de chaleur unidimensionnelle à travers un élément de volume
d’une sphère.
On considère un mince élément de coque sphérique d'épaisseur dr dans une sphère. La
densité du cylindre est ρ, la chaleur spécifique est C et la longueur est L. la surface de
sphère à la direction normale de transfert de chaleur à n'importe quel endroit est où r
est la valeur du rayon à cet endroit est A=4πR2où r est la valeur du rayon à cet en
droit. Notez que la surface de transfert de chaleur A dépend de r dans ce cas, et donc
elle varie selon l'emplacement. Lorsqu'il n'y a pas de génération de chaleur .Un bilan
énergétique sur cet élément de coque cylindrique mince d'épaisseur ∆r pendant un
petit intervalle de temps ∆t peut être exprimé comme:

En utilisant A= 4πr2 ou lien de A=2ΠrL.

L‟équation de la conduction de chaleur transitoire unidimensionnelle pour une sphère

*conductivité thermique variable :

1 𝜕 𝜕𝑇 𝜕𝑇
𝑟2. 𝐾 .
+ 𝑒𝑔𝑒𝑛 = 𝜌. 𝐶. … … … … … … … … … … … … … … . . . I. 26
𝑟 2 𝜕𝑟 𝜕𝑟 𝜕𝑡

*conductivité thermique constante :

.
𝑒𝑔𝑒𝑛
1 𝜕 𝜕𝑇 1 𝜕𝑇
𝑟2 + = . … … … … … … … … … … … … … … … … … I. 27
𝑟2 𝜕𝑟 𝜕𝑟 𝑘 𝛼 𝜕𝑡

𝒌
𝛂=
𝜌. 𝐶

8
*Régime parementant
.
𝑒𝑔𝑒𝑛
1 𝜕 𝜕𝑇
𝑟2. 𝐾 + = 0 … … … … … … … … … … … … … … … … … I. 28
𝑟 2 𝜕𝑟 𝜕𝑟 𝑘
∂T
=0
∂t
*Régime transitoire pas de génération de chaleur (𝒆.𝒈𝒆𝒏 = 𝟎)
1 𝜕 𝜕𝑇 1 𝜕𝑇
𝑟2 = . … … … … … … … … … … … … … … … … … … … … I. 29
𝑟 2 𝜕𝑟 𝜕𝑟 𝛼 𝜕𝑡

*Régime parementant pas de génération de chaleur

∂T .
=, 𝑒𝑔𝑒𝑛 =0
∂t
𝜕 2 𝜕𝑇
𝑟 = 0 … … … … … … … … … … … … … … … … … … … … … … I. 30
𝜕𝑟 𝜕𝑟
𝜕𝑇 𝑑𝑇
Ou 𝑟 2 +2 =0
𝜕𝑟 𝑑𝑟

I.4.Equation général de la conduction de chaleur

Dans ce paragraphe nous développons l‟équation générale gouvernante dans tel
système en coordonnées rectangulaire cylindrique et sphérique :

*Coordonnées rectangulaires (cartésiennes)

∂ ∂T ∂ ∂T ∂ ∂T ∂T
k + k + k + e.gen = ρ. C. … … … … … … … … … … . . I. 31
∂x ∂x ∂y ∂y ∂z ∂z ∂t
*Conductivité thermique constante
.
𝜕 2 𝑇 ∂2T ∂2T 𝑒𝑔𝑒𝑛 1
+ + + =𝛼 .𝜕𝑇 ……………………………………………... I. 32
𝜕2𝑥 ∂2y ∂2z 𝐾 𝜕𝑡

𝒌
𝛂 = 𝜌.𝐶 Équation de fourrier

*Régime permanent
𝜕𝑇
. =0
𝜕𝑡
.
𝜕2𝑇 𝜕2𝑇 𝜕2𝑇 𝑒𝑔𝑒𝑛
+ + + = 0… … … … … … … … … … … … … … … … … … … . I. 33
𝜕2𝑥 𝜕2𝑦 𝜕2𝑧 𝐾
Equation de poisson
*Régime transitoire pas de génération de chaleur
𝜕2𝑇 𝜕2𝑇 𝜕2𝑇 1 𝜕𝑇
+ + = . …………………………………….. I. 34
𝜕2𝑥 𝜕2𝑦 𝜕2𝑧 𝛼 𝜕𝑡

9
Équation de diffusion
*Régime permanent pas génération de chaleur :

𝜕2𝑇 𝜕2𝑇 𝜕2𝑇

+ + = 0 … … … … … … … … … … … … … … … … … … … … I. 35
𝜕2𝑥 𝜕2𝑦 𝜕2𝑧

Equation de la place

*Coordonnées cylindriques
1 ∂ ∂T 1 ∂ ∂T ∂ ∂T ∂T
2
k.r + 2 k + k .r + e.gen = ρ. C. … … … … I. 36
r ∂r ∂r r ∂Φ ∂Φ ∂Z ∂Z ∂t
*Coordonnées sphériques
1 𝜕 𝜕𝑇 1 𝜕 𝜕𝑇 1 𝜕 𝜕𝑇
𝑘 . 𝑟2 + . 𝑘 + . 𝑘 sin 𝜃 .
+𝑒𝑔𝑒𝑛 =
𝑟2 𝜕𝑟 𝜕𝑟 𝑟 2 .sin 2 𝜃 𝜕𝛷 𝜕𝛷 𝑟 2 .sin 2 𝜃 𝜕𝜃 𝜕𝑟
𝜕𝑇
𝜌. 𝐶. … … … … … … … … … … … … … … … … … … … … … … … … … I. 37
𝜕𝑡
Exemple (I.3)

Un récipient sphérique est soumis à une température spécifiée sur la surface intérieure
et à une convection sur la surface extérieure. La formulation mathématique, la
variation de la température et le taux de transfert de chaleur doivent être déterminés
pour un transfert de chaleur unidimensionnel régulier.
Hypothèses
1- La conduction thermique est stable et unidimensionnelle car il n'y a pas de
changement avec le temps et il y a une symétrie thermique autour du point médian.
2- La conductivité thermique est constante.
3- Il n'y a pas de génération de chaleur
Les propriétés de la conductivité thermique est donnée comme étant k = 30 W/m. °C.
Notant que le transfert de chaleur est unidimensionnel dans la direction radiale r, la
formulation mathématique de ce problème peut être exprimée comme suit:

10
 d  2 dT  k T1
r 0
dr  dr 
r1 T
r2
et T (r1 )  T1  0C h

dT (r2 )
k  h[T (r2 )  T ]
dr

(b) L'intégration de l'équation différentielle une fois par rapport à r donne

dT
r2  C1
dr

En divisant les deux côtés de l'équation ci-dessus par r pour l'amener à une forme
facilement intégrable puis en intégrant,
dT C1

dr r 2

C1
T (r )    C2
r

Où C1 et C2 sont des constantes arbitraires. L'application des conditions aux limites

donne
C1
r = r1: T (r1 )    C2  T1
r1

C1  C 
r = r2:  k 2
 h  1  C2  T 
r2  r2 

Résoudre pour donne simultanément

r2 (T1  T ) C1 T1  T r2
C1  and C2  T1   T1 
r k r1 r k r1
1 2  1 2 
r1 hr2 r1 hr2

En se substituant à la solution générale, la variation de température est déterminée

comme étant

C1 C  1 1 T1  T  r2 r2 
T (r )    T1  1  C1     T1      T1
r r1  r1 r  r k  r1 r 
1 2 
r1 hr2
(0  25)C  2.1 2.1 
     0C  29.63(1.05  2.1 / r )
2.1 30 W/m C  2 r 
1 
2 (18 W/m2  C)( 2.1 m)

(c) Le taux de conduction thermique à travers le mur est

11
 dT C r (T  T )
Q  kA  k (4r 2 ) 21  4kC1  4k 2 1
dx r r k
1 2 
r1 hr2
(2.1 m)(0  25)C
 4 (30 W/m  C)  23,460 W
2.1 30 W/m  C
1 
2 (18 W/m 2  C)(2.1 m)
I.3 Conditions aux limites
Une des caractéristiques des équations de type elliptique c‟est qu‟ils nécessitent des
conditions aux limites aux frontières du domaine d‟étude. Nous présentons en ce qui
suit les conditions aux limites spécifiques à l‟équation de l‟énergie. Il existe 03 types
de conditions aux limites.

I.3.1 Conditions de Direchlet ‘Température imposée’

Figure I.4 Conditions aux limites du 1er type.

La valeur de la grandeur à déterminée est connue aux différentes limites (frontières)

de la géométrie du domaine d‟étude. Cette condition aux limite est facile à
programmé, par contre il est difficile de maintenir tout une surface à une température
fixe au laboratoire.
I.3.2 Flux imposé ‘condition de Neuman’

Aux frontières la densité du flux est connue (voir figure I.5) lorsque la paroi est isolée
le flux thermique dans ce cas est nul. Cette condition est très facile à mettre en œuvre
au laboratoire.

12
 Figure I.5 Conditions aux limites du 2ertype.

I.3.3 Condition aux limites du 3eme type.’ Condition de Cauchy

Lorsque les frontières du domaine d‟étude sont en contact avec un fluide en

mouvement dont on connait le coefficient convectif et la température, la condition aux
limites imposée dans ce cas est la continuité des flux convectif et diffusif à l‟interface
(frontière solide/fluide).

Figure I.6 Conditions aux limites du 3er type.

La condition à l‟interface limite s‟écrit comme suit :

Ax=0 flux convectif = flux diffusif
L‟expression de la densité du flux convectif est donnée par la loi de refroidissement
de Newton :

𝑄 𝑜𝑢 𝛷0 = 𝑕0 𝑇0 −𝑇∞ … … … … … … … … … … … … … … … … I. 38

L‟expression du flux diffusif est donnée par la loi de Fourrier :

𝜕𝑇
𝑄 𝑜𝑢 𝛷0 = -k … … … … … … … … … … … … … … … … … … I. 39
𝜕𝑥

13
A l‟interface

𝜕𝑇
𝛷0 = -k =𝑕0 𝑇0 −𝑇∞ … … … … … … … … … … … … … … … … … … I. 40
𝜕𝑥

Noter qu‟à partir de la relation I. 40 on détermine le nombre de Nusselt comme suit :

𝑕 0 𝐿𝐶
𝑁𝑢 = … … … … … … … … … … … … … … … … … … … … … … … … . I. 41
𝑘

Avec :
LC : longueur caractéristique.
Donc en réécrivant l‟équation I. 41 , on aura :
𝑑𝑇
𝑕0 −
𝑑𝑥
= … … … … … … … … … … … … … … … … … … … … … … I. 42
𝐾 𝑇0 −𝑇∞
Multipliant les deux membres de l‟équation I. 42 par la longueur caractéristique LC
𝑑𝑇
𝑕 0 𝐿𝐶 − 𝐿
𝑑𝑥 𝐶
𝑁𝑢 = = … … … … … … … … … … … … … … … … … … … … … I. 43
𝑘 𝑇0 −𝑇∞

14
 Chapitre II

CLASSIFICATION DES EQUATIONS AUX DERIVEES

PARTIELLES.
II.1Définition :

*Une équation aux dérivées partielles (EDP) est une relation faisant intervenir les
variables indépendantes𝑥1 ,𝑥2 … … … … .𝑥1𝑛 , la fonction 𝑓 et ses dérivées partielles.
Par exemple, si 𝑓 est une fonction de deux variables, une EDP peut s‟écrire par la
relation :

𝜕𝑓 𝜕𝑓 𝜕 2 𝑓 𝜕 2 𝑓 𝜕2𝑓 𝜕2𝑓 𝜕3𝑓 𝜕3𝑓 𝜕3𝑓 𝜕3𝑓

F 𝑥, 𝑦, , , 2
, 2
, , , 3
, 3
, 2
, , . . … … . . II. 1
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦 𝜕𝑥𝜕𝑦 𝜕𝑦𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥𝜕 𝑦 𝜕𝑥 2 𝜕𝑦

* On appelle ordre de l‟EDP l‟ordre le plus élevé des dérivées partielles intervenant
dans l‟EDP, par exemple :

𝜕3𝑓 𝜕2𝑓 𝜕2𝑓 𝜕𝑓

+3 +𝑥 + +𝑓 + 𝐶 = 0 est d‟ordre 3 … … … … … II. 2
𝜕𝑋 2 𝜕𝑦 𝜕𝑥 2 𝜕𝑦 2 𝜕𝑥

2 2
𝜕2𝑓 𝜕2𝑓 𝜕2𝑓
− + − 𝑐 = 0 Est d‟ordre 2..… … … … … … … … II. 3
𝜕𝑋 2 𝜕𝑦 2 𝜕𝑥𝜕𝑦

L‟EDP est dite linéaire si 𝐹 est linéaire par rapport à ses arguments 𝑓 et ses dérivées
partielles, et si les coefficients qui les lient ne dépendent que de (𝑥, 𝑦) ; sinon elle est
non linéaire. Par exemple, l‟EDP du second ordre :

𝜕2𝑓 𝜕2𝑓 𝜕2𝑓 𝜕𝑓 𝜕𝑓

𝑎1 2
+𝑎2 2
+ 𝑎3 + 𝑎4 + 𝑎5 + 𝑎6 𝑓 + 𝑎7 = 0 … … ….… II. 4
𝜕𝑋 𝜕𝑦 𝜕𝑥𝜕𝑦 𝜕𝑥 𝜕𝑦

Est linéaire si les 𝑎i ne dépendent que de (𝑥, 𝑦).

II.2. Classification mathématique des EDP linéaires du second ordre (cas de deux
variables indépendantes) :

𝜕2𝛷 𝜕2𝛷 𝜕2𝛷 𝜕𝛷 𝜕2𝛷

𝐴 +𝐵 +𝐶 +D + 𝐸 + 𝐸𝛷 + 𝐺 𝑥, 𝑦 = 0 … … … II. 5
𝜕𝑥2 𝜕𝑥𝑦 𝜕2𝑧 𝜕𝑥 𝜕𝑦

Généralement les équations ou dérivées partielles sont classe on trois catégories

appelles : elliptique, parabolique, hyperbolique.

15
pour illustrer cette classification nous considérons l‟équation au dérivées partielles du
deuxième ordre au deux variables in dépendants x, y donner par Forsyth et Asonen
1967. ici nous supposons une équation linéaire qui est A, B, C, D, E et F sont fonction
de (x, y) mais ne dépend pas de la variable 𝞥
La variable 𝞥 = Température, pression, densité ou vitesse.
L‟équation de (1) dépend de A, B, C, D, E, F.
L‟équation au dérivée partielle de (1) un point (x,y) est appel :
Lorsque la quantité Δ = (b2 -4ac) < 0 l‟équation (2) est dite du type elliptique.
*Lorsque la quantité Δ = (b2 -4ac) = 0 l‟équation (2) est dite du type parabolique.
*Lorsque la quantité Δ = (b2 -4ac) >0 l‟équation (2) est dite du type hyperbolique
Cette classification est purement mathématique, une particularité des équations de
type elliptique c‟est qu‟elles nécessitent des conditions aux limites à toutes les
frontières du domaine d‟étude.
II.3. Classification mathématique dans le cas général (n variables
indépendantes) :
Si 𝑓 est une fonction de 𝑛 variables indépendantes, les EDP linéaires du second ordre
Sont du type :
𝑛 𝜕2𝑓 𝑛 𝜕𝑓
𝑖=1 𝑎𝑖 𝑥 1 …..𝑥 𝑛 𝜕𝑥 2 + 𝑖=1 𝑏𝑖 𝑥 1 …..𝑥 𝑛 𝜕𝑥 +c 𝑥1 … . . 𝑥𝑛 𝑓 + 𝑑 𝑥1 … . . 𝑥𝑛 =
𝑖 𝑖
0 … … . … … … … … … … … … … … … … … … … … … … … … … … … … … … … II. 6

* Si tous les 𝑎i sont non nuls et de même signe, l‟EDP est de type elliptique.
* Si tous les 𝑎i sont non nuls et sont ; à une exception près, de même signe, l‟EDP est
de type hyperbolique.
* Si un seul des 𝑎i est nul (noté 𝑎𝑖0 et tous les autres de même signe et si 𝑏i0 est non
nul, l‟EDP est de type parabolique.
Les fonctions 𝑎i et 𝑏i étant dépendantes des variables (𝑥1, … . .n), la classification est
évidemment
Fonction du point (𝑥1, … . .n) considéré. Une EDP peut donc être de différents types
suivant les points considérés : on dit qu‟elle est de type mixte.
Exemples (II.1) : soient (𝑥, 𝑦) une fonction de deux variables et (𝑥, 𝑦, 𝑡) une fonction
de trois variables.

𝜕2𝑇 𝜕2𝑇
2
+ = 0 Est une EDP elliptique.
𝜕𝑥 𝜕𝑦 2

16
𝜕2𝑔 𝜕2𝑔 𝜕2𝑔
= + Est une EDP hyperbolique
𝜕𝑡 2 𝜕𝑥 2 𝜕𝑦 2

𝜕𝑔 𝜕2𝑔 𝜕2𝑔
= 2
+ Est une EDP parabolique
𝜕𝑡 𝜕𝑥 𝜕𝑦 2

𝜕2𝑓 𝜕2𝑓
𝑥 2
+ =0 Elliptique pour x> 0
𝜕𝑥 𝜕𝑦 2

Hyperbolique x< 0
Parabolique x=0
II.4. Classification physique des EDP :
De nombreux phénomènes physiques se rangent dans l‟une des classes suivantes :
*Les problèmes d‟équilibre étudient l‟état stationnaire d‟un phénomène (champ,
chaleur……) dans un domaine borné ou non. Ils sont gouvernés par l‟EDP elliptiques.
*Les problèmes de valeurs propres sont en général des extensions des problèmes
d‟équilibre dans lesquels les valeurs critiques de certains paramètres doivent être
déterminées. C‟est le cas par exemple de la résonance des circuits électriques.
*Les problèmes d‟évolution étudient l‟évolution avec le temps d‟un phénomène
(champ, chaleur, vibration,….) à partir d‟un état initial donné. Ils sont gouvernés par
des EDP hyperboliques ou des EDP paraboliques.

II.4.1.‐ Equation de la chaleur :

La conduction thermique à l‟intérieur d‟un domaine D bidimensionnel provoque un
changement de la température (𝑡, 𝑥, 𝑦), qui régi, en l‟absence de source de chaleur par
l‟EDP :
𝜕𝑇 𝜕 𝜕𝑇 𝜕 𝜕𝑇
𝜌𝑐 = 𝑘 + … … … … … … … … … … … … … … … II. 7
𝜕𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑥

Où 𝑘, 𝜌, 𝑐 sont respectivement la conductivité thermique, la masse volumique et la

chaleur spécifique du solide constituant le domaine D.
Lorsque k dépend seulement de la position (𝑥, 𝑦), l‟EDP est linéaire; si k dépend de la
température T, l‟EDP est non linéaire.
Dans la majorité des cas rencontrés, on considère k comme constante et l‟équation de
la chaleur peut être sous la forme:
𝜕𝑇 𝑘 𝜕2𝑇 𝜕2𝑇
= 2
+ = α ∆𝑇 … … … … … … … … … … … … … … … … … . II. 8
𝜕𝑥 𝜌𝑐 𝜕𝑋 𝜕𝑌 2

𝜕2𝑇 𝜕2𝑇
∆𝑇 = 2
+ désigne le la placien de T
𝜕𝑋 𝜕𝑌 2

17
 𝑘
α= est le coefficient de T
𝜌𝑐

De manière générale, si T dépend de n variables d‟espace (𝑥1 , … . . 𝑥𝑛 ) on a :

𝜕𝑇 2
𝑖=𝑛 𝜕 𝑇
=α 𝑖=1 𝜕𝑥 2
𝜕𝑡 𝑖

Tous les problèmes de diffusion sont régis par ce type d‟équations.

Exemples(II.2):

∂ ∂T ∂ ∂T ∂ ∂T ∂T
k + k + k + e.gen = ρ. C. ……………… (1.31)
∂x ∂x ∂y ∂y ∂z ∂z ∂t

∂ ∂T
-2-D bidimensionnelle (x,y) k =0
∂z ∂z
∂T
- Régime permennt ∂t = 0

-Pas génération de chaleur e.gen = 0

Propriétés constante (k=cte)

L‟équation de la place qui est elliptique (1.31)

𝜕2𝑇 𝜕2𝑇
+ = 0 K= cte
𝜕2𝑥 𝜕2𝑦
𝜕2𝛷 𝜕2𝛷 𝜕2𝛷 𝜕𝛷 𝜕2𝛷
𝐴 2
+𝐵 +𝐶 +D +𝐸 + 𝐸𝛷 + 𝐺 𝑥, 𝑦 = 0
𝜕𝑋 𝜕𝑥𝑦 𝜕2𝑧 𝜕𝑥 𝜕𝑦

Δ= (b2 -4.a.c) =-4 < 0, A=1, B=0, C=1 l‟équation elliptique

L‟équation de la conductivité de chaleur avec la génération de chaleur
𝜕2𝑇 𝜕2𝑇 𝒆.𝒈𝒆𝒏
+ + = 0.
𝜕2𝑥 𝜕2𝑦 𝑘
Laquelle est une équation elliptique
La caractéristique de l‟équation elliptique, laquelle est exige la spécifique des
conditions au limite appropriées sur toutes les limites.
* les équations de la conduction thermique en régime variable (transitoire) sans
génération de chaleur avec K= cet.
𝜕2𝑇 1 𝜕𝑇
= . Équation parabolique unidimensionnelle (1-D)
𝜕2𝑥 𝛼 𝜕𝑡
𝜕2𝛷 𝜕2𝛷 𝜕2𝛷 𝜕𝛷 𝜕2𝛷
𝐴 2
+𝐵 +𝐶 +D +𝐸 + 𝐸𝛷 + 𝐺 𝑥, 𝑦 = 0
𝜕𝑋 𝜕𝑥𝑦 𝜕2𝑧 𝜕𝑥 𝜕𝑦

A=0 B=0 C = 0 ∆=0

L‟équation d‟ordre deux 2eme Ordre

18
𝜕2 𝛷 1 𝜕 2 𝛷𝑇
= . … … … … … … … … … … … … … … … … … … … … … … I. 5
𝜕 2 𝑥 𝐶 2 𝜕𝑡
On t et le temps x est la variable d‟espace, et c est la variable de la vitesse de la
propagation de l‟onde
1 1
A = 1, B = 0, C = ∆= >0
𝐶2 𝐶2
Est l‟équation hyperbolique

 l‟équation de la conduction de chaleur (2D) régime transitoire, avec la

génération de chaleur et propriété constante K=Cst.

𝜕 2 𝑇 𝜕 2 𝑇 𝒆.𝒈𝒆𝒏 1 𝜕𝑇
+ + = .
𝜕 2 𝑥 𝜕𝑦 2 𝑘 𝛼 𝜕𝑡
Est une équation parabolique dans le temps

II.5.Résolution des systèmes d’équations linières

Dans cette partie, nous intéressant à la résolution de l‟équation algébrique linière qui la forme
générale suivant :

𝑎11 𝑥1 𝑎12 𝑥2 ⋯ 𝑎1𝑛 𝑥𝑛 =𝑏1

𝑎21 𝑥1 𝑎22 𝑥2 ⋯ 𝑎21 𝑥𝑛 = 𝑏2
………………………… II. 11
⋮ ⋮ ⋯ ⋮
𝑎11 𝑥1 𝑎12 𝑥2 ⋯ 𝑎𝑛𝑛 𝑥𝑛 = 𝑏3
Où les coefficients 𝑎11 𝑎12 𝑎𝑛𝑛 , sont des coefficients constants et n le nombre
d‟équations.

Pour de petit nombre d‟équation n ≤ 3

L‟équation linière (parfois non linéaire) peut être résolue par des techniques simples
(Méthode cramer)

*cependant pour quatre équations ou plus la résolution deviennent ordeuse et les

ordinateurs doivent être utilisés.

La méthode d‟élimination de Gauss (Méthode directe), algorithme de thomas

(Méthode directe) et la méthode de Gauss-seidel (Méthode itérative) seront présentés
dans ce chapitre ainsi leurs algorithme.

Notation matricielle :

19
 𝑎11 𝑎12 𝑎13 ⋯ 𝑎1𝑚
𝑎 𝑎22 𝑎23 𝑎2𝑚
𝐴 = 21 ⋯
⋮ ⋮ ⋮ ⋮
𝑎𝑛1 𝑎𝑛2 𝑎𝑛3 ⋯ 𝑎𝑛𝑚
La matrice dans la figure (II.1)a n rangées et m colonnes ayant une dimension de n
par m (ou n x m) elle est appelé n x m

Les matrices où n=1 𝐵 = 𝑏1 𝑏2 ⋯ 𝑏𝑚 est appelées vecteurs rangées

 les matrices où m=1

𝐶1
𝐶
𝐶= 2
⋮
𝐶𝑛
 les matrices où n=m sont appelées matrices carrées par exemple une matrice u
par u

𝑎12 𝑎12 𝑎13 𝑎14

𝑎 𝑎22 𝑎23 𝑎24
𝐴 = 𝑎21 𝑎32 𝑎33 𝑎34
31
𝑎41 𝑎42 𝑎43 𝑎44
*une matrice symétrique

5 1 1
𝐴 = 1 3 7
2 7 8
𝑎𝑖𝑗 = 𝑎𝑗𝑖

𝑎𝑖𝑗 Les éléments de la matrice, i désigne la lecture des éléments dans la direction
horizontale, et j la lecture des éléments dans la direction verticale.
 une matrice triangulaire supérieure

𝑎11 𝑎12 𝑎13

𝐴 = 0 𝑎22 𝑎23
0 0 𝑎33
 une matrice triangulaire inférieure
𝑎11 0 0
𝐴 = 𝑎21 𝑎22 0
𝑎31 𝑎31 𝑎33

20
  une matricé d’identité
1 0 0
𝐴 = 0 1 0
0 0 1
 une matricé tri-diagonale à la forme suivante:
𝑎12 𝑎12 0 0
𝑎 𝑎22 𝑎23 0
𝐴 = 21 𝑎33 𝑎34
0 𝑎32
0 0 𝑎43 𝑎44

II.5.1. Méthode de Gauss-seidel :

Quand le nombre d‟équation est très large une méthode itérative préfère à la méthode
directe de solution.

-la Méthode itérative converger bien pour un système d‟équation ayant une diagonale
dominante (pas de zéro) parmi les Méthode itérative ,on a la méthode de (gauss-
seidel),cette Méthode est très utilisée par apport à autre méthodes itératives.

*supposer que nous avons un système d‟équation de n équations

𝐴 𝑋 = 𝐵
*prenons n=3, un système d‟équation de 3X3

*Si les éléments de la diagonale sont tous nuls

𝑎11 𝑎12 𝑎13 𝑥1 𝑏1

A= 𝑎21 𝑎22 𝑎23 = 𝑥2 , 𝐵 = 𝑏2
𝑎31 𝑎32 𝑎33 𝑥3 𝑏3
 le système d‟équation (n=3)

𝑎11 𝑥1+𝑎12 𝑥2 + 𝑎13 𝑥3 = 𝑏1

𝑎21 𝑥1+𝑎22 𝑥2 + 𝑎23 𝑥3 = 𝑏2 … … … … … … … … … … … … . II. 12
𝑎31 𝑥1+𝑎32 𝑥2 + 𝑎33 𝑥3 =𝑏3
La première équation peut être résolue pour 𝑥1

La deuxième équation peut être résolue pour 𝑥21

21
La troisième équation peut être résolue pour 𝑥3

b 1 −a 12 x 2 −a 13 x 3
=> x1 = ……………………………. II. 13. a
a 11

b2 − a21 x1 − a23 x3
=> x2 = … … … … … … … … … … . II. 13. b
a22
b3 − a31 x1 − a32 x2
=> x3 = … … … … … … … … … … . . II. 13. c
a33
Maintenant nous peuvent commencer la procédure de solution on choisit les valeurs
estimées𝑥1 , 𝑥2 , 𝑥3 .

Une méthode simple pour obtenir les estimations initiales et de supposer que des
valeurs 𝑥2 , 𝑥3 sont nulles ces zéros sont remplacés de (II.13-a) pour calculer 𝑥1 =
𝑏1
.alors nous remplaçons cette valeur de 𝑥1 dans l‟équation de (II.13-b) pour calculer
𝑎 11

la nouvelle valeur de𝑥2 .

La procédure est répétée pour l‟équation (II.13.c) pour calculer une nouvelle
estimation pour𝑥3 .

Ainsi nous retournons a la première équation est répété la procédure jusqu‟à votre
solution converge assis près des valeurs critères.

 la convergence peut être vérifiée en utilisant le critère suivant :

𝑗 𝑗 −1
𝑥 𝑖 −𝑥 𝑖
ɛ𝑖 = 𝑗 <ɛ𝑖
𝑥𝑖

pour tout i(𝑥𝑖 ,i=1,,,,,n)

j et j-1 sont respectivement la présente et la précédente itérations.

Exemple (II.3)

Utiliser la méthode de Gauss-seidel pour obtenir la solution du système d‟équation

suivant :

3𝑥1 − 0.1𝑥2 − 0.2𝑥3 = 7.85

0.1𝑥1 + 7𝑥2 − 0.3𝑥3 = -19.3… … … … … … … … … … … … … … … … . II. 14

22
0.3𝑥1 − 0.2𝑥2 + 10𝑥3 = 71.4

Sachant que la solution exacte est 𝑥1 = 3, 𝑥2 = 2.5, et 𝑥3 = 7

Estimer 𝑥2 = 0 et 𝑥3 = 0
7.85−0.1 x 2 −0.2 x 3
 x1 = …………………………………..(II. 14. a)
3
−19.3−0.1 𝑥 1 +0.3 𝑥 3
 𝑥2 = ……………………………...…(II. 14. b)
7

71.4−0.3𝑥 1 −0.2𝑥 2
 𝑥3 = … … … … … … … … … … … … … (II. 14. c)
10
7.85+0+0
L‟équation(II.14.a) 𝑥1 = = 2.6166
3

Cette valeur avec 𝑥3 =0 est remplacée à l‟équation (II.14.b) pour calculer

−19.3−0.1(2.6166)+0
𝑥2 = =-2.7945
7

La première itération est terminée, en remplaçant la valeur de 𝑥1 𝑒𝑡 𝑥2 dans l‟équation

(II.14.c) pour obtenir 𝑥3

71.4−0.3(2.6166 )−0.2(−2.7945)
𝑥3 = =7.0056.
10

Pour la 2emeitération la même

7.85 − 0.1 −2.7945 + 0.2(7.0056)

ɛ𝑖 = 0.31% 𝑥1 = = 2.9905
3
−19.3 − 0.1(2.9905) + 0.3(7.0056)
ɛ𝑖 = 0.015% 𝑥2 = = −2.49962
7
71.4 − 0.3 2.990557 + 0.2 −2.49962
ɛ𝑖 = 0.0012 𝑥3 =
10
= 7.0002
La méthode est, par conséquence, convergée vers la solution exacte des itérations en
plus peuvent être appliquées pour améliorer la solution.

23
II.5.2Amélioration de la convergence en utilisant la relaxation :

La relaxation une légère modification de la méthode de Gauss-seidel et elle est

désignée pour augmenter la convergence, donnée par la relation suivante :

𝑥𝑖𝑁𝑒 =λ𝑥𝑖𝑁𝑒 + (1 −λ)𝑥𝑖𝑜𝑙𝑑

λ: facteur de relaxation, 0< λ<2
Si 0< λ<1 :=> sous relaxation
Si 0< λ<2 :=> sur relaxation
Si 𝝀 = 𝟏 la méthode est dite Gauss- Seidel

II.5.3Méthode de Jacobi

Exemple : soit le système suivant :

*le système d‟équation (n=3)

𝑎11 𝑥1+𝑎12 𝑥2 + 𝑎13 𝑥3 =𝑏1

𝑎21 𝑥1+𝑎22 𝑥2 + 𝑎23 𝑥3 =𝑏2 … … … … … … … … … … … … … …(II. 15)
𝑎31 𝑥1+𝑎32 𝑥2 + 𝑎33 𝑥3 =𝑏3
La première équation peut être résolue pour 𝑥1

La deuxième équation peut être résolue pour 𝑥2

La troisième équation peut être résolue pour 𝑥3

𝑏1 −𝑎 12 𝑥 2 −𝑎 13 𝑥 3
 𝑥1 = … … … … … … … … … … … … …(II. 15. a)
𝑎 11

𝑏2 −𝑎 21 𝑥 1 −𝑎 23 𝑥 3
=> 𝑥2 = … … … … … … … … … … … … … ..(II. 15. b)
𝑎 22

𝑏3 −𝑎 31 𝑥 1 −𝑎 32 𝑥 2
=> 𝑥3 = … … … … … … … … … … … … … ..(II. 15. c)
𝑎 33

Maintenant nous peuvent commencer la procédure de solution on choisit les valeurs

estimer 𝑥1 , 𝑥2 , 𝑥3 .

Une méthode simple pour obtenir les estimations initiales et de supposer que des
valeurs 𝑥2 , 𝑥3 sont nulles ces zéros sont remplacés de (15-a) pour calculer 𝑥1 =

24
𝑏1
.alors nous remplaçons cette valeur de 𝑥1 dans l‟équation de (2.15-b) pour calculer
𝑎 11

la nouvelle valeur de 𝑥2 .

La procédure est répétée pour l‟équation (15.c) pour calculer une nouvelle estimation
pour 𝑥3 .

A la 1ère itération on donne des valeurs initiales à𝑥1 , 𝑥2 , 𝑥3 , soit (𝑥1 = 𝑥2 = 𝑥3 =

0), on obtient 𝑥1 , 𝑥2 , 𝑥3 ensuite on passe à la 2ème itération (itération 2), on calcule

𝑥1 , 𝑥2 , 𝑥3 à partir de (𝑥1 , 𝑥2 , 𝑥3 ) de 1ère itération et ainsi de suite jusqu‟à la

convergence :
II.5.4. Méthode de résolution directe :
Le problème essentiel de la méthode directe est la place mémoire. Cette méthode
présente des avantages et des inconvénients par rapport à la méthode itérative :
*Avantage :
-Le temps de calcul en général plus petit pour une même précision ;
-Dans certains cas la méthode itérative peut ne pas converger.
*Inconvénients :
-Occupation mémoire importante ;
-Risque d‟erreur d‟arrondis importante si certains pivots sont trop petit ;
-La méthode directe n‟est pas applicable aux équations non linéaires.

II.5.4.1 Méthode d’élimination de Gauss :

Soit le système donné sous cette forme
𝑎11 𝑥1 𝑎12 𝑥2 ⋯ 𝑎1𝑛 𝑥𝑛 =𝑏1
𝑎21 𝑥1 𝑎22 𝑥2 ⋯ 𝑎21 𝑥𝑛 = 𝑏2
… … … … … … … … … … … … …(II. 16)
⋮ ⋮ ⋯ ⋮
𝑎11 𝑥1 𝑎12 𝑥2 ⋯ 𝑎𝑛𝑛 𝑥𝑛 = 𝑏3

𝑎11 𝑎12 𝑎13 ⋯ 𝑎1𝑚 𝑥1 𝑏1

𝑎 𝑎22 𝑎23 𝑎2𝑚 𝑥2 𝑏2
𝐴 = 21 ⋯ ⋮ =
⋮ ⋮ ⋮ ⋮ ⋮
𝑎𝑛1 𝑎𝑛2 𝑎𝑛3 ⋯ 𝑎𝑛𝑚 𝑥𝑛 𝑏𝑛

𝐴. 𝑋 = 𝑏

25
La méthode consiste à transformer la matrice A en une matrice triangulaire supérieure
S, le vecteur b subit les mêmes opérations et devient b’
𝑆. 𝑋 =𝑏‟
𝑎11 𝑎12 𝑎13 ⋯ 𝑎1𝑛 𝑥1 𝑏1
′ ′ ′ 𝑥2 𝑏2′
0 𝑎22 𝑎23 ⋯ 𝑎2𝑛
0 0 ′
𝑎33 ⋯ ′
𝑎3𝑛 ⋮ = 𝑏3′
⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
0 0 0 ⋯ ′
𝑎𝑛𝑛 𝑥𝑛 𝑏𝑛′
Exemple (II.4) : Résoudre le système suivant :

3𝑥1 − 𝑥2 + 2𝑥3 = 12 … … (1)

𝑥1 + 2𝑥2 + 3𝑥3 = 11 … … (2)
2𝑥1 − 2𝑥2 − 𝑥3 = 2 … … . . (3)

3 −1 2 𝑥1
1 2 3 𝑥2
2 −2 −1 𝑥3
Etape 1 : Elimination de 𝑥i des équations (2) et (3)
𝑎 21 1
eq2 − . eq1= eq2− . eq1
𝑎 11 3

𝑎 31 1
eq2− . eq1= eq3− . eq1
𝑎 11 3

On obtient le système suivant :

3𝑥1 − 𝑥2 + 2𝑥3 = 12 (1)

0. 𝑥1 + 2.33𝑥2 + 2.33𝑥3 = 7 (2′ )
0𝑥1 − 1.33𝑥2 − 2.33𝑥3 = −6 (3′ )

3 −1 2 𝑥1 12
0 2.33 2.33 𝑥2 = 7
𝑥
0 −1.33 −2.33 3 −6
Etape 2 : Elimination de 𝑥3 des équations (3’)
′
𝑎 32 (−1.33)
eq3’− ′ . eq2’= eq3’− . eq2’
𝑎 22 2.33

3𝑥1 − 𝑥2 + 2𝑥3 = 12 (1)

0. 𝑥1 + 2.33𝑥2 + 2.33𝑥3 = 7 (2′ )
0𝑥1 − 0𝑥2 − 𝑥3 = −2 (3′′ )

26
 3 −1 2 𝑥1 12
0 2.33 2.33 𝑥2 = 7
0 0 −1 𝑥3 −2
Etape 3 : substitution en arrière
De l‟équation (3′′ ) on obtient
x3=2

De (2′ ) on obtient :

2,33𝑥2 + 2,33𝑥3 = 7

x2=1
De (𝟏) on obtient :
𝑥1=(12+ 𝑥2 - 2𝑥3)/3
𝑥1 = 3

II.5.5.2Algorithme de Thomas (Méthode direct)

Considérer un système de N équation algébrique le coefficient de la matrice tri
diagonal est les suivants

𝑏1 𝑐1 0 𝑂 0 0 𝑇1 𝑑1
𝑎1 𝑏2 𝑐2 0 0 0 𝑇2 𝑑2
0 𝑎3 𝑏3 𝐶3 0 0 𝑇3 𝑑
= 3
⋮ ⋮ ⋱ ⋱ ⋱ 𝑂 ⋮ ⋮
0 0 0 ⋱ 𝑏𝑛 −1 𝑐𝑛 −1 𝑇𝑛 −1 𝑑𝑛−1
0 0 0 0 𝑎𝑛 𝑏𝑛 𝑇𝑁 𝑑𝑁
1ere étape :

La première équation (Eq(1)) est choisie comme‟ pivot „

𝑎 𝑎2
Multipliée x 2 . (Eq(1) x ) est soustrait à partir de la deuxième équation (Eq(2))
𝑏1 𝑏1
pour éliminer a2 de l‟équation (Eq(2))

𝑎2
Nouveaua2 = a2 – xb1=0
𝑏1

𝑎2
Nouveaub2 = b2 – xC1=0
𝑏1

𝑎2
Nouveaud2 = d2 – xd1=0
𝑏1

27
2ere étape :

La deuxième équation Eq (2) modifiée

Est choisie comme pivot, une approche similaire est suivie pour éliminer, a3

Nouveaua3= 0

𝑎3
Nouveaub3 = b3 – xC2=0
𝑏2

𝑎3
Nouveaud3= d3 – xd2=0
𝑏2

3ere étape :

La procédure est continue jusqu‟à aN est éliminé (aN=0)

De la dernière équation Eq(N)

En générale :
𝑎𝑖
Remplacer bi par bi – x Ci-1=0 pour i=2,3,…N
𝑏 𝑖−1

𝑎𝑖
Remplacer di par di – x di-1=0 pour i=2,3,…N
𝑏 𝑖−1

les inconnues Ti(𝑇1 , 𝑇2 , 𝑇3 ……𝑇𝑁 ) sont calculés par substitution, en arriérée

commençant de la dernière équation :

𝑑𝑁
𝑇𝑁 =
𝑏𝑁
𝑇 𝑑 𝐶 ∗𝑇 , i=N-1, N-2 ………
𝑖= 𝑖− 𝑖 𝑖−1
𝑏𝑖

Exemple (II.5):considère de système suivant une matrice à coefficient tri diagonale

Résoudre le système avec la méthode de thomas.

−1 1 0 0 𝑇1 −40
1 −2 1 0 𝑇2 −30
=
0 1 −2 1 𝑇3 −30
0 0 1 −2 𝑇4 −30

28
 −𝑇1 + 𝑇2 =
𝑇1 −𝑇2 𝑇3
= −30
𝑇2 2𝑇3
𝑇4
𝑇3 2𝑇4
Appliquons les équations (2.14a) et (2.14b)

𝑏1 𝑐1 0 𝑇1 𝑑1
𝑎2 𝑏2 𝑐2 𝑇2 𝑑2
=
0 𝑎3 𝑏3 𝑇3 𝑑3
0 0 𝑎4 𝑇4 𝑑4
Remplacer bi par :

𝑎𝑖
b i- 𝑐𝑖−1
𝑏 𝑖−1

𝑏1= − 1
1
𝑏2= − 2 − 𝑋1 = −1
−1
1
𝑏3= − 2 − 𝑋1 = −1
−1
1
𝑏4= − 2 − 𝑋1 = −1
−1
𝑎𝑖
Remplacer di par : di- 𝑑𝑖−1 , i=2, 3,4
𝑏 𝑖−1

𝑑1= − 1
1
𝑑2= −30 − 𝑋 − 40 = −70
−1
1
𝑏3= − 30 − 𝑋 − 70 = −100
−1
1
𝑏4= −30 − 𝑋 − 100 = −130
−1
Appliquons l‟Eq 2.20 pour calculer :

𝑇1 , 𝑇2 , 𝑇3 , , 𝑇4
𝑑4 −130
𝑇4 = = =130
𝑏4 −1

𝑑3− 𝐶3 ∗ 𝑇4 −100 − 1 ∗ 130

𝑇3 = = = 230
𝑏3 −1

29
 𝑑2−𝐶2 ∗ 𝑇3 −70 − 1 ∗ 230
𝑇2 = = = 300
𝑏3 −1
𝑑1− 𝐶1 ∗ 𝑇2 −40 − 1 ∗ 300
𝑇1 = = = 340
𝑏3 −1

30
III. Méthode des différences finies
C‟est une méthode d‟approximation des équations. On cherche une solution exacte à
𝜕 𝜕 𝜕2
partir de la discrétisation des opérateurs différentiels ( , , )sur un maillage.
𝜕𝑥 𝜕𝑦 𝜕𝑥 2

III.1 Approximation des opérateurs différentiels.

Le développement autour du point i (dans le maillage) d‟une grandeur quelconque donne

1 𝑑𝛷 1 𝑑 2𝛷 1 𝑑 3𝛷
𝛷𝑖−1 = 𝛷𝑖 −
1ǃ
∆x + ∆x 2
𝑑𝑥 𝑖 2ǃ
+ ∆x 3
𝑑𝑥 2 𝑖 3ǃ 𝑑𝑥 32 𝑖
… … … . …(III. 1)

1 𝑑𝛷 1 𝑑 2𝛷 1 𝑑 3𝛷
.𝛷𝑖∓1 = 𝛷𝑖 +
1ǃ
∆x + ∆x 2
𝑑𝑥 𝑖 2ǃ
+ ∆x 3
𝑑𝑥 2 𝑖 3ǃ 𝑑𝑥 3 𝑖
… … … …(III. 2)

Sachant que ∆x=hi+1-hi= cts… … … … … … … … … … … … … ….(III. 3)

𝑑𝛷 𝛷 𝑖 −𝛷 𝑖−1
= +……… … … … … … … … … … … … … … … … …(III. 5)
𝑑𝑥 𝑖 ∆𝑥

Ce qui signifie que la dérivée d‟ordre 1, au point i, est approchée par différences
finies régressives d‟ordre 1.
En retenant les premiers deux termes du développement de la relation (3.2) on obtient
𝑑𝛷 𝛷 𝑖+1 −𝛷 𝑖
= +……… … … … … … … … … … … … … … … ….(III. 6)
𝑑𝑥 𝑖 ∆𝑥

Ce qui signifie que la dérivée d‟ordre 1, au point i, est approchée par différences finies
progressives d‟ordre 1.
En soustrayant la relation (III. 1) de la relation (III. 2) on obtient l‟approximation par
différences finies centrées d‟ordre 2

𝑑𝛷 𝛷 𝑖+1 −𝛷 𝑖−1
= +…………… … … … … … … … … … … … … …(III. 7)
𝑑𝑥 𝑖 ∆𝑥
En additionnant les relations (III. 1) et (III. 2) on obtient l‟approximation de la dérivée
de deuxième ordre par différences finies centrées d‟ordre 2.

𝑑 2𝛷 𝛷 𝑖+1 −2𝛷 𝑖+ 𝛷 𝑖−1

= … … … … … … … … … … … … … … … … … … ….(III. 8)
𝑑𝑥 2 𝑖 ∆x 2

III.2 Maillage

Positionner des nœuds où seront stockées des différentes grandeurs à étudier dans le
domaine d‟étude.

31
III.2.1 Maillage structuré et non-structuré

Ces maillages sont présentés dans la figure suivante, Les avantages et les
inconvénients des deux maillages sont :

Figure III.1 Maillage structuré et non-structuré

Le maillage structuré est simple à mettre en œuvre, par contre, il est limité pour des
géométries simples. Le maillage non-structuré est très compliqué à mettre en œuvre
(utilisation des logiciels commerciaux comme Gambit, Catia..). Il est utilisé pour les
géométries complexes comme le corps humain (étude de la portance et de la trainée
générées par les nageurs), écoulement autour des ailes d‟avion, aubes…
Nous utilisons en ce qui suit des maillages structurés.
III.3 Génération du maillage structuré

III.3.1 Maillage régulier

C‟est un maillage où la distance entre les différents nœuds est constante. Il est facile à
mettre en œuvre.
III.3.1.1 Procédure

On fixe le nombre de nœuds n (cas unidimensionnel) pour le cas bidimensionnel, on

fixe n et m.
Calculer la distance entre les nouds par la relation suivante
𝐿
∆𝑥 = … … … … … … … … … … … … … … … … … … … … … … … … … . (III. 9)
𝑛−1
On utilise ce type de maillage lorsque la variable à déterminer dans le domaine d‟étude
varie faiblement.

32
III.3.2 Maillage irrégulier
Recommander si la variable varie fortement dans le domaine d‟étude, dans ce cas, il est
préférable de raffiner le maillage dans certaines zones ou les gradients varient fortement
afin de capter au mieux l‟information.

III.3.2.1 Procédure
Fixer le nombre de nœuds selon la direction où la variable à déterminer varie fortement
Calculer le premier pas selon la somme d‟une suite géométrique définit comme suit :
1−𝑟
∆𝑥 1 = 𝐿 … … … … … … … … … … … … … … … … … … … . . (III. 10)
1−𝑟 𝑛 −1
Les autres pas sont déterminés par l‟expression suivante :
∆𝑥 𝑖 + 1 = ∆𝑥 𝑖 . 𝑟 𝑛−1 … … … … … … … … … … … … … … … … … (III. 11)
Avec
r: Raison de la suite géométrique (généralement 0.8<r<1.2) pour obtenir une solution
stable.
L : longueur du domaine à discrétisé.
III.4 Résume
Pour résoudre les équations aux dérivées partielles par la méthode des différences
finies, on procède comme suit :
1- Mailler le domaine d‟étude
2- Discrétiser les équations ainsi que les conditions aux limites
3- Constitue le système d‟équation globale
4- Résoudre ce système
III.5 Etude numérique de la conduction unidimensionnel stationnaire

III.5.1 Exemple 1

Soit une barre de longueur L. les températures aux bouts de celle-ci sont T(0) = 10°C
𝐿 𝐿 3𝐿
et T(L) = 50°C. Déterminer les températures aux points x=4 , 2 𝑒𝑡 4

Solution
1/ maillage

33
 Figure III.2 Géométrie du problème.
Le nombre de nœuds est n= 5
2/ Modèle mathématique.
L’équation de la pure conduction en 1-D stationnaire est :
d2 T
= 0
dx 2
T 0 = 10°C
T L = 50°C
Le modèle mathématique comporte les équations qui régissent le phénomène
physique ainsi que les conditions aux limites imposées.
3/ discrétisation:
L‟opérateur différentiel du deuxième ordre est évalué par le développement en série
de Taylor, on obtient :
𝑇𝑖−1 − 2𝑇𝑖 + 𝑇𝑖+1
∆x 2
Noter que cette équation est valable uniquement aux nœuds internes.
Construction du système d‟équation

𝑖= 2 𝑇1 − 2𝑇2 + 𝑇3 = 0
𝑖= 3 𝑇2 − 2𝑇3 + 𝑇4 = 0
𝑖= 4 𝑇3 − 2𝑇4 + 𝑇5 = 0

La discrétisation des conditions aux limites donne :

𝑇1 = 10°𝐶
𝑇5 = 50°𝐶

L‟écriture du système d‟équation pour tous les nœuds donne :

34
 𝑇1 = 10
𝑇1 − 2𝑇2 + 𝑇3 = 0
𝑇2 − 2𝑇3 + 𝑇4 = 0
𝑇3 − 2𝑇4 + 𝑇5 = 0
𝑇5 = 50

Mise sous forme matricielle elle s‟écrit :

1 0 0 0 0 𝑇1 10
1 −2 1 0 0 𝑇2 0
0 1 −2 1 0 𝑇3 = 0
0 0 1 −2 1 𝑇4 0
0 0 0 0 1 𝑇5 50
La matrice associée aux systèmes est tri-diagonales, on utilise l‟algorithme de Thomas
[1] pour obtenir la solution.

On obtient finalement :
𝑇1 10
𝑇2 0
𝑇3 = 0
𝑇4 0
𝑇5 50
III.5.2 Exemple 2

Figure III.3 : Géométrie du problème.

Solution :
Le modèle mathématique du présent exemple est :

35
 𝑑2 𝑇
𝒌 + 𝒒= 𝟎
𝑑𝑥 2
𝑇 0 = 𝟏𝟎𝟎°𝑪
𝑇 𝐿 = 𝟐𝟎𝟎°𝑪

L‟équation de la pure conduction en 1-D stationnaire avec source interne. La

discrétisation du modèle mathématique donne :

𝑇𝑖−1 −2𝑇𝑖+ 𝑇𝑖+1

k +𝑞𝑖 = 0
∆x 2

Comme le chauffage est homogène qi est constant dans tout le domaine d‟étude, Le
système d‟équations obtenu pour tous les nœuds est :

𝑇1 = 100
𝑞2 ∆𝑥 2
𝑇1 − 2𝑇2 + 𝑇3 = −
𝑘
𝑞 ∆𝑥 2
𝑇2 − 2𝑇3 + 𝑇4 = − 2
𝑘
𝑞2 ∆𝑥 2
𝑇3 − 2𝑇4 + 𝑇5 = −
𝑘
𝑇5 = 200
Noter que
K
q2 = q3 = q4 = 1000
m3

III.5.3 Exemple 3
Résoudre l‟exemple 2 pour le cas où la face « Ouest » est maintenue à 100° C et la
face « Est » est isolée.
Solution.
𝑑2 𝑇
+ 𝒒= 𝟎
𝑑𝑥 2
𝑇 0 = 𝟏𝟎𝟎°𝑪
à 𝒙 = 𝑳 𝒑𝒂𝒓𝒐𝒊 𝒊𝒔𝒐𝒍é𝒆 𝑸𝑳 = 𝟎
A la face Est, la paroi est isolée, on introduit la définition de la densité du flux par la
loi de Fourrier
𝑑𝑇
𝑄𝐿 = =0
𝑑𝑥

36
 𝑑𝑇
On évalue le gradient par un développement en série de Taylor d‟ordre 1 à gauche
𝑑𝑥

𝑑𝑇 𝑇𝑖 − 𝑇𝑖−1
=
𝑑𝑥 ∆𝑥
A x = L c'est-à-dire au dernier nœud, dans notre cas n = 5

𝑇5 − 𝑇4
∆𝑥
𝑇5 = 𝑇4
𝑇1 = 100
𝑞2 ∆𝑥 2
𝑇1 − 2𝑇2 + 𝑇3 = −
𝑘
𝑞 ∆𝑥 2
𝑇2 − 2𝑇3 + 𝑇4 = − 2
𝑘
𝑞 ∆𝑥 2
𝑇3 − 2𝑇4 + 𝑇5 = − 2
𝑘
𝑇5 = 𝑇4
Sous forme matricielle, on aura :

100
𝑞 2 ∆𝑥 2
1 0 0 0 0 𝑇1 𝑘
1 −2 1 0 0 𝑇2 𝑞 2 ∆𝑥 2
0 1 −2 1 0 𝑇3 = 𝑘
0 0 1 −2 1 𝑇4 𝑞 2 ∆𝑥 2
0 0 0 −1 1 𝑇5 𝑘
𝑞 2 ∆𝑥 2
𝑘
III.5.4 Exemple 4
Une plaque d'uranium est soumise à une isolation d'un côté et à une convection de
l'autre côté. La formulation de différences finies de ce problème doit être obtenue, et
les températures nodales dans des conditions stables doivent être déterminées.
Hypothèses
1- Le transfert de chaleur à travers le mur est régulier car il n'y a aucune indication de
changement avec le temps.
2- Le transfert de chaleur est unidimensionnel car la plaque est grande par rapport à
son épaisseur.
3- La conductivité thermique est constante.

37
4- Le transfert de chaleur par rayonnement est négligeable.
Propriétés La conductivité thermique est donnée comme étant k = 28 W / m. ° C.

Le nombre de nœuds est spécifié pour être M = 6. Alors l'espacement nodal ∆X

devient
L 0.05 m
x    0.01 m
M 1 6 -1
Ce problème implique 6 températures nodales inconnues, et donc nous devons avoir 6
équations pour les déterminer de manière unique. Le nœud 0 est sur une frontière
isolée, et nous pouvons donc le traiter comme une note intérieure en utilisant le
concept d'image miroir. Les nœuds 1, 2, 3 et 4 sont des nœuds intérieurs, et donc pour
eux, nous pouvons utiliser la relation de différence finie générale exprimée comme
Tm1  2Tm  Tm1 g m
  0, pour m= 0, 1, 2, 3, et 4
x 2
k
Enfin, l'équation aux différences finies pour le nœud 5 sur la surface droite soumise à
convection est obtenue en appliquant un bilan énergétique sur l'élément demi-volume
autour du nœud 5 et en prenant la direction de tous les transferts de chaleur vers le
nœud considéré:

T1  2T0  T1 g
nœud 0 (surface gauche - isolé) :  0
x 2 k
T0  2T1  T2 g g
Nœud 1 (interieur ) :  0
x 2
k isolé
x
T1  2T2  T3 g h, T
nœud 2 (interieur ) :   0      
x 2 k 0 1 2
3 4 5
T2  2T3  T4 g
nœud 3 (interieur ) :  0
x 2 k
T3  2T4  T5 g
nœud 4 (interieur ) :  0
x 2 k
T  T5
nœud 5 ( surface droit - convection) : h(T  T5 )  k 4  g (x / 2)  0
x
Où x  0.01 m, g  6 10 W/m , k  28 W/m C, h  60 W/m  C, and T  30C.
5 3 2

Ce système de 6 équations avec six températures inconnues constitue la formulation

de différence finie du problème.
(b) Les 6 températures nodales dans des conditions stables sont déterminées en
résolvant les 6 équations ci-dessus simultanément avec un résolveur d'équations à
T0 =556.8C, T1 =555.7C, T2 =552.5C, T3 =547.1C, T4 =539.6C,
etT5=530.0C.

38
 CHAPITRE IV

DISCRETISATION DE L’EQUATION DE LA CONDUCTION DE

CHALEUR (2-D) ENREGIME PERMANENT PAR LA METHODE
DES DIFFERENCES FINIES

IV. Discrétisation de l’équation de la conduction de chaleur (2-D) en

Régime permanent par la méthode des différences finies
Considère une région rectangulaire dans laquelle la conduction de chaleur est
important dans la direction x,y
maintenant, divisons le plan (x-y) de la région en maille rectangulaire ayant des
nœuds espaces de ∆x et ∆y dans les directions x,y ,et considérons ∆z=1 dans la
direction Z.
Notre objectif est de chaleur les temperature aux nœuds m = 1, 2, 3……..M-1
suivant x et N=1, 2, 3……..N -1 suivant y les coordonnée du nœud (m, n) sont
simplement x = m.∆x et y = n.∆y

d
yn
P
W w e E Δ
Y

d
ys
s

Δ
X
d d
x
FigIV.1 : élément de volume xe (m, n) pour la conduction (2-D)
d’un nœud général
w
en coordonnées cartésiennes.

40
Maintenant considérons un élément de volume de dimensions ∆x.∆y.1 centré ou tour
du nœud (m, n) dans une région dont laquelle la chaleur est générée a un taux de e
(m/m3) et la conductivité thermique K est constante (voir FigIV.2)

L‟équation d‟énergie sur l‟élément de volume de la fig IV.2 peut être exprimée
comme, en régime permanent.

𝑡𝑜𝑢𝑥 𝑑𝑒 𝑔𝑒𝑛𝑒𝑟𝑎𝑡𝑖𝑜𝑛
𝑡𝑜𝑢𝑥 𝑑𝑒 𝑙𝑎 𝑐𝑜𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛
𝑑𝑒 𝑐𝑕𝑎𝑙𝑒𝑢𝑟 𝑎 𝑡𝑜𝑢𝑥 𝑑𝑒
𝑑𝑒 𝑐𝑕𝑎𝑙𝑒𝑢𝑟 𝑎𝑢𝑥 ′
𝑙 𝑖𝑛𝑡𝑒𝑟𝑖𝑒𝑢𝑟 𝑐𝑕𝑎𝑛𝑔𝑒𝑚𝑒𝑛𝑡
𝑠𝑢𝑟𝑓𝑎𝑐𝑒 + =
𝑑𝑒 𝑙′𝑒𝑙𝑒𝑚𝑒𝑛𝑡 𝑑𝑒 𝑙 ′ 𝑒𝑛𝑒𝑟𝑔𝑖𝑒 𝑑𝑒
(𝑙𝑒𝑓𝑡, 𝑟𝑒𝑔𝑕𝑡, 𝑏𝑜𝑡𝑡𝑜𝑚, 𝑡𝑜𝑝
𝑙′𝑒𝑙𝑚𝑒𝑛𝑡

∆E 𝑒𝑙𝑒𝑚𝑒𝑛𝑡
Qcond, left+ Qcond,right Qcond,top + Qcond,bottom+ Egen,elem= … … (IV. 1)
∆𝑡
Remarque: les températures entre les nœuds adjacents doivent variés linéairement, la
surface d‟échange dans la direction x, Ax= ∆y.1=∆y
Et surface d‟échange suivant y, Ay= ∆x.1=∆x
𝑇𝑚 −1,𝑛 − 𝑇𝑚 ,𝑛 (𝑇𝑚 +1,𝑛 − 𝑇𝑚 ,𝑛 ) (𝑇𝑚 ,𝑛 −1 − 𝑇𝑚 ,𝑛 )
k∆y + k∆y + k∆x +
∆x ∆x ∆x
(𝑇𝑚 +1,𝑛 − 𝑇𝑚 ,𝑛 )
k∆y + em,n(∆x. ∆y)=0… … … … … … … … … … … … (IV. 2)
∆x
Divisant l‟eq(4-2) par ∆x. ∆y, on obtient :

(𝑇𝑚 −1,𝑛 −2𝑇𝑚 ,𝑛 + 𝑇𝑚 +1,𝑛 ) (𝑇𝑚 −1,𝑛 − 2𝑇𝑚 ,𝑛 + 𝑇𝑚 +1,𝑛 ) 𝑒𝑚 ,𝑛

+ +
∆x 2 ∆x 2 𝑘
= 0 … … … … … … … … … … … … … … … … … … … … … . (IV. 3)
Pour M=1,2 ,3…… M-1
N=1, 2,3……… M-1
Dans le cas ∆x= ∆y=L, l‟eq(4-3) devient :
𝑒 𝑚 ,𝑛 𝐿 2
𝑇𝑚 −1,𝑛 − 4 𝑇𝑚 ,𝑛 + 𝑇𝑚 +1,𝑛 +𝑇𝑚 −1,𝑛 + 𝑇𝑚 +1,𝑛 + =0… … … … . . (IV. 4)
𝑘
Elle peut être exprimée dans cette forme
𝑒 𝑛𝑜𝑒𝑢𝑑 𝐿 2
𝑇𝑙𝑒𝑓𝑡 + 𝑇𝑟𝑖𝑔 𝑕𝑡 + 𝑇𝑏𝑜𝑡𝑡𝑜𝑚 +𝑇𝑡𝑜𝑝 − 4 𝑇𝑛𝑜𝑒𝑢𝑑 + =0… … … … . . (IV. 6)
𝑘
Exemple IV.1 (conduction de chaleur 2-D, régime permanent dans une barre en L)

Considère le transfert de chaleur en régime permanent dans un corps solide L dont la

section et données dans la figure Exemple (IV.1) le transfert de chaleur dans la
direction Z est négligeable et donc le transfert de chaleur est bidimensionnelle a

41
 conductivité thermique du corps 15/m.°C. Et la chaleur générée dans corps e =
2X106/m3

-la surface gauche du corps et isolé et la surface en bas est maintenue à une
température uniforme de 90°C la totalité de la surface en haut est soumis à la
température d‟air T∞=215°C avec un coefficient de transfert de chaleur convectif et la
surface droite est soumise au flux de chaleur q.R=500/m2

Le nombre de nœuds total égale 15 avec ∆x=∆y=1.2cm

Six nœuds sont en bas de la surface et ainsi leur température est température connue

*obtenir les équations de la différence finis et calculer les températures nodales dans
chaque nœud en sachant les équations obtenir.
Élément de volume
1 2 3 Convection
Δy

Isolée
4 5 6 7 8 9
Δy

10 11 12 13 14 15
Δx Δx Δx Δx Δx
90°C
Solution exemple IV.1 : ∆𝑥
h. 𝑇 − 𝑇1 Convection
Supposition 2 ∞
1 2
1-Régime permanent et 2-D

2- K-Cte ∆𝑦 𝑇2 − 𝑇1
Δy/2

k.
2 ∆𝑥
3- rayonnement est négligeable
Isolé
Nœud 1

Δx/2
∆𝑥 𝑇4 − 𝑇1
k.
2 ∆𝑦
4

42
 ∆𝑦 𝑇2 − 𝑇1 ∆𝑥 𝑇4 − 𝑇1 ∆𝑥 ∆𝑥 ∆𝑦
k. + 0 + k. +h. 𝑇∞ − 𝑇1 + 𝑒1. . =0
2 ∆𝑥 2 ∆𝑦 2 2 2

𝑕𝑙 𝑕𝑙 𝑒1. 𝑙 2
− 2+ 𝑇 + 𝑇2 + 𝑇4 = − 𝑇∞ − =0
𝑘 1 𝑘 2𝑘

Nœud 2
h . ∆𝑥 𝑇∞ − 𝑇2
Convection

1 2 3

∆𝑦 𝑇1 − 𝑇2 ∆𝑦 𝑇3 − 𝑇2
Δy/2

k. k.
2 ∆𝑦 2 ∆𝑥

Δx

5 𝑇5 − 𝑇2
k. ∆𝑥
∆𝑦
∆𝑦 𝑇1 − 𝑇2 ∆𝑦 𝑇3 − 𝑇2 𝑇3 − 𝑇2 ∆𝑦
k. 2 + k. + k. ∆𝑥 + h . ∆𝑥 𝑇∞ − 𝑇2 + 𝑒2. ∆𝑥. =0
∆𝑥 2 ∆𝑥 ∆𝑦 2

2𝑕𝑙 2𝑕𝑙 𝑒2. 𝑙 2

𝑇1 − 4 + 𝑇2 + 𝑇3 − 2𝑇5 = − 𝑇 − =0
𝑘 𝑘 ∞ 𝑘

Nœud 3
∆𝑥
𝑕 (𝑇 − 𝑇3 )
2 ∞
2 3
∆𝑦 𝑇∞ − 𝑇3
𝑘 ( ) ∆𝑦
2 ∆𝑥 𝑕 (𝑇 − 𝑇3 )
2 ∞

∆𝑥 𝑇6 − 𝑇3
6 𝑘 ( )
2 ∆𝑦

∆𝑦 𝑇2 − 𝑇3 ∆𝑥 𝑇6 − 𝑇3 ∆𝑥 ∆𝑦
k. + k. + h . (∆𝑥 + ∆𝑦) 𝑇∞ − 𝑇3 + 𝑒3. . =0
2 ∆𝑥 2 ∆𝑦 2 2

2𝑕𝑙 2𝑕𝑙 𝑒3. 𝑙 2

𝑇1 − 2 + 𝑇3 + 𝑇6 + 𝑇4 = − 𝑇 − =0
𝑘 𝑘 ∞ 2𝑘

43
Remarque : le point 4 se situe sur la limite isolée et peut être traité comme nœud
interne, en remplaçant l‟isolation par un miroir ceci met une image réfléchie au nœud
5 du côté gauche du nœud 4

Miroir
Nœud 4 :

5 5

𝑒4. 𝑙 2
𝑇5 + 𝑇5 + 𝑇1 + 𝑇10 − 4 𝑇4 − = 0………………………………(4)
𝑘

Nœud 5 :
𝑒5. 𝑙 2
𝑇4 + 𝑇6 + 𝑇11 + 𝑇2 − 4 𝑇5 − = 0………………………………..(5)
𝑘
2

4 5 6
3
Nœud 6 :
∆𝑥 𝑇7 − 𝑇6
𝑘 ( )
2 ∆𝑦 11
∆𝑦
𝑕 (𝑇 − 𝑇6 )
𝑇5 − 𝑇6 2 ∞
𝑘 ∆𝑦( ) ∆𝑦
∆𝑥 𝑕 (𝑇 − 𝑇6 ) 7
5 2 ∞

6
∆𝑦 𝑇7 − 𝑇6
𝑘 ( )
2 ∆𝑥

𝑇12 − 𝑇6
𝑘 ∆𝑥( )
12 ∆𝑦

44
 𝑇5 − 𝑇6 𝑇12 − 𝑇6 ∆𝑥 𝑇3 − 𝑇6 ∆𝑥
𝑘. ∆𝑦 + k. ∆𝑥 + k. +h. 𝑇 − 𝑇6
∆𝑥 ∆𝑦 2 ∆𝑦 2 ∞
∆𝑦 ∆𝑦 𝑇7 − 𝑇6 3
+h. 𝑇∞ − 𝑇6 + k. + 𝑒6. . ∆𝑦∆𝑥 = 0
2 2 ∆𝑥 4
2𝑕𝑙 2𝑕𝑙 𝑒6. 𝑙 2
𝑇3 − 2𝑇5 − 6 + 𝑇6 + 𝑇4 = −180 − 𝑇∞ − 3 = 0……………(6)
𝑘 𝑘 2𝑘

Nœud 7
h . ∆𝑥 𝑇∞ − 𝑇7

6 Δy/2 7 8

∆𝑦 𝑇8 − 𝑇7
∆𝑦 𝑇6 − 𝑇7 k.
𝑘. 2 ∆𝑥
2 ∆𝑥
Δx
𝑇13 − 𝑇7
13 k. ∆𝑥
∆𝑦
∆𝑦 𝑇6 − 𝑇7 ∆𝑦 𝑇8 − 𝑇7 𝑇13 − 𝑇7
𝑘. + k. + k. ∆𝑥 + h . ∆𝑥 𝑇∞ − 𝑇7
2 ∆𝑥 2 ∆𝑥 ∆𝑦
∆𝑦
+ 𝑒7. ∆𝑥 = 0
2
2𝑕𝑙 2𝑕𝑙 𝑒7. 𝑙 2
𝑇6 − 4 + 𝑇7 + 𝑇8 = −180 − 𝑇∞ − = 0………………(7)
𝑘 𝑘 𝑘
Le nœud 8 identique au nœud 7 la formulation de différence finie peut être obtenue à
partir du nœud 7
Nœud 8
2hl 2hl e.8 l2
T7 − 4 + T8 + T9 = −180 − T − =0
k k ∞ k

Nœud 9 ∆𝑥
𝑕 (𝑇 − 𝑇9 )
2 ∞
∆x=∆y=l
2 9
𝑇15 = 90°𝐶
∆𝑦 𝑇∞ − 𝑇9
𝑘 ( ) ∆𝑦
2 ∆𝑥 𝑞𝑅. 𝑥
2

∆𝑥 𝑇15 − 𝑇9
15 𝑘 ( )
2 ∆𝑦
45
 2hl q.R l 2hl e.9 l2
T8 − 2 + T9 = −90 − − T − =0
k k k ∞ k

Le système de 9 équations pour détermination les températures nodales inconnus

−2.064 𝑇1 + 𝑇2 +𝑇4 = -11.2 ……………………………………….......(1)

𝑇1 + 4.128 𝑇2 + 𝑇3 + 2𝑇5 = -22.4……………………………………..(2)

𝑇2 − 2.128 𝑇3 + 𝑇6 = -12.8 …………………………………………...(3)

𝑇1 − 4𝑇4 + 2𝑇5 = -109.2…………………………………………….....(4)

𝑇2 + 𝑇4 − 4𝑇5 + 𝑇6 = -109.2………………………………………......(5)

𝑇3 + 2𝑇5 − 6.128𝑇6 + 𝑇7 = -212.0……………………………………(6)

𝑇6 − 4.128𝑇7 + 𝑇8 = -202.4…………………………………………...(9)

𝑇7 − 4.128𝑇8 + 𝑇9 = -202.4…………………………………………..(10)

𝑇8 − 2.064𝑇9 = -105.2…………………………………………………………(11)
En utilisant la méthode de Gauss-Seidel au thomas, la solution de système donne

T1 = 112.1°CT2 = 110.8°CT3 = 106.6°CT4 = 109.4°C

T5 = 108.1°CT6 = 103.2°CT7 = 97.3°CT8 = 96.3°C
T8 = 97.6°C .

Exemple IV.2
Un corps long et solide est soumis à un transfert thermique bidimensionnel régulier.
Les températures nodales inconnues et le taux de perte de chaleur de la surface
inférieure à travers une section de 1 m de long doivent être déterminés.

Hypothèses
260 305 350
1. Le transfert de chaleur à travers le corps doit être régulier 200°C 
g
et bidimensionnel. Isolé
3 290 2

2. La chaleur est générée uniformément dans le corps. 5 cm
240 1 325
3. Le transfert de chaleur par rayonnement est négligeable 
0
Convection
h, T
46
Propriétés : La conductivité thermique est donnée comme étant k = 45 W / m. ° C.
Analyse L'espacement nodal est donné ∆X = ∆Y = L = 0,05 m, et la forme générale
de différence finie d'un nœud intérieur pour une conduction thermique
bidimensionnelle stable est exprimée comme
g node l 2
Tleft  Ttop  Tright  Tbottom  4Tnode  0
k
Ou
g nodel 2 g 0 l 2 (8 10 6 W/m3 )(0.05 m) 2
   93.5C
k k 214 W/m C
Les équations aux différences finies pour les nœuds limites sont obtenues en
appliquant un bilan énergétique sur les éléments volumiques et en prenant la direction
de tous les transferts de chaleur vers le nœud considéré:

l 240  T1 290  T1 l 325  T1 e. 0 l 2

Nœud 1 ( convection) : k  kl k  hl (T  T1 )  0
2 l l 2 l 2k
e. 0 l 2
Nœud 2 (interieur ) : 350  290  325  290 - 4 T2  0
k
e. 0 l 2
Nœud 3 (interieur ) : 260  290  240  200 - 4T3  0
k

Ou k  45 W/m. C, h  50 W/m 2 .C, e.  8  10 6 W/m 3 , T  20C

Substitution T1 = 280.9°C, T2 = 397.1°C, T3 = 330.8°C,

(b) Le taux de perte de chaleur de la surface inférieure à travers une section de 1 m de
longueur est :

Q   Q
m
element, m   hA
m
surface, m (Tm  T )

 h(l / 2)( 200  T )  hl (240  T )  hl (T1  T )  h(l / 2)(325  T )

 (50 W/m2  C)(0.05 m  1 m)[(200 - 20)/2  (240 - 20)  (280.9 - 20)  (325 - 20)/2]C  1808 W

47
 Chapitre V

DISCRETISATION DE L’EQUATION DE LA CONDUCTION DE

CHALEUR EN REGIME TRANSITOIRE PAR LA METHODE DE
LA DIFFERENCE FINIE

V.1 Discrétisation de l’équation de la conduction de chaleur en

régime transitoire par la méthode de la différence finie
Dans ce chapitre, nous allons utiliser la méthode de la différence finie pour résoudre
le problème transitoire, dans ce problème les températures changent avec le temps
aussi bien que la position, et donc la solution des différance finies du problème
T
transitoire exige un discrétisation en temps et en espace (Fig V.1)

i+
Tm-1 Tm i+1 Tm+1 i+1
1 ΔT i+1

i
Tm-1 i Tm i Tm+1 i

1
Δx Δx

0 1 2 m-1 m M+
1

Fig V.1 formulation de la différence du problème transitoire en temps et en espace

Dans les problèmes transitoires l‟exposant (i) est utilisé comme compteur
d‟incrément de temps correspond à la condition initiale spécifique en général „i‟ à ti= i
∆t
(∆t incrément du temps, le pas de temps)

𝑇𝑚𝑖 : Représentela température aux nœuds m à l‟incrément de temps

la chaleur transférée de changement

dans l′ élméntde volume la chaleur générée de
à dans l′ élmént de l′ énergie l′ élément
partir des toutes ses + volume = de
surface à durant volume
durant ∆t durant
∆t ∆t

48
 .
On 𝑡𝑜𝑢𝑡𝑒𝑠 𝑄 . +∆t . Egen . élément = ∆Eélément … … … … … … … . . V. 1
𝑙𝑒𝑠𝑠𝑢𝑟𝑓𝑒𝑐𝑒

Avec
∆Eélément = Et+∆t − Et = mcp Tt+∆t − Tt = ρ. dx. A. C. Tt+∆t Tt . V. 2
ρ est la masse volumique, Cp est la chaleur spécifique de l‟élément en divisant
l‟équation V. 1 par ∆t

∆E élément ρ.V élément .C P ∆t

𝑡𝑜𝑢𝑡𝑒𝑠 𝑄. + ∆t . Egen
.
. élément = = … V. 3
𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑒𝑐𝑒 ∆t ∆t

Pour n‟importe quels nœuds dans l‟élément de volume l‟équation est :

. 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑡𝑜𝑢𝑡𝑒𝑠 𝑄 . + ∆t . Egen . élément = ρ. Vélément . CP … … V. 4
𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑒𝑐𝑒 ∆t

Ou𝑇𝑚𝑖+1 𝑒𝑡 𝑇𝑚𝑖 représentela température du nœud m aux temps ti = i∆tet ti+1 =(i+1)∆tet
𝑇𝑚𝑖+1 −𝑇𝑚𝑖 représente le changement de température du nœud m durant l‟intervalle de
temps ∆t entre les incrément i et i+1 (Fig V.2)

∆E= ρ. v. CP . ∆t = ρ. v. CP . 𝑇𝑚𝑖+1 −𝑇𝑚𝑖 … … … … … … … . … … … … … V. 5

Elément de volume

Nœud m’

Fig. V.2 de changement en régime dans l‟élément de volume durant un intervalle de

temps ∆t

𝑇𝑚𝑖+1 = 𝑇𝑚𝑖 Pas de changement de temperature en fonction de temps

 Méthode explicite :

. 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑡𝑜𝑢𝑡𝑒𝑠 𝑄 . + Egen . élément = ρ. Vélément . CP … … … … … V. 6
𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑎𝑐𝑒 ∆t

 Méthode implicite :
i+1 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑡𝑜𝑢𝑡𝑒𝑠 𝑄 .𝑖+1 + Egen . elment = ρ. Vélément . CP … … … … … V. 7
𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑒𝑐𝑒 ∆t

49
 V.2.Conduction de chaleur transitoire dans une paroi plane
Considère la conduction de chaleur transitoire dans une paroi plane en régime
transitoire d‟une paroi plane d‟épaisseur L avec génération de chaleur e.(x, t) ,qui
doit varier avec le temps et la position, et la conductivité thermique k est constante et
𝐿
un maillage de deux dimensions ∆x = 𝑀 et des nœuds à,1,2……,M dans la direction x

voir figure (Fig V.3)

Noter que Vélément = A. ∆x, pour un nœud interne la formulation des différance finies
en régime transitoire peut être exprimée en utilisant l‟équation (V.3) comme.
Elément de volume …….

A Paroi plane ėm A
Tm i+1
(𝑇𝑚 +1 − 𝑇𝑚 )
(𝑇𝑚 −1 − 𝑇𝑚 ) 𝐾. 𝐴.
𝐾. 𝐴. ∆𝑥
∆𝑥 Tm i

Δx Δx

0 1 2 m-1 r m m+1 M-1 M

Fig V.3 les nœuds et l‟élément de volume pour la formulation de la différence finie
transitoire de la conduction 1-D dans une paroi plane.

𝑇𝑚𝑖+1 −𝑇𝑚𝑖 𝑇𝑚𝑖+1 −𝑇𝑚𝑖 𝑇𝑚𝑖+1 −𝑇𝑚𝑖

𝐾. 𝐴. + 𝐾. 𝐴. + e.gen élément . A. ∆x = ρ. CP A. ∆x … V. 8
∆t ∆t ∆t
.∆x
L‟Eq 5.6 par l‟équation V.5 devient
k.A
.
𝑒𝑚 . ∆𝑥 2 ∆𝑥 2
𝑇𝑚 −1 − 2𝑇𝑚 + 𝑇𝑚 +1 + = 𝑇𝑚𝑖+1 −𝑇𝑚𝑖 … … … … … … … . V. 9
𝑘 𝛼∆t

𝑘
𝛼= Diffusivité thermique de la paroi 𝑚2 /𝑠
ρ.C P

𝛼 ∆t
𝜏= … … … … … … … … … … … … … … … … … … … … … … … … … V. 10
∆𝑥 2
.
𝑒𝑚 . ∆𝑥 2 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑇𝑚 −1 − 2𝑇𝑚 + 𝑇𝑚 +1 + = … … … … … … … … … … … V. 11
𝑘 𝜏

50
V.3.1.Formulation explicite des différences finies :
.
𝑒𝑖 𝑚 . ∆𝑥 2 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑇𝑚𝑖 −1 − 2𝑇𝑚𝑖 + 𝑖
𝑇𝑚 +1 + = … … … … … … … … … … V. 12
𝑘 2
𝑖.
𝑒𝑚 . ∆𝑥 2
𝑇𝑚𝑖+1 =𝜏 𝑇𝑚𝑖 −1 +𝑇𝑚𝑖 +1 + 1 − 2𝜏 𝑇𝑚𝑖 +𝜏 … … … … … . V. 13
𝑘
m =1, 2,3…..M-1

V.3.2Formulation Implicite des différences finies :

∆𝑥 2 𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑇𝑚𝑖+1
−1 − 𝑇𝑚𝑖+1 + 𝑇𝑚𝑖+1
+1 + 𝑖+1.
𝑒𝑚 . . = … … … … … … … . V. 14
𝑘 𝜏
∆𝑥 2
𝜏𝑇𝑚𝑖+1
−1 − 1 − 2𝜏 𝑇𝑚𝑖+1 + 𝜏𝑇𝑚𝑖+1
+1 +𝜏 𝑖+1.
𝑒𝑚 . . + 𝑇𝑚𝑖 = 0 … … . V. 15
𝑘
Par exemple, la formulation de la conduction aux limites (cas de la conduction) au
nœud 0 (voir Fig.V.4)

Pour le cas de méthode explicite peut être exprimé comme

𝑕. 𝐴 𝑇∞𝑖 − 𝑇0𝑖 𝑐𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑜𝑛

𝑒0𝑖. . 𝑇1𝑖 −𝑇0𝑖

q0 𝐾. 𝐴. ∆x h, T
x

    
0
1 2 3 4

Fig V.4 schéma explicite par la conduction aux limites (cas de la convection) côte
gauche d‟une paroi plane

𝑇1𝑖 −𝑇0𝑖 ∆𝑥 ∆𝑥 𝑇0𝑖+1 −𝑇0𝑖

𝐾. 𝐴. + 𝑕. 𝐴 𝑇∞𝑖 − 𝑇0𝑖 + 𝑒𝑚
𝑖.
. 𝐴. = ρ. CP 𝐴. … … . . V. 16
∆x 2 2 ∆t

𝑕∆𝑥 𝑕∆𝑥 𝜏. 𝑒𝑖.0 ∆𝑥2.

𝑇0𝑖+1 = 1 + 2. 𝜏 − 2. 𝜏 𝑇0𝑖 + 2. 𝜏. 𝑇1𝑖 + 2. 𝜏 . 𝑇∞ + … … . . V. 17
𝑘 𝑘 k
V.3. 1 Critères de la stabilité pour la méthode explicite
La méthode explicite est facile à utiliser mais elle souffre d‟un défaut indésirable qui
sévèrement son utilisation, la méthode explicite n‟est pas conditionnellement stable.

51
*Si ∆t n‟est pas assis petit, la solution obtenues par la méthode explicite peuvent
osciller et divergent de la solution réelle.

-pour entier tels divergence, la valeur de ∆t dont être inférieur à une certaine limite
établie par le critère de stabilité

𝛼∆𝑡 1
𝜏= ≤
∆𝑥 2 2
Exemple V.1 :

Considérer une plaque large d'uranium initialement à température uniforme est

soumise à une isolation d'un côté et à une convection de l'autre côté. La formulation
des différences finies transitoires de ce problème doit être obtenue, et les températures
nodales après 5 min et dans des conditions stables doivent être déterminées

Hypothèses

1- Le transfert de chaleur est unidimensionnel car la plaque

Insulate g
est grande par rapport à son épaisseur. h, T
d x
2- La conductivité thermique est constante. 
0 1 2 3 4
3- Le transfert de chaleur par rayonnement est négligeable.

Propriétés La conductivité et la diffusivité sont données pour être k = 28 W/m°C et

  12.5  106 m2 / s .

L'espacement nodal est donné ∆x = 0,02 m. Ensuite, le nombre de nœuds devient

M  L / x  1 = 0,08 / 0,02 + 1 = 5. Ce problème implique 5 températures nodales
inconnues, et donc nous devons avoir 5 équations. Le nœud 0 est sur une frontière
isolée, et nous pouvons donc le traiter comme une note intérieure en utilisant le
concept d'image miroir. Les nœuds 1, 2 et 3 sont des nœuds intérieurs, et donc pour
eux, nous pouvons utiliser la relation de différence finie explicite générale exprimée
comme.

a) la Méthode explicite

emi x 2 Tmi 1  Tmi

Tmi 1  2Tmi  Tmi 1  
k 
𝑖.
𝑒𝑚 2
EqV. 13 : 𝑇𝑚𝑖+1 =𝜏 𝑇𝑚𝑖 −1 +𝑇𝑚𝑖 +1 + 1 − 2𝜏 𝑇𝑚𝑖 + 𝜏 . ∆𝑥 … … . V. 13
𝑘

52
L'équation aux différences finies pour le nœud 4 sur la surface droite soumise à
convection est obtenue en appliquant un bilan énergétique sur l'élément demi-volume
autour du nœud 4 et en prenant la direction de tous les transferts de chaleur vers le
nœud considéré:

e . x 2
Noeud 0 (isolé) : T0i 1   (T1i  T1i )  (1  2 )T0i  
k
e . x 2
Noeud 1 (intérieur) : T1i 1   (T0i  T2i )  (1  2 )T1i  
k
e x 2
.
Noeud 2 (intérieur) : T2   (T1  T3 )  (1  2 )T2  
i 1 i i i

k
e x 2
.
Noeud 3 (intérieur) : T3   (T2  T4 )  (1  2 )T3  
i 1 i i i

k
T3  T4
i i
. x x T4i 1  T4i
Noeud 4 (intérieur) : h(T  T4 )  k
i
e  C
x 2 2 t

𝑇4𝑖 −𝑇3𝑖 𝑖 𝑖 𝑖.
∆𝑥 ∆𝑥 𝑇4𝑖+1 −𝑇4𝑖
𝐾. 𝐴. + 𝑕. 𝐴 𝑇∞ − 𝑇4 + 𝑒4. 𝐴. = ρ. CP 𝐴.
∆t 2 2 ∆t
En divisant cette équation par k.A/∆𝑥. 2

∆𝑥 𝑖 ∆𝑥 𝑇4𝑖+1 −𝑇4𝑖
2𝑕. 𝑇 − 𝑇4𝑖 + 2 𝑇3𝑖 − 𝑇4𝑖 + 𝑒2.𝑖. =
2 ∞ 𝑘 𝜏
Ou x  0.02 m, g 0  10 6 W/m3 , k  28 W/m C, h  35 W/m2  C, T  20C , et   12.5 10 6
m2/s.

La limite supérieure du pas de temps ∆t est déterminée à partir des critères de stabilité
qui nécessitent que tous les coefficients primaires soient supérieurs ou égaux à zéro.
Le coefficient le plus petit dans ce cas, et donc les critères de stabilité pour ce
problème peuvent être exprimés comme.

Le critère de stabilité

𝑕∆𝑥 1 𝛼∆𝑡
1 + 2. 𝜏 − 2. 𝜏 ≥0→ 𝜏≤ → 𝜏=
𝑘 2 1+
𝑕∆𝑥 ∆𝑥 2
𝑘

(0.02 m) 2
t   15.6 s
2(12.5 10 6 m 2 /s)[1  (35 W/m2 .C)( 0.02 m) /(28 W/m.C)]

53
 ∆𝑥 2
∆t ≤ 𝑕 ∆𝑥 ≤ 15.55
2𝛼 1+
𝑘

Par conséquent, tout pas de temps inférieur à 15,5 s peut être utilisé pour résoudre ce
problème. Pour plus de commodité, choisissons le pas de temps ∆t = 15 s. Ensuite, le
nombre de Fourier maillé devient

𝛼∆𝑡
En prenant ∆t =15s ⇒ 𝜏 = = 0.45875
∆𝑥 2

Remarque : n‟importe quel ∆t≤ 15.5𝑠 peut être choisir pour résoudre ce problèmes

En substituant cette valeur et d'autres quantités données, les températures nodales

après 5x60 / 15 = 20 pas de temps (5 min) sont déterminées comme étant

Après 5 min: T0 =228.9C, T1 =228.4C, T2 =226.8C, T3 =224.0C, et T4

=219.9 C
(b) Le temps nécessaire pour un fonctionnement transitoire à établir est déterminé en
augmentant le nombre de pas de temps jusqu'à ce que les températures nodales ne
changent plus. Dans ce cas, un fonctionnement stable est établi en ---- min, et les
températures nodales dans des conditions stables sont déterminées comme étant

T0 =2420C, T1 =2413C, T2 =2391C, T3 =2356C, et T4 =2306 C

Discussion
La solution stable peut être vérifiée indépendamment en obtenant la formulation de
différence finie stable et en résolvant les équations résultantes simultanément.

Time [s] T1 [C] T2 [C] T3 [C] T4 [C] T5 [C] Rangée

0 100 100 100 100 100 1
15 106.7 106.7 106.7 106.7 104.8 2
30 113.4 113.4 113.4 112.5 111.3 3
45 120.1 120.1 119.7 119 117 4
60 126.8 126.6 126.3 125.1 123.3 5
75 133.3 133.2 132.6 131.5 129.2 6
90 139.9 139.6 139.1 137.6 135.5 7
105 146.4 146.2 145.4 144 141.5 8
120 152.9 152.6 151.8 150.2 147.7 9
135 159.3 159.1 158.1 156.5 153.7 10
… … … … … … …
… … … … … … …
3465 1217 1213 1203 1185 1160 232
3480 1220 1216 1206 1188 1163 233
3495 1223 1220 1209 1192 1167 234
3510 1227 1223 1213 1195 1170 235

54
 3525 1230 1227 1216 1198 1173 236
3540 1234 1230 1219 1201 1176 237
3555 1237 1233 1223 1205 1179 238
3570 1240 1237 1226 1208 1183 239
3585 1244 1240 1229 1211 1186 240
3600 1247 1243 1233 1214 1189 241
Tableau V.1
Variation de 𝑇1 , 𝑇2 , 𝑇3 , 𝑇4 ,𝑇5 , 𝑇6 , 𝑇7 et 𝑇8 avec le temps obtenues par la
méthode explicit
V.4.Conduction thermique transitoire 2-D

Elément de volume
m, n+1
n+1

Δy

m-1, n m, n m+1, n
n

Δy

m, n-1
n-1
Δx Δx
m-1 m M+1

Fig V.5.L‟élément de volume d‟un nœud interne (m, n) pour la conduction thermique
transitoire en coordonnés cartésiennes.

Considérer une région rectangulaire dans laquelle la conduction thermique est

importante dans ∆Z=1 dans la direction Z.

La chaleur peut être généré dans le milieu à un taux de 𝑒 . (x, y, t) la quelle peut varier
avec le temps et la position avec la conductivité thermique k du milieu supposer
constant.

diviser le plan x,y de la région en maille rectangulaire de nœuds (espace) ayant des
nœuds uniforme ∆x et ∆y, et considérer un nœud intérieur général(m,n) dont sec
cordonnée sont x = m.∆x et y = m.∆y (voir fig V.5).

Noter que l‟élément de volume centré au tour de l‟intérieur général

V élément=∆x. ∆y.1=∆x. ∆y

55
La formulation des différance finies en régime transitoire pour en nœuds intérieure
général sur la base de l‟équation

𝑇𝑚𝑖+1 −𝑇𝑚𝑖
𝑡𝑜𝑢𝑡𝑒𝑠 𝑄. + Egen
.
. élément = ρ. Vélément . CP
𝑙𝑒𝑠 𝑠𝑢𝑟𝑓𝑒𝑐𝑒 ∆t

𝑇𝑚 −1,𝑛 −𝑇𝑚 ,𝑛 𝑇𝑚 +1,𝑛 −𝑇𝑚 ,𝑛 𝑇𝑚 ,𝑛−1 −𝑇𝑚 ,𝑛

𝐾. ∆𝑦. + 𝐾. ∆𝑦. + 𝐾. ∆𝑥.
∆t ∆t ∆t
𝑇𝑚 ,𝑛+1 −𝑇𝑚 ,𝑛 𝑇𝑚𝑖+1 𝑖
,𝑛 −𝑇𝑚 ,𝑛
+ 𝐾. ∆𝑥. +e.gen élément . ∆x. ∆y = ρ. CP A. ∆x
∆t ∆t
Si la ∆x.∆y=L.

𝑒 𝑚 ,𝑛 𝐿 2 𝑇𝑚𝑖+1 𝑖
,𝑛 −𝑇𝑚 ,𝑛
𝑇𝑚 −1,𝑛 − 4 𝑇𝑚 ,𝑛 + 𝑇𝑚 +1,𝑛 +𝑇𝑚 −1,𝑛 + 𝑇𝑚 +1,𝑛 + =
𝑘 τ

Avec:

𝑘
𝛼= (𝐷𝑖𝑓𝑓𝑢𝑠𝑖𝑣𝑖𝑡é 𝑡𝑕𝑒𝑟𝑚𝑖𝑞𝑢𝑒)
ρ. CP
𝛼 ∆t
𝜏= (𝑛𝑜𝑚𝑏𝑟𝑒 𝑑𝑒 𝐹 𝑜𝑢𝑟𝑖𝑒𝑟)
∆𝑥 2

 Méthode explicite :
.𝑖
𝑖 𝑖 𝑖 𝑖 𝑖 𝑒𝑛𝑜𝑒𝑢𝑑 .𝐿2 𝑇𝑚𝑖+1 𝑖
,𝑛 −𝑇𝑚 ,𝑛
𝑇𝑙𝑒𝑓𝑡 + 𝑇𝑟𝑖𝑔 𝑕𝑡 +𝑇𝑏𝑜𝑡𝑡𝑜𝑚 + 𝑇𝑡𝑜𝑝 − 4𝑇𝑛𝑜𝑒𝑢𝑑𝑒 + =
𝑘 2

.𝑖
𝑖+1 𝑖 𝑖 𝑖 𝑖 𝑖 τ 𝑒𝑛𝑜𝑒𝑢𝑑 .𝐿2
𝑇𝑛𝑜𝑒𝑢𝑑𝑒 = 𝜏 . (𝑇𝑙𝑒𝑓𝑡 + 𝑇𝑟𝑖𝑔 𝑕𝑡 +𝑇𝑏𝑜𝑡𝑡𝑜𝑚 + 𝑇𝑡𝑜𝑝 ) + (1-4 𝜏)𝑇𝑛𝑜𝑒𝑢𝑑𝑒 + 𝑘

m=1, 2, 3,……M-1 valable pour les nœuds internes.

n=1, 2, 3,……N-1.

Critère de stabilité :

1-4.𝜏 ≥ 0
𝛼 ∆t 1
𝜏= ≤
∆𝑥 2 𝛼
Remarque : la méthode implicite est inconditionnellement stable, n‟importe quelle
valeur de ∆t peut être utilisé à la méthode implicite.

Méthode implicite :

En exprimant l‟Eq 5.19 à l‟incrément de temps i+1 le côté gauche au bien i, on

obtient la formulation implicite

56
 .𝑖+1
𝑖+1 𝑖+1 𝑖+1 𝑖+1 𝑖+1 𝑒𝑛𝑜𝑒𝑢𝑑 .𝐿2 𝑇𝑚𝑖+1 𝑖
,𝑛 −𝑇𝑚 ,𝑛
𝑇𝑙𝑒𝑓𝑡 + 𝑇𝑟𝑖𝑔 𝑕𝑡 +𝑇𝑏𝑜𝑡𝑡𝑜𝑚 + 𝑇𝑡𝑜𝑝 − 4𝑇𝑛𝑜𝑒𝑢𝑑𝑒 + =
𝑘 τ

Exemple V.2 :

La conduction de chaleur à travers une barre longue est solide en forme de L avec des
conditions aux limites spécifiées est considérée. La température au coin supérieur
(noeud # 3) du corps après 2, 5 et 30 min doit être déterminée avec la méthode des
différences finies explicites transitoires.

Hypothèses
h, T
1- Le transfert de chaleur à travers le corps est considéré 1 2
 3

Comme transitoire et bidimensionnel. qL isolé

4
 5 6 7 8
2- La conductivité thermique est constante.

3- La génération de chaleur est uniforme.

140
Propriétés : La conductivité et la diffusivité sont données pour être k = 15 W/m°C et
.  106 m2 / s .
  32

Une analyse

L'espacement nodal est donné ∆x = ∆y = l = 0,015 m. Les équations aux différences

finies explicites sont déterminées sur la base du bilan énergétique pour le cas
transitoire exprimé en

T i 1  Tmi
 Q i  G element
All sides
i
 Velement C m
t

Les quantités h, 𝑇∞ et 𝑒 . ne changent pas avec le temps, et donc nous n'avons pas
besoin d'utiliser l'exposant i pour elles. De plus, les expressions du bilan énergétique
peuvent être simplifiées en utilisant les définitions de la diffusivité thermique et le
nombre de Fourier maillé sans dimension   t / l 2 où x  y  l . On note que tous

les nœuds sont des nœuds limites sauf le nœud 5 qui est un nœud intérieur. Par
conséquent, nous devrons nous fier aux bilans énergétiques pour obtenir les équations
aux différences finies. En utilisant les bilans énergétiques, les équations aux
différences finies pour chacun des 8 nœuds sont obtenues comme suit:

57
 l l l T2i  T1i l T4i  T1i . l
2
l 2 T1i 1  T1i
Noeud1: q L  h (T  T1 )  k
i
k e  C
2 2 2 l 2 l 4 4 t
l T1i  T2i l T3i  T2i T i  T2i l2 l 2 T i 1  T2i
Noeud2: hl (T  T2i )  k k  kl 5  e.   C 2
2 l 2 l l 2 2 t
l T2i  T3i l T6i  T3i l2 l 2 T i 1  T3i
Noeud3: hl (T  T3i )  k k  e.   C 3
2 l 2 l 4 4 t
 hl   hl e . 3l 2 
(Il peut être réorganisé comme T3
i 1
 1  4  4 T3i  2  T4i  T6i  2 T  )
 k  k 2k 

l T1i  T4i l 140  T4i T5i  T4i . l

2
l 2 T4i 1  T4i
Noeud4: q L l  k k  kl e  C
2 l 2 l l 2 2 t
 e .l 2 
Noeud 5 (intérieur): T5
i 1
 1  4 T5i    T2i  T4i  T6i  140  
 k 
Noeud6:
l T3i  T6i T5i  T6i 140  T6i l T7i  T6i . 3l
2
3l 2 T6i 1  T6i
hl (T  T )  k
i
 kl  kl k e  C
t
6
2 l l l 2 l 4 4

l T6i  T7i l T8i  T7i 140  T7i l2 l 2 T i 1  T7i

Noeud7 hl (T  T7i )  k k  kl  e.   C 7
2 l 2 l l 2 2 t
l l T7i  T8i l 140  T8i l2 l 2 T i 1  T8i
Noeud8: h (T  T8i )  k k  e.   C 8
2 2 l 2 l 4 4 t
Ou e.  2  10 7 W/m 3 , q L  8000 W/m 2 , l = 0.015 m, k =15 W/mC, h = 80

W/m2C, et T = 25C.

La limite supérieure du pas de temps ∆t est déterminée à partir des critères de stabilité
qui nécessitent que le coefficient de Tmi dans Tmi1 l'expression (le coefficient primaire)
soit supérieur ou égal à zéro pour tous les nœuds. Le plus petit coefficient primaire
dans les 8 équations ci-dessus est le coefficient de T3i dans T3i 1 l'expression car il est
exposé à la plupart des convections par unité de volume (cela peut être vérifié), et
donc les critères de stabilité pour ce problème peuvent être exprimés comme :

hl 1 l2
1  4  4 0    t 
k 4(1  hl / k ) 4 (1  hl / k )
Depuis   t / l 2 En substituant les quantités données, la valeur maximale admissible
du pas de temps est déterminée comme étant

(0.015 m) 2
t   16.3 s
4(3.2 10 6 m 2 /s)[1  (80 W/m 2 .C)(0.015 m) /(15 W/m. C)]

58
Par conséquent, tout pas de temps inférieur à 16,3 s peut être utilisé pour résoudre ce
problème. Pour plus de commodité, nous ne choisissons que le pas de temps soit ∆t =
15 s. Ensuite, le nombre de Fourier maillé devient

t (3.2  10 6 m 2 /s) (15 s)

   0.2133 (Pourt = 15 s)
l2 (0.015 m) 2
En utilisant la condition initiale spécifiée comme solution au temps t = 0 (pour
i = 0), le balayage des 9 équations ci-dessus donnera la solution à des intervalles de 15
s. À l'aide d'un ordinateur, la solution au nœud du coin supérieur (nœud 3) est
déterminée à 441, 520 et 529 ° C à 2, 5 et 30 min, respectivement. On peut montrer
que la solution d'état stationnaire au nœud 3 est 531 ° C.

Time T1 [C] T2 [C] T3 [C] T4 [C] T5 [C] T6 [C] T7 [C] T8 [C] Rangée
[s]
0 140 140 140 140 140 140 140 140 1
15 203.5 200.1 196.1 207.4 204 201.4 200.1 200.1 2
30 265 259.7 252.4 258.2 253.7 243.7 232.7 232.5 3
45 319 312.7 300.3 299.9 293.5 275.7 252.4 250.1 4
60 365.5 357.4 340.3 334.6 326.4 300.7 265.2 260.4 5
75 404.6 394.9 373.2 363.6 353.5 320.6 274.1 267 6
90 437.4 426.1 400.3 387.8 375.9 336.7 280.8 271.6 7
105 464.7 451.9 422.5 407.9 394.5 349.9 286 275 8
120 487.4 473.3 440.9 424.5 409.8 360.7 290.1 277.5 9
135 506.2 491 456.1 438.4 422.5 369.6 293.4 279.6 10
… … … … … … … … … …
… … … … … … … … … …
1650 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 111
1665 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 112
1680 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 113
1695 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 114
1710 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 115
1725 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 116
1740 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 117
1755 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 118
1770 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 119
1785 596.3 575.7 528.5 504.6 483.1 411.9 308.8 288.9 120
Tableau V.2
Variation de 𝑇1 , 𝑇2 , 𝑇3 , 𝑇4 , 𝑇5 , 𝑇6 , 𝑇7 et 𝑇8 avec le temps obtenues par la
méthode explicit

59
Bibliographie
1. H.K. Versteeg, W Malalasekera, An introduction to computational fluid dynamics:
the finite volume method, Longman (1998).

2. S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington
D.C (1980),

3. D.R. Croft, D.G. Lilley, Heat transfer calculations using finite difference equations,
Elsevier science and technologie (1977)

4. W.J. Minckowicz, E.M. Sparrow, J.Y. Murthy, Handbook of numerical heat

transfer, Wiley (2006).

5. C.A.J. Fletcher, Computational techniques for fluid dynamics 1, Springer (1988).

6. M. Boumahrat et J. Gourdin, „Méthodes Numériques Appliquées‟, Office des

Publications Universitaires (OPU), 1983.
7. M. Necati Ozisik, “Finite Difference Methds in Heat Transfer”; Mechanical and
Aerospace Engineering Department North Carolina State University

60