Mna 28 12 2022

Edition 2022 L.
MESSAOUDI
METHODES NUMERIQUES
APPLIQUEES
Tome I : Cours
Université MB Batna 2
Maître de Conféren es.
Département de Mé anique.
Fa ulté de Te hnologie.
Université Mustafa Benboulaïd Batna 2.
- -
Table des matières
Liste des gures v

Liste des tableaux vii
1 Notions sur les équations aux dérivées partielles 1
1.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Dénitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.3 Classi ation des EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3.1 Cara téristiques des EDP . . . . . . . . . . . . . . . . . . . . . 4
1.4 Diérents types de onditions aux limites (C.L) . . . . . . . . . . . . . 5

er
1.4.1 Conditions de Diri hlet (ou 1 type) . . . . . . . . . . . . . . . 5
ème
1.4.2 Conditions de Neumann (ou 2 type) . . . . . . . . . . . . . 6
ème
1.4.3 Conditions mixtes (ou 3 type) . . . . . . . . . . . . . . . . . 6
1.5 Exer i es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2 Méthodes Analytiques 7
2.1 Introdu tion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Prin ipe de superposition . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.3 Méthode de séparation des variables . . . . . . . . . . . . . . . . . . . . 8
2.3.1 Equation de diusion ave onditions de Diri hlet . . . . . . . . 9
2.3.2 Equation de Lapla e en oordonnées re tangulaires . . . . . . . 12
2.3.3 Equation d'onde . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.4 Exer i es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3 Méthodes des diéren es nies 19

3.1 Introdu tion sur les méthodes numériques . . . . . . . . . . . . . . . . 19
3.2 Prin ipe et bases de la MDF . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3 Formule de Taylor en oordonnées re tangulaires . . . . . . . . . . . . . 25
3.4 S hémas de dis rétisation . . . . . . . . . . . . . . . . . . . . . . . . . . 25
i
TABLE DES MATIÈRES
3.5 Notions de stabilité, onsistan e et onvergen e . . . . . . . . . . . . . 25
3.5.1 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.5.2 Consistan e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.5.3 Convergen e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.6 Méthodes de dis rétisation . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.6.1 Méthode expli ite . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.6.2 Méthode impli ite . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.6.3 Méthode de Crank-Ni holson . . . . . . . . . . . . . . . . . . . 28
3.7 Appli ation de la méthode expli ite à l'équation de diusion . . . . . . 29
3.7.1 Consistan e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.7.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.8 Appli ation de la méthode impli ite à l'équation de diusion . . . . . . 31
3.8.1 Consistan e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.8.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.9 Formulation matri ielle du s héma expli ite appliqué à l'équation de
diusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.10 Formulation matri ielle du s héma impli ite appliqué à l'équation de
diusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.11 Exemple résolu par wxMaxima . . . . . . . . . . . . . . . . . . . . . . 36

ème
3.12 Cas des onditions aux limites du 2 type . . . . . . . . . . . . . . . 38
3.12.1 Cas de l'approximation entrée . . . . . . . . . . . . . . . . . . 39
3.12.2 Cas de l'approximation dé entrée d'ordre 1 . . . . . . . . . . . . 41
3.12.3 Cas de l'approximation dé entrée d'ordre 2 . . . . . . . . . . . . 43
3.13 Equations non linéaires . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.13.1 Dis rétisation par un s héma expli ite . . . . . . . . . . . . . . 45
3.13.2 Dis rétisation par un s héma impli ite . . . . . . . . . . . . . . 46
3.14 Problèmes multidimensionnels . . . . . . . . . . . . . . . . . . . . . . . 46
3.14.1 Dis rétisation par un s héma expli ite . . . . . . . . . . . . . . 46
3.14.2 Dis rétisation par un s héma impli ite . . . . . . . . . . . . . . 47
3.15 Equation de Lapla e . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.15.1 Dis rétisation par la formule à 5 points (Runge 1908) . . . . . . 48
3.15.2 Conditions limites de Diri hlet variables . . . . . . . . . . . . . 52
3.15.3 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . 53
3.15.4 Conditions aux limites de Neumann . . . . . . . . . . . . . . . . 54
3.15.5 Equation de Poisson . . . . . . . . . . . . . . . . . . . . . . . . 59
3.15.6 Méthode d'ordre élevé . . . . . . . . . . . . . . . . . . . . . . . 60
3.15.7 Cas du domaine non re tangulaire . . . . . . . . . . . . . . . . . 62
3.15.8 Equation non-linéaire et problème tridimensionnel . . . . . . . . 64
3.16 Equation d'onde . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
- ii -
TABLE DES MATIÈRES
3.16.1 Résolution de l'équation d'onde . . . . . . . . . . . . . . . . . . 64
3.16.1.1 Dis rétisation . . . . . . . . . . . . . . . . . . . . . . . 65
3.16.1.2 Stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.16.1.3 Autres s hémas de dis rétisation . . . . . . . . . . . . 66
3.16.2 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . 67
3.17 Exer i es
68
4 Méthode des Volumes Finis 81

4.1 Equations de onservation . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.2 Equation de diusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.2.1 Problème de diusion unidimenionnel . . . . . . . . . . . . . . . 85
4.2.2 Exemple 1 : Equation de diusion 1D . . . . . . . . . . . . . . . 87
4.2.3 Exemple 2 : Equation de diusion 1D ave terme sour e . . . . 89
4.2.4 Exemple 3 : Equation de diusion 1D ave onve tion à la surfa e
et bord isolé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.2.5 Problème de diusion bidimensionnel . . . . . . . . . . . . . . . 93
4.2.6 Problème tridimensionnel . . . . . . . . . . . . . . . . . . . . . 99
4.2.7 Problèmes symétriques . . . . . . . . . . . . . . . . . . . . . . . 100
4.2.8 Règles de base à satisfaire . . . . . . . . . . . . . . . . . . . . . 100
4.3 Problème de onve tion-diusion . . . . . . . . . . . . . . . . . . . . . 101
4.3.1 S héma entré . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.3.2 S héma dé entré (Upwind) . . . . . . . . . . . . . . . . . . . . . 103
4.3.3 S héma exponentiel . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.3.4 S héma Hybride . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.3.5 S héma Power Law . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.3.6 Forme générale . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.3.7 Exemple 1 : S héma entré . . . . . . . . . . . . . . . . . . . . . 105
4.3.8 Exemple 2 : S héma Upwind . . . . . . . . . . . . . . . . . . . . 111
4.4 Equation de diusion 1D instationnaire . . . . . . . . . . . . . . . . . . 115
4.4.1 S héma expli ite . . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.4.2 S héma impli ite . . . . . . . . . . . . . . . . . . . . . . . . . . 116
4.4.3 S héma de Crank-Ni holson . . . . . . . . . . . . . . . . . . . . 116
4.5 Exer i es . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
Bibliographie 121
A Formulation à 9 points 125
B Approximation des dérivées aux interfa es 129
- iii -
TABLE DES MATIÈRES
- iv -
Liste des gures
1.1 Etapes de la simulation numérique. . . . . . . . . . . . . . . . . . . . . 2
1.2 Domaine d'étude pour les problèmes elliptiques. . . . . . . . . . . . . . . 4
1.3 Domaine d'étude pour les problèmes paraboliques et hyperboliques. . . . . 5
2.1 Problème original. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.2 Problème modié. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.3 Dis tribution de la température en fon tion du temps et de l'espa e. . . . 12
2.4 Equation de Lapla e dans un re tangle. . . . . . . . . . . . . . . . . . . 13
2.5 Distribution de la température en fon tion de x et de y. . . . . . . . . . 15
2.6 Variation du dépla ement en fon tion du temps et de l'espa e. . . . . . . 17
3.1 Passage du domaine ontinu au domaine dis ret. . . . . . . . . . . . . . 20
3.2 Approximation par séries de Taylor de la fon tion ex . . . . . . . . . . . 23
3.3 Diérents types de s hémas de dis rétisation. . . . . . . . . . . . . . . . 25
3.4 S héma expli ite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.5 S héma impli ite. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.6 S héma de Crank-Ni holson. . . . . . . . . . . . . . . . . . . . . . . . . 28
3.7 Dis rétisation du domaine pour l'équation de la haleur. . . . . . . . . . 33
3.8 Dis rétisation du domaine spatio-temporel. . . . . . . . . . . . . . . . . 37
3.9 Conditions aux limites de type Neumann. . . . . . . . . . . . . . . . . . 38
3.10 Dis rétisation du domaine pour l'équation de Lapla e. . . . . . . . . . . 48
3.11 S héma ellulaire de la formulation à 5 points : (a) β = 1, (b) β6= 1. . . . 49
3.12 Equation de Lapla e. Exemple 1. . . . . . . . . . . . . . . . . . . . . . . 49
3.13 Equation de Lapla e. Exemple 3. . . . . . . . . . . . . . . . . . . . . . . 51
3.14 Equation de Lapla e. (C.L) de Neumann. . . . . . . . . . . . . . . . . . 54
3.15 Résolution de l'équation de Poisson. Répartition des isothermes. . . . . 60
3.16 Représentation s hématique de la formulation à 9 points (β = 1). . . . . 61
3.17 Formulation à 9 points. Exemple 1. . . . . . . . . . . . . . . . . . . . . 61
3.18 Domaine physique non re tangulaire. . . . . . . . . . . . . . . . . . . . 63
v
LISTE DES FIGURES
3.19 Transformation du domaine physique. . . . . . . . . . . . . . . . . . . . 63
3.20 Exo-3-26. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
3.21 Exo-3-31. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
3.22 Exo-3-34. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.23 Exo-3-36. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.24 Exo-3-37. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
3.25 Exo-3-38. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.26 Exo-3-39. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
3.27 Exo-3-45. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
3.28 Exo-3-46. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
4.1 Equations de onservation. . . . . . . . . . . . . . . . . . . . . . . . . . 82

4.2 Disposition du (v. ) 1D. . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.3 Exemple 1 : Equation de diusion 1D. . . . . . . . . . . . . . . . . . . . 87
4.4 Exemple 2 : Equation de diusion 1D ave terme sour e. . . . . . . . . 89
4.5 Exemple 3 : Equation de diusion 1D, onve tion en surfa e et bord isolé. 91
4.6 Disposition du (v. ) 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.7 Exemple 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
4.8 Problème ondu to- onve tif 1D. . . . . . . . . . . . . . . . . . . . . . . 101
4.9 Exemple 1 : S héma entré. . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.10 Exemple1- Cas 1 : u=0,1 m/s et ∆x=0,2 m. . . . . . . . . . . . . . . . 108
B.1 Approximation de la dérivée aux interfa es. . . . . . . . . . . . . . . . . 130
- vi -
Liste des tableaux
3.1 Approximation entrée en O(h2 ). . . . . . . . . . . . . . . . . . . . . . . 21
3.2 Approximation dé entrée avant du 1er ordre O(h). . . . . . . . . . . . . . 22
3.3 Approximation dé entrée arrière du 1er ordre O(h) . . . . . . . . . . . . 22
3.4 Approximation dé entrée avant du 2nd ordre O(h2 ). . . . . . . . . . . . . 22
3.5 Approximation dé entrée arrière du 2nd ordre O(h2 ). . . . . . . . . . . . 22
3.6 Inuen e négative de la simple pré ision. . . . . . . . . . . . . . . . . . 23
3.7 Inuen e positive de la double pré ision. . . . . . . . . . . . . . . . . . . 24
3.8 Inuen e de l'approximation. . . . . . . . . . . . . . . . . . . . . . . . . 24
4.1 Fon tion à prendre pour A (|Pe |). . . . . . . . . . . . . . . . . . . . . . 105
vii
LISTE DES TABLEAUX
- viii -
Chapitre 1
Notions sur les équations aux dérivées
partielles
1.1 Introdu tion

Les s ientiques et les ingénieurs utilisent plusieurs te hniques pour la résolution des
problèmes de hamp (diusion de la haleur, propagation d'ondes ...et ). Ces te hniques
peuvent être expérimentales, analytiques ou numériques.
Les méthodes expérimentales sont très hers, prennent beau oup de temps et dans
ertains as, elles sont hasardeuses et même dangereuses. Elles ne permettent pas sou-
vent une grande exibilité des paramètres de variation.
La plupart des méthodes analytiques ne s'appliquent que dans des as limités.
Pour des problèmes relatifs à des systèmes de forme géométrique omplexe ou à des
milieux à ara téristiques non uniformes ou non isotropes, qui est le as de la plupart
des problèmes ren ontrés en pratique, il est né essaire de faire appel aux méthodes
numériques.
Les problèmes ren ontrés dans le domaine des s ien es de l'ingénieur sont sou-
vent représentés (ou modélisés) par des équations aux dérivées partielles (EDP) qui
modélisent les phénomènes physiques présents (é oulement de uides, transfert de ha-
leur, vibration de stru tures, propagation d'ondes, hamp éle tromagnétique ...et ).
Les étapes depuis la modélisation jusqu'à la programmation sont s hématisées sur la
gure (Fig.1.1).
1.2 Dénitions
Une équation aux dérivées partielles (EDP) est une équation faisant intervenir
une fon tion in onnue de plusieurs variables ainsi que ertaines de ses dérivées
partielles.
1
Ch1.NOTIONS SUR LES ÉQUATIONS AUX DÉRIVÉES PARTIELLES.
Fig. 1.1: Etapes de la simulation numérique.
- 2 -
1.3. Classi ation des EDP.
Une équation diérentielle ordinaire (EDO) ontient seulement des dérivées par
rapport à une seule variable.
On appelle ordre d'une EDP, l'ordre de la plus grande dérivée présente dans
l'équation.
Exemples :
nd
1. L'équation i- ontre est une EDO du 2 degré : α f ′′ (x) + β f ′ (x) + γ f (x) = λ.
nd
2. L'équation i-dessous est une EDP du 2 degré :
∂2φ ∂2φ ∂2φ ∂φ ∂φ

a 2
+ b + c 2
+d +e +f φ=g (1.1)
∂x ∂x ∂y ∂y ∂x ∂y
Elle peut aussi s'é rire sous la forme :
a φxx + b φxy + c φyy + d φx + e φy + f φ = g ⇐⇒ Lφ = g
Une EDP est linéaire si l'équation est linéaire par rapport aux dérivées par-
tielles de la fon tion in onnue ( -à-d qu'il n'existe pas de produit des variables
dépendantes et/ou produit de ses dérivées).
Une EDP est non linéaire si l'équation ontient un produit des variables dépen-
dantes et/ou un produit de ses dérivées.
Exemple s:
1- L'équation i-dessous est une EDP linéaire du se ond degré.
∂φ(x,t) ∂ 2 φ(x,t) ∂φ(x,t)

∂t
+ 12 k1 x2 ∂x2
+ k2 x ∂x
− k3 φ(x, t) = 0
2- L'équation i-dessous est une EDP non linéaire du premier degré.
∂φ(x,t) ∂φ(x,t)
∂t
+ φ(x, t) ∂x
= g(x, t)
On appelle problème aux limites une EDP munie de onditions aux limites (C.L)
sur la totalité de la frontière du domaine sur lequel elle est posée.
Un problème est bien posé si pour toute donnée (se ond membre, domaine,
données aux bords, . . . ), il admet une solution unique et si ette solution dépend
ontinûment de la donnée.
1.3 Classi ation des EDP

Considérons les équations aux dérivées partielles du se ond degré de type (1.1).
Leur lassi ation se fait de la manière suivante :
- 3 -
L'EDP est elliptique si : ∆ = b2 − 4 a c < 0

L'EDP est parabolique si : ∆ = b2 − 4 a c = 0
L'EDP est hyperbolique si : ∆ = b2 − 4 a c > 0
Les termes elliptique, parabolique et hyperbolique sont dérivés de l'équation qua-
dratique :
a x2 + b xy + c y 2 + d x + e y + f = 0
qui représente une hyperbole, une parabole ou une ellipse selon la valeur de ∆.
Exemples :
∂2u ∂2u
Equation d'onde :
∂t2
= c2 ∂x2
a = c2 , b = 0, c = −1 ⇒ ∆ = 4c2 don 'est une EDP de type hyperbolique.
∂u ∂2u
Equation de diusion ( haleur ) :
∂t
= ∂x2
a = 1, b = 0, c = 0 ⇒ ∆ = 0 don 'est une EDP de type parabolique.
∂2u ∂2u
L'équation :
∂x2
+ ∂y 2
+λu = f
a = 1, b = 0, c = 1 ⇒ ∆ = −4 don 'est une EDP de type elliptique.
2
∂2u
- L'équation de Lapla e :
∂x2
+ ∂∂yu2 =0 dans e as u est dite harmonique.
∂2u 2
- L'équation de Poisson :
∂x2
+ ∂∂yu2 =f
∂2u ∂2u
Equation de Tri omi :
∂x2
+x ∂y 2
=0
a = 1, b = 0, c = x ⇒ ∆ = −4x don 'est une EDP de type mixte. Elle est elliptique
si x > 0 et hyperbolique si x < 0.
1.3.1 Cara téristiques des EDP

Les EDP elliptiques sont asso iées ave des phénomènes stationnaires. Elles
modélisent souvent un problème intérieur (fermé ) (Fig.1.2) et alors la solution
est souvent bornée (exp : les équations de Lapla e et de Poisson).
(C.L)
Domaine fermé
Fig. 1.2: Domaine d'étude pour les problèmes elliptiques.
Les EDP hyperboliques apparaissent dans les problèmes de propagation. La
solution est souvent ouverte et elle avan e alors indéniment vers l'extérieur
à partir d'une ondition initiale satisfaisant toujours les onditions aux limites
(C.L) spé iées (Fig.1.3).
- 4 -
1.4. Diérents types de onditions aux limites (C.L).
Propagation de la solution
(C.L)
(C.L)
Domaine ouvert
(C.I)
Fig. 1.3: Domaine d'étude pour les problèmes paraboliques et hyperboliques.
Les EDP paraboliques sont généralement asso iées ave des problèmes dans les-
quels la quantité qui nous intéresse varie lentement omparée aux mouvements
aléatoires qui produisent les variations. Comme l'EDP hyperbolique, le domaine
de la solution de l'EDP parabolique est souvent ouvert. Les onditions initiale
(C.I) et limites (C.L) typiquement asso iées ave ette équation ressemblent à
eux des problèmes hyperboliques à l'ex eption de la seule ondition initiale à

er
t=0 né essaire puisque l'équation est du 1 ordre en temps. Aussi, les équa-
tions paraboliques et hyperboliques sont résolues par des te hniques similaires
alors que l'équation elliptique est souvent plus di ile et requiert des te hniques
diérentes.
1.4 Diérents types de onditions aux limites (C.L)

1.4.1 Conditions de Diri hlet (ou 1er type)
Dans e as, la valeur de la variable dépendante est spé iée le long de la frontière
du domaine.
φ(r) = 0, r sur S ondition homogène.
φ(r) = p(r), r sur S ondition non homogène.
r étant la variable indiquant la position et S la frontière du domaine.
- 5 -
1.4.2 Conditions de Neumann (ou 2ème type)

Dans e as, le gradient normal de la variable dépendante est spé ié le long de la
frontière du domaine.
∂φ(r)
∂n
= 0, r sur S ondition homogène.
∂φ(r)
∂n
= q(r), r sur S ondition non homogène.
∂φ(r)
∂n
étant la dérivée normale de φ le long de la frontière S du domaine.
1.4.3 Conditions mixtes (ou 3ème type)

Dans e as, une ombinaison des deux C.L pré édentes spé iée le long de la fron-
tière du domaine.
∂φ(r)
∂n
+ h(r) φ(r) = 0, r sur S ondition homogène.
∂φ(r)
∂n
+ h(r) φ(r) = w(r), r sur S ondition non homogène.
1.5 Exer i es
1-01 : Classer les EDP suivantes (parabolique, hyperbolique ou elliptique) :
∂u
1.
∂t
+ 2 ∂u
∂x
= 0.
∂2φ ∂2φ
2. (1 − M) ∂x2
+ ∂y 2
=0 (étudier en fon tion de M ).
∂2u 2
3.
∂x2
+ x ∂∂yu2 = 0.
∂T ∂2T ∂2T
4.
∂t
= ∂x2
+ ∂y 2
.
1-02 : Classer les EDP suivantes :

2 2 ∂2u
1.
∂ u
4 ∂∂yu2 − 4 ∂x∂y + ∂x2
+ 3 u = 0.
∂u ∂2u ∂2u
2.
∂x
+ ∂t2
+ ∂x2
= f (t).
∂2φ ∂ φ 2 ∂2φ
3. x2 ∂x2
− 2 xy ∂x∂y + y2 ∂y 2
+ x ∂φ
∂x
+y ∂φ
∂y
= 0.
∂2T ∂ T 2 2
4.
∂x2
+ 6 ∂x∂y − 16 ∂∂yT2 = 0.
- 6 -
Chapitre 2
Méthodes Analytiques
2.1 Introdu tion

La solution la plus satisfaisante d'un problème de hamp est ertainement la solu-
tion mathématique exa te. Cependant, dans beau oup de as pratiques, ette solution
analytique ne peut être obtenue et on a alors re ours aux méthodes numériques qui sont
approximatives. La solution analytique, si elle existe, est intéressante an de s'assurer
de la solution obtenue numériquement.
Parmi les méthodes analytiques les plus utilisées, on distingue :
Méthode de séparation des variables (MSV) ou méthode de Fourier.
Méthodes de perturbation.
Méthode des séries de puissan e.
Méthode analogiques.
Méthode des transformées de Lapla e . . . et .
Nous nous intéressons uniquement à la MSV qui est très puissante et d'un emploi
plus général que les autres méthodes qui sont spé iques à ertains type de problèmes
physiques.
Avant de détailler la MSV, nous allons, en premier lieu, voir un prin ipe très im-
portant :
2.2 Prin ipe de superposition

Si haque membre de l'ensemble des fon tions φn , n = 1, 2, · · · , N est solution
de l'EDP : Lφ = 0 (L : Opérateur diérentiel) ave des (C.L) données, alors la
ombinaison linéaire :
N
X
φ = φ0 + an φn
n=1
7
Ch2.MÉTHODES ANALYTIQUES.
satisfait aussi l'équation L φ = g.
Théorème d'uni ité : Ce théorème garanti que la solution obtenue de l'EDP ave
des (C.L) données est la seule possible.
2.3 Méthode de séparation des variables

Le prin ipe de ette méthode est de her her les solutions d'une EDP sous la forme :
φ(x, y, z, t) = X(x) . Y (y) . Z(z) . T (t) (2.1)
Pour voir si la MSV indépendantes peut s'appliquer à un problème physique donnée,
nous devons onsidérer l'EDP dé rivant le problème, la forme du domaine de la solution
et les (C.L).
Par exemple, pour appliquer ette méthode à un problème faisant intervenir deux
variables x et y, trois hoses doivent être onsidérées :
1. L'opérateur diérentiel L doit être séparable, -à-d. il doit être une fon tion de
φ(x, y) .
2. Toutes les (C.L) et (C.I) doivent être prises sur des surfa es à oordonnées
onstantes -à-d. x = C te , y = C te .
3. L'opérateur linéaire dénissant la (C.L) à x = C te (ou y = C te ) ne doit pas
ontenir des dérivées partielles en φ par rapport à y (ou x) et ses oe ients
doivent être indépendants de y (ou de x).
Exemples :
∂2φ ∂2φ ∂2φ

Lφ = ∂x2
+ ∂x ∂y
+ ∂y 2
viole la règle (1).
Si la région de la solution n'est pas un re tangle à fa es parallèles aux axes x et
y alors la règle (2) est violée.

∂φ
Ave une (C.L) sur une partie ou x=0 et
∂x
=0 sur une autre partie, la règle
(3) est violée.
Nous allons maintenant appliquer la MSV aux EDP les plus ren ontrés telles que
l'équation de diusion de la haleur, l'équation de Lapla e et l'équation d'onde en
respe tant les étapes suivantes :
1. Séparer les variables indépendantes.
2. Trouver les solutions parti ulières aux équations séparées satisfaisant les (C.L).
3. Combiner linéairement es solutions pour obtenir la solution générale (prin ipe
de superposition).
4. Satisfaire la (C.I) pour déterminer les oe ients restant.
- 8 -
2.3. Méthode de séparation des variables.
2.3.1 Equation de diusion ave onditions de Diri hlet

Soit une barre de longueur L et de très faible se tion. Considérons le transfert de
haleur entre ses extrémités maintenues à des températures xes θ1 , la température
initiale étant θ0 (Fig.2.1).
∂2θ ∂θ
α = 0≤x≤ L, t≥0 (2.2)
∂x2 ∂t



θ = θ0 t=0 et 0≤x≤L

θ = θ1 t>0 et x=0


θ = θ

t>0 et x=L
1
θ1
θ1
0 L x
θ0
Fig. 2.1: Problème original.
Pour fa iliter la résolution du système (2.2), il est re ommandé d'homogénéiser les
(C.L), pour ela, on pose :
θ − θ1 α x
T = , τ= t, X=
θ0 − θ1 L2 L
T, τ et X sont alors des variables adimensionnelles (sans unités).
Le système (2.2) devient alors (Fig.2.2) :
∂2T ∂T
2
= 0≤X ≤1, τ ≥0 (2.3)
∂X ∂τ



 T =1 τ =0 et 0≤X≤1

T =0 τ >0 et X =0


T = 0

τ >0 et X =1
- 9 -
T=0
T=0
0 1 X
T=1
Fig. 2.2: Problème modié.
Appliquons maintenant la MSV, on pose alors : T (X, τ ) = f (X) . g(τ )

le système (2.3) devient alors :
f ′′ (X) . g(τ ) = f (X) . g ′(τ )

 f (0) = 0
 f (1) = 0
En séparant les variables :
f ′′ (X) g ′ (τ )
=
f (X) g(τ )
et omme f ne dépend que de X et g ne dépend que de τ , e i ne peut être satisfait
que si les deux rapports sont égaux à une même onstante :
f ′′ (X) g ′ (τ )
= = C te
f (X) g(τ )
Pour des onsidérations physiques ainsi que pour éviter de trouver une solution
triviale au problème, ette onstante doit être égale un nombre négatif, soit −λ2 . Nous
obtenons alors deux équations diérentielles ordinaires (EDO) :

f ′′ + λ2 f = 0
g ′ + λ2 g = 0
La solution de la première équation est donnée par les deux formes :
f (X) = A ei λ X + B e−i λ X = C sin(λ X) + D cos(λ X)

f (0) = 0 =⇒ D = 0
- 10 -
f (1) = 0 =⇒ C sin(λ) = 0 =⇒ sin(λ) = 0 =⇒ λ = n π , n = 1, 2, · · · ∞.
Nous avons trouvé un ensemble inni de valeurs dis rètes de λ pour lesquelles
l'équation n'a pas de solutions triviales ; e i nous donne :
fn (X) = Cn sin(n π X)
La solution de la deuxième équation est donnée par la méthode des oe ients
indéterminés :
2 2 π2 τ
g(τ ) = K e− λ τ
= K e− n et par suite :
2 π2 τ
gn (τ ) = Kn e− n
En remplaçant les deux fon tions :
2 π2 τ
Tn (X, τ ) = fn (X) . gn(τ ) = Cn sin(n π X) . e− n
et en appliquant le théorème de superposition, nous aurons la solution sous la forme :
∞ ∞
2 π2
X X
T (X, τ ) = Tn (X, τ ) = Cn sin(n π X) . e−n τ
n=1 n=1
Maintenant déterminons la onstante en utilisant la ondition initiale :
∞
P
T (X, 0) = Cn sin(n π X) = 1
n=1
2
RL R1 2
Cn = L
T (X, 0) sin(n π X) dX = 2 sin(n π X) dX = nπ
[1 − (−1)n ]
0 0
La solution est triviale pour tout n paire, on ne prendra don que les valeurs im-
paires, pour ela on pose :
n = 2k + 1, k = 0, 1, 2, · · · ∞. Dans e as :
4
Ck =
(2k + 1) π
et par suite, la solution du problème est :
∞
4 X 1
sin [(2k + 1) π X] . e−[(2k+1) π τ ]
2 2
T (X, τ ) = (2.4)
π k=0 (2k + 1)
et en revenant aux variables dimensionnelles :
∞

(2k+1)2 π 2

4(θ0 − θ1 ) X 1 (2k + 1) π − α
L2
t
θ(x, t) = θ1 + sin x .e (2.5)
π k=0
(2k + 1) L
- 11 -
Les distributions de la température en fon tion du temps et de l'espa e sont repré-
sentées sur la gure (Fig.2.3).
Fig. 2.3: Dis tribution de la température en fon tion du temps et de l'espa e.
2.3.2 Equation de Lapla e en oordonnées re tangulaires
Soit une plaque re tangulaire soumise à des températures onstantes omme il est
indiqué sur la gure (Fig.2.4).
∂2T ∂2T
+ =0 0 ≤ x ≤ a ,0 ≤ y ≤ b (2.6)
∂x2 ∂y 2

T (0, y) = 0 , T (x, 0) = 0
T (a, y) = 0 , T (x, b) = T0
- 12 -
y
T0
T=0
T=0
0 a x
T=0
Fig. 2.4: Equation de Lapla e dans un re tangle.
Appliquons maintenant la MSV, on pose alors : T (x, y) = f (x) . g(y)

f ′′ (x) . g(y) = −f (x) . g ′′ (y)



 f (0) = 0

f (a) = 0


 g(0) = 0

f ′′ (x) g ′′ (y)
=−
f (x) g(y)
et omme f ne dépend que de x et g ne dépend que de y , e i ne peut être satisfait
f ′′ (x) g ′′ (y)
=− = C te = −λ2
f (x) g(y)
Nous obtenons alors deux équations diérentielles ordinaires (EDO) :

f ′′ + λ2 f = 0
g ′′ − λ2 g = 0
La solution de la première équation est donnée par les deux formes :
f (x) = A ei λ x + B e−i λ x = C sin(λ x) + D cos(λ x)

f (0) = 0 =⇒ D = 0
f (a) = 0 =⇒ C sin(λ a) = 0 =⇒ sin(λ a) = 0
- 13 -
nπ
=⇒ λ a = n π =⇒ λ = a
, n = 1, 2, · · · ∞
nπ
fn (x) = Cn sin( x)
a
La solution de la deuxième équation est donnée par les deux formes :
g(y) = A eλ y + B e−λ y = C sh(λ y) + D ch(λ y)

g(0) = 0 =⇒ D = 0 et par suite :
nπ
gn (y) = Cn sh( y)
a
Tn (x, y) = fn (x) . gn (y) = Cn sin( naπ x) . sh( naπ y)
∞ ∞
X X nπ nπ
T (x, y) = Tn (x, y) = Cn sin( x) . sh( y)
n=1 n=1
a a
Maintenant déterminons la onstante en utilisant la dernière ondition aux limites :
∞
Cn sin( naπ x) . sh( naπ b) = T0
P
T (x, b) =
n=1
Ra Ra 2 T0
Cn . sh( naπ b) = 2
a
T (x, b) sin( naπ x) dx = 2
a
T0 sin( naπ x) dx = nπ
[1 − (−1)n ]
0 0
La solution est triviale pour tout n paire, on ne prendra don que les valeur impaires,
pour ela on pose :
n = 2k − 1, k = 1, 2, · · · ∞. Dans e as :
4 T0 1
Ck = . h i
(2k − 1) π sh (2k−1) π b
a
et par suite, la solution du problème est :
h i h i
(2k−1) π
4 T0
∞
X sin a
x . sh (2k−1) a
π
y
T (x, y) = h i (2.7)
π k=1 (2k − 1) sh (2k−1) π
b
a
La distribution de la température en fon tion de x et de y est représentée sur la
- 14 -
gure (Fig.2.5).
Fig. 2.5: Distribution de la température en fon tion de x et de y.
2.3.3 Equation d'onde

Soit une orde maintenue par es deux extrémités et dénie par le problème suivant :
∂2u 2
2 ∂ u
=c 0≤x≤L , t≥0 (2.8)
∂t2 ∂x2

u(0, t) = 0 , u(x, 0) = 2 sin( Lπ x)
∂u
u(L, t) = 0 , ∂t
(x, o) = −sin( 2Lπ x)
Appliquons maintenant la MSV, on pose alors : u(x, t) = f (x) . g(t)

f (x) . g ′′(t) = c2 f ′′ (x) . g(t)

 f (0) = 0
 f (L) = 0
f ′′ (x) 1 g ′′ (t)
= 2
f (x) c g(t)
et omme f ne dépend que de x et g ne dépend que de t, e i ne peut être satisfait
- 15 -
f ′′ (x) 1 g ′′ (t)
= 2 = C te = −λ2
f (x) c g(t)
Nous obtenons alors deux équations diérentielles ordinaires (EDO) :

f ′′ + λ2 f = 0
g ′′ + c2 λ2 g = 0
La solution de la première équation est données par :
f (x) = A sin(λ x) + B cos(λ x)

f (0) = 0 =⇒ B = 0
f (L) = 0 =⇒ A sin(λ L) = 0 =⇒ sin(λ L) = 0
nπ
=⇒ λ L = n π =⇒ λ = L
, n = 1, 2, · · · ∞.
nπ
fn (x) = An sin( x)
L
La solution de la deuxième équation est donnée par :
g(t) = C sin(c λ t) + D cos(c λ t)
un (x, t) = fn (x) . gn (t) = an sin( nLπ x) . sin( c nL π t) + bn sin( nLπ x) . cos( c nL π t)
∞ ∞ h
X X nπ cnπ nπ cnπ i
u(x, t) = un (x, t) = an sin( x) . sin( t) + bn sin( x) . cos( t)
n=1 n=1
L L L L
Maintenant, déterminons les onstantes en utilisant les ondition initiales :

P∞
u(x, 0) = n=1 bn sin( nLπ x) = 2 sin( Lπ x)
P∞
=⇒ 2 sin( Lπ x) = b1 sin( Lπ x) + n=2 bn sin( nLπ x)
=⇒ b1 = 2 et bn = 0 pour n > 1.
- 16 -
2.4. Exer i es.
∂u
P∞ cnπ

∂t
(x, o) = n=1 an L
sin( nLπ x) = −sin( 2Lπ x)
P∞
=⇒ −sin( 2Lπ x) = a1 cπ
L
sin( Lπ x) + a2 2cπ
L
sin( 2Lπ x) + n=3 an cnπ
L
sin( nLπ x)
=⇒ a2 = − 2 Lπ c et an = 0 pour n 6= 2
et par suite, la solution nale du problème est :
π πc L 2π 2πc
u(x, t) = 2 sin( x) . cos( t) − sin( x) . sin( t) (2.9)
L L 2πc L L
La représentation graphique du dépla ement en fon tion du temps et de l'espa e
est représenté sur la gure (Fig.2.6).
Fig. 2.6: Variation du dépla ement en fon tion du temps et de l'espa e.
2.4 Exer i es
2-01 : En utilisant la MSV, déterminer la solution du problème suivant :
∂T ∂2T
∂t
− ∂x2
=0 0≤x≤1, t>0
T (0, t) = 0 , T (1, t) = 0 , T (x, 0) = 100
- 17 -
∂2T ∂2T
∂x2
+ ∂y 2
=0 0≤x≤1, 0≤y≤2
T (0, y) = 0 , T (1, y) = 0 , T (x, 0) = 0 , T (x, 2) = 1/2
∂2T ∂2T
∂x2
+ ∂y 2
=0 0≤x≤2, 0≤y≤1
T (0, y) = 0 , T (2, y) = 0 , T (x, 0) = 0 , T (x, 1) = 2
2-04 : En utilisant la MSV, déterminer la solution du problème suivant pour les deux
as de f (x) :
∂2u 2
2 ∂ u
= c 0≤x≤L , t≥0
∂t2 ∂x2

u(0, t) = 0 , u(x, 0) = 3 sin( Lπ x)
u(L, t) = 0 ∂u
, ∂t
(x, o) = f (x)
1. f (x) = 0.
2. f (x) = 3 sin( Lπ x) − sin( 2Lπ x)
- 18 -
Chapitre 3
Méthodes des diéren es nies
3.1 Introdu tion sur les méthodes numériques

En général, les méthodes numériques donnent une solution approximative ave une
pré ision susante pour les problèmes ren ontrés en pratique dans les domaines de
transferts thermiques, de mé anique des uides , d'éle tromagnétisme, . . .et . Parmi
es méthodes, on distingue d'une manière générale :
Méthode des diéren es nies (MDF).
Méthode des éléments nis (MEF).
Méthode des volumes nis (MVF).
Méthodes spe trales (MS).
Il existe d'autres méthodes spé iques à haque domaine (méthode des lignes, méthode
de Monté Carlo, . . .et ).
Les méthodes MDF et MVF sont basées sur une dis rétisation du domaine et le
rempla ement de l'opérateur diérentielle par un quotient diérentiel.
La méthode MEF est basée sur la formulation variationnelle.
La méthode MS est basée sur la re her he d'une solution appro hée sous forme
d'un développement sur une ertaine famille de fon tions (Fourier, polynmes,
splines, . . .et ). Ces méthodes sont souvent oûteuses mais pré ises.
3.2 Prin ipe et bases de la MDF

L'obje tif de la méthode des diéren es nies (MDF) est de transformer une équa-
tion ontinue valable sur un domaine ontinu (ou domaine physique ) en un système
à N équations à N in onnues asso iées à un domaine dis ret (ou domaine de al ul )
appelé maillage.(Fig. 3.1).
L'approximation des dérivées d'une fon tion f (x) par diéren es nies est basée sur
le développement en séries de Taylor de ette fon tion de la manière suivante :
19
Ch3.MÉTHODES DES DIFFÉRENCES FINIES.
P
Y ij
(C.L)
∆y
1111111
0000000
0000000
1111111
0000000
1111111
Domaine
i ∆x X
Fig. 3.1: Passage du domaine ontinu au domaine dis ret.
h ′ h2 ′′ h3 ′′′ h4 (4)
f (x + h) = f (x) + f (x) + f (x) + f (x) + f (x) + · · · (3.1)
1! 2! 3! 4!
h ′ h2 ′′ h3 ′′′ h4 (4)
f (x − h) = f (x) − f (x) + f (x) − f (x) + f (x) − · · · (3.2)
1! 2! 3! 4!
2h ′ (2h)2 ′′ (2h)3 ′′′ (2h)4 (4)

f (x + 2h) = f (x) + f (x) + f (x) + f (x) + f (x) + · · · (3.3)
1! 2! 3! 4!
2h ′ (2h)2 ′′ (2h)3 ′′′ (2h)4 (4)

f (x − 2h) = f (x) − f (x) + f (x) − f (x) + f (x) − · · · (3.4)
1! 2! 3! 4!
h étant le pas spatial (∆x, ∆y ou ∆z ) ou temporel ∆t.
Nous pouvons aussi faire sortir les sommes et les diéren es des séries :
h4 (4)
f (x + h) + f (x − h) = 2 f (x) + h2 f ′′ (x) + f (x) + · · · (3.5)
12
h3 ′′′
f (x + h) − f (x − h) = 2 h f ′(x) + f (x) + · · · (3.6)
3
4h4 (4)
f (x + 2h) + f (x − 2h) = 2 f (x) + 4h2 f ′′ (x) + f (x) + · · · (3.7)
3
8h3 ′′′
f (x + 2h) − f (x − 2h) = 4h f ′(x) + f (x) + · · · (3.8)
3
Notons que les sommes ontiennent seulement les dérivées paires par ontre, les
- 20 -
3.2. Prin ipe et bases de la MDF.
diéren es ontiennent seulement les dérivées impaires. Le nombre d'équations impli-
quées et le nombre de termes à onserver dans haque équation dépend de la dérivée
et du degré de pré ision désiré.
L'équation (3.6) nous donne :
f (x+h)−f (x−h) h2
f ′ (x) = 2h
− 6
f ′′′ (x) + · · ·
Ou en ore :
f (x + h) − f (x − h)
f ′ (x) ≃ + O(h2 ) (3.9)
2h
qui est une approximation entrée du se ond ordre de la première dérivée de f (x).
2 2
O(h ) est appelée Erreur de tron ature en h .
Similairement, l'équation (3.5) nous donne :
f (x+h)−2 f (x)+f (x−h) h2

f ′′ (x) = h2
− 12
f (4) (x) + · · ·
ou en ore :
f (x + h) − 2 f (x) + f (x − h)
f ′′ (x) ≃ + O(h2 ) (3.10)
h2
qui est une approximation entrée du se ond ordre de la dérivée se onde de f (x).
Nous pouvons aussi obtenir de la même manière et en utilisant les équations (3.1)-
(3.8) les dérivées f ′′′ (x) et f (4) (x) omme il est indiqué dans les tableaux (Tab.3.1-3.5).
f (x − 2h) f (x − h) f (x) f (x + h) f (x + 2h)

′
2h f ( x) -1 0 +1
h2 f ′′ (x) +1 -2 +1
2h3 f ′′′ (x) -1 2 0 -2 +1
h4 f (4) (x) +1 -4 +6 -4 +1
Tab. 3.1: Approximation entrée en O(h2 ).
L'équation (3.1) nous donne :
f (x+h)−f (x)
f ′ (x) = h
− h2 f ′′ (x) + · · ·
Ou en ore :
f (x + h) − f (x)
f ′ (x) ≃ + O(h) (3.11)
h
- 21 -
er
qui est une approximation dé entrée avant du 1 ordre de la première dérivée de f (x).
er
De même, l'équation (3.2) nous donne l'approximation dé entrée arrière du 1
ordre de la première dérivée de f (x) :
f (x) − f (x − h)
f ′ (x) ≃ + O(h) (3.12)
h
f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h)

h f ′ (x) -1 +1
h2 f ′′ (x) +1 -2 +1
h3 f ′′′ (x) -1 +3 -3 +1
h4 f (4) (x) +1 -4 +6 -4 +1
Tab. 3.2: Approximation dé entrée avant du 1er ordre O(h).
f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

′
h f (x) -1 +1
h2 f ′′ (x) +1 -2 +1
h3 f ′′′ (x) -1 +3 -3 +1
h4 f (4) (x) +1 -4 +6 -4 +1
Tab. 3.3: Approximation dé entrée arrière du 1er ordre O(h) .
Ces approximations pouvant être utilisées par exemple aux limites du domaine où
nd
l'approximation entré ne peut être appliquée. Les approximations dé entrées du 2
ordre sont données dans les tableaux (Tab.3.4 & 3.5) :
f (x) f (x + h) f (x + 2h) f (x + 3h) f (x + 4h) f (x + 5h)

′
2h f (x) -3 +4 -1
h2 f ′′ (x) +2 -5 +4 -1
2h3 f ′′′ (x) -5 +18 -24 +14 -3
h4 f (4) (x) +3 -14 +26 -24 +11 -2
Tab. 3.4: Approximation dé entrée avant du 2nd ordre O(h2 ).
f (x − 5h) f (x − 4h) f (x − 3h) f (x − 2h) f (x − h) f (x)

′
2h f (x) +1 -4 +3
h2 f ′′ (x) -1 +4 -5 +2
2h3 f ′′′ (x) +3 -14 +24 -18 +5
h4 f (4) (x) -2 +11 -24 +26 -14 +3
Tab. 3.5: Approximation dé entrée arrière du 2nd ordre O(h2 ).
- 22 -
3.2. Prin ipe et bases de la MDF.
Re ommandations :
Utiliser une double pré ision dans les al uls.
nd
Utiliser le 2 ordre O(h2 ) ou plus.
Exemples :
1- La ourbe (Fig.3.2) nous montre l'eet du nombre de termes à onsidérer durant
l'approximation de la fon tion ex par les séries de Taylor.
450
Exp(x)
7 termes
400 6 termes
5 termes
4 termes
350
300
250
200
150
100
50
0
3 4 5 6
x
Fig. 3.2: Approximation par séries de Taylor de la fon tion ex.
2- Soit à al uler la dérivée se onde de la fon tion f (x) = e−x au point x = 1.
Ave 6 dé imales : f ′′ (1) = f (1) = 0.367879.

L'approximation de ette dérivée est donnée par la relation (3.10).
h f (x + h) f (x − h) f ′′ (x)
0.64 0.193980 0.697676 0.380611
0.32 0.267135 0.506617 0.371038
0.16 0.313486 0.431711 0.368699
0.08 0.339596 0.398519 0.368214
0.04 0.353455 0.382893 0.368480
0.02 0.360595 0.375311 0.370098
0.01 0.364219 0.371577 0.376706
0.005 0.366045 0.369723 0.403174
0.0025 0.366961 0.368800 0.509054
0.00125 0.367420 0.368340 0.932579
0.001 0.367512 0.368248 1.250222
Tab. 3.6: Inuen e négative de la simple pré ision.
- 23 -
Ave 9 dé imales : f ′′ (1) = f (1) = 0.367879441.
h f (x + h) f (x − h) f ′′ (x)
0.64 0.193980042 0.697676326 0.380609097
0.32 0.267135302 0.506616992 0.371029414
0.16 0.313486181 0.431710523 0.368664921
0.08 0.339595526 0.398519041 0.368075685
0.04 0.353454682 0.382892886 0.367928494
0.02 0.360594940 0.375311099 0.367891704
0.01 0.364218980 0.371576691 0.367882507
0.005 0.366044635 0.369723445 0.367880208
0.0025 0.366960891 0.368800290 0.367879633
0.00125 0.367419879 0.368339578 0.367879489
0.001 0.367511746 0.368247505 0.367879472
Tab. 3.7: Inuen e positive de la double pré ision.
Nous remarquons, d'après es deux tests, l'eet néfaste du hoix de la simple pré i-
sion par rapport à la double pré ision qui onverge bien vers la solution exa te quand
le pas h diminue.
3- Soit à al uler la dérivée de la fon tion f (x) = Ln(x) au point x = 1. Prenons
h = 0.1 et omparons les approximations entrée (3.9), dé entrée avant (3.11) et
dé entrée arrière (3.12).
Approximation f ′ (1) erreur relative (%)
Centrée 1.003353477 0.34

Dé entrée avant 0.953101798 4.69
Dé entrée arrière 1.053605157 5.36
Tab. 3.8: Inuen e de l'approximation.
Nous remarquons que l'approximation entrée est plus pré ise que elles dé entrées
et que les erreurs al ulées sont pro hes de elles prévues par la théorie :

Pour les approximations dé entrées : E = h2 f ′′ (1) = 0.1
2
.(−1) = 0.05.

h2 ′′′ (0.1)2
Pour l'approximation entrée : E = 6 f (1) = 6
.2 = 0.00333.
- 24 -
3.3. Formule de Taylor en oordonnées re tangulaires.
3.3 Formule de Taylor en oordonnées re tangulaires

Etant donnée que le prin ipe de la dis rétisation est basée sur la formule de Taylor, il
est alors utile de la rappeler dans le as d'une fon tion f (x, y) au point (x+∆x, y+∆y) :
2
∂ ∂ 1 ∂ ∂
f (x+∆x, y+∆y) = f (x, y)+ ∆x + ∆y f (x, y)+ ∆x + ∆y f (x, y)+
∂x ∂y 2! ∂x ∂y
n−1 n
1 ∂ ∂ 1 ∂ ∂
+···+ ∆x + ∆y f (x, y) + ∆x + ∆y f (x, y)
(n − 1)! ∂x ∂y n! ∂x ∂y
(3.13)
3.4 S hémas de dis rétisation

Avant de rempla er les dérivées partielles par les diéren es nies, il faut d'abord
hoisir l'un des s hémas de dis rétisation suivants (Fig.3.3) :
Schéma décentré arrière

i−2 i−1 i (ou retardé)
Schéma centré
i−1 i i+1
Schéma décentré avant

i i+1 i+2 (ou avancé)
Fig. 3.3: Diérents types de s hémas de dis rétisation.
3.5 Notions de stabilité, onsistan e et onvergen e

3.5.1 Stabilité
La notion de stabilité est relative à la solution de l'équation aux diéren es, on dira
que le s héma est stable si, quelque soient ∆x et ∆y (du même ordre), la solution du
problème dis rétisé reste bornée quand ∆t tends vers 0.
Pour étudier la stabilité, on utilise souvent l'analyse de Von-Neumann qui onsiste
à étudier toute solution de la forme : ψ(t) ej β x (j 2 = −1).

Si on substitue ette forme dans l'équation aux diéren es on obtient une relation
entre ψ(t) et ψ(t + ∆t) .
- 25 -
Théorème : Le s héma est stable ( ondition susante) si le fa teur d'ampli ation
ξ vérie :

ξ = ψ(t + ∆t) ≤ 1

(3.14)
ψ(t)
Remarque : L'analyse de la stabilité par Von-Neumann est valable uniquement
pour les problèmes linéaires à oe ient onstants ar e sont les seuls problèmes
pour lesquels nous pouvons avoir une solution analytique. D'autres méthodes
existent pour l'analyse de la stabilité (matri ielle, perturbation dis rète ).
3.5.2 Consistan e
Le terme de onsistan e appliqué à une ertaine pro édure de onstru tion de s hé-
mas aux diéren es signie que ette pro édure appro he la solution de l'EDP onsi-
dérée et non la solution d'une autre EDP. La qualité de ette onsistan e s'appelle la
pré ision du s héma.
Théorème : Soient les opérateurs D dis ret et L ontinu et soit ε un paramètre
aussi petit que l'on veut.
Si |Df − Lf | ≤ ε quand ∆x et ∆t −→ 0 alors la méthode est onsistante.
3.5.3 Convergen e
Le s héma est onvergent si la solution du problème dis rétisé tend vers la solution
du problème aux dérivées partielles, lorsque ∆x et ∆y tendent simultanément vers 0.
Théorème de Lax : Pour une EDP linéaire, si le ritère de onsistan e est vérié,
une ondition né essaire et susante de onvergen e et que le s héma soit stable.
3.6 Méthodes de dis rétisation

3.6.1 Méthode expli ite
L'équation expli ite est obtenue en é rivant le terme spatial à l'instant n où la
solution est onnue. L'in onvénient prin ipal de ette méthode (Fig.3.4) est qu'elle né-
essite de hoisir un pas de temps ∆t susamment petit, sinon la solution de l'équation

dis rétisée devient instable.
- 26 -
3.6. Méthodes de dis rétisation.
n+1
n Solution connue
jusqu’en n
i−1 i i+1
Fig. 3.4: S héma expli ite.
Exemple : La dis rétisation de l'équation de la haleur nous donne :
∂T ∂2T Tin+1 − Tin n

Ti−1 − 2 Tin + Ti+1
n
= −→ = (3.15)
∂t ∂x2 ∆t ∆x2
∆t
∈ 0, 12

Cette méthode onverge pour λ= ∆x2
ave une erreur de dis rétisation
en O(∆t + ∆x2 ).
1
Pour λ= 6
, on démontre que l'erreur est en O(∆t2 + ∆x2 ) .
A ause de e risque d'instabilité, il est préférable, lorsque ela est possible,
d'utiliser l'une des méthodes impli ites.
Il existe aussi d'autres méthodes de dis rétisations expli ites telles que elles de Saul'Yev,
Dufort-Frankel et Barakat et Clark ...et .
3.6.2 Méthode impli ite

On obtient une équation impli ite en é rivant le terme spatial à l'instant (n + 1) où
la solution n'est pas onnue (Fig.3.5).
n+1
n Solution connue
jusqu’en n
i−1 i i+1
Fig. 3.5: S héma impli ite.
Exemple :
La dis rétisation de l'équation de la haleur nous donne :
∂T ∂2T Tin+1 − Tin T n+1 − 2 Tin+1 + Ti+1

n+1
= −→ = i−1 (3.16)
∂t ∂x2 ∆t ∆x2
- 27 -
Un avantage essentiel de ette méthode est qu'elle est universellement stable. La
seule limitation sur ∆t est elle qui maintient les erreurs de tron ature dans des limites
a eptables.
∆t
Cette méthode onverge quelque soit λ= ∆x2
= C te .
Il existe aussi la méthode impli ite de Leasonen en O(∆t + ∆x2 ) et qui est in on-
ditionnellement stable.
3.6.3 Méthode de Crank-Ni holson

Elle onsiste à é rire le se ond membre de l'équation dis rétisée omme la demi-
somme des se onds membres des méthodes impli ite et expli ite.
- Exemple : La dis rétisation de l'équation de la haleur nous donne :

n+1
Tin+1 − Tin n
− 2 Tin + Ti+1
n
1 Ti−1 − 2 Tin+1 + Ti+1
n+1

1 Ti−1
= + (3.17)
∆t 2 ∆x2 2 ∆x2
Le premier ro het de l'équation étant onnu, on voit que e s héma est en fait
un s héma impli ite. L'avantage prin ipal de ette équation est que pour une valeur
donnée de ∆x, l'erreur de tron ature sur le terme en ∆t est nettement plus petite que
dans les méthodes impli ite et expli ite. Elle exprime en eet la dérivée par rapport
∆t
au temps t à l'aide de diéren es entrées de pas
2
au lieu de la faire à l'aide de
diéren es à droite ou à gau he de pas ∆t (Fig.3.6).
n+1
n+1/2
n
i−1 i i+1
Fig. 3.6: S héma de Crank-Ni holson.
Cette méthode onverge ave une erreur de dis rétisation en O(∆t + ∆x2 ).
Pour λ= √1 , on démontre que l'erreur est en O(∆x6 ).
20
Ce s héma est in onditionnellement stable.
- 28 -
3.7. Appli ation de la méthode expli ite à l'équation de diusion.
Une forme plus générale de l'équation (3.17) s'é rit :
n+1
Tin+1 − Tin Ti−1 − 2 Tin+1 + Ti+1
n+1
n
Ti−1 − 2 Tin + Ti+1
n

= (1 − θ) +θ (3.18)
∆t ∆x2 ∆x2
ave 0 ≤ θ ≤ 1. Cette méthode est omplètement expli ite pour θ=0 et partiellement
impli ite pour 0 < θ < 1. Pour θ = 1, elle devient omplètement impli ite et pour θ=
1/2, on obtient la méthode de Crank-Ni holson qui présente en générale des solutions
raisonnablement stables et pré ises.
3.7 Appli ation de la méthode expli ite à l'équation

de diusion
L'équation de la haleur dis rétisée ave une s héma expli ite (3.15) peut se mettre
sous la forme :
Tin+1 = λ Ti−1
n
+ (1 − 2λ) Tin + λ Ti+1
n
(3.19)
Etudions maintenant la onsistan e et la stabilité de e s héma :
3.7.1 Consistan e
Les approximations des dérivées première et se onde sont :
∂T Tin+1 − Tin
≃ + O(∆t)
∂t ∆t
n
∂2T Ti−1 − 2 Tin + Ti+1
n
≃ + O(∆x2 )
∂x2 ∆x2
∂T ∂2T
f= ∂t
− ∂x2
Appliquons les opérateurs dis ret D et ontinu L à l'équation f :
h n+1 n n −2 T n +T n
i
T −T Ti−1 ∂2T
|Df − Lf | = i ∆t i + O(∆t) − i i+1
+ O(∆x2 ) − ∂T −

∆x2 ∂t ∂x2
|Df − Lf | = |O(∆t) + O(∆x2 )| = O(∆t + ∆x2 )
- 29 -
Si ∆x−→ 0 et ∆t −→ 0 alors O(∆t + ∆x2 ) −→ 0 =⇒ le s héma est onsistant
ave l'EDP.
3.7.2 Stabilité
Pour étudier la stabilité, utilisons l'analyse de Von-Neumann :
On suppose que la solution est de la forme : T (x, t) = ψ(t) ej β x , on a alors :
Tin = ψ(t) ej β x n
Ti−1 = ψ(t) ej β (x−∆x)
Tin+1 = ψ(t + ∆t) ej β x n

Ti+1 = ψ(t) ej β (x+∆x)
En remplaçant es termes dans l'équation (3.19) et en simpliant par ej β x , nous
aurons :
ψ(t + ∆t) = λ ψ(t) e−j β ∆x + (1 − 2 λ) ψ(t) + λ ψ(t) ej β ∆x
ψ(t+∆t)
= λ ej β ∆x + e−j β ∆x + (1 − 2 λ) = 2λ cos(β ∆x) + (1 − 2λ)

ψ(t)
D'où le fa teur d'ampli ation :
ξ = 1 − 2 λ (1 − cos(β ∆x)) (3.20)
|ξ| ≤ 1 =⇒ −1 ≤ ξ ≤ +1
 


1 − 2 λ (1 − cos(β ∆x)) ≥ −1 

 1 − 4 λ sin2( β ∆x
2
) ≥ −1
 
⇐⇒

 

1 − 2 λ (1 − cos(β ∆x)) ≤ 1
 1 − 4 λ sin2( β ∆x ) ≤ 1

2
∆t
Etant donné que λ= ∆x2
>0 alors la deuxième équation est toujours vériée ∀ λ.
La première équation nous donne la ondition :
∆t 1
λ= 2
≤
∆x 2
qui est la ondition de stabilité de e s héma. On dira alors qu'il est onditionnellement
stable.
C'est à e niveau que e pose la di ulté ar le hoix du pas de temps est stri tement
∆x2
lié au pas d'espa e et doit don obligatoirement satisfaire la ondition : ∆t ≤ 2
.
- 30 -
3.8. Appli ation de la méthode impli ite à l'équation de diusion .
3.8 Appli ation de la méthode impli ite à l'équation

de diusion
L'équation de la haleur dis rétisée ave une s héma impli ite (3.16) peut se mettre
sous la forme :
Tin = −λ Ti−1
n+1
+ (1 + 2λ) Tin+1 − λ Ti+1
n+1
(3.21)
Etudions maintenant la onsistan e et la stabilité de e s héma :
3.8.1 Consistan e
Les approximations des dérivées première et se onde sont :
∂T T n+1 − Tin
≃ i + O(∆t)
∂t ∆t
n+1
∂2T Ti−1 − 2 Tin+1 + Ti+1
n+1
≃ + O(∆x2 )
∂x2 ∆x2
∂T ∂2T
f= ∂t
− ∂x2
Appliquons les opérateurs dis ret D et ontinu L à l'équation f :
n+1
−2 Tin+1 +Ti+1
n+1
h n+1 n i
T −T Ti−1 2 ∂T ∂2T
|Df − Lf | = i ∆t i + O(∆t) − + O(∆x ) − −

∆x 2 ∂t ∂x2
|Df − Lf | = |O(∆t) + O(∆x2 )| = O(∆t + ∆x2 )
Si ∆x−→ 0 et ∆t −→ 0 alors O(∆t + ∆x2 ) −→ 0 =⇒ le s héma est onsistant ave
l'EDP.
3.8.2 Stabilité
On suppose que la solution est de la forme : T (x, t) = ψ(t) ej β x , on a alors :
n+1
Tin = ψ(t) ej β x Ti−1 = ψ(t + ∆t) ej β (x−∆x)
Tin+1 = ψ(t + ∆t) ej β x n+1

Ti+1 = ψ(t + ∆t) ej β (x+∆x)
aurons :
- 31 -
ψ(t) = −λ ψ(t + ∆t) e−j β ∆x + (1 + 2 λ) ψ(t + ∆t) − λ ψ(t + ∆t) ej β ∆x
ψ(t) = −2λ cos(β ∆x) ψ(t + ∆t) + (1 + 2 λ) ψ(t + ∆t)
ψ(t)
ψ(t+∆t)
= 1 + 2 λ (1 − cos(β ∆x)) = 1 + 4 λ sin2( β ∆x
2
)
D'où le fa teur d'ampli ation :
1
ξ= (3.22)
1 + 4 λ sin2( β ∆x
2
)
∀ ∆x et ∀ ∆t , nous avons toujours ∀λ > 0 : |ξ| ≤ 1.
On dira alors que e s héma est in onditionnellement stable.
3.9 Formulation matri ielle du s héma expli ite ap-

pliqué à l'équation de diusion
Reprenons notre problème dis rétisé par la méthode expli ite (3.19) ave les (C.L)
dénies dans le problème ontinu (2.3) :
Tin+1 = λ Ti−1
n
+ (1 − 2λ) Tin + λ Ti+1
n
n = 0, 1, 2, · · · , nmax , i = 1, 2, · · · , imax
La dis rétisation du milieu est faite selon la gure (Fig.3.7) i-dessous : haque ligne
représente le problème unidimensionnel (1D) au pas temps orrespondant.
Dans e as, nous avons l'abs isse et l'ordonnée de haque n÷ud :
xi = (i − 1) ∆x , tn = n ∆t
- Conditions aux limites (C.L) :

T n = α = 0
1
n = 0, 1, 2, · · · , nmax
T n = β = 0
imax
- 32 -
3.9. Formulation matri ielle du s héma expli ite appliqué à l'équation de diusion.
n
t Ti
nmax
∆t
n+1
n
β
n−1
α
3
2
1
0
∆x i max x
1 2 3 i−1 i i+1
γ
Fig. 3.7: Dis rétisation du domaine pour l'équation de la haleur.
- Conditions initiales (C.I) :
Ti0 = γ = 1 , i = 2, 3, · · · , imax − 1
Pour tout point intérieur au domaine Tin , on utilise l'équation dis rétisée ave :
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
Dis rétisons maintenant l'équation aux pas de temps 0 et 1 avant d'aboutir à la
formulation matri ielle nale :
n=0 : Ti1 = λ Ti−1

0
+ (1 − 2λ) Ti0 + λ Ti+1
0
T21 = λ T10 + (1 − 2λ) T20 + λ T30
T31 = λ T20 + (1 − 2λ) T30 + λ T40

. . . .
.. .. .. ..
Ti1max −1 = λ Ti0max −2 + (1 − 2λ) Ti0max −1 + λ Ti0max
D'après ette première étape, les Ti1 sont tous onnues.
n=1 : Ti2 = λ Ti−1

1
+ (1 − 2λ) Ti1 + λ Ti+1
1
- 33 -
T22 = λ T11 + (1 − 2λ) T21 + λ T31
T32 = λ T21 + (1 − 2λ) T31 + λ T41

. . . .
. . . .
. . . .
Ti2max −1 = λ Ti1max −2 + (1 − 2λ) Ti1max −1 + λ Ti1max
A ette étape, les Ti2 sont tous onnues et ainsi de suite jusqu'à nmax .
A partir de es détails, il est maintenant fa ile de passer à la formulation matri ielle :
 n+1   n
T2 − λα (1 − 2λ) λ 0 ··· 0
.. .
.
   

 T3 


 λ (1 − 2λ) λ . . 



. .. .. ..

. =   Ti 
   
 .  0 . . . 0 
..
   
..
Timax −2

. λ (1 − 2λ) λ
   
   .  
Timax −1 − λβ 0 ··· 0 λ (1 − 2λ)
(3.23)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
Nous obtenons ainsi un système d'équations fa ile à résoudre pour haque pas de
temps ar haque équations ontient une seule in onnue. Les résultats sont inje tés dans
le système au pas de temps suivant et ainsi de suite. . . . Ce système de N équations à
N in onnues peut se mettre sous la forme :
A.T = B
B est le ve teur olonne orrespondant aux N températures in onnues, T est un
ve teur olonne orrespondant aux termes non homogènes des équations aux diéren es
(températures onnues). La matri e A = (aij ) des oe ients des températures dans les
équations aux diéren es est une matri e (N × N) présentant de nombreux zéros. C'est
une matri e bande et, parti ulièrement pour e problème, 'est une matri e tridiagonale.
- 34 -
3.10. Formulation matri ielle du s héma impli ite appliqué à l'équation de diusion.
3.10 Formulation matri ielle du s héma impli ite ap-

pliqué à l'équation de diusion
Reprenons notre problème dis rétisé par la méthode impli ite (3.21) :
Tin = −λ Ti−1
n+1
+ (1 + 2λ) Tin+1 − λ Ti+1
n+1
n = 0, 1, 2, · · · , nmax , i = 1, 2, · · · , imax
Dis rétisons maintenant l'équation aux pas de temps 0 et 1 :
n=0 : Ti0 = −λ Ti−1

1
+ (1 + 2λ) Ti1 − λ Ti+1
1
T20 = −λ T11 + (1 + 2λ) T21 - λ T31
T30 = −λ T21 + (1 + 2λ) T31 - λ T41

. . . .
.. .. .. ..
Ti0max −1 = −λ Ti1max −2 + (1 + 2λ) Ti1max −1 - λ Ti1max
On remarque alors que haque équation ontient plusieurs in onnues. En remplaçant
les valeurs onnues à partir des (C.L) et (C.I), nous aurons le système suivant dont la
solution sera utilisée à l'étape suivante :
    1
γ + λα (1 + 2λ) −λ 0 ··· 0
   .. .. 
γ −λ (1 + 2λ) −λ . .
    
    
. .. .. ..

.. = 0 0   Ti 
   
 . . . 
.
   
..


 γ  
  .. . −λ (1 + 2λ) −λ

 
γ + λβ 0 ··· 0 −λ (1 + 2λ)
i = 2, 3, · · · , imax − 1
n=1 : Ti1 = −λ Ti−1

2
+ (1 + 2λ) Ti2 − λ Ti+1
2
- 35 -
T21 = −λ T12 + (1 + 2λ) T22 - λ T32
T31 = −λ T22 + (1 + 2λ) T32 - λ T42

. . . .
. . . .
. . . .
Ti1max −1 = −λ Ti2max −2 + (1 + 2λ) Ti2max −1 - λ Ti2max
Nous pouvons maintenant passer à la formulation matri ielle :
 n   n+1
T2 + λα (1 + 2λ) −λ 0 ··· 0
.. .
.
   

 T3 


 −λ (1 + 2λ) −λ . . 



.. .. .. ..

= 0 0   Ti 
   
 .   . . . 
..
   
..


 Timax −2 


 . . −λ (1 + 2λ) −λ

 
Timax −1 + λβ 0 ··· 0 −λ (1 + 2λ)
(3.24)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
A haque étape, on résout alors un système de N équations à N in onnues. A part
la première et la dernière équations qui ontiennent ha une 2 in onnues, toutes les
autres équations en ontiennent 3. Le système d'équations obtenu est don dépendant
ontrairement à elui du as expli ite.
3.11 Exemple résolu par wxMaxima

Il s'agit de déterminer la distribution instantanée de la température à travers une
barre de longueur L maintenue à es extrémités par une température nulle et initiale-
ment à la température unité (problème déni par (2.3)).
Pour simplier, nous allons prendre les pas de temps et d'espa e omme suit :
5 2 ∆t
∆t = 1000
, ∆x = 10
=⇒ λ= ∆x2
= 0.125
Choisissons aussi omme nombre de divisions maximal spatial : imax = 11 et temporel :
nmax = 15.
- 36 -
3.11. Exemple résolu par wxMaxima.
Nous aurons alors la dis rétisation suivante (Fig.3.8) du domaine spatio-temporel :
t
15
14
9
13 T7
12
11
10
valeure connue
9
valeure inconnue
8
7
6
5
4
3 1
T10
2
1
0
1 2 3 4 5 6 7 8 9 10 11 x
Fig. 3.8: Dis rétisation du domaine spatio-temporel.
E rire un programme ave le logi iel de al ul symbolique wxMaxima pour résoudre
e problème et tra er les ourbes né essaires de la distribution de température à travers
ette barre.
Ce programme doit être stru turé de la manière suivante :
Initialisation.
Données.
Cal uls intermédiaires.
Conditions aux limites et initiale.
Bou le prin ipale.
Résultats numériques.
Tra é des ourbes.
Changer les valeurs de ∆t et ∆x pour voir leurs eets sur la stabilité de la solution.
- 37 -
3.12 Cas des onditions aux limites du 2ème type

Dans le as d'un ux imposé aux bornes du domaine, des onsidérations spé iales
doivent être prises puisque dans e as, la dis rétisation omprends les points situés aux
limites du domaine (Fig.3.9). Une manière simple de traiter e problème est d'utiliser
un s héma dé entré omme approximation aux bornes.
Points fictifs extérieurs au domaine
x=0 x=L
0 1 2 3 i−1 i i+1 i max −1 i max i max +1
Fig. 3.9: Conditions aux limites de type Neumann.
- Conditions aux limites (C.L) :
Dans le domaine ontinu :
∂T
∂T

∂x x=0
=α et/ou
∂x x=L
=β
e qui orrespond dans le domaine dis rétisé à :
∂T ∂T

∂x i=1
=α et/ou
∂x i=imax
=β
En onsidérant par exemple le s héma expli ite -à-d que les variables sont évaluées
au temps onnu n, nous aurons d'après les diérentes possibilités des tableaux (3.1. à
3.5) :
- Au premier ordre :
- A gau he, on hoisi un s héma dé entré avant :
T2n − T1n
α= ∆x
⇒ T1n = T2n − α ∆x .
- A droite on hoisi un s héma dé entré arrière :
Tinmax − Tinmax −1
β= ∆x
⇒ Tinmax = β ∆x + Tinmax −1 .
- Au se ond ordre :
- A gau he, on hoisi un s héma dé entré avant :
−T3n + 4 T2n − 3 T1n 1

α= 2 ∆x
⇒ T1n = 3
(4 T2n − T3n − 2 α ∆x) .
- 38 -
3.12. Cas des onditions aux limites du 2ème type.
ou un s héma entré en faisant intervenir des points tifs :
T2n − T0n
α= 2 ∆x
⇒ T0n = T2n − 2 α ∆x .
- A droite on hoisi un s héma dé entré arrière :

3 Tinmax − 4 Tinmax −1 + Tinmax −2 1
Tinmax = 4 Tinmax −1 − Tinmax −2 + 2 β ∆x

β= 2 ∆x
⇒ 3
.
ou un s héma entré en faisant intervenir des points tifs :
Tinmax +1 − Tinmax −1
β= 2 ∆x
⇒ Tinmax +1 = Tinmax −1 + 2 β ∆x .
3.12.1 Cas de l'approximation entrée

En appliquant par exemple le s héma entré su essivement aux s hémas expli ite
et impli ite de l'équation de la haleur, nous aurons :
- S héma expli ite :
Tin+1 = λ Ti−1
n
+ (1 − 2λ) Tin + λ Ti+1
n
Ce s héma nous indique que la matri e est tridiagonale et que pour les noeuds in-
ternes les oe ients de ette matri e sont : λ, (1 − 2 λ) et λ. Il reste à déterminer
eux de la première et la dernière ligne en fon tion des (C.L).
- A gau he : T1n+1 = λ T0n + (1 − 2λ) T1n + λ T2n

T1n+1 = λ (T2n − 2 α ∆x) + (1 − 2λ) T1n + λ T2n
⇒ T1n+1 + 2 λ α ∆x = (1 − 2λ) T1n + 2 λ T2n
- A droite : Tin+1
max
= λ Tinmax −1 + (1 − 2λ) Tinmax + λ Tinmax +1
Tin+1
max
= λ Tinmax −1 + (1 − 2λ) Tinmax + λ (Tinmax −1 + 2 β ∆x)
⇒ Tin+1
max
− 2 λ β ∆x = 2 λ Tinmax −1 + (1 − 2λ) Tinmax
D'où la forme matri ielle :
 n+1   n
T1 + 2λα∆x (1 − 2λ) 2λ 0 ··· 0
.. .
.
   

 T2 


 λ (1 − 2λ) λ . . 



.. .. .. ..

=   Ti 
   
 .  0 . . . 0 
..
   
..
Timax −1

. λ (1 − 2λ) λ
   
   .  
Timax − 2λβ∆x 0 ··· 0 2λ (1 − 2λ)
(3.25)
- 39 -
n = 0, 1, 2, · · · , nmax , i = 1, 2, 3, · · · , imax
On remarque que la dimension de la matri e a augmenté par rapport à elle déterminée
ave la (C.L) de Diri hlet (3.23).
- S héma impli ite :
Tin = −λ Ti−1
n+1
+ (1 + 2λ) Tin+1 − λ Ti+1
n+1
Ce s héma nous indique que la matri e est tridiagonale et que pour les noeuds in-
ternes les oe ients de ette matri e sont : −λ, (1 + 2 λ) et −λ. Il reste à déterminer
- A gau he : T1n = −λ T0n+1 + (1 + 2λ) T1n+1 − λ T2n+1
T1n = −λ (T2n+1 − 2 α ∆x) + (1 + 2λ) T1n+1 − λ T2n+1
⇒ T1n − 2 λ α ∆x = (1 + 2λ) T1n+1 − 2 λ T2n+1
- A droite : Tinmax = −λ Tin+1

max −1
+ (1 + 2λ) Tin+1
max
− λ Tin+1
max +1
Tinmax = −λ Tin+1
max −1
+ (1 + 2λ) Tin+1
max
− λ (Tin+1
max −1
+ 2 β ∆x)
⇒ Tinmax + 2 λ β ∆x = −2 λ Tin+1
max −1
+ (1 + 2λ) Tin+1
max
 n   n+1
T1 − 2λα∆x (1 + 2λ) −2λ 0 ··· 0
.. .
.
   

 T2 


 −λ (1 + 2λ) −λ . . 



. .. .. ..

.  =   Ti 
   
 . 0 . . . 0 
.
   
..
Timax −1 .

. −λ (1 + 2λ) −λ
   
   .  
Timax + 2λβ∆x 0 ··· 0 −2λ (1 + 2λ)
(3.26)
n = 0, 1, 2, · · · , nmax , i = 1, 2, 3, · · · , imax
- 40 -
On remarque, dans e as aussi, que la dimension de la matri e a augmenté par
rapport à elle déterminée ave la (C.L) de Diri hlet (3.24).
Puisque la matri e est tridiagonale, alors elle peut être fa ilement résolue ave la
fa torisation LU ( ou Cholesky ).
3.12.2 Cas de l'approximation dé entrée d'ordre 1

Considérons maintenant, toujours pour le as de l'équation de la haleur dis réti-
sée par les s hémas expli ite et impli ite d'Euler, le as de l'approximation dé entrée
d'ordre 1 utilisée pour dis rétiser les onditions aux limites de type Neumann (CLN) :
Tin+1 = λ Ti−1
n
+ (1 − 2λ) Tin + λ Ti+1
n
Ce s héma nous indique que la matri e est tridiagonale et que pour les n÷uds in-
on utilise dans e as, l'approximation dé entrée avant, nous aurons alors :
T2n+1 = λ (T2n − α ∆x) + (1 − 2λ) T2n + λ T3n
⇒ T2n+1 + λ α ∆x = (1 − λ) T2n + λ T3n
- A droite : Tin+1
max −1
= λ Tinmax −2 + (1 − 2λ) Tinmax −1 + λ Tinmax
on utilise dans e as, l'approximation dé entrée arrière, nous aurons alors :
Tin+1
max −1
= λ Tinmax −2 + (1 − 2λ) Tinmax −1 + λ (Tinmax −1 + β ∆x)
⇒ Tin+1
max −1
− λ β ∆x = λ Tinmax −2 + (1 − λ) Tinmax −1
- 41 -
 n+1   n
T2 + λα∆x (1 − λ) λ 0 ··· 0
.. .
..
   

 T3 


 λ (1 − 2λ) λ . 



. .. .. ..

. =   Ti 
   
 .  0 . . . 0 
.
   
..
Timax −2 ..

. λ (1 − 2λ) λ
    
   
Timax −1 − λβ∆x 0 ··· 0 λ (1 − λ)
(3.27)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
On remarque que la dimension de la matri e est inférieure à elle déterminée ave l'ap-
proximation entrée.
Tin = −λ Ti−1
n+1
+ (1 + 2λ) Tin+1 − λ Ti+1
n+1
T2n = −λ (T2n+1 − α ∆x) + (1 + 2λ) T2n+1 − λ T3n+1
⇒ T2n − λ α ∆x = (1 + λ) T2n+1 − λ T3n+1
- A droite : Tinmax −1 = −λ Tin+1

max −2
+ (1 + 2λ) Tin+1
max −1
− λ Tin+1
max
Tinmax −1 = −λ Tin+1
max −2
+ (1 + 2λ) Tin+1
max −1
− λ (Tin+1
max −1
+ β ∆x)
⇒ Tinmax −1 + λ β ∆x = −λ Tin+1
max −2
+ (1 + λ) Tin+1
max −1
- 42 -
 n   n+1
T2 − λα∆x (1 + λ) −λ 0 ··· 0
.. .
..
   

 T3 


 −λ (1 + 2λ) −λ . 



. .. .. ..

.  =   Ti 
   
 . 0 . . . 0 
.
   
..
Timax −2 ..

. −λ (1 + 2λ) −λ
    
   
Timax −1 + λβ∆x 0 ··· 0 −λ (1 + λ)
(3.28)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
On remarque, dans e as aussi, que la dimension de la matri e est inférieure à elle
déterminée ave l'approximation entrée.
3.12.3 Cas de l'approximation dé entrée d'ordre 2

Considérons maintenant, toujours pour le as de l'équation de la haleur dis réti-
sée par les s hémas expli ite et impli ite d'Euler, le as de l'approximation dé entrée
d'ordre 2 utilisée pour dis rétiser les onditions aux limites de type Neumann (CLN) :
Tin+1 = λ Ti−1
n
+ (1 − 2λ) Tin + λ Ti+1
n
T2n+1 = λ
3
(4 T2n − T3n − 2 α ∆x) + (1 − 2λ) T2n + λ T3n
⇒ T2n+1 + 32 λ α ∆x = (1 − 32 λ) T2n + 32 λ T3n
- A droite : Tin+1
max −1
= λ Tinmax −2 + (1 − 2λ) Tinmax −1 + λ Tinmax
- 43 -
Tin+1
max −1
= λ Tinmax −2 + (1 − 2λ) Tinmax −1 + λ
3
(4 Tinmax −1 − Tinmax −2 + 2β ∆x)
⇒ Tin+1
max −1
− 23 λ β ∆x = 32 λ Tinmax −2 + (1 − 23 λ) Tinmax −1
n+1 n
(1 − 23 λ) 2
  
T2 + 32 λα∆x 3
λ 0 ··· 0
   .. .. 

 T3 


 λ (1 − 2λ) λ . . 



.. .. .. ..

= 0 0   Ti 
   
 .  . . . 
.
   
..
Timax −2

   .. . λ (1 − 2λ) λ
 
   
2
Timax −1 − 3
λβ∆x 0 ··· 0 2
λ (1 − 32 λ)
3
(3.29)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
On remarque que la dimension de la matri e est inférieure à elle déterminée ave l'ap-
proximation entrée.
Tin = −λ Ti−1
n+1
+ (1 + 2λ) Tin+1 − λ Ti+1
n+1
T2n = − λ3 (4 T2n+1 − T3n+1 − 2 α ∆x) + (1 + 2λ) T2n+1 − λ T3n+1
⇒ T2n − 23 λ α ∆x = (1 + 32 λ) T2n+1 − 23 λ T3n+1
- A droite : Tinmax −1 = −λ Tin+1

max −2
+ (1 + 2λ) Tin+1
max −1
− λ Tin+1
max
- 44 -
3.13. Equations non linéaires.
Tinmax −1 = −λ Tin+1
max −2
+ (1 + 2λ) Tin+1
max −1
− λ
3
(4 Tin+1
max −1
− Tin+1
max −2
+ 2 β ∆x)
⇒ Tinmax −1 + 32 λ β ∆x = − 23 λ Tin+1
max −2
+ (1 + 32 λ) Tin+1
max −1
n n+1
(1 + 23 λ) − 23 λ
  
T2 − 23 λα∆x 0 ··· 0
.. .
..
   

 T3 


 −λ (1 + 2λ) −λ . 



. .. .. ..

..  =   Ti 
   
 0 . . . 0 
.
   
..
Timax −2 .

. −λ (1 + 2λ) −λ
   
   .  
2
Timax −1 + 3
λβ∆x 0 ··· 0 2
−3λ (1 + 32 λ)
(3.30)
n = 0, 1, 2, · · · , nmax , i = 2, 3, · · · , imax − 1
On remarque, dans e as aussi, que la dimension de la matri e est inférieure à elle
déterminée ave l'approximation entrée
3.13 Equations non linéaires

Considérons l'équation non linéaire suivante :
∂T ∂T ∂2T
+ u(T ) . = α(T ) . 2 (3.31)
∂t ∂x ∂x
où la vitesse de onve tion u(T ) et le oe ient de diusion (diusivité thermique)
α(T ) dépendent de la température T .
3.13.1 Dis rétisation par un s héma expli ite

Tin+1 − Tin T n − Ti−1
n
T n − 2 Tin + Ti+1
n
+ uni . i+1 = αin . i−1 (3.32)
∆t 2∆x ∆x2
Les oe ients non-linéaires sont simplement évalués au n÷ud de base (i, n) où Tin
est onnue, don uni et αin aussi. L'équation (3.32) peut être fa ilement résolue pour
Tin+1 ar les oe ients non-linéaires ne posent au une ompli ation numérique. C'est
le as typique de tous les s hémas expli ites.
- 45 -
3.13.2 Dis rétisation par un s héma impli ite
n+1 n+1 n+1 n+1 n+1

Tin+1 − Tin n+1 Ti+1 − Ti−1 n+1 Ti−1 − 2 Ti + Ti+1
+ ui . = αi . (3.33)
∆t 2∆x ∆x2
Dans e as, les oe ients non-linéaires posent de sérieux problèmes d'ordre nu-
mérique. un+1
i et αin+1 dépendent de Tin+1 qui est in onnu. L'équation (3.33) appliquée
à haque n÷ud du maillage donne un système non-linéaire d'équations ouplées. Ce
système peut être résolu en retardant (lagging) les oe ients non-linéaires ( -à-d on
pose : un+1
i = uni et αin+1 = αin ) et par itérations en utilisant la méthode de Newton ou
par linéarisation temporelle.
3.14 Problèmes multidimensionnels

Considérons l'équation de diusion linéaire bidimensionnelle :
∂2T ∂2T

∂T
=α + (3.34)
∂t ∂x2 ∂y 2
3.14.1 Dis rétisation par un s héma expli ite

La dis rétisation de l'équation (3.34) par un s héma expli ite donne :
n+1 n n n n n n n
Ti,j − Ti,j Ti−1,j − 2 Ti,j + Ti+1,j Ti,j−1 − 2 Ti,j + Ti,j+1

=α + (3.35)
∆t ∆x2 ∆y 2
On peut montrer que e s héma est onditionnellement stable ave la ondition :

1 1 1
α ∆t 2
+ ≤ (3.36)
∆x ∆y 2 2
α ∆t 1
Si ∆x = ∆y alors ette ondition de stabilité sera λ= ∆x2
≤ 4
qui est en ore plus
1
restri tive que le as unidimensionnel (λ ≤ ), e qui fait que ette méthode devient
2
moins pratique.
n+1
L'équation (3.35) peut être résolue dire tement pour Ti,j et au une omplexité
numérique supplémentaire ne sera rajouté à ause du terme supplémentaire. De même

∂2T
pour le as 3D où on rajoute en ore le terme dis rétisé orrespondant à sans au un
∂z 2
problème d'ordre numérique. Ce i est le as typique de tous les s hémas expli ites.
On montre que e s héma est onditionnellement stable ave la ondition :

1 1 1 1
α ∆t + + ≤ (3.37)
∆x2 ∆y 2 ∆z 2 2
- 46 -
3.15. Equation de Lapla e.
α ∆t 1
Si ∆x = ∆y = ∆z alors ette ondition de stabilité sera λ= ∆x2
≤ 6
qui est en ore
plus restri tive que le as bidimensionnel.
3.14.2 Dis rétisation par un s héma impli ite
La dis rétisation de l'équation (3.34) par un s héma impli ite donne :
" #
n+1 n n+1 n+1 n+1 n+1 n+1 n+1
Ti,j − Ti,j Ti−1,j − 2 Ti,j + Ti+1,j Ti,j−1 − 2 Ti,j + Ti,j+1
=α + (3.38)
∆t ∆x2 ∆y 2
Ce s héma est in onditionnellement stable. En appliquant l'équation (3.38) en
haque n÷ud du maillage, on obtient une matri e bande pentadiagonale (5) qui né es-
site beau oup d'eorts de al ul.
La méthode SOR peut être appliquée pour les problèmes 2D, mais ette appro he est
assez lourde pour les n÷uds 3D. Les méthodes ADI et AFI Approximation-Fa torization-
Impli ite peuvent être utilisées pour réduire la matri e bande en deux (ou 3 pour les
problèmes 3D) systèmes de matri es tridiagonales qui peuvent être résolus su essive-
ment par l'algorithme de Thomas.
3.15 Equation de Lapla e
Dans le domaine de l'énergétique, l'équation de Lapla e est souvent ren ontré en
mé anique des uides dans les é oulements potentiels (irrotationnels) et en transfert
thermique dans la ondu tion dans un solide par exemple.
Soit par exemple à résoudre un problème régis par l'équation de Lapla e et modé-
lisant le transfert thermique à travers une plaque re tangulaire. Cette équation s'é rit
sous la forme :
∂2T ∂2T
+ =0 (3.39)
∂x2 ∂y 2
Nous passons du domaine ontinu au domaine dis rêt en divisant le premier en un
ensemble de mailles aussi petites que l'on peut (Fig.3.10).
- 47 -
Y (C.L) y T i,j
j max
b
(C.L) ∆y
j+1
(C.L) j
j−1
4
3
2
x
1
∆x
0 a X 1 2 3 i−1 i i+1 i max
(C.L)
Domaine continu Domaine discrétisé
Fig. 3.10: Dis rétisation du domaine pour l'équation de Lapla e.
3.15.1 Dis rétisation par la formule à 5 points (Runge 1908)

L'équation (3.39) est dis rétisée par l' approximation entrée d'ordre 2 en espa e
2 2
O(∆x + ∆y ) de la manière suivante :
Ti−1,j − 2 Ti,j + Ti+1,j Ti,j−1 − 2 Ti,j + Ti,j+1

+ =0 (3.40)
∆x2 ∆y 2
∆x
On posant : β= ∆y
nous pouvons al uler les températures nodales par :
Ti−1,j + Ti+1,j + β 2 (Ti,j−1 + Ti,j+1)

Ti,j = (3.41)
2 (1 + β 2 )
et l'équation dis rétisée s'é rit :
−2 (1 + β 2 ) Ti,j + Ti−1,j + Ti+1,j + β 2 (Ti,j−1 + Ti,j+1 ) = 0 (3.42)
L'étude de la onsistan e de e s héma s'ee tue de la même manière que pour
l'équation de diusion. La onvergen e sera réalisée quand la diéren e Eij entre la
solution al ulée Tij et la solution exa te T ij tend vers 0 quand ∆x et ∆y tendent vers
0 ( -à-d en rendant la maillage de plus en plus n).
Pour simplier, on prend souvent ∆x = ∆y et don β = 1. L'équation (3.41) se
réduit à :
Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1
Ti,j = (3.43)
4
−4 Ti,j + Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1 = 0 (3.44)
- 48 -
La représentation s hématique de e s héma est (Fig.3.11) :
β2
(a) 1 (b)
1 −4 1 1 −2(1 + β 2 ) 1
PSfrag repla ements

1
β2
Fig. 3.11: S héma ellulaire de la formulation à 5 points : (a) β = 1, (b) β6= 1.
Exemple 1 :
Soit à déterminer les températures nodales dans la plaque arrée soumise aux ondi-
tions aux limites de Diri hlet indiquées sur la gure (Fig.3.12).
y 100
4
T2,4 T3,4 T4,4 50
3
75 T2,3 T3,3 T4,3
2
T2,2 T3,2 T4,2
1
1 2 3 4 5 x
0
Fig. 3.12: Equation de Lapla e. Exemple 1.
Prenons 4 divisions selon X, d'où un nombre de n÷uds : Nx = 5.

Prenons 4 divisions selon Y, d'où un nombre de n÷uds : Ny = 5.
Nombre d'équations du système à résoudre : N = (Nx − 2) × (Ny − 2) = 9.
Puisque ∆x = ∆y , alors β=1 et on applique alors la formule (3.43) à haque n÷ud
interne du maillage, nous aurons le système de 9 équations à 9 in onnues suivant :
- 49 -



−4 T2,2 + T1,2 + T3,2 + T2,1 + T2,3 =0


−4 T3,2 + T2,2 + T4,2 + T3,1 + T3,3 =0






−4 T4,2 + T3,2 + T5,2 + T4,1 + T4,3 =0

 

T2,1 = T3,1 = T4,1 =0

 
 
−4 T2,3 + T1,3 + T3,3 + T2,2 + T2,4 =0

 

 

 T = T

= T1,4 = 75
1,2 1,3
−4 T3,3 + T2,3 + T4,3 + T3,2 + T3,4 = 0 ave les (C.L) :
 T = T = T5,4 = 50
  5,2 5,3
−4 T4,3 + T3,3 + T5,3 + T4,2 + T4,4 =0

 

 

 T = T

= T4,5 = 100
 2,5 3,5
−4 T2,4 + T1,4 + T3,4 + T2,3 + T2,5 =0






−4 T3,4 + T2,4 + T4,4 + T3,3 + T3,5 =0






−4 T + T + T + T + T

=0
4,4 3,4 5,4 4,3 4,5
La forme matri ielle A.T = B de e système s'é rit de la manière suivante :
     
−4 1 0 1 0 0 0 0 0 T2,2 −75
 1 −4 1 0 1 0 0 0 0   T3,2   0 
     
     
 0 1 −4 0 0 1 0 0 0   T4,2   −50 
     
 1 0 0 −4 1 0 1 0 0   T2,3   −75 
     
     
 0 1 0 1 −4 1 0 1 0  .  T3,3  =  0 
     
 0 0 1 0 1 −4 0 0 1   T4,3   −50 
     
     
 0 0 0 1 0 0 −4 1 0   T2,4   −175 
     
 0 0 0 0 1 0 1 −4 1 T −100 
     
  3,4  
0 0 0 0 0 1 0 1 −4 T4,4 −150
La solution numérique est déterminée en utilisant le logi iel de al ul wxMaxima :
T2,2 = 42.85, T3,2 = 33.26, T4,2 = 33.93, T2,3 = 63.17, T3,3 = 56.25,
T4,3 = 52.45, T2,4 = 78.57, T3,4 = 76.11, T4,4 = 69.64.
Exemple 2 :
E rire la forme matri ielle A.T = B du système à résoudre pour le problème de la
gure (Fig.3.12) en utilisant le s héma donné par la relation (3.41) et en supposant les
(C.L) in onnues .
En suivant la même démar he que pour l'exemple 1 et en posant γ = −2 (1 + β 2 ),
- 50 -
nous aurons :
     
γ 1 0 β2 0 0 0 0 0 T2,2 −T1,2 − β 2 T2,1
2
 1 γ 1 0 β 0 0 0 0   T3,2 −β 2 T3,1
     
  
     
 0 1 γ 0 0 β 2 0 0 0   T4,2   −T5,2 − β 2 T4,1 
     
 2
 β 0 0 γ 1 0 β 2 0 0   T2,3 −T1,3
    
  
     
 0 β 2 0 1 γ 1 0 β 2 0  .  T3,3 = 0 
     
 0 0 β 2 0 1 γ 0 0 β 2   T4,3 −T5,3
     
  
     
 0 0 0 β 2 0 0 γ 1 0   T2,4   −T1,4 − β 2 T2,5 
     
 0 0 0 0 β 2 0 1 γ 1   T3,4 −β 2 T3,5
     
  
0 0 0 0 0 β2 0 1 γ T4,4 2
−T5,4 − β T4,5
La matri e A que nous avons obtenu est dite matri e bande symétrique.
Exemple 3 :
E rire la forme matri ielle A.T = B pour le problème de la gure (Fig.3.13) et en
utilisant le s héma donné par la relation (3.41).
y
T1,7 T2,7 T3,7 T4,7
T1,6 T2,6 T3,6 T4,6
T1,5 T2,5 T3,5 T4,5
T1,4 T2,4 T3,4 T4,4
T1,3 T2,3 T3,3 T4,3
T1,2 T2,2 T3,2 T4,2
x
T1,1 T2,1 T3,1 T4,1
Fig. 3.13: Equation de Lapla e. Exemple 3.
En suivant la même démar he que pour l'exemple 1 et en posant γ = −2 (1 + β 2 ),
- 51 -
nous aurons :
     
γ 1 β2 0 0 0 0 0 0 0 T2,2 −T1,2 − β 2 T2,1
β2 0 0   −T4,2 − β 2 T3,1
     

 1 γ 0 0 0 0 0  
  T3,2  



 β2 0 γ 1 β2 0 0 0 0 0  
  T2,3  
  −T1,3 

0 β2 γ 0 β2
     

 1 0 0 0 0  
  T3,3  
  −T4,3 

 0 0 β2 0 γ 1 β2 0 0 0   T2,4   −T1,4 
. =
     
0 β2 0 β2
 

 0 0 1 γ 0 0  
  T3,4  
  −T4,4 


 0 0 0 0 β2 0 γ 1 β2 0  
  T2,5  
  −T1,5 

0 β2 1 0 β2
     

 0 0 0 0 γ  
  T3,5  
  −T4,5 


 0 0 0 0 0 0 β2 0 γ 1  
  T2,6   −T − β 2 T
  1,6 2,7


0 0 0 0 0 0 0 β2 1 γ T3,6 2
−T4,6 − β T3,7
3.15.2 Conditions limites de Diri hlet variables

Dans e as, l'un des otés de la plaque (ou bien tous) est soumis à une tempé-
rature variable (CLDV) régie par une fon tion quel onque (polynomiale, sinusoïdale,
logarithmique, ...et ).
Exemple 4 :
Reprenons le même exemple 1 de la gure (Fig.3.12) mais en supposant que la
température du oté supérieur n'est pas onstante et qu'elle est donnée par la relation
T (x) = 2.104 x2 . Les dimensions de la plaque étant de (12 × 12) cm2, e qui nous donne
un pas axial ∆x = 3 cm. Dans e as les nouvelles (C.L) seront :



 T2,1 = T3,1 = T4,1 =0


T = T

= T1,4 = 75
1,2 1,3
(C.L) :


 T5,2 = T5,3 = T5,4 = 50


T = 18,

T3,5 = 72, T4,5 = 162
2,5
ave la (CLDV) al ulée omme suit :
Ti,5 = T (xi ) = 2.104 x2i où xi = (i − 1) ∆x ave (i = 2, 3, 4)

e qui nous donne :

2
T = 2.104 ∆x2 = 2.104 (3.10−2 ) = 18
 2,5


2 2
T3,5 = 2.104 (2 ∆x2 ) = 2.104 (6.10−2 ) = 72

T = 2.104 (3 ∆x2 )2 = 2.104 (9.10−2 )2

= 162

4,5
- 52 -
et la nouvelle forme matri ielle de e problème s'é rira :
     
−4 1 0 1 0 0 0 0 0 T2,2 −75
1 −4 1 0 1 0 0 0 0 T3,2 0
     
     
     

 0 1 −4 0 0 1 0 0 0  
  T4,2  
  −50 

1 0 0 −4 1 0 1 0 0 T2,3 −75
     
     
     
 0 1 0 1 −4 1 0 1 0 . T3,3 = 0 
     
0 0 1 0 1 −4 0 0 1 T4,3 −50
     
     
     

 0 0 0 1 0 0 −4 1 0  
  T2,4  
  −93 

0 0 0 0 1 0 1 −4 1 T3,4 −72
     
     
0 0 0 0 0 1 0 1 −4 T4,4 −212
3.15.3 Méthodes de résolution

Les méthodes utilisées pour la résolution de e genre de systèmes sont les méthodes
dire tes et indire tes (ou itératives ).
Méthodes dire tes : Elimination de Gauss, Algorithme de Thomas, Cholesky,
...et . Ces méthodes onsomment beau oup de temps, les erreurs d'arrondi aug-
mentent ave la taille du système et deviennent atastrophiques !. Elles sourent
beau oup du problème de sto kage des éléments nuls qui o upent inutilement
la mémoire, elle onviennent don pour les systèmes de faible taille.
Méthodes itératives : Ja obi, Gauss-Siedel, SOR, ADI, ...et . Dans es mé-
thodes, les erreurs d'arrondi sont orrigées à haque itération et elles utilisent
uniquement les éléments non nuls don pas de problème de mémoire et de e fait
onviennent bien aux systèmes larges. Les méthodes re ommandées sont dans
l'ordre : SOR et ADI en premier, Gauss-Siedel en se ond et enn Ja obi qui
onverge plus lentement.
Remarques :
La méthode d'élimination de Gauss né essite un nombre de multipli ation égal
1
à : Nm = 3
N 3 + N 2 − 31 N soit 321 multipli ations pour notre exemple où
N = 9. Si e hire est doublé (N

= 18) alors Nm = 2262.
Pour un maillage de (5×7) n÷uds, -à-d 15 in onnues, la matri e orrespondante
ontient 225 éléments dont seulement 59 diérents de zéro, soit 26 %.
Pour un maillage de (9 × 13) n÷uds, -à-d 77 in onnues, la matri e orrespon-

dante ontient 5929 éléments dont seulement 169 diérents de zéro, soit 2.8 %.
- 53 -
3.15.4 Conditions aux limites de Neumann

Considérons que la plaque présente par exemple un plan de symétrie verti al (Fig.3.14).
Don au un transfert de haleur ne s'ee tue à travers e plan et il n'y aura don au un
∂T

gradient de température -à-d :
∂x AA′
= 0. Nous aurons alors d'après ette gure :
A
5
j=1
i=1 2 3 4
A’
Fig. 3.14: Equation de Lapla e. (C.L) de Neumann.

∂T
=0 (j = 2, 3, 4) (3.45)
∂x 3,j
E rivons la formulation à 5 points pour tous les points appartenant à l'axe AA′
-à-d de oordonnées (3,j ), d'après (3.41) nous aurons :
2 (1 + β 2 ) T3,j = T2,j + T4,j + β 2 (T3,j−1 + T3,j+1) (j = 2, 3, 4) (3.46)
Remarquons que les T4,j sont en dehors du domaine de al ul. Pour les al uler,
on utilisera l'équation (3.45) en la dis rétisant par un s héma entré de la manière
suivante :
∂T T4,j − T2,j
α= = + O(∆x2 )
∂x 3,j
2 ∆x
d'où l'on tire :
T4,j = T2,j + 2 α ∆x (3.47)
En remplaçant (3.47) dans (3.46) nous aurons les 3 équations :
2 (1 + β 2 ) T3,j = 2 T2,j + β 2 (T3,j−1 + T3,j+1) + 2 α ∆x (j = 2, 3, 4) (3.48)
- 54 -
Dans notre as parti ulier α = 0, alors la relation (3.48) devient :
2 (1 + β 2 ) T3,j = 2 T2,j + β 2 (T3,j−1 + T3,j+1 ) (j = 2, 3, 4) (3.49)
Généralisons maintenant l'équation (3.48) dans le as où la ondition de Neumann
(CLN) est pres rite sur l'un des otés d'une plaque re tangulaire, les autres otés
étant soumis à des onditions de Diri hlet (CLD). Cette (CLN) sera dis rétisée par
des s hémas d'ordre 1 et 2. Cette ondition sera dis rétisée i-dessous par des s hémas
d'ordre 1 et 2.
a- Dis rétisation par un s héma entré
CLN à gau he : d'après (3.41) nous aurons :
2 (1 + β 2 ) T1,j = T0,j + T2,j + β 2 (T1,j−1 + T1,j+1 )
T2,j −T0,j

∂T
α= ∂x 1,j
= 2 ∆x
=⇒ T0,j = T2,j − 2 α ∆x
2 (1 + β 2 ) T1,j = 2 T2,j + β 2 (T1,j−1 + T1,j+1) − 2 α ∆x (j = 2, 3, · · · jmax − 1) (3.50)
CLN à droite : d'après (3.41) nous aurons :
2 (1 + β 2 ) Timax ,j = Timax −1,j + Timax +1,j + β 2 (Timax ,j−1 + Timax ,j+1)
∂T Timax +1,j −Timax −1,j

α= ∂x imax ,j
= 2 ∆x
=⇒ Timax +1,j = Timax −1,j + 2 α ∆x
2 (1+β 2) Timax ,j = 2 Timax −1,j +β 2 (Timax ,j−1 +Timax ,j+1)+2 α ∆x (j = 2, 3, · · · jmax −1)
(3.51)
CLN en bas : d'après (3.41) nous aurons :
2 (1 + β 2 ) Ti,1 = Ti−1,1 + Ti+1,1 + β 2 (Ti,0 + Ti,2 )

∂T Ti,2 −Ti,0
α= ∂y
= 2 ∆y
=⇒ Ti,0 = Ti,2 − 2 α ∆y
i,1
2 (1 + β 2 ) Ti,1 = Ti−1,1 + Ti+1,1 + 2 β 2 (Ti,2 − α ∆y) (i = 2, 3, · · · imax − 1) (3.52)
CLN en haut : d'après (3.41) nous aurons :
2 (1 + β 2 ) Ti,jmax = Ti−1,jmax + Ti+1,jmax + β 2 (Ti,jmax−1 + Ti,jmax+1 )
- 55 -

∂T Ti,jmax +1 −Ti,jmax −1
α= ∂y
= 2 ∆y
=⇒ Ti,jmax +1 = Ti,jmax −1 + 2 α ∆y
i,jmax
2 (1 + β 2) Ti,jmax = Ti−1,jmax + Ti+1,jmax + 2 β 2 (Ti,jmax −1 + α ∆y) (i = 2, 3, · · · imax − 1)

(3.53)
b- Dis rétisation par des s hémas dé entrés d'ordre 1
2 (1 + β 2 ) T2,j = T1,j + T3,j + β 2 (T2,j−1 + T2,j+1)
∂T T2,j −T1,j

α= ∂x 1,j
= ∆x
=⇒ T1,j = T2,j − α ∆x
(1 + 2 β 2) T2,j = T3,j + β 2 (T2,j−1 + T2,j+1 ) − α ∆x (j = 2, 3, · · · jmax − 1) (3.54)
2 (1 + β 2 ) Timax −1,j = Timax −2,j + Timax ,j + β 2 (Timax −1,j−1 + Timax −1,j+1)
Timax ,j −Timax −1,j

∂T
α= ∂x imax ,j
= ∆x
=⇒ Timax ,j = Timax −1,j + α ∆x
(1+2 β 2) Timax −1,j = Timax −2,j +β 2 (Timax −1,j−1+Timax −1,j+1)+α ∆x (j = 2, 3, · · · jmax −1)
(3.55)
2 (1 + β 2 ) Ti,2 = Ti−1,2 + Ti+1,2 + β 2 (Ti,1 + Ti,3 )

∂T Ti,2 −Ti,1
α= ∂y
= ∆y
=⇒ Ti,1 = Ti,2 − α ∆y
i,1
(2 + β 2 ) Ti,2 = Ti−1,2 + Ti+1,2 + β 2 (Ti,3 − α ∆y) (i = 2, 3, · · · imax − 1) (3.56)
2 (1 + β 2 ) Ti,jmax −1 = Ti−1,jmax −1 + Ti+1,jmax −1 + β 2 (Ti,jmax −2 + Ti,jmax )

∂T Ti,jmax −Ti,jmax −1
α= ∂y
= ∆y
=⇒ Ti,jmax = Ti,jmax −1 + α ∆y
i,jmax
(2+β 2) Ti,jmax −1 = Ti−1,jmax −1 +Ti+1,jmax −1 +β 2 (Ti,jmax −2 +α ∆y) (i = 2, 3, · · · imax −1)

(3.57)
- 56 -
- Dis rétisation par des s hémas dé entrés d'ordre 2
2 (1 + β 2 ) T2,j = T1,j + T3,j + β 2 (T2,j−1 + T2,j+1 )
−T3,j +4 T2,j −3 T1,j

∂T 1
α= ∂x 1,j
= 2 ∆x
=⇒ T1,j = 3
(4 T2,j − T3,j − 2 α ∆x)
2 (1 + 3 β 2) T2,j = 2 T3,j + 3 β 2 (T2,j−1 + T2,j+1) − 2 α ∆x (j = 2, 3, · · · jmax − 1) (3.58)
2 (1 + β 2 ) Timax −1,j = Timax −2,j + Timax ,j + β 2 (Timax −1,j−1 + Timax −1,j+1)
∂T 3 Timax ,j −4 Timax −1,j +Timax −2,j

α= ∂x imax ,j
= 2 ∆x
1
=⇒ Timax ,j = 3
(4 Timax −1,j − Timax −2,j + 2 α ∆x)
2 (1+3 β 2) Timax −1,j = 2 Timax −2,j +3 β 2 (Timax −1,j−1 +Timax −1,j+1)+2 α ∆x (j = 2, 3, · · · jmax −1)
(3.59)
2 (1 + β 2 ) Ti,2 = Ti−1,2 + Ti+1,2 + β 2 (Ti,1 + Ti,3 )

∂T −Ti,3 +4 Ti,2 −3 Ti,1 1
α= ∂y
= 2 ∆y
=⇒ Ti,1 = 3
(4 Ti,2 − Ti,3 − 2 α ∆y)
i,1
2 (3 + β 2 ) Ti,2 = 3 (Ti−1,2 + Ti+1,2 ) + 2 β 2 (Ti,3 − α ∆y) (i = 2, 3, · · · imax − 1) (3.60)
2 (1 + β 2 ) Ti,jmax −1 = Ti−1,jmax −1 + Ti+1,jmax −1 + β 2 (Ti,jmax −2 + Ti,jmax )

∂T 3 Ti,jmax −4 Ti,jmax −1 +Ti,jmax −2
α= ∂y
= 2 ∆y
i,jmax
1
=⇒ Ti,jmax = 3
(4 Ti,jmax −1 − Ti,jmax −2 + 2 α ∆y)
2 (3+β 2) Ti,jmax −1 = 3 (Ti−1,jmax −1 + Ti+1,jmax −1 )+2 β 2 (Ti,jmax −2 +α ∆y) (i = 2, 3, · · · imax −1)

(3.61)
- 57 -
Exemple :
Soit une plaque arrée de 15 cm de oté et de faible épaisseur de telle sorte qu'on
puisse négliger le transfert de haleur dans ette dire tion. La plaque est soumise à
une température nulle sur les otés gau he et bas, à 100 °C par le haut et isolée sur le
oté droit. On prendra Nx = Ny = 4 et on dis rétisera la CLN par une approximation
entrée.
La formulation à 5 points ave β=1 (puisque ∆x = ∆y ) appliquée aux 4 n÷uds
intérieurs nous donne quatre équations mais ave 6 in onnues, les deux équations sup-
plémentaires sont obtenues par l'équation (3.51) appliquée aux n÷uds internes du bords
droit. Nous aurons don dans l'ordre :



−4 T2,2 + T3,2 + T2,3 = 0


−4 T3,2 + T4,2 + T2,2 + T3,3 = 0





−4 T + 2 T + T = 0

4,2 3,2 4,3


−4 T2,3 + T3,3 + 100 + T2,2 = 0


−4 T3,3 + T4,3 + T2,3 + 100 + T3,2 = 0






−4 T + 2 T + T + 100 = 0

4,3 3,3 4,2
     
−4 1 0 1 0 0 T2,2 0
     
 1 −4 1 0 1 0   T3,2   0 
     
 0 2 −4 0 0 1   T   0 
  4,2  
. =
 
 
 1 0 0 −4 1 0   T2,3   −100 
     
 0 1 0 1 −4 1   T   −100 
   3,3   
0 0 1 0 2 −4 T4,3 −100
T2,2 = 17.37, T3,2 = 25.75, T4,2 = 28.08, T2,3 = 43.73, T3,3 = 57.57, T4,3 = 60.81.
Exer i e :
Que devient ette forme matri ielle si on dis rétise la CLN par des approximations
dé entrées d'ordre 1 et 2?
Solution :
- Pour l'approximation dé entrée arrière d'ordre 1 :
- 58 -
     
−4 1 1 0 T2,2 0
     

 1 −3 0 1  
. T3,2  
= 0 


 1 0 −4 1  
  T2,3  
  −100 

0 1 1 −3 T3,3 −100
T2,2 = 16.84, T3,2 = 24.21, T2,3 = 43.15, T3,3 = 55.79.
- Pour l'approximation dé entrée arrière d'ordre 2 :
     
−4 1 1 0 T2,2 0
     
 2 −8 0 3   T3,2   0 
.
  T  =  −100
   
 1 0 −4 1 
   2,3   
0 3 2 −8 T3,3 −300
T2,2 = 17.56, T3,2 = 26.27, T2,3 = 43.97, T3,3 = 58.34.
Remarque :
Si nous utilisons, dans et exemple, des approximations dé entrées avant, nous
aurons alors 2 équations supplémentaires à ajouter et les matri es obtenues seront
d'ordre 6 omme pour le as de l'approximation entrée.
3.15.5 Equation de Poisson

L'équation de Poisson est l'équation non homogène de l'équation de Lapla e. On
la ren ontre dans les problèmes de diusion de la masse, de diusion de la haleur
( ondu tion), d'é oulement de uides in ompressibles, ...et .
Considérons par exemple notre plaque re tangulaire de ondu tivité thermique k

.
ave une sour e de haleur Q (résistan e par exemple) au milieu. L'EDP représentant
e problème est :
∂2T ∂2T Q̇
2
+ 2
=− (3.62)
∂x ∂y k
En dis rétisant ette équation par un s héma à 5 points, on aura :
Q˙i,j
−2 (1 + β 2 ) Ti,j + Ti−1,j + Ti+1,j + β 2 (Ti,j−1 + Ti,j+1) + ∆x2 ( )=0 (3.63)
k
- 59 -
Si β=1 -à-d ∆x = ∆y alors (3.63) devient :
Q˙i,j
Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1 + ∆x2 ( k
)
Ti,j = (3.64)
4
Q˙i,j
−4 Ti,j + Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1 + ∆x2 ( )=0 (3.65)
k
L'équation de Poisson peut se résoudre de la même manière que l'équation de La-
pla e et tout e que nous avons dit pour ette dernière est appli able à l'équation de
Poisson.
Exemple :
Soit une plaque re tangulaire (k = 0.4 J/cm.s.°C) de 1.5 cm de hauteur de 1 cm
de largeur et de faible épaisseur de telle sorte que l'on puisse négliger le transfert de
haleur dans ette dire tion. La plaque est soumise à une température nulle sur tous
es otés et hauée en son milieu par une résistan e développant un ux de haleur
de 400 J/cm3 .s. En utilisant le logi iel wxMaxima et en é rivant un programme ave
un maillage de 20 × 30, nous obtenons la répartition des isothermes à l'intérieur de la
plaque (Fig.3.15).
Fig. 3.15: Résolution de l'équation de Poisson. Répartition des isothermes.
3.15.6 Méthode d'ordre élevé

Parmi les méthodes d'ordre élevés, on distingue prin ipalement la formulation à 9
2 2
points (Annexe A) pré ise au se ond ordre (O(∆x + ∆y )) :
2 (5 − β 2 )
Ti−1,j−1 + Ti+1,j−1 + Ti−1,j+1 + Ti+1,j+1 + (Ti−1,j + Ti+1,j ) + . . .
1 + β2
- 60 -
2 (5β 2 − 1)
... (Ti,j−1 + Ti,j+1 ) − 20 Ti,j = 0 (3.66)
1 + β2
4
Si β=1 -à-d ∆x = ∆y alors (3.66) devient pré ise au quatrième ordre (O(∆x ))
et sa représentation s hématique est donné par la gure (Fig.3.16).
−20 Ti,j + 4 (Ti−1,j + Ti+1,j + Ti,j−1 + Ti,j+1) + Ti−1,j−1 + Ti+1,j−1 + Ti−1,j+1 + Ti+1,j+1 = 0
(3.67)
1 4 1
4 −20 4
1 4 1
Fig. 3.16: Représentation s hématique de la formulation à 9 points (β = 1).
Exemple 1 :
Soit à déterminer la forme matri ielle du problème de la gure (Fig.3.17).
y 100
3 20
T2,3 T3,3
60
2
T2,2 T3,2
j=1 x
i=1 2 3 4
0
Fig. 3.17: Formulation à 9 points. Exemple 1.
Nombre de n÷uds selon X : Nx = 4.

Nombre de n÷uds selon Y : Ny = 4.
Nombre d'équations du système à résoudre : N = (Nx − 2) × (Ny − 2) = 4.
En appliquant la formule (3.67) à haque n÷ud interne du maillage, nous aurons le
système d'équations suivant :
- 61 -



−20 T2,2 + T1,1 + T3,1 + T1,3 + T3,3 + 4 (T1,2 + T3,2 + T2,1 + T2,3 ) = 0


−20 T + T + T + T + T + 4 (T + T + T + T ) = 0

3,2 2,1 4,1 2,3 4,3 2,2 4,2 3,1 3,3


−20 T2,3 + T1,2 + T3,2 + T1,4 + T3,4 + 4 (T1,3 + T3,3 + T2,2 + T2,4 ) = 0


−20 T + T + T + T + T + 4 (T + T + T + T ) = 0

3,3 2,2 4,2 2,4 4,4 2,3 4,3 3,2 3,4
 


 T2,1 = T3,1 = 0 

 T1,1 = 0.5 (0 + 60) = 30

 

T = T

= 60 T

= 0.5 (0 + 20) = 10
1,2 1,3 4,1
ave les (C.L) : et aux oins :


 T4,2 = T4,3 = 20 

 T1,4 = 0.5 (60 + 100) = 80

 

T = T

= 100 T

= 0.5 (100 + 20) = 60
2,4 3,4 4,4
     
−20 4 4 1 T2,2 −330
     

 4 −20 1 4   T3,2
.
  −110 
= 

 4 1 −20 4   T
  2,3
  −880 
  
1 4 4 −20 T3,3 −660
T2,2 = 37.14, T3,2 = 26.66, T2,3 = 63.33, T3,3 = 52.86.
Exemple 2 :
Déterminer la forme matri ielle du problème de la gure (3.12) et omparer les
résultats entre les deux s hémas.
Exer i e :
1. Pour le problème de l'exemple 1, omparer les résultats obtenus ave les s hémas
à 5 points et à 9 points ave la solution exa te. Dresser un tableau ave les
pour entages relatifs à haque solution.
2. Comparer les pré isions des s hémas à 9 points pour β=1 et pour β 6= 1.
3.15.7 Cas du domaine non re tangulaire

Toutes les méthodes qui ont été présentées jusqu'i i ne on ernent que les domaines
physiques re tangulaires et les maillages re tangulaires. Plusieurs simpli ations ré-
sultent de e hoix :
Les n÷uds du maillage oïn ident ave les limites du domaine physique, don
les (C.L) peuvent être appliquées.
Le maillage est uniforme et orthogonal, don l'approximation des dérivées par-
tielles par diéren es nies est pré ise.
- 62 -
L'espa e du maillage adja ent aux limites du domaine est uniforme et orthogo-
nal.
x x
Domaine physique rectangulaire Domaine physique en 1/4 de cercle

Maillage rectangulaire Maillage rectangulaire
Fig. 3.18: Domaine physique non re tangulaire.
- L'espa e du maillage adja ent à la limite du quart de er le n'est pas uniforme.
Evidemment, une nouvelle appro he en diéren es nies est né essaire.
Il existe plusieurs appro hes pour modéliser les domaines physiques non re tangu-
laires :
1. Limite physique approximative : raner au maximum le maillage an de faire
oïn ider les n÷uds ave les limites du domaine physique.
Cette appro he n'est pas re ommandée.
2. Changement de oordonnées ( ylindriques ou sphériques).
3. Approximation par diéren es nies non uniforme.
4. Transformation des espa es (Fig.3.19).
y η
Transformation
J −1
x ξ
Domaine physique Domaine de calcul
Fig. 3.19: Transformation du domaine physique.
- 63 -
3.15.8 Equation non-linéaire et problème tridimensionnel

- Les équations de Lapla e et de Poisson déjà onsidérées sont des EDP linéaires.
Par onséquent, les équations aux diéren es qui en dé oulent sont linéaires.
Quand une EDP non linéaire est résolue par diéren es nies, il en résulte un
système d'équations aux diéren es non linéaire. On utilise alors deux pro édures pour
sa résolution :
Méthode itérative : → linéarisation → S.O.R. (en 2D et 3D ).

Méthode de Newton : → perturbation → système linéaire → Ja obi, Gauss-
Siedel, S.O.R.
- Les problèmes tridimensionnels (3D ) peuvent être résolus de la même manière que
les problèmes 2D . La seule di ulté est liée à la grande dimension (taille) du système
obtenu en 3D.
- La meilleur solution pour les EDP elliptiques 2D ou 3D linéaires ou non est
l'utilisation de la méthode multi-maillage (Multigrid Method : Brandt 1977).
3.16 Equation d'onde

L'équation d'onde est ren ontrée dans les problèmes d'ingénierie sous deux formes :
Dans les domaines de vibrations et de hamp a oustique :
∂2u 2
2 ∂ u
− c =0 (3.68)
∂t2 ∂x2
u : amplitude de l'onde ; c : vitesse de propagation de l'onde.
Dans les domaines de mé anique des uides et de transferts thermiques :
∂T ∂T
+u =0 (3.69)
∂t ∂x
u : vitesse de onve tion.
3.16.1 Résolution de l'équation d'onde

La résolution de l'équation (3.68) peut se faire en la transformant en un système
d'équations ouplées de type (3.69) :

∂f ∂g
 ∂t + c ∂x =0



(3.70)


 ∂g + c ∂f

=0
∂t ∂x
- 64 -
3.16. Equation d'onde.
Soit le problème mathématique ave ses (C.L) déni par :
∂2u 2
2 ∂ u
= c 0≤x≤π , t≥0 (3.71)
∂t2 ∂x2

u(0, t) = 0 , u(x, 0) = f (x)
u(π, t) = 0 ∂u
, ∂t
(x, o) = g(x)
3.16.1.1 Dis rétisation

L'équation (3.71) est dis rétisée par exemple par un s héma entré d'ordre 2 en
temps et en espa e de la manière suivante :
uin−1 − 2 uni + un+1

i
n n n
2 ui−1 − 2 ui + ui+1
= c
∆t2 ∆x2
∆t 2
λ = c2 ∆x = C te , nous aurons l'équation dis rétisée

en posant : suivante :
un+1 = 2 (1 − λ) uni + λ uni−1 + uni+1 − uin−1

i (3.72)
C'est un s héma à trois étapes temporelles (n − 1, n, n + 1), nous devons don

−1
al uler ui an d'initialiser la pro édure de al ul. Pour ela, utilisons la (C.I) de
nd
Neumann et exprimons là par un s héma aux diéren es entré du 2 ordre :
∂u u1i − u−1
i
(x, o) = = g(xi ) =⇒ u−1 1
i = ui − 2 ∆t g(xi )
∂t 2 ∆t
E rivons maintenant l'équation (3.72) au pas de temps 0 :
u1i u0i u0i−1 u0i+1 u−1

= 2 (1 − λ) +λ + − i
En utilisant la (C.I) u0i = f (xi ) et en remplaçant u−1
i par sa valeur, nous aurons :
λ
u1i = (1 − λ) f (xi ) + [f (xi−1 ) + f (xi+1 )] + ∆t g(xi ) (3.73)
2
Maintenant, nous pouvons al uler toutes les valeurs de l'équation (3.72) en utilisant
les (C.L) suivantes : un1 = 0 et unimax = 0.
3.16.1.2 Stabilité
On suppose que la solution est de la forme : u(x, t) = ψ(t) ej β x , on a alors :
uni = ψ(t) ej β x un+1

i = ψ(t + ∆t) ej β x uin−1 = ψ(t − ∆t) ej β x
uni−1 = ψ(t) ej β (x−∆x) uni+1 = ψ(t) ej β (x+∆x)
- 65 -
Le fa teur d'ampli ation ξ est donné par :
ψ(t + ∆t) ψ(t)

ξ= = (3.74)
ψ(t) ψ(t − ∆t)
aurons :
ψ(t + ∆t) = 2 (1 − λ) ψ(t) + λ e−j β ∆x + ej β ∆x ψ(t) − ψ(t − ∆t)

1
ξ = 2 (1 − λ) + 2λ cos(β ∆x) − ξ
⇐⇒ ξ 2 − 2 [1 − λ (1 − cos(β ∆x))] ξ + 1 = 0
ξ 2 − 2 1 − 2 λ sin2( β ∆x

2
) ξ+1=0
On posant : k = 1 − 2 λ sin2( β ∆x
2
) nous aurons l'équation :
√
ξ2 − 2 k ξ + 1 = 0 dont les solutions sont : ξ1,2 = k ± k2 − 1
Ce fa teur d'ampli ation devant satisfaire la ondition : |ξ| ≤ 1, nous aurons alors
la onditions de stabilité suivante : λ ≤ 1.
3.16.1.3 Autres s hémas de dis rétisation

a- S héma expli ite d'Euler :
un+1 − uni un − uni

i
+ c i+1 =0 (3.75)
∆t ∆x
un+1 − uni un − uni−1
i
+ c i+1 =0 (3.76)
∆t 2 ∆x
Les deux s hémas (3.75) et (3.76) sont pré is d'ordre 1 (ordre minimum dans l'équa-
tion dis rétisée) et sont tous les deux in onditionnellement instables.
b- S héma expli ite Upwind :
un+1
i − uni uni − uni−1
+c =0 (c > 0) (3.77)
∆t ∆x
∆t
Le s héma (3.77) et pré is d'ordre 1 et est onditionnellement stable : 0 ≤ c ∆x ≤ 1.
- 66 -
3.16. Equation d'onde.
Il peut aussi s'é rire :

un+1 ∆t
= uni − c ∆x uni − uni−1

 (c > 0)
 i


(3.78)


un+1
 ∆t
= uni − c ∆x uni+1 − ui n

(c < 0)
i
- S héma expli ite de Lax :
un+1 1
uni+1 + uni−1

i − un − uni−1
2
+ c i+1 =0 (3.79)
∆t 2 ∆x
Ce s héma est pré is d'ordre 1 en temps et d'ordre 2 en espa e et il est ondition-
∆t
nellement stable : c ≤ 1.
∆x
Il existe aussi le s hémas de Lax-Wendro pré is d'ordre 2 et le s hémas de Lax-
Wendro à deux étapes pré is aussi d'ordre 2 mais spé ialement pour les équations
non-linéaires.
d- S héma expli ite de Leapfrog :
un+1 − uin−1 un − uni−1

i
+ c i+1 =0 (3.80)
2 ∆t 2 ∆x
Ce s héma est pé is d'ordre 2 en temps et en 2 espa e et il est onditionnellement
∆t
stable : c ≤ 1.
∆x
e- S héma impli ite d'Euler :
un+1 − uni un+1 − un+1

i
+ c i+1 i−1
=0 (3.81)
∆t 2 ∆x
Ce s héma est pré is d'ordre 1 et il est in onditionnellement stable.
Il existe aussi plusieurs autres s hémas de dis rétisation tels que le s héma expli ite
de Ma -Cormak (prédi teur- orre teur) utilisée pour les EDP non-linéaires ren ontrées
en mé anique des uides par exemple ; le s héma Upwind du 2nd ordre, le s héma
trapézoïdal, le s héma de Rusanov et elui de Runge-Kutta.
3.16.2 Méthodes de résolution

Comme les matri es trouvées pour e genre de problèmes sont souvent reuses, alors
les méthodes de résolution utilisée sont dans l'ordre dé roissant :
S.O.R.(meilleur), Gauss-siedel (préférée), Ja obi (plus simple) et enn elle d'éli-
mination de Gauss.
- 67 -
3.17 Exer i es
3-01 : Parmi les s hémas i-dessous, dire lequel est expli ite ou impli ite :
1. Tin+1 − Tin + λ (Tin − Ti−1

n
) = 0.
2. Tin+1 − Tin + λ (Ti+1
n n+1
− Ti−1 ) = 0.
3. Ti+1,j + Ti−1,j + Ti,j+1 + Ti,j−1 − 4 Ti,j = 0.
4. Tin+1 − 2 (1 − λ) Tin − λ (Ti+1
n n
− Ti−1 ) + Tin−1 = 0.
5. un+1
i − 2 uni + un+1
i−1 = 0.
6. 2 ui,j + ui+1,j − ui−1,j = 0.

7. Tin+1 − Tin−1 + λ (Ti+1
n n
− Ti−1 ) = 0.
3-02 : Dans l'équation de la haleur, dis rétiser :
1. Le terme temporel par une approximation dé entrée avant d'ordre 2.
2. Le terme temporel par une approximation dé entrée arrière d'ordre 1.
3. Le terme temporel par une approximation entrée d'ordre 2.
4. Le terme spatial par une approximation dé entrée arrière d'ordre 2.
3-03 : La dis rétisation de l'équation d'onde donne :
un+1 −2 un n−1
i +ui
un n n
i−1 −2 ui +ui+1
i
∆t2
= c2 ∆x2
1. Ce s héma est-il expli ite ou impli ite ? Justier.
∆t2
2. Déterminer le fa teur d'ampli ation de e s héma (λ = c2 ∆x2
).
3-04 : Dis rétiser l'équation de la haleur par :
1. Une approximation entrée d'ordre 2 en temps et en espa e.
2. Une approximation dé entrée avant d'ordre 1 en temps et d'ordre 2 en espa e.
3. Une approximation dé entrée arrière d'ordre 1 en temps et d'ordre 2 en espa e.
4. Une approximation entrée d'ordre 2 en temps et dé entrée avant d'ordre 1 en
espa e.
- 68 -
3.17. Exer i es
.
3-05 : Soit le s héma de dis rétisation de l'équation de la haleur :
Tin+1 = Tin + ∆t n+1

− 2 Tin+1 + Ti+1
n+1
1 2 n
− 2 Tin + Ti+1
n

∆x2 3
Ti−1 + 3
Ti−1
∆t
2. Déterminer le fa teur d'ampli ation de e s héma (λ = ).
∆x2
3-06 : Considérons l'équation de Lapla e en deux dimensions.
En dis rétisant le terme en x par une approximation dé entrée arrière d'ordre 1,

∆x
le terme en y par une approximation dé entrée avant d'ordre 2 et en posant β = ∆y
,
donner le s héma de dis rétisation (équation) et dessiner les ellules orrespondantes
pour β quel onque et pour β = 2.
3-07 : Considérons l'équation de la haleur unidimensionnelle.
1. Classer ette EDP en justiant votre réponse.
2. Cette équation est dis rétisée par le s héma de Crank-Ni holson. E rire l'équa-
tion dis rétisée et dire si e s héma est expli ite ou impli ite.
∆t
3. Déterminer son fa teur d'ampli ation (on posera λ= ∆x2
).
3-08 : Considérons l'équation de Lapla e appliquée à une plaque re tangulaire.
∆x
1. En posant β= ∆y
, retrouver le s héma à 5 points.
2. Cette plaque est soumise à 100◦ C en bas, 20◦ C à gau he, 40◦ C à droite et 10◦ C
en haut. Elle est dis rétisée par 4 divisions selon X et 5 selon Y.
3. Dessiner les ellules orrespondantes à e s héma pour β = 2.
4. E rire les équations né essaires à la résolution de e problème.
5. Déterminer la forme matri ielle du problème A.T = B.
3-09 : La dis rétisation l'équation de la haleur par la méthode de Ri hardson
donne :
Tin+1 −Tin−1 n −2 T n +T n
Ti−1
2 ∆t
= i
∆x2
i+1
1. Ce s héma est-il impli ite ou expli ite ? Justier.
∆t
2. Déterminer le fa teur d'ampli ation de e s héma (λ = ).
∆x2
- 69 -
3. Etudier sa stabilité.
Pour rendre e s héma stable, on utilise le s héma de Dufort-Frankel :
n −(T n+1 −T n−1 )+T n

Tin+1 −Tin−1 Ti−1
2 ∆t
= i
∆x2
i i+1
1. Ce s héma est-il impli ite ou expli ite ?
2. Déterminer son fa teur d'ampli ation.
3-10 : Dis rétiser l'équation d'onde par :
1. Une approximation entrée d'ordre 2 en temps et en espa e.
2. Une approximation dé entrée arrière d'ordre 1 en temps et dé entrée avant d'ordre
2 en espa e.
En dis rétisant le terme en x par une approximation dé entrée arrière d'ordre 2,

∆x2
le terme en y par une approximation entrée et en posant β= ∆y 2
, donner le s héma
de dis rétisation (équation) et dessiner les ellules orrespondantes pour β=1 et pour
β = 2.
3-12 : Soit le s héma de dis rétisation suivant :
Tin+1 −Tin n −2 T n +T n
Ti−1
∆t
=2 i
∆x2
i+1
∆t
2. Etudier sa stabilité (on posera : λ= ∆x2
).
Soit une plaque re tangulaire dis rétisée par 3 divisions selon X et 4 selon Y. En
∆x2
appliquant le s héma à 5 points et en posant β= ∆y 2
, déterminer la forme matri ielle
du problème A.T = B (les onditions de Diri hlet étant onnues).
3-14 : Soient les deux EDP suivantes :
∂2u ∂2u ∂u ∂2u

∂t2
− ∂x2
=0 et
∂t
− ∂x2
=0
- 70 -
3.17. Exer i es
.
1. Dis rétiser l'équation de la haleur par une approximation dé entrée arrière
d'ordre 1 en temps et d'ordre 2 en espa e.
2. Dis rétiser l'équation d'onde par une approximation entrée en temps et dé en-
trée avant d'ordre 2 en espa e.
3-15 : Soit l'EDP i-dessous :
∂2T ∂T ∂T ∂2T
∂x2
+ ∂x
+ ∂y
+ ∂y 2
=0
∆x
En dis rétisant tous les termes par des approximations entrées et en posant β= ∆y
,
donner le s héma de dis rétisation (équation) et dessiner les ellules orrespondantes
pour β et ∆x quel onques ainsi que pour β=1 et ∆x = 1.

Cette équation est appliquée à une plaque re tangulaire soumise à 100◦ C en bas,
20◦ C à gau he, 40◦ C à droite et 10◦ C en haut. Cette plaque est dis rétisée par 5
divisions selon X et 3 selon Y.

E rire les équations né essaires à la résolution de e problème (β = 1 et ∆x = 1).
Déterminer la forme matri ielle du problème A.T = B .
3-16 : Soit l'équation de la haleur dis rétisée par les trois approximations suivantes :
Tin+1 −Tin−1 T n −2 T n +T n
1.
2 ∆t
− i−1 ∆xi2 i+1 = 0.
T n −Tin−1 Tin −2 Ti+1
n +T n
2. i
∆t
− ∆x2
i+2
= 0.
n n−1 n−2 n n +5 T n −2 T n
Ti−3 −4 Ti−2
3 Ti −4 Ti +Ti
3.
2 ∆t
+ ∆x2
i−1 i
= 0.
Donner, pour haque équation, le type d'approximation hoisie pour haque terme.
3-17 : La dis rétisation de l'équation de onve tion
∂u
∂t
+ c ∂u
∂x
=0
par le s héma de Leapfrog donne :
un+1 −un−1 un n
i+1 −ui−1
i
2 ∆t
i
+c 2 ∆x
=0
2. Quelle est sa pré ision ? Justier.
3. Déterminer son fa teur d'ampli ation.
- 71 -
3-18 : Dis rétiser l'équation de la haleur 2D par un s héma expli ite entré d'ordre 2
en temps et en espa e, ensuite par un s héma impli ite dé entré arrière d'ordre 1 en
temps et dé entré avant d'ordre 1 en espa e.
3-19 : Soit l'équation de la haleur 1D dénie par :
∂T ∂2T
∂t
− ∂x2
=0 0≤ x≤L, t> 0
∂T
T (0, t) = α , ∂x
(L, t) =β, T (x, 0) = γ
En onsidérant le s héma expli ite d'Euler et en dis rétisant la (C.L) par une ap-
proximation dé entrée arrière d'ordre 1, déterminer la forme matri ielle de e problème.
3-20 : Déterminer la forme matri ielle pour l'équation de Lapla e appliquée à une
plaque (Nx = Ny = 4) en utilisant le s héma à 5 points et en dis rétisant la (C.L) de
Neumann par une approximation entrée. On donne :

∂T
Tg = 10°C , Td = 30°C , Tb = 10°C , ∂y
= 0.
h
plaque (Nx = Ny = 4) en utilisant le s héma à 9 points et en dis rétisant les (C.L) de
Neumann par des approximations dé entrées d'ordre 1. On donne :

∂T ∂T

Tb = 100°C , Th = 40°C , ∂x g
= 0, ∂x d
= 0.
3-22 : Dis rétiser l'équation de la haleur 2D par un s héma impli ite entré d'ordre
2 en temps et en espa e, ensuite par un s héma expli ite dé entré avant d'ordre 1 en
temps et dé entré arrière d'ordre 1 en espa e.
∂T ∂2T
∂t
− ∂x2
=0 0≤ x≤L, t> 0
∂T
∂x
(0, t) =α, T (L, t) = β , T (x, 0) = γ
proximation dé entrée avant d'ordre 1, déterminer la forme matri ielle de e problème.
Neumann par des approximations entrées. On donne :

∂T

∂T
Tb = 100°C , Th = 40°C , ∂x g
= 0, ∂x d
= 0.
- 72 -
3.17. Exer i es
.
∂T ∂2T
∂t
− ∂x2
=0 0≤ x≤L, t> 0
∂T
∂x
(0, t) =α, T (L, t) = β , T (x, 0) = γ
proximation dé entrée avant d'ordre 2, déterminer la forme matri ielle de e problème.
3-26 : Obtenir la forme matri ielle pour l'équation de Lapla e appliquée à la plaque
s hématisée en dessous en utilisant le s héma à 9 points et en dis rétisant la (C.L) de
Neumann par :
1. une approximations dé entrée d'ordre 1.
2. une approximations dé entrée d'ordre 2.
3. une approximations entrée.
On donne :
∂T

T1 = 40°C , T2 = 20°C , T3 = 100°C , α = ∂x g
= 10°C/cm, ∆x = 0.2 cm.
T3
α 3
T2
2
1
1 2 3 4 5
T1
Fig. 3.20: Exo-3-26.

∂T ∂T

Tb = 100°C , Th = 40°C , ∂x g
= 10, ∂x d
= 20, ∆x = 2.
∂T ∂2T
∂t
− ∂x2
=0 0≤ x≤L, t> 0
- 73 -
∂T ∂T
∂x
(0, t) =α, ∂x
(L, t) =β, T (x, 0) = γ
En onsidérant le s héma expli ite d'Euler et en dis rétisant les (C.L) par des ap-
proximations dé entrées d'ordre 1, déterminer la forme matri ielle de e problème.

∂T ∂T

Tb = 100°C , Th = 40°C , ∂x g
= 10, ∂x d
= 20, ∆x = 1.
∂T ∂2T
∂t
− ∂x2
=0 0≤ x≤L, t> 0
∂T ∂T
∂x
(0, t) =α, ∂x
(L, t) =β, T (x, 0) = γ
En onsidérant su essivement les s hémas impli ite et expli ite d'Euler et en dis-
rétisant les (C.L) par des approximations dé entrées d'ordre 2, déterminer les formes
matri ielles de es problèmes.
s hématisée en dessous en utilisant le s héma à 5 points et en dis rétisant la (C.L) de
Neumann par une approximations entrée. On donne :
∂T

T1 = 40°C , T2 = 20°C , T3 = 100°C , α = ∂x g
= 0.
T3
α 3
T2
2
1
1 2 3 4 5
T1
Fig. 3.21: Exo-3-31.
- 74 -
3.17. Exer i es
.
3-32 : Soit une barre, de longueur L et de très faible se tion, soumise aux onditions
α et β respe tivement à ses extrémités gau he et droite. Ce phénomène est régit par
l'équation suivante :
∂T 2
∂t
− 3 ∂∂xT2 = 0
2. Dis rétiser le terme temporel par une approximation dé entrée avant d'ordre 1
et le terme spatial par une approximation entrée au temps n.

∆t
3. En posant λ= ∆x2
, é rire l'équation dis rétisée et justier si le s héma obtenu
est expli ite ou impli ite.
4. Etudier la stabilité de e s héma en utilisant le ritère de Von-Neumann.
5. En supposant que α et β sont des onditions de Diri hlet, déterminer la forme
matri ielle du problème (pré iser les bornes pour i et n).

∂T
6. En suppose maintenant que α est une ondition de Neumann ( α = ∂x x=0
). En
dis rétisant ette ondition par une approximation dé entrée avant d'ordre 1,
déterminer la nouvelle forme matri ielle du problème (pré iser les bornes pour
i et n).
3-33 : Soit l'EDP suivante :
∂2T ∂2T
∂x2
+ 2 ∂T
∂x
+ 2 ∂T
∂y
+ ∂y 2
=0
2. Dis rétiser tous les termes par des approximations entrées.
∆x
3. En posant β= ∆y
, é rire l'équation dis rétisée.
4. Dessiner les ellules orrespondantes pour β et ∆x quel onques ainsi que pour
β=2 et ∆x = 1/2.
5. Cette équation est appliquée à une plaque re tangulaire soumise à 100◦ C en bas,
◦ ◦ ◦
20 C à gau he, 40 C à droite et 10 C en haut. En dis rétisant ette plaque par
4 divisions selon X et 3 selon Y :
(a) E rire les équations né essaires à la résolution de e problème (β = 2 et
∆x = 1/2).
(b) Déterminer la forme matri ielle du problème A.T = B .

∂T
6. En suppose maintenant une ondition de Neumann ( α = = 10). En dis-
∂x g
rétisant ette ondition par une approximation entrée, déterminer la nouvelle
forme matri ielle du problème.
- 75 -
s hématisée en dessous en utilisant le s héma à 5 points et en dis rétisant les (C.L) de
Neumann par une approximations dé entrées d'ordre 1. On donne :
T1 = 100°C , T2 = 20°C , T3 = 10°C , T4 = 40°C ,

∂T ∂T
α= ∂x d
= 2 y + 10 °C/cm, β= ∂y
= 3 x + 10 °C/cm, ∆x = 2 cm, ∆y = 2 cm.
h
T3
α
5
T2 3
T4
1
1 2 3 4 5 6
T1
Fig. 3.22: Exo-3-34.
3-35 : Même exer i e que 3-34 mais ave le s héma à 9 points au lieu du s héma à 5
points.
3-36 : En appliquant le s héma de dis rétisation à 5 points à l'équation de Lapla e pour
la plaque s hématisée i-dessous, déterminer la matri e A et le ve teur B du système
obtenu : A . T = B.
- 76 -
3.17. Exer i es
.
100
4 60
20 3
1
1 2 3 4 5 6 7
Fig. 3.23: Exo-3-36.
3-37 : En appliquant le s héma de dis rétisation à 5 points à l'équation de Lapla e pour
la plaque s hématisée i-dessous, déterminer la matri e A et le ve teur B du système
obtenu : A . T = B.
5
20
4 60
60 3
1
1 2 3 4 5 6 7
40
Fig. 3.24: Exo-3-37.
s hématisée en dessous en utilisant le s héma à 5 points et en dis rétisant les (C.L) de
Neumann par une approximations dé entrées d'ordre 1. On donne :
T1 = 10°C , T2 = 50°C , T3 = 30°C , T4 = 100°C ,

∂T ∂T
α= ∂x d
= 5°C/cm, β= ∂y
= 20°C/cm, ∆x = 2 cm, ∆y = 1 cm.
h
- 77 -
T3
α
5
T2 3
T4
2
1
1 2 3 4 5
T1
Fig. 3.25: Exo-3-38.
i-dessous en utilisant le s héma à 9 ave β = 1. On donne :
T1 = 10°C , T2 = 100°C , T3 = 2.104 x2 , T4 = 20°C , T5 = 10°C , T6 = 104 y 2 ,

∆x = 2 cm.
T3
5 T4
T5
T2 3
T6
2
1
1 2 3 4 5
T1
Fig. 3.26: Exo-3-39.
3-40 : Même problème que 3-39 ave les données suivantes :
T1 = 5.105 x2 , T2 = 3.105 y 2 , T3 = 60°C, T4 = 30°C , T5 = 10°C , T6 = 100°C,

∆x = 3 cm.
∂2T ∂T
3-41 : Quelle est la nature de l'EDP suivante : ∂x2
− ∂y
= 0.
- 78 -
3.17. Exer i es
.
Cette EDP est dis rétisée par les trois s hémas suivants :
2 Ti,j −5 Ti+1,j +4 Ti+2,j −Ti+3,j 4T −3 T −T
1.
∆x2
− i,j+1 2 ∆yi,j i,j+2 = 0.
Ti,j −2 Ti−1,j +Ti−2,j T −T
2.
∆x2
+ i,j ∆yi,j+1 = 0.
2 Ti,j −Ti+1,j −Ti−1,j T −Ti,j
3.
∆x2
− i,j−1 ∆y
= 0.
Donner, pour haque équation, le type d'approximation hoisie pour haque terme.
3-42 : Soit le s héma de dis rétisation suivant :

n+1
Tin+1 −Tin n −2 T n +T n −2 Tin+1 +Ti+1
n+1
h i h i
Ti−1 Ti−1
∆t
= (1 − θ) i
∆x2
i+1
+θ ∆x2
où θ est une onstante omprise entre 0 et 1.

∆t
1. Déterminer son fa teur d'ampli ation (prendre : λ= ∆x2
).
2. Etudier sa stabilité pour les as où θ = 0, θ = 1/2 et θ = 1.
3-43 : Soit le s héma de dis rétisation de l'équation de la haleur :
Tin+1 = Tin + λ n+1

− 2 Tin+1 + Ti+1
n+1 λ n
− 2 Tin + Ti+1
n

3
Ti−1 + 3
Ti−1
1. Déterminer son fa teur d'ampli ation en appliquant l'analyse de Von-Neumann.
2. Etudier sa stabilité sa hant que λ > 0.
3-44 : En utilisant le s héma à 9 points, déterminer la forme matri ielle A . T = B pour

2
l'équation de Lapla e appliquée à une plaque re tangulaire de (24 × 18) cm dis rétisée
par 4 divisions selon x et 3 selon y. Les deux onditions de Neumann seront dis rétisée
par des approximations dé entrées d'ordre 1.
On donne : Les températures en °C et les gradients en °C/ m.

∂T ∂T

Tb = 200 x + 10, Tg = 300 y + 20, Th = ∂y
= 5, Td = ∂x
= 10.
2 ∆x2
(54 x 12 cm ) i-dessous en utilisant le s héma à 5 ave β= ∆y 2
. On donne :
T0 = 10°C , T1 = 50°C , T2 = 20°C , T3 = 103 x2 , T4 = 2.103 y .
T3
T2
T4
T0
T1
Fig. 3.27: Exo-3-45.
- 79 -
i-dessous en utilisant le s héma à 5 ave β = 1. On donne :
T1 = 10°C , T2 = 100°C , T3 = 2.104 x2 , T4 = 20°C , T5 = 10°C , T6 = 104 y 2 ,

∆x = 2 cm.
T5
5
T4
4
T3
T6
3
T2 2
1
1 2 3 4 5
T1
Fig. 3.28: Exo-3-46.
3-47 : Considérer l'équation de Lapla e en deux dimensions et obtenir le s héma à 5

∆x
points en fon tion de β= ∆y
.
Obtenir la forme matri ielle en appliquant e s héma à la plaque de l'exer i e 3-45
pour β = 1. On donne :
T0 = 10°C , T1 = 50°C , T2 = 20°C , T3 = 40◦ C , T4 = 60◦ C .
- 80 -
Chapitre 4
Méthode des Volumes Finis
Dans le but de simuler les é oulements de uides, de transferts thermiques ainsi
que d'autres phénomènes physiques relatifs, il est né essaire de dé rire la physique
asso iée par des termes mathématiques. La grande partie des phénomènes qui nous
intéresse sont gouvernés par des prin ipes de onservation et sont régis par les EDP
exprimant es prin ipes. Par exemple, les équations de quantité de mouvement ex-
priment la onservation de la quantité de mouvement, l'équation d'énergie exprime la
onservation de l'énergie totale, ...et .
Dans e paragraphe, nous allons faire ressortir l'équation de onservation typique
et examiner es propriétés mathématiques.
4.1 Equations de onservation

Les équations typiques qui dé rivent la onservation de la masse, de la quantité de
mouvement, d'énergie ou des espè es himiques sont é rites en fon tion de quantités
spé iques -à-d des quantités exprimées par unité de masse. Considérons la quantité
spé ique φ qui peut être soit la quantité de mouvement par unité de masse, soit
l'énergie par unité de masse ou autre. Considérons un volume de ontrle (v.c) de
dimensions ∆x, ∆y et ∆z (Fig.4.1).
Nous voulons exprimer les variations de la quantité spé ique φ dans le volume
de ontrle (v.c) en fon tion du temps. Supposons que φ est gouvernée par le prin-
ipe de onservation qui stipule que : L'a umulation de φ dans le (v.c) pendant le
temps ∆t est égale au ux net de φ dans le (v.c) augmenté de la génération net de φ
à l'intérieur du (v.c). En traduisant e prin ipe nous aurons sous forme mathématique :
(ρ φ ∆v)t+∆t − (ρ φ ∆v)t = (Jx − Jx+∆x ) ∆y ∆z ∆t + (Jy − Jy+∆y ) ∆x ∆z ∆t . . .

. . . + (Jz − Jz+∆z ) ∆x ∆y ∆t + S ∆v ∆t
81
Ch4.MÉTHODE DES VOLUMES FINIS.
1111111
0000000 111111
000000
Y Jy+∆y
Jz
0000000
1111111
0000000
1111111 000000
111111
0000000
1111111 000000
111111
0000000
1111111 000000
111111
0000000
1111111 000000
111111
0000000
1111111 ∆x
000000
111111
000000
111111
0000000
1111111 000000
111111
Jx+∆x
0000000
1111111 000000
111111
Jx o
0000000
1111111
S X
0000000
1111111 000000
111111
∆y
0000000
1111111 000000
111111
000000
111111
∆z
0000000
1111111 000000
111111
Z 0000000
1111111 000000
111111
Jz+∆ z Jy
Fig. 4.1: Equations de onservation.
S représente le terme sour e (terme de génération) de φ par unité de volume et par
unité de temps.
Pour les phénomènes physiques qui nous intéresse, la quantité φ est transportée
par deux mé anismes primaires : la diusion due à la ollision des molé ules et la
onve tion due au mouvement du uide. Dans beau oup de as, les ux diusif et
onve tif s'é rivent par exemple selon l'axe X :

∂φ
Flux diusif : Jdif f = −Γ ∂x
Flux onve tif : Jconv = ρ u φ
Flux net : Jx = (ρ u φ − Γ ∂φ )
∂x x
et Jx+∆x = (ρ u φ − Γ ∂φ )
∂x x+∆x
où Γ est une ara téristique du uide, (ρ u)x est le ux massique à travers la fa e x du
→
− →
− →
− →
−
(v.c) et u la omposante du hamp de vitesse donné par : q =u i +v j +w k.
En divisant l'équation de onservation par (∆v ∆t), on obtient :
(ρ φ)t+∆t − (ρ φ)t Jx − Jx+∆x Jy − Jy+∆y Jz − Jz+∆z

= + + +S
∆t ∆x ∆y ∆z
et en prenant les limites quand ∆x, ∆y ,∆z et ∆t tendent vers 0, on aura :
∂(ρ φ) ∂Jx ∂Jy ∂Jz

=− − − +S
∂t ∂x ∂y ∂z
ou en ore en développant :
∂ ∂ ∂ ∂ ∂ ∂φ ∂ ∂φ ∂ ∂φ
(ρ φ) + (ρ u φ) + (ρ v φ) + (ρ w φ) = (Γ ) + (Γ ) + (Γ )+S
∂t ∂x ∂y ∂z ∂x ∂x ∂y ∂y ∂z ∂z
(4.1)
- 82 -
4.1. Equations de onservation.
et sous forme ve torielle ompa te :
∂ →
− →
− →
−
(ρ φ) + ∇ . (ρ −
→
q φ) = ∇ . (Γ ∇ φ) + S (4.2)
∂t
L'équation (4.2) est la forme onservative des équations gouvernantes pour les é ou-
lements de uides, les transferts thermiques et de masse ainsi que d'autres équations
de transport.
L'équation (4.1) peut aussi s'é rire :
∂ ∂Jx ∂Jy ∂Jz

(ρ φ) + + + =S (4.3)
∂t ∂x ∂y ∂z
En l'absen e du terme sour e et en mode stationnaire, la divergen e du ux est
nulle :
− −
→ →
∇. J =0 (4.4)
où :
−
→ →
− →
− →
−
J = Jx i + Jy j + Jz k
Considérons maintenant ertains as spé iques de l'équation de onservation de φ.

1- Equation de onservation de la masse : En posant φ = 1 et S = 0 l'équation
(4.2) devient :
∂ →
−
(ρ) + ∇ . (ρ →
−
q)=0 (4.5)
∂t
2- Equation d'énergie : Soit h l'enthalpie spé ique, k la ondu tivité thermique.
k
En posant :φ = h, Γ = et S = Sh l'équation (4.2) devient :
Cp
∂ −
→ → k −
− →
(ρ h) + ∇ . (ρ −
→
q h) = ∇ . ( ∇ h) + Sh (4.6)
∂t Cp
3- Equation de quantité de mouvement : En posant : φ = u, Γ = µ et S = − ∂P
∂x
+ Sx
l'équation (4.2) devient :
∂ →
− →
− →
− ∂P
(ρ u) + ∇ . (ρ −
→
q u) = ∇ . (µ ∇ u) − + Sx (4.7)
∂t ∂x
Cette équation est elle de quantité de mouvement pour un uide Newtonien se-
lon l'axe X , Sx ontient les parties du tenseur des ontraintes qui n'apparaissent pas
∂P
dire tement dans le terme de diusion et
∂x
est le gradient de pression selon l'axe X.
De même, nous pouvons aussi é rire l'équation de transport des espè es himiques
d'un mélange, ...et .
- 83 -
L'équation (4.2) est appelée équation générale de transport s alaire et elle
peut aussi s'é rire de la manière suivante :
∂ −−→
(ρ φ) + div . (ρ →
−
q φ) = div . (Γ grad φ) + S (4.8)
∂t
Exemples :
1- Equation de diusion (Elliptique) :
∂ −
→ →
− −
→ →
−
obtenue pour
∂t
=0 et q = 0 =⇒ ∇ . (Γ ∇ φ) + S = 0.
selon l'axe X (1D) on aura : ∂x ∂
(Γ ∂φ
∂x
) + S = 0.
k
si φ = T , Γ = on obtient l'équation de la haleur unidimensionnelle stationnaire :
Cp
∂ k ∂T
(
∂x Cp ∂x
) + S = 0 et si k = C te et Cp = C te alors ette équation devient :
∂2T S
α 2
+ =0
∂x ρ
2- Equation de la haleur (parabolique) :
∂ ∂
même problème que l'exemple 1 mais instationnaire :
∂t
(ρ T ) = ( k ∂T )
∂x Cp ∂x
+ S.
si ρ = C te , k = C te et Cp = C te alors ette équation devient :
∂T ∂2T S
=α 2 +
∂t ∂x ρ
3- Equation d'onde (Hyperbolique) :
obtenue pour S = 0, Γ = 0, φ = u et
−
→
q = C te = c.
où u et c sont respe tivement l'amplitude et la élérité de l'onde.
∂ ∂
∂t
(ρ u) + ∂x (ρ c u) = 0
si ρ = C te alors ette équation devient :
∂u ∂u
+c =0
∂t ∂x
Remarque :
Dans la majorité des problèmes ren ontrés en ingénierie, l'équation (4.8) peut avoir
un omportement mixte : le terme de diusion tend à la rendre elliptique, les termes
instationnaire et onve tif tendant à la rendre parabolique ou hyperbolique.
4.2 Equation de diusion

Dans e paragraphe, nous allons tourner notre attention vers un pro essus important
qu'on appelle diusion.
- 84 -
4.2. Equation de diusion.
Les opérateurs de diusion sont ommun en transferts de haleur, de masse et de
quantité de mouvement et peuvent aussi être utilisés pour modéliser l'éle trostatique,
le rayonnement ainsi que d'autre phénomènes physiques.
Nous onsidérons la dis rétisation et la solution de l'équation s alaire de transport
pour les problèmes de diusion stationnaires et instationnaires.
4.2.1 Problème de diusion unidimenionnel

La méthode des volumes nis (MVF) ou en ore méthode des volumes de ontrles
(v.c) divise le domaine d'étude en un nombre ni de ellules ou (v.c) à travers lesquels
la onservation de φ est respe tée.
Examinons le pro essus de dis rétisation à travers l'équation de diusion 1D ave
terme sour e :
d dφ
(Γ ) + S = 0 (4.9)
dx dx
Considérons un maillage unidimensionnel ave des mailles ( ellules) omme indiqué
sur la gure (Fig.4.2). Les valeurs dis rètes de φ sont sto kées aux entres des mailles
et représentées par W, P et E. Les fa es de la maille sont notées w et e.
δx
Limites du (v.c)
W P E
PSfrag repla ements w e x
δxw δxe
Fig. 4.2: Disposition du (v. ) 1D.
Considérons la maille asso iée ave le n÷ud P. Commençons par intégrer l'équation
(4.9) à travers le (v.c), nous aurons :
d
(Γ dφ d
(Γ dφ
R R R
dx dx
) + S dv = 0 =⇒ dx dx
) dv + S dv = 0 =⇒
(vc) (vc) (vc)
Re Re
d
dx
(Γ dφ
dx
) A dx + S A dx = 0 =⇒
w w
Ze
dφ dφ
ΓA − ΓA + S A dx = 0 (4.10)
dx e dx w
w
Nous allons maintenant supposer un prol de φ dé rivant sa variation entre les
n÷uds entraux (entre W et P et entre P et E ). Si nous supposons une variation
- 85 -
linéaire de φ, nous pouvons alors exprimer les dérivées aux interfa es du (v.c) (Annexe
B) et l'équation (4.10) devient :
Γe Ae (φE − φP ) Γw Aw (φP − φW )
− + S̄ A δx = 0 (4.11)
δxe δxw
ΓP + ΓE ΓP + ΓW
où : Γe = 2
et Γw = 2
S̄ étant la moyenne de S dans le (v.c).
Dans les situations pratiques, le terme sour e peut être une fon tion de la variable
dépendante φ. Dans es as, la (MVF) utilise une approximation linéaire de e terme
de la forme :
S̄ ∆v = Su + SP φP (4.12)
On dit que le terme sour e est linéarisé (nous verrons par la suite sa forme générale).
En substituant (4.12) dans (4.11) et en regroupant les termes, nous aurons :

Γw Aw Γe Ae Γw Aw Γe Ae
+ − SP φP = φW + φ E + Su (4.13)
δxw δxe δxw δxe
Remarque :
On note que l'équation (4.13) n'est pas totalement exa te à ause de l'approxima-
tion linéaire de φ entre les noeuds.
L'équation (4.13) peut s'é rire sous la forme suivante :
aP φP = aW φW + aE φE + Su (4.14)
ave :
Γw Aw Γe Ae
aW = ; aE = ; aP = aW + aE − SP
δxw δxe
Des équations identiques à (4.14) peuvent être établies pour l'ensemble des mailles
du domaine ex epté eux des mailles situées sur les limites du domaine. Ces dernières
devront être traitées diéremment en fon tion du type de C.L. (Diri hlet, Neumann ou
mixte ). Ce i qui nous onduira à un système d'équations algébriques dont le nombre
est égal à elui des noeuds entraux et qui pourra être résolu par exemple par une
méthode itérative.
Remarques :
On note les points suivants pour le pro essus de dis rétisation :
- 86 -
1. On ommen e par é rire le prin ipe de onservation à travers la maille ( ellule
ou v.c) qui nous onduit aux valeurs de φ qui satisfons ette onservation. De
e fait, la onservation est satisfaite pour haque maille du maillage.
2. La onservation ne garantie pas la pré ision. La solution pour φ peut être non
pré ise mais onservative.
− Γ dφ

3. La quantité
dx e
est le ux diusif à travers la fa e e. L'équilibre dans la
maille est é rit en fon tion des ux à travers ses fa es. Le gradient de φ doit
don être évalué sur les fa es de la maille.
4. Les prols supposés pour φ et S ne sont pas né essairement identiques.
4.2.2 Exemple 1 : Equation de diusion 1D

AB de lon-
Cal uler la répartition de la température stationnaire à travers la barre
2
gueur L = 0, 5 m, de se tion A = 0, 01 m , de ondu tivité thermique k = 1000 W/m°C
et soumise à ses extrémités aux températures TA = 100°C et TB = 500°C (Fig.4.3).
Solution :
d
L'équation à dis rétiser est la suivante :
dx
(k dT
dx
) =0
Divisons notre barre en 5 parties égales, don δx = 0, 1 m.
Nous avons d'après les données : kw = ke = k , δxw = δxe = δx et Aw = Ae = A.
δx
1 2 3 4 5
TA TB
PSfrag repla ements
δx/2 δx δx δx δx δx/2
Fig. 4.3: Exemple 1 : Equation de diusion 1D.
Noeuds internes : d'après l'équation (4.14) on aura : aP φP = aW φW + aE φE ave :
kA kA 2kA
aW = δx
; aE = δx
; aP = δx
=⇒ SP = (aW + aE ) − aP = 0; Su = 0
Cette équation s'applique pour les n÷uds 2, 3 et 4. Les n÷uds 1 et 5 sont des n÷uds
limites et devront être traités diéremment.
- 87 -
N÷ud 1 : En intégrant l'équation à travers le (v.c) entourant le n÷ud 1 on aura
d'après (4.11) :
k A (TE −TP ) k A (TP −TA )

δx
− δx/2
= 0 =⇒ 3 kδxA TP = 0 TW + kA
δx
TE + 2 kδxA TA
et en omparant ave (4.14) on aura :
kA
aW = 0 ; aE = δx
; aP = 3 kδxA =⇒ SP = −2 kδxA ; Su = 2 kδxA TA
i i P =1 et E = 2.
d'après (4.11) :
k A (TB −TP ) k A (TP −TW )

δx/2
− δx
= 0 =⇒ 3 kδxA TP = kA
δx
TW + 0 TE + 2 kδxA TB
kA
aW = δx
; aE = 0 ; aP = 3 kδxA =⇒ SP = −2 kδxA ; Su = 2 kδxA TB
i i P =5 et W = 4.
kA 1000.0,01
δx
= 0,1
= 100.
Nous aurons alors le système algébrique suivant :



300 T1 = 100 T2 + 200 TA


200 T2 = 100 T3 + 100 T1





200 T3 = 100 T4 + 100 T2


200 T4 = 100 T5 + 100 T3






300 T

= 100 T4 + 200 TB
5
     
−300 100 0 0 0 T1 −2.104
100 −200 100 0 0 T2 0
     
     
     
 0 100 −200 100 .
0 T3 = 0 
     
0 0 100 −200 100   T4 0
     
   
0 0 0 100 −300 T5 −105
La solution numérique est déterminée en utilisant le logi iel de al ul Maple :
T1 = 140 °C; T2 = 220 °C; T3 = 300 °C; T4 = 380 °C; T5 = 460 °C.
- 88 -
Exer i e :
Résoudre e problème ave Maple et tra er la solution en la omparant ave la
solution exa te.
4.2.3 Exemple 2 : Equation de diusion 1D ave terme sour e

Cal uler la répartition de la température stationnaire à travers la plaque (Fig.4.4)
de longueur L = 2 cm, de se tion A = 1 m2 , de ondu tivité thermique k = 0, 5 W/m°C

et soumise à ses extrémités aux températures TA = 100°C et TB = 200°C et à une sour e
de haleur q = 1000 kW/m3 .
111111
000000
000000
111111
z
y
000000
111111
x
000000
111111
000000
111111
A
000000
111111
000000
111111
T
000000
111111
B
TA
000000
111111
L
000000
111111
Fig. 4.4: Exemple 2 : Equation de diusion 1D ave terme sour e.
Solution :
d
dx
(k dT
dx
) +q =0
N÷uds internes : d'après l'équation (4.14) on aura : aP TP = aW TW + aE TE + Su
kA kA
ave : aW = δx
; aE = δx
; aP = 2 kδxA =⇒ SP = 0; Su = q A δx
Cette équation s'applique pour les n÷uds 2, 3 et 4. Les n÷uds 1 et 5 sont des n÷uds
limites et devront être traités diéremment.
d'après (4.11) :
k A (TE −TP )
δx
− k A (TP −TA )
δx/2
+ q A δx = 0=⇒ 3 kδxA TP = 0 TW + kδxA TE + 2 kδxA TA + q A δx
- 89 -
kA
aW = 0 ; aE = δx
; aP = 3 kδxA =⇒ SP = −2 kδxA ; Su = 2 kδxA TA + q A δx
i i P =1 et E = 2.
d'après (4.11) :
k A (TB −TP )
δx/2
− k A (Tδx
P −TW )
+ q A δx = 0=⇒ 3 kδxA TP = kA
δx
TW + 0 TE + 2 kδxA TB + q A δx
kA
aW = δx
; aE = 0 ; aP = 3 kδxA =⇒ SP = −2 kδxA ; Su = 2 kδxA TB + q A δx
i i P =5 et W = 4.
kA 0,5.1
δx
= 0,004
= 125 et q A δx = 106 .1.0, 004 = 4000.
Nous aurons alors le système algébrique suivant :




375 T1 = 125 T2 + 4000 + 250 TA


250 T2 = 125 T1 + 125 T3 + 4000





250 T3 = 125 T2 + 125 T4 + 4000


250 T4 = 125 T3 + 125 T5 + 4000






375 T

= 125 T4 + 4000 + 250 TB
5
     
−375 125 0 0 0 T1 −29.103
 125 −250 125 0 0   T2   −4000 
     
     
 0 125 −250 125 0  .  T3  =  −4000 
     
 0 0 125 −250 125 T −4000 
     
  4  
0 0 0 125 −375 T5 −54.103
T1 = 150 °C; T2 = 218 °C; T3 = 254 °C; T4 = 258 °C; T5 = 230 °C.
Exer i e :
solution exa te.
- 90 -
4.2.4 Exemple 3 : Equation de diusion 1D ave onve tion à

la surfa e et bord isolé
Cal uler la répartition de la température stationnaire d'une ailette ylindrique de
se tion uniforme A, de longueur L = 1m et ondu tivité thermique k, exposée à
un ux de haleur de température ambiante T∞ = 20°C, soumise à une température
TB = 100°C sur son bord gau he et isolée sur son bord droit (Fig.4.5).
1111111
0000000
0000000
1111111
0000000
1111111
0000000
1111111 Too A
0000000
1111111
0000000
1111111
T B bord isolé
0000000
1111111
0000000
1111111
L
0000000
1111111
0000000
1111111
0000000
1111111
Fig. 4.5: Exemple 3 : Equation de diusion 1D, onve tion en surfa e et bord isolé.
Solution :

d
dx
(k A dT
dx
) − h P (T − T∞ ) = 0
h étant le oe ient de transfert onve tif et P le périmètre de l'ailette.
hP
En la divisant par kA et en posant n2 = kA
ette équation devient :
d dT
( )
dx dx
− n2 (T − T∞ ) = 0
Le 1er terme de ette équation est traité de la même manière que les exemples 1 et 2.
Le 2nd terme est évalué en supposant que le terme sour e est onstant pour haque
(v.c), don :
R d dT R 2
( ) dv −
dx dx
n (TP − T∞ ) dv = 0
(vc) (vc)
=⇒ A dT − A dT − [n2 (TP − T∞ )] A δx = 0

dx e dx w
TE −TP TP −TW
− [n2 (TP − T∞ )] δx = 0

=⇒ δx
− δx
2 1 1
+ n2 δx TP = TE + n2 δx T∞

=⇒ δx δx
TW + δx
- 91 -
1 1 2
aW = δx
; aE = δx
; aP = δx
+ n2 δx =⇒ SP = −n2 δx ; Su = n2 δx T∞
N÷ud 1 : La fa e ouest du (v.c) qui entoure le n÷ud 1 est soumise à une C.L de
type Diri hlet et don elle est traitée de la même manière que pour l'exemple 1.
h i
TE −TP TP −TB
− [n2 (TP − T∞ )] δx = 0

δx
− δx/2
3 1 2
+ n2 δx TP = 0 TW + TE + n2 δx T∞ +

=⇒ δx δx δx
TB
1 3 2 2
aW = 0 ; aE = δx
; aP = δx
+ n2 δx =⇒ SP = −n2 δx − δx
; Su = n2 δx T∞ + δx
TB
i i P =1 et E = 2.
N÷ud 5 : La fa e est du (v.c) qui entoure le n÷ud 5 est soumise à une C.L de type
Neumann, 'est à dire à un ux qui est dans e as nul (fa e isolée), don :
TP −TW
− [n2 (TP − T∞ )] δx = 0

0− δx
1 1
+ n2 δx TP = TW + 0 TE + n2 δx T∞

=⇒ δx δx
1 1
aW = δx
; aE = 0 ; aP = δx
+ n2 δx =⇒ SP = −n2 δx ; Su = n2 δx T∞
i i P =5 et W = 4.
1
En prenant
δx
=5 et n2 = 25, nous aurons alors le système algébrique suivant :



20 T1 = 5 T2 + 5 T∞ + 10 TB


15 T2 = 5 T1 + 5 T3 + 5 T∞





15 T3 = 5 T2 + 5 T4 + 5 T∞


15 T4 = 5 T3 + 5 T5 + 5 T∞






10 T

= 5 T4 + 5 T∞
5
- 92 -
     
−20 5 0 0 0 T1 −1100
 5 −15 5 0 0   T2   −100 
     
     
 0 5 −15 5 0  .  T3  =  −100 
     
 0 0 5 −15 5   T4   −100 
     
0 0 0 5 −10 T5 −100
T1 = 64, 22 °C; T2 = 36, 91 °C; T3 = 26, 50 °C; T4 = 22, 60 °C; T5 = 21, 30 °C.
Exer i e :
solution exa te.
Comparer ensuite les solutions ave 5 (v.c) et ave 10 (v.c).
4.2.5 Problème de diusion bidimensionnel

Considérons l'équation de diusion stationnaire 2D en oordonnées re tangulaires
ave terme sour e :

∂ ∂φ ∂ ∂φ
(Γ )+ (Γ )+S =0 (4.15)
∂x ∂x ∂y ∂y
Intégrons ette équation à travers le (v.c) entourant le point P (Fig.4.6), nous
aurons :
δx
y
N (v.c)
δyn n
PSfrag repla ements W P E δy

w e x
δys δxw s δxe

S
Fig. 4.6: Disposition du (v. ) 2D.

R h i
∂
∂x
(Γ ∂φ
∂x
) + ∂
∂y
(Γ ∂φ
∂y
) + S dv = 0
(vc)
- 93 -
∂
(Γ ∂φ ∂
(Γ ∂φ
R R R
=⇒ ∂x ∂x
) dv + ∂y ∂y
) dv + S dv = 0
(vc) (vc) (vc)
Re Rn
=⇒ d
dx
(Γ dφ
dx
) A dx + d
dy
(Γ dφ
dy
) A dy + S̄ ∆v = 0
w s
e n
=⇒ Γ A dφ + Γ A dφ + S̄ δx δy δz = 0
dx w dy
s

dφ dφ dφ dφ

=⇒ Γe Ae dx e
− Γw Aw dx w
+ Γn An dy
− Γs As dy
+ S̄ δx δy δz = 0
n s
Γe Ae (φE −φP ) Γw Aw (φP −φW ) Γn An (φN −φP ) Γs As (φP −φS )

=⇒ δxe
− δxw
+ δyn
− δys
+ Su + SP φ P = 0
et en réarrangeant les termes, nous aurons :

Γw Aw Γe Ae Γs As Γn An Γw Aw Γe Ae
+ + + − SP φP = φW + φE · · ·
δxw δxe δys δyn δxw δxe

Γs As Γn An
···+ φS + φ N + Su (4.16)
δys δyn
Cette équation peut aussi se mettre sous la forme :
aP φP = aW φW + aE φE + aS φS + aN φN + Su (4.17)
ave :
Γw Aw Γe Ae Γs As Γn An
aW = ; aE = ; aS = ; aN = ; aP = aW +aE +aS +aN −SP
δxw δxe δys δyn
ou en ore sous forme plus ompa te :

P
 a φ = nb anb φnb + Su
 P P


(4.18)


a = P a − S

P nb nb P
ave : nb = {W, E, S, N}.
Remarques :
- Les équations que nous venons d'établir s'appliquent uniquement aux n÷uds in-
ternes. Les extrémités doivent toujours être traitées à part.
-
P P anb
Si SP = 0 alors on aura : aP = anb =⇒ aP
= 1.
nb nb
- 94 -
-
P
anb φnb
Si Su = 0 alors on aura : φP = nb
aP
. Puisque φP est la somme pondérée
des valeurs de voisins, elle est toujours bornée par es valeurs.
- Si S 6= 0, φP n'est pas né essairement bornée par les valeurs de ses voisins et peut
sortir de es limites et 'est parfaitement logique. Ce dépassement est di té par
les valeurs de Su et SP par rapport aux anb .
- Les aires des fa ettes sont données par : Ae = Aw = δy.δz et As = An = δx.δz .
δz étant l'épaisseur de la plaque.
- L'équation dis rétisée exprime l'équilibre entre les ux et le terme sour e. De e
fait, la onservation à travers haque (v.c) est garantie. Cependant, la onserva-
tion totale dans le domaine de al ul n'est pas garantie du fait que le transfert
diusif d'une maille entre dans la maille suivante. Par exemple, en é rivant
l'équilibre pour E nous devons assurer que le ux utilisé dans la fa e e est bien
J~e et qu'il soit dis rétisé exa tement omme indiqué plus haut.
- Les oe ients aP et anb sont tous de même signe. Dans notre as, ils sont tous
positifs et e i a un sens physique. Si par exemple la température en E augmente,

on s'attend à e que elle de P augmente et ne diminue pas.
- Dans l'équation linéarisée du terme sour e, nous avons imposé que SP soit négatif.
Ce i aussi a un sens physique. Si, par exemple, S est une sour e de haleur,
nous ne voudrions pas avoir la situation où quant T augmente, S augmente
indéniment. Nous ontrlons alors ette situation ave le s héma numérique en
imposant que SP soit maintenu négatif.
Exemple :
2
Soit une plaque bidimensionnelle (0, 3x0, 4 m ) d'épaisseur 1 cm et de ondu tivité
thermique 1000 W/m°C. La fa e Ouest reçoit un ux onstant égal à 500 kW/m2 et les
otés Est et Sud sont isolés (Fig.4.7). La fa e Nord est maintenue à une température
de 100°C. En utilisant un maillage uniforme ave ∆x = ∆y = 0, 1 m, al uler la
distribution de la température dans ette plaque. En négligeant le transfert thermique
selon l'axe Z , l'équation utilisée sera don l'équation (4.15) sans terme sour e :
∂ ∂T ∂ ∂T
(k )+ (k )=0
∂x ∂x ∂y ∂y
- 95 -
T0 =100 °C
y
10 11 12
q = 500 kW/m2
7 8 9
q=0
4 5 6
1 2 3
x
q=0
Fig. 4.7: Exemple 2D.
Solution :
- Noeuds internes (5 et 8) :
Aw = Ae = ∆y.e = 0, 1.0, 01 = 10−3 m2 .
As = An = ∆x.e = 0, 1.0, 01 = 10−3 m2 .
Aw = Ae = As = An = A et kw = ke = ks = kn = k .
kA 103 .10−3
aW = aE = aS = aN = ∆x
= 0.1
= 10.
P
Su = 0 ; aP = 40 =⇒ SP = 0 et aP TP = nb anb Tnb d'où les équations :

40 T = 10 T4 + 10 T6 + 10 T2 + 10 T8 (5)
5
40 T = 10 T7 + 10 T9 + 10 T5 + 10 T11 (8)
8
- N÷ud 11 :
aN = 0; aW = aE = aS = 10;
aP = 50 =⇒ SP = −20; Su = 2 k∆x
A
T0 = 2000 d'où l'équation :
50 T11 = 10 T10 + 10 T12 + 10 T8 + 2000 (11)
- N÷ud 2 :
aS = 0; aE = aW = aN = 10; aP = 30 =⇒ SP = 0; Su = 0 d'où l'équation :
30 T2 = 10 T1 + 10 T3 + 10 T5 (2)
- N÷uds (4 et 7) :
aW = 0; aE = aS = aN = 10;
- 96 -
aP = 30 =⇒ SP = 0; Su = q A = 500 d'où les équations :

30 T = 10 T1 + 10 T5 + 10 T7 + 500 (4)
4
30 T = 10 T4 + 10 T8 + 10 T10 + 500 (7)
7
- N÷uds (6 et 9) :
aE = 0; aW = aS = aN = 10;
aP = 30 =⇒ SP = 0; Su = 0 d'où les équations :

30 T = 10 T3 + 10 T5 + 10 T9 (6)
6
30 T = 10 T6 + 10 T8 + 10 T12 (9)
9
- N÷ud 3 :
aE = aS = 0; aW = aN = 10; aP = 20 =⇒ SP = 0; Su = 0
d'où l'équation :
20 T3 = 10 T2 + 10 T6 (3)
- N÷ud 10 :
aW = aN = 0; aE = aS = 10; aP = 40 =⇒ SP = −20;
Su = 2 k∆x
A
T0 + q A = 2500 d'où l'équation :
40 T10 = 10 T7 + 10 T11 + 2500 (10)
- N÷ud 1 :
aW = aS = 0; aE = aN = 10; aP = 20 =⇒ SP = 0; Su = 500 d'où
l'équation :
20 T1 = 10 T2 + 10 T4 + 500 (1)
- N÷ud 12 :
aE = aN = 0; aW = aS = 10; aP = 40 =⇒ SP = −20; Su = 2000 d'où
l'équation :
40 T12 = 10 T9 + 10 T11 + 2000 (12)
En mettant les équations dans l'ordre (1) à (12), nous aurons la forme matri ielle
- 97 -
A.T = B de e système :
     
−20 10 0 10 0 0 0 0 0 0 0 0 T1 −500
     
 10 −30 10 0 10 0 0 0
  0 0 0 0 T2   0 
     
 0 10 −20 0 0 10 0 0
  0 0 0 0 T3   0 
     
     
 10 0 0 −30 10 0 10 0
  0 0 0 0 T4   −500 
     
 0 10 0 10 −40 10 0 10 0 0 0 0   T5   0 
     
     
 0 0 10 0 10 −30 0 0 10 0 0 0   T6   0 
 . = 
 0 0 0 10 0 0 −30 10 0 10 0 0   T7   −500 
     
     
 0 0 0 0 10 0 10 −40 10 0 10 0   T8   0 
     
 0 0 0 0 0 10 0 10 −30 0 0 10   T9   0 
     
     
 0 0 0 0 0 0 10 0 0 −40 10 0   T10   −2500 
     
 0 0 0 0 0 0 0 10 0 10 −50 10   T11   −2000 
     
0 0 0 0 0 0 0 0 10 0 10 −40 T12 −2000
En utilisant le logi iel de al ul Maple, nous aurons la solution numérique :
   
T1 260, 036
   

 T2  
  227, 798 


 T3  
  212, 164 

   

 T4  
  242, 274 


 T5  
  211, 195 

   

 T6  
= 196, 529 


 T7  
  205, 591 

   

 T8  
  178, 178 


 T9  
  166, 229 

   

 T10  
  146, 322 


 T11  
  129, 696 

T12 123, 981
Remarque :
La ompréhension de e problème 2D dépend fortement de elles des problèmes 1D
étudiés dans les exemples (1 à 3) et les résultats trouvés dans es derniers peuvent
dire tement être appliqués pour e problème.
- 98 -
4.2.6 Problème tridimensionnel

Considérons l'équation de diusion stationnaire 3D en oordonnées re tangulaires
ave terme sour e :
∂ ∂φ ∂ ∂φ ∂ ∂φ
(Γ )+ (Γ ) + (Γ )+S =0 (4.19)
∂x ∂x ∂y ∂y ∂z ∂z
L'intégration de ette équation à travers le (v.c) englobant le n÷uds interne P ne
pose au une di ulté par rapport au problème 2D et nous aurons alors :

Γe Ae dφ Aw dφ An dφ As dφ dφ dφ

− Γw + Γn
− Γs

+ Γt At − Γb Ab + S̄ δx δy δz = 0
dx e dx w dy dy dz t dz b
n s
Les indi es b et t étant su essivement pour Bottom (bas) et Top (Haut).
En approximant les dérivées et en réarrangeant les termes, nous aurons :

Γw Aw Γe Ae Γs As Γn An Γb Ab Γt At Γw Aw Γe Ae
+ + + + + − SP φP = φW + φE · · ·
δxw δxe δys δyn δzb δzt δxw δxe

Γs As Γn An Γb Ab Γt At
···+ φS + φN + φB + φ T + Su (4.20)
δys δyn δzb δzt
Cette équation peut aussi se mettre sous la forme :
aP φP = aW φW + aE φE + aS φS + aN φN + aB φB + aT φT + Su (4.21)
ave :
Γw Aw Γe Ae Γs As Γn An Γb Ab Γt At
aW = ; aE = ; aS = ; aN = ; aB = ; aT = ;
δxw δxe δys δyn δzb δzt
aP = aW + aE + aS + aN + aB + aT − SP
et sous forme ompa te :

P
 a φ = nb anb φnb + Su
 P P


(4.22)


a = P a − S

P nb nb P
ave : nb = {W, E, S, N, B, T }.
- 99 -
4.2.7 Problèmes symétriques
Pas en ore implémenté !
4.2.8 Règles de base à satisfaire
La MVF est basée sur les règles suivantes :
Règle N°1 : Règle de l'uni ité des ux aux interfa es :
Les ux sortant d'un (v.c) par une fa e doit être le même que elui entrant par la
fa e du (v.c) voisin.
Règle N°2 : Règle des oe ients positifs :
Les oe ients aP et anb doivent toujours être positifs.
Règle N°3 : Règle de la pente négative :
Pour le terme sour e, nous devons toujours avoir SP ≤ 0 .
Règle N°4 : Règle de la somme des oe ients :

P
Nous avons dans le as d'absen e du terme sour e : aP = anb
nb
P
et d'une manière générale, lorsqu'il y a un terme sour e : aP ≥ anb .
nb
- 100 -
4.3. Problème de onve tion-diusion.
4.3 Problème de onve tion-diusion

Considérons le problème ondu to- onve tif unidimensionnel (Fig.4.8) dont l'équa-
tion peut être déduite de l'équation générale de transport (4.8) :
δx
uw ue
PSfrag repla ements
W P E
w e x
δxw δxe
Fig. 4.8: Problème ondu to- onve tif 1D.

d d dφ
(ρ u φ) = Γ (4.23)
dx dx dx
Puisqu'il y a mouvement alors l'équation de ontinuité est toujours satisfaite :
d
(ρ u) = 0 (4.24)
dx
L'intégration des équations (4.23) et (4.24) à travers le (v.c) nous donne :

= Γ A dφ − Γ A dφ



 (ρ u A φ)e − (ρ u A φ)w dx e dx w



(ρ u A) − (ρ u A)

=0
e w
En supposant un prol linéaire pour le terme diusif nous aurons les équations
simpliées suivantes :

Fe φe − Fw φw = De (φE − φP ) − Dw (φP − φW )



(4.25)


F − F

=0
e w
ave :
Γw Aw Γe Ae
Fw = ρw uw Aw ; Fe = ρe ue Ae ; Dw = ; De =
δxw δxe
F = ρuA : orrespond au débit massique et peut être positif ou négatif.
ΓA k
D= δx
: orrespond à la diusion de haleur ou de matière (Γ = Cp
, µ, . . .) et doit
être positif.
- 101 -
Le problème posé i i est la détermination des variables φ aux niveaux des interfa es,
'est-à-dire : φw et φe ?
Nous allons donner les diérents s hémas utilisés dans la littérature an d'approxi-
mer le terme onve tif aux interfa es du (v.c).
4.3.1 S héma entré

L'approximation du terme onve tif par un s héma entré est donnée par :

1
 φ = (φW + φP )
 w 2


(4.26)


φ 1
= (φP + φE )

e 2
En remplaçant e s héma dans (4.25), nous aurons :
Fe Fw
2
(φP + φE ) − 2
(φW + φP ) = De (φE − φP ) − Dw (φP − φW )
Fw Fe Fw Fe

=⇒ Dw − 2
+ De + 2
φ P = Dw + 2
φ W + De − 2
φE
Fw Fe Fw Fe

=⇒ Dw + 2
+ De − 2
+ (Fe − Fw ) φP = Dw + 2
φ W + De − 2
φE
qui toujours s'é rit sous la forme :
X
aP φP = anb φnb + Su = aW φW + aE φE + Su
nb
ave :
Fw Fe
aW = Dw + ; aE = De − ; aP = aW +aE +(Fe − Fw ) ; Sp = Fw −Fe ; Su = 0
2 2
In onvénients :
|F |
Ce s héma est limité à des faibles valeurs du rapport (sinon on perd la dominan e
D
de la diagonale de la matri e) :
|F | ρ u δx
= Pe = < 2
D Γ
|F |
représente soit le nombre de Pe let en transfert thermique soit le nombre de
D
Reynolds en é oulement de uide. Il représente don l'importan e de la onve tion par
rapport à la diusion.
- 102 -
Pour montrer et in onvénient, examinons l'exemple simple suivant :
Supposons que Fw = Fe = 4 et Dw = De = 1 , e qui nous donne un nombre de
Pe let égal à 4 et les oe ients aW = 3, aE = −1 et aP = 2.

3.100−200
Si φW = 100 et φE = 200 alors nous aurons : φP = 2
= 50 !
3.200−100
Si φW = 200 et φE = 100 alors nous aurons : φP = 2
= 250 !
Ces valeurs sont en dehors de l'intervalle dans lequel la solution devrait y être -à-d
[100, 200]. Dans et exemple, la règle des oe ients positifs est violée.
4.3.2 S héma dé entré (Upwind)

Ce s héma propose que la valeur à l'interfa e soit prise égale à elle de l'amont :
 
 φ = φW si Fw > 0  φ = φP si Fe > 0
 w  e

 

et (4.27)

 

φ

= φP si Fw < 0 φ

= φE si Fe < 0
w e
Soit [a, b] = max (a, b) le s héma Upwind se résume alors à :

 F φ = [Fw , 0] φW − [−Fw , 0] φP
 w w


(4.28)


F φ

= [Fe , 0] φP − [−Fe , 0] φE
e e
En remplaçant e s héma dans (4.25), l'équation dis rétisée s'é rit alors :
aP φP = aW φW + aE φE + Su
ave :
aW = Dw +[Fw , 0] ; aE = De +[−Fe , 0] ; aP = Dw +[−Fw , 0]+De +[Fe , 0] ; Sp = 0; Su = 0
In onvénients :
|F |
Ce s héma est non approprié pour les faibles valeurs du rapport
D
= |Pe |.
|F | dφ
La diusion est surestimée pour les grandes valeurs de
D
(normalement
dx
→ 0) .
- 103 -
4.3.3 S héma exponentiel

L'intégration de l'équation (4.23) pour 0≤x≤L ave les onditions aux limites
φ0 = φ(x = 0) et φL = φ(x = L) nous donne :
F
φ − φ0 exp( D x) − 1 exp(Pe Lx ) − 1
= F
= (4.29)
φL − φ0 exp( D )−1 exp(Pe ) − 1
ρuL
ave : Pe = Γ
F
Nous pouvons alors représenter les variations de φ(x) pour diérentes valeurs de D
.
F
- Pour
D
= 0 =⇒ u = 0 don pas de terme onve tif (problème de diusion pur),
d2 φ
=⇒ dx2
= 0 =⇒ φ a un prol linéaire.
F
- Pour les valeurs élevées et positives de , les valeurs de φ sont pro hes de φ0 .
D
F
- Pour les valeurs élevées et négatives de , les valeurs de φ sont pro hes de φL .
D
L'intégration de l'équation (4.23) sur le (vc) nous donne :

φP − φE φW − φP
Fe φP + − Fw φW + =0 (4.30)
exp(Pee ) − 1 exp(Pew ) − 1
L'équation dis rétisée s'é rit alors :
ave :
Fw
Fw exp( D ) Fe
aW = Fw
w
; aE = Fe
; aP = aW + aE + (Fe − Fw )
exp( D w
) −1 exp( D e
) − 1
In onvénients :
dφ L
- Lorsque |Pe | prend des valeurs élevées alors
dx
=0 à x= 2
don pas de diusion,
alors qu'ave un s héma Upwind, un prol linéaire est utilisé pour la diusion. Don
la diusion est surestimée.
- Le s héma paraît approprié mais n'est pas assez utilisé ar d'une part, le al ul
des exponentielles revient très her numériquement et, d'autre part, le s héma n'est
pas exa t pour les situations 2D et 3D ave terme sour e.
4.3.4 S héma Hybride

C'est un s héma ayant la qualité du s héma exponentiel mais ave des expressions
simples dont le al ul n'est pas her (expressions algébriques).
- Ce s héma est identique au s héma entré si : |Pe | ≤ 2.

- Ce s héma est identique au s héma dé entré (Upwind) si : |Pe | > 2.
- 104 -
In onvénients :
L'erreur est maximale au voisinage de |Pe | ≃ 2.
4.3.5 S héma Power Law

La loi en puissan e ou polynomiale améliore l'enveloppe de la variation exponen-
tielle. Ce s héma est re ommandé pour les situations de onve tion-diusion.
4.3.6 Forme générale

On ombine les diérents s hémas en une seule forme ompa te :
ave :
aW = Dw A (|Pew |) + [Fw , 0] ; aE = De A (|Pee |) + [−Fe , 0] ; aP = aW + aE − SP
où A (|Pe |) qui prendra une fon tion selon le s héma hoisi :
S héma A (|Pe |)
Centré 1 − 0.5 |Pe |
Upwind 1
Hybride [0, 1 − 0.5 |Pe |]
|Pe |
Exponentiel
exp(|Pe |)−1
0, (1 − 0.1 |Pe |)5

Power Law
Tab. 4.1: Fon tion à prendre pour A (|Pe |).
4.3.7 Exemple 1 : S héma entré

La propriété φ est transportée par onve tion-diusion à travers le domaine 1D
(Fig.4.9). L'équation gouvernante est l'équation (4.23) ave les onditions aux limites
φ0 = φ(x = 0) = 1 et φL = φ(x = L) = 0. En utilisant 5 (v.c) égaux et le s héma
entré pour la onve tion et la diusion al uler la distribution de φ(x) dans les as
suivants :
a)- u = 0, 1 m/s ; b)- u = 2, 5 m/s ;

et omparer le résultat ave la solution analytique :
φ − φ0 exp(ρ u x/Γ) − 1
=
φL − φ0 exp(ρ u x/Γ) − 1
- 105 -
)- Re al uler la solution pour u = 2, 5 m/s ave 20 (v.c) et omparer le résultat
ave la solution analytique.
On donne : L = 1 m ; ρ = 1, 0 kg/m3 ; Γ = 0, 1 kg/m.s ; ∆x = 0, 2 m ; A = 1 m2 ;
δx
u u
PSfrag repla ements
φ0 = 1 1 2 3 4 5 φL = 0
W P E
δx/2 δx δx δx δx δx/2
Fig. 4.9: Exemple 1 : S héma entré.
Solution :
1er as :
Fw = Fe = F = ρ u A = 1.0, 1.1 = 0, 1.
ΓA 0,1.1 F
Dw = De = D = ∆x
= 0,2
= 0, 5. Don :
D
= 0, 2.
- N÷uds internes (2, 3 et 4) :
aW = Dw + F2w = D + F2 = 0, 5 + 0, 05 = 0, 55.
aE = De − F2e = D − F2 = 0, 5 − 0, 05 = 0, 45.
F F
aP = D + 2
+D− 2
= 2 D = 2 . 0, 5 = 1 =⇒ SP = 0, 55 + 0, 45 − 1 = 0.
Su = 0 .
En appliquant : aP φP = aW φW + aE φE + Su , nous aurons les équations :

φ = 0, 55 φ1 + 0, 45 φ3 (2)
 2


φ3 = 0, 55 φ2 + 0, 45 φ4 (3)


φ = 0, 55 φ + 0, 45 φ

(4)
4 3 5
- N÷ud gau he (1) :
Fe φe − Fw φw = De (φE − φP ) − Dw (φP − φW )
φP +φE

=⇒ F 2
− F φ0 = D (φE − φP ) − 2 D (φP − φ0 )
F F

=⇒ 3 D + 2
φP = 0 φW + D − 2
φE + (2 D + F ) φ0
F
+ (2 D + F ) φP = 0 φW + D − F2 φE + (2 D + F ) φ0

=⇒ D− 2
- 106 -
qui s'é rit sous la forme :
aP φP = aW φW + aE φE + Su et aP = aW + aE − SP
ave :
F F
aW = 0; aE = D − ; aP = 3 D + ; SP = − (2 D + F ) ; Su = (2 D + F ) φ0
2 2
A.N : aW = 0; aE = 0, 45; aP = 1, 55; SP = −1, 1; Su = 1, 1.

D'où l'équation :
1, 55 φ1 = 0, 45 φ2 + 1, 1 (1)
- N÷ud droit (5) :
φW +φP

=⇒ F φL − F 2
= 2 D (φL − φP ) − D (φP − φW )
F F

=⇒ 3 D − 2
φP = D + 2
φW + 0 φE + (2 D − F ) φL
F
+ (2 D − F ) φP = D + F2 φW + 0 φE + (2 D − F ) φL

=⇒ D+ 2
qui s'é rit sous la forme :
aP φP = aW φW + aE φE + Su et aP = aW + aE − SP
ave :
F F
aW = D + ; aE = 0; aP = 3 D − ; SP = − (2 D − F ) ; Su = (2 D − F ) φL
2 2
A.N : aW = 0, 55; aE = 0; aP = 1, 45; SP = −0, 9; Su = 0 .

D'où l'équation :
1, 45 φ5 = 0, 55 φ4 (5)
En mettant les équations (1) à (5) sous la forme matri ielle A.T = B , nous aurons
- 107 -
le système :
     
−1, 55 0, 45 0 0 0 φ1 −1, 10
 0, 55 −1, 00 0, 45 0 0   φ2   0
     

     
 0 0, 55 −1, 00 0, 45 0  .  φ3  =  0 
     
0 0 0, 55 −1, 00 0, 45   φ4   0
     
 
0 0 0 0, 55 −1, 45 φ5 0
   
φ1 0, 9421
 φ2   0, 8006
   

   
 φ3  =  0, 6276 
   
 φ4   0, 4163
   

φ5 0, 1579
- Comparaison :
La omparaison entre les résultats numériques et exa tes est mentionnée dans le
tableau et la gure i-dessous :
Noeud Abs isse (m) Solution Numérique Solution Exa te Erreur relative (%)
1 0,1 0,9421 0,9389 -0,341
2 0,3 0,8006 0,7963 -0,540
3 0,5 0,6276 0,6222 -0,868
4 0,7 0,4163 0,4098 -1,586
5 0,9 0,1579 0,1502 -5,126
Fig. 4.10: Exemple1- Cas 1 : u=0,1 m/s et ∆x=0,2 m.
- 108 -
2ème as :
F 2,5
Dans e as :
D
= 0,5
= 5.
De la même manière que le 1er as, nous aurons le système suivant :
     
−2, 75 −0, 75 0 0 0 φ1 −3, 50
 1, 75 −1, 00 −0, 75 0 0   φ2   0
     

     
 0 1, 75 −1, 00 −0, 75 0  .  φ3  =  0 
     
0 0 1, 75 −1, 00 −0, 75   φ4   0
     
 
0 0 0 1, 75 −0, 25 φ5 0
La solution numérique est donnée par le logi iel de al ul Maple :
   
φ1 1, 0356
 φ2   0, 8694
   

   
 φ3  =  1, 2573 
   
 4   0, 3521
φ
   

φ5 2, 4644
- Comparaison :
1 0,1 1,036 1,0000 -3,60
2 0,3 0,8694 1,0000 13,60
3 0,5 1,2570 1,0000 -25,70
4 0,7 0,3521 0,9994 64,77
5 0,9 2,2464 0,9179 -168,40
- 109 -
On remarque bien, d'après le tableau et la gure i-dessus, que la solution numérique
ne orrespond pas du tout à la solution analytique à ause du non respe t de la ondition

|F |
D
< 2. On dit que la solution diverge.
3ème as :
F 2,5
Dans e as : ∆x = 0, 05 et
D
= 2
= 1, 25.
De la même manière que le 2ème as, nous aurons les résultats suivants :
Noeuds aW aE SP aP Su
1 0 0, 75 −6, 50 7, 25 6, 50
2 − 19 3, 25 0, 75 0 4, 00 0
20 3, 25 0 −1, 50 4, 75 0
- Comparaison :
- 110 -
1 0,025 1,0000 1,0000 0
2-11 1,0000 1,0000 0
12 0,575 1,0000 1,0000 0
13 0,625 1,0000 0,9999 -0,01
14 0,675 0,9999 0,9997 -0,02
15 0,725 0,9998 0,9990 -0,08
16 0,775 0,9989 0,9964 -0,25
17 0,825 0,9954 0,9875 -0,80
18 0,875 0,9800 0,9558 -2,53
19 0,925 0,9135 0,8474 -7,80
20 0,975 0,6250 0,4621 -35,25
F
Le ranement du maillage de 5 à 20 (v.c) a permis de réduire le rapport
D
de 5
ème
(dans le 2 as) à 1, 25 et don on onstate bien que la solution onverge ave une
erreur relative a eptable sauf pour le dernier n÷ud. Ce dernier problème pourra être
réglé en augmentant susamment le nombre de n÷uds.
4.3.8 Exemple 2 : S héma Upwind

La variable φ est transportée par onve tion-diusion à travers le domaine 1D
(Fig.4.9) de longueur L. ave les onditions aux limites φ0 = φ(x = 0) = 1 et
φL = φ(x = L) = 0. En utilisant 5 (v.c) égaux et le s héma Upwind pour le terme
onve tif, al uler la distribution de φ(x) dans les deux as suivants.
- 111 -
ΓA
1. Fw = Fe = F = ρ u A = 0, 1 et Dw = De = D = δx
= 0, 5.
ΓA
2. Fw = Fe = F = ρ u A = 2, 5 et Dw = De = D = δx
= 0, 5.
Solution :
Puisque F >0 alors le s héma Upwind nous donne : φw = φW et φe = φP .
1er as : F = 0, 1
- N÷uds internes (2, 3 et 4) :
D'après l'équation (4.25) nous auvons :
F φe − F φw = D (φE − φP ) − D (φP − φW )
et en remplaçant φw et φe nous aurons :
F φp − F φW = D (φE − φP ) − D (φP − φW )
=⇒ (2 D + F ) φP = (D + F )φW + D φE
En omparant ave : aP φP = aW φW + aE φE + Su , nous aurons les oe ients :
aW = D + F = 0, 5 + 0, 1 = 0, 6.
aE = D = 0, 5. Su = 0 .
aP = 2 D + F = 2 . 0, 5 + 0, 1 = 1, 1 =⇒ SP = (D + F ) + D − (2 D + F ) = 0.
et nalement les équations :

1, 1 φ2 = 0, 6 φ1 + 0, 5 φ3 (2)



1, 1 φ3 = 0, 6 φ2 + 0, 5 φ4 (3)


1, 1 φ = 0, 6 φ + 0, 5 φ

(4)
4 3 5
- N÷ud gau he (1) :
=⇒ F (φP ) − F φ0 = D (φE − φP ) − 2 D (φP − φ0 )
=⇒ (3 D + F ) φP = 0 φW + D φE + (2 D + F ) φ0
- 112 -
aW = 0.
aE = D = 0, 5. Su = (2 D + F ) φ0 = 1, 1.
aP = 3 D + F = 1, 6 =⇒ SP = D − (3 D + F ) = −(2 D + F ) = −1, 1.
et nalement l'équation :
1, 6 φ1 = 0, 5 φ2 + 1, 1 (1)
- N÷ud droit (5) :
=⇒ F (φL ) − F φW = 2 D (φL − φP ) − D (φP − φW )
=⇒ 3 DφP = (D + F )φW + 0 φE + (2 D − F ) φL
aW = D + F = 0, 6.
aE = 0. Su = (2 D − F ) φL = 0.
aP = 3 D = 1, 5 =⇒ SP = (D + F ) + 0 − 3 D = −(2 D − F ) = −0, 9.
et nalement l'équation :
1, 5 φ5 = 0, 6 φ4 (5)
En mettant les équations (1) à (5) sous la forme matri ielle A.T = B , nous aurons
le système :
     
−1, 6 0, 5 0 0 0 φ1 −1, 10
 0, 6 −1, 1 0, 5 0 0   φ2   0
     

     
 0 0, 6 −1, 1 0, 5 0  .  φ3  =  0 
     
 0 0 0, 6 −1, 1 0, 5 φ 0
     
  4   
0 0 0 0, 6 −1, 5 φ5 0
   
φ1 0, 9348
 φ2   0, 7914
   

   
 φ3  =  0, 6194 
   
 φ4   0, 4129
   

φ5 0, 1652
- 113 -
- Comparaison :
N÷ud Abs isse (m) Solution Numérique Solution Exa te Erreur relative (%)
1 0,1 0.9348 0.9389 -0,4367
2 0,3 0.7914 0.7963 -0,6153
3 0,5 0.6194 0.6222 -0,4500
4 0,7 0.4129 0.4098 0,7565
5 0,9 0.1652 0.1502 9,9870
2ème as : F = 2, 5
De la même manière que le 1er as, nous aurons les résultats suivants :
N÷ud aW aE SP aP Su
1 0,0 0,5 -3,5 4,0 3,5
2 3,0 0,5 0,0 3,5 0,0
3 3,0 0,5 0,0 3,5 0,0
4 3,0 0,5 0,0 3,5 0,0
5 3,0 0,0 1,5 1,5 0,0
- Comparaison :
N÷ud Abs isse (m) Solution Numérique Solution Exa te Erreur relative (%)
1 0,1 1,0000
2 0,3 1,0000
3 0,5 1,0000
4 0,7 0,9994
5 0,9 0,9179
- Remarque :
Nous remarquons que pour le même nombre de n÷uds, le s héma entré n'est pas
apable de produire une solution pour le même problème alors que le s héma Upwind
donne une solution plus réelle malgré qu'elle reste éloignée de la solution exa te surtout
vers la limite droite.
- 114 -
4.4. Equation de diusion 1D instationnaire.
4.4 Equation de diusion 1D instationnaire

L'équation de diusion instationnaire unidimensionnelle (ou équation de la haleur)
est donnée par :

∂T ∂ ∂T
ρ Cp = (k )+S (4.31)
∂t ∂x ∂x
L'intégration de ette équation à travers le (v.c) et dans l'intervalle de temps
[t, t + ∆t] donne :
t+∆t t+∆t t+∆t

∂T ∂
(k ∂T
R R R R R R
ρ Cp ∂t
dv dt = ∂x ∂x
) dv dt + S dv dt
t (vc) t (vc) t (vc)
Re
t+∆t t+∆t t+∆t
ρ Cp dT (k A dT dT
R R R
dt
dt dv = dx
) e − (k A dx
) w dt + S ∆v dt
w t t t
t+∆t
R h i t+∆t
(TE −TP ) (TP −TW )
ρ Cp (TP − TP0 ) ∆v =
R
(ke Ae δxe
− kw Aw δxw
dt + S ∆v dt
t t
ave la notation TP0 qui orrespond au temps t et TP au temps t + ∆t.

En posant :
t+∆t
Z
TP dt = θ TP + (1 − θ) TP0 ∆t

It = (4.32)
nous aurons alors les diérents s hémas temporels de dis rétisation :
s héma expli ite pour θ = 0,

s héma omplètement impli ite pour θ = 1,
1
s héma impli ite de Crank-Ni holson pour θ = .
2
(TP −TP0 )
h i
ke kw
ρ Cp ∆t
δx = θ ( δx e
(TE − TP ) − δxw
(TP − TW )
h i
ke kw
+(1 − θ) ( δx e
(TE0 − TP0 ) − δxw
(TP0 − TW
0
) + S δx

δx kw ke kw 0 ke
ρ Cp ∆t
+θ δxw
+ δxe
TP = δxw
[θ TW + (1 − θ) TW ]+ δxe
[θ TE + (1 − θ) TE0 ]
h i
δx kw ke
+ ρ Cp ∆t
− (1 − θ) δx w
− (1 − θ) δxe
TP0 + S δx
qui nalement peut se mettre sous la forme :
0
+aE θ TE + (1 − θ) TE0 + a0P − (1 − θ) aW − (1 − θ) aE TP0 +b

aP TP = aW θ TW + (1 − θ) TW
ave :
- 115 -
δx kw ke
aP = θ (aW + aE ) + a0P ; a0P = ρ Cp ; aW = ; aE = ; b = S δx
∆t δxw δxe
4.4.1 S héma expli ite

Pour e s héma nous avons : θ=0 et en posant : b = Su + SP TP0 nous aurons :
0
+ aE TE0 + a0P − (aW + aE − SP ) TP0 + Su

aP TP = aW TW
où :
δx kw ke
aP = a0P ; a0P = ρ Cp ; aW = ; aE =
∆t δxw δxe
4.4.2 S héma impli ite

Pour e s héma nous avons : θ=1 et don nous aurons :
aP TP = aW TW + aE TE + a0P TP0 + Su
où :
δx kw ke
aP = a0P + aW + aW − SP ; a0P = ρ Cp ; aW = ; aE =
∆t δxw δxe
4.4.3 S héma de Crank-Ni holson

1
Pour e s héma nous avons : θ= 2
et don nous aurons :
TE − TE0 0
TW − TW h
0 aW aE i 0
aP TP = aW + aE + aP − − TP + b
2 2 2 2
où :
1 1 δx kw ke 1
aP = (aW + aW )+a0P − SP ; a0P = ρ Cp ; aW = ; aE = ; b = Su + SP TP0
2 2 ∆t δxw δxe 2
- 116 -
4.5. Exer i es.
4.5 Exer i es
4-01 : Soit un volume de ontrle (v.c) de très faibles dimensions ∆x, ∆y et ∆z . En
onsidérant que l'a umulation de la variable s alaire φ dans le (v.c) pendant un temps
∆t est égale au ux net de φ dans le (v.c) augmentée de la génération de φ à l'intérieur
du (v.c), retrouver l'équation générale de transport s alaire de la variable φ.

1. Identier haque terme de ette équation.
2. Simplier ette équation pour le as d'un problème ondu to- onve tif unidimen-
sionnel, stationnaire et sans terme sour e.
3. En utilisant la MVF, intégrer ette équation à travers un (v.c) onvenablement
hoisi (faire un s héma).
4. En utilisant les s hémas linéaire et entré su essivement pour approximer les
termes diusif et onve tif, é rire la forme générale de dis rétisation en expli i-
tant tous les oe ients (anb , ap , Sp et Su ).

5. Que représente physiquement le nombre de Pe let.
6. Quelle est l'in onvénient du s héma entré utilisé pour approximer le terme
onve tif ? Justier votre réponse et donner un exemple numérique.
4-02 : A- Soit un domaine unidimensionnel onstitué d'une barre horizontale de lon-
gueur L soumise à une température onstante TA à son bout gau he et TB à son bout
droit. Cette barre est onstituée de deux tronçons (L1 , A1 , k1 ) et (L2 , A2 , k2 ). En
dis rétisant l'équation de diusion par volumes nis en 7 (v.c) égaux, déterminer la
forme matri ielle du problème. On donne :

L1 = 0, 7 m ; L2 = 0, 7 m ; TA = 300 °C ; TB = 600 °C ;
k1 = 100 W/m°C ; k2 = 200 W/m°C ; A1 = 0, 01 m2 ; A2 = 0, 02 m2.
B- Si TB = TA , A2 = A1 et k2 = k1 , déterminer la nouvelle forme matri ielle du
problème et al uler les températures aux diérents n÷uds.
4-03 : Considérons le transfert ondu tif stationnaire et bidimensionnel dans une
plaque arrée sans terme sour e. Cette plaque de té L a une épaisseur e et une
ondu tivité thermique k. Les fa es Nord et Ouest sont maintenues à une température
onstante T0 et les fa es Est et Sud sont soumises à un ux de haleur onstant q.
En utilisant une dis rétisation par volumes nis sur un maillage re tangulaire (Nx x
Ny ), déterminer les équations dis rétisées des n÷uds suivants : 1, 2, 4, 5, 6, 8, 13, 14
et 16. On donne :
L = 1, 6 m ; e = 2 cm ; T0 = 100 °C ; q = 105 W/m2 ;

δx = 0, 4 m ; δy = 0, 4 m ; k = 500 W/m°C .
- 117 -
plaque re tangulaire sans terme sour e. Cette plaque de largeur L et de hauteur H, a
une épaisseur e et une ondu tivité thermique k1 dans la première moitié de sa largeur
et k1 dans la se onde moitié. La fa e Nord est soumise à un ux de haleur onstant
q, la fa e Ouest est maintenue à une température onstante T0 et les fa es Est et Sud
sont isolées.
Ny ), déterminer les équations dis rétisées des n÷uds suivants : 1, 2, 3, 4, 5, 6, 7, 8,
17, 18, 19 et 20. On donne :
L = 1, 6 m ; e = 2 cm ; T0 = 40 °C ; q = 100 kW/m2 ; H = 1, 0 m ;
δx = 0, 4 m ; δy = 0, 2 m ; k1 = 500 W/m°C ; k2 = 200 W/m°C .
plaque re tangulaire sans terme sour e. Cette plaque de largeur L et de hauteur H, a
une épaisseur e et une ondu tivité thermique k1 dans la première moitié de sa largeur
et k1 dans la se onde moitié. La fa e Nord est soumise à un ux de haleur onstant
q, la fa e Est est maintenue à une température onstante T0 et les fa es Ouest et Sud
sont isolées.
Ny ), déterminer les équations dis rétisées des n÷uds suivants : 1, 2, 3, 4, 5, 6, 7, 8,
17, 18, 19 et 20. On donne :
L = 1, 6 m ; e = 2 cm ; T0 = 40 °C ; q = 100 kW/m2 ; H = 1, 0 m ;
δx = 0, 4 m ; δy = 0, 2 m ; k1 = 500 W/m°C ; k2 = 200 W/m°C .
plaque re tangulaire ave terme sour e. Cette plaque de largeur L et de hauteur H,

a une épaisseur e et une ondu tivité thermique k. La fa e Nord est soumise à un ux
de haleur onstant q, la fa e Est est maintenue à une température onstante T0 et les
fa es Ouest et Sud sont isolées.
Ny ), déterminer les équations dis rétisées des n÷uds suivants : 1, 2, 4, 5, 6, 8, 17, 18
et 20. On donne :
L = 1, 6 m ; e = 2 cm ; H = 1, 0 m ; q = 100 kW/m2 ; T0 = 40 °C ;
3
δx = 0, 4 m ; δy = 0, 2 m ; k = 500 W/m°C ; S = 5000 W/m .
plaque re tangulaire sans terme sour e. Cette plaque de largeur L et de hauteur H,
- 118 -
4.5. Exer i es.
a une épaisseur e et une ondu tivité thermique k. Les fa es Est et Ouest sont mainte-
nues à une température onstante T0 et les fa es Sud et Nord sont soumises à un ux
de haleur onstant q.
Ny ), déterminer la forme matri ielle du problème. On donne :
L = 1, 0 m ; e = 0, 1 cm ; H = 0, 4 m ; T0 = 100 °C ;
δx = 0, 2 m ; δy = 0, 1 m ; k = 1000 W/m°C ; q = 600 kW/m2 .
4-08 : Considérons le transfert ondu tif stationnaire et unidimensionnel dans une
barre sans terme sour e. Cette barre de longueur L, de se tion A, de ondu tivité ther-
mique k est soumise à une température onstante TA à son bout droit et à un ux de
haleur onstant q à son bout gau he.
En utilisant une dis rétisation par volumes nis en 3 (v.c) égaux, déterminer la forme
matri ielle du problème et al uler les températures aux diérents n÷uds. On donne :
L = 0, 15 m ; A = 0, 01 m2 ; k = 1000 W/m°C ; TA = 50 °C ; q = 50 kW/m2 .
4-09 : La variable φ est transportée par onve tion-diusion à travers le domaine uni-
dimensionnel de longueur L ave les (C.L) suivantes :
φ(x = 0) = φ0 = 2 et φ(x = L) = φL = 1
En utilisant 5 (v.c) égaux et le s héma entré pour le terme onve tif, déterminer la
Fw = Fe = F = 0, 1 et Dw = De = D = 0, 5.
φ(x = 0) = φ0 = 1 et φ′ (x = L) = φ′L = 0
En utilisant 3 (v.c) égaux et le s héma Upwind pour le terme onve tif, déterminer la
forme matri ielle du problème et al uler les valeurs de φ aux diérents n÷uds. On
donne :
ΓA
Fw = Fe = F = ρ u A = 0, 2 et Dw = De = D = δx
= 0, 4.
φ′ (x = 0) = φ′0 = 0 et φ(x = L) = φL = 1
- 119 -
forme matri ielle du problème et al uler les valeurs de φ aux diérents n÷uds. On
donne :
ΓA
Fw = Fe = F = ρ u A = −0, 2 et Dw = De = D = δx
= 0, 4.
4-12 : Un large four industriel est supporté par une olonne en brique de forme arrée de
té L, de faible épaisseur e et d ondu tivité thermique k. Cette olonne est soumise à
une température onstante T0 sur les fa es Est, Ouest et Nord. La fa e Sud est refroidie
par un ourant d'air à une température T∞ et un oe ient de transfert onve tif h. En
onsidérant l'équation i-dessous et en dis rétisant la olonne par volumes nis en 4
(v.c) égaux, déterminer le système d'équations algébriques à résoudre pour e problème.
On donne :
L = 1m; e = 10 mm ; k = 100 W/m°K ;

T0 = 500 °K ; T∞ = 300 °K ; h = 10 W/m2°K .
∂ ∂T ∂ ∂T hP
( )+ ( ) − n2 (T − T∞ ) = 0 avec n2 =
∂x ∂x ∂y ∂y k Ac
P et Ac sont respe tivement le périmètre et l'aire de la surfa e onve tée.
4-13 : Considérons le transfert ondu tif stationnaire et unidimensionnel dans une
barre sans terme sour e. Cette barre de longueur L, de se tion A, de ondu tivité ther-
mique k est soumise à une température onstante TA à son bout gau he et TB à son
bout droit.
En utilisant une dis rétisation par volumes nis en 3 (v.c) égaux, déterminer la forme
matri ielle du problème et al uler les températures aux diérents n÷uds. On donne :
L = 0, 15 m ; A = 0, 01 m2 ; k = 1000 W/m°C ; TA = 10 °C ; TB = 50 °C .
φ(x = 0) = φ0 = 1 et φ(x = L) = φL = 0

ΓA
Fw = Fe = F = ρ u A = 0, 1 et Dw = De = D = δx
= 0, 5.
- 120 -
Bibliographie
[1℄ Joe D. Homann, Numeri al Methods for Engineers and S ientists , 2nd Ed.,
M Graw-Hill, New York, 2001.
[2℄ Chapra, S. C., and Canale, R. P., Numeri al Methods for Engineers , 3rd Ed.,
M Graw-Hill, New York, 1998.
[3℄ Aslak Tveito and Ragnar Winther, Introdu tion to Partial Dierential Equations :
A Computational , Springer-Verlag, New York, 1998.
[4℄ Y. Pin hover and J. Rubinstein, An Introdu tion to Partial Dierential Equa-
tions , Cambridge University Press, 2005.
[5℄ Ri hard hyberman, Elementary Applied Partial Dierential Equations , 2nd. Ed.,
Prenti e Hall, USA, 1983.
[6℄ Klauss A. Homann, Steve T. Chiang, Computational Fluid Dynami s , 4th Ed.,
Engineering Edu ation System, USA, 2000.
[7℄ William F. Ames, Numeri al Methods for Partial Dierential Equations , 2nd
Ed., A ademi Press, USA, 1977.
[8℄ H. F. Weinberger, A First Cours in Partial Dierential Equations with Complexe
Variables and Transform Methods , Dover Publi ations, New York, 1965.
[9℄ David L. Powers, Boundary Value Problems and Partial Dierential Equations ,
5th Ed., Elsevier A ademi Press, USA, 2006.
[10℄ David R. Croft, David G. Lilley, Heat transfer al ulations using Finite Dieren e
Equations , Applied S ien e Publishers, England, 1977.
[11℄ Henry J. Ri ardo, A Modern Introdu tion to Dierential Equations , 2nd Ed.,
Elsevier, Canada, 2009.
[12℄ Jayathi Y. Murthy, Numeri al Methods in Heat, Mass, and Momentum Transfer ,
Purdue University, 2002.
[13℄ K. W. Morton, David Mayers, Numeri al Solution of Partial Dierential Equa-
tions , 2nd Ed., Cambridge, 2005.
121
BIBLIOGRAPHIE
[14℄ Peter J. Collins, Dierential and Integral Equations , Oxford University Press,
New York, 2006.
[15℄ C. Pozrikidis, Fluid Dynami s : Theory, Computation, and Numeri al Simula-
tion , 2nd Ed., Springer, New York, 2009.
[16℄ Aslak Tveito, Ragnar Winther, Introdu tion to Partial Dierential Equations : A
Computational Approa h , Springer-Verlag, New York, 1998.
[17℄ John C. Tannehill, Dale A. Anderson, Ri hard H. Plet her, Computational Fluid
Me hani s And Heat Transfer , 2nd Ed., Taylor&Fran is, USA, 1997.
[18℄ H. K. Versteeg, W. Malalasekera, An Introdu tion to Computational Fluid Dyna-
mi s. The Finite Volume Method , 1st Ed., Longman s ienti &Te hni al, USA,
1995.
[19℄ Suhas V. Patankar, Numeri al Heat Transfer and Fluid Flow , Hemisphere Pu-
blishing, USA, 1980.
- 122 -
Annexes
123
Annexe A
Formulation à 9 points
Démonstration de la formulation à 9 points utilisée pour l'équation de Lapla e.
h2 h3
f (x + h) = f (x) + h f ′(x) + 2
f ′′ (x) + 6
f ′′′ (x) + O(h4 )
h2 h3
f (x − h) = f (x) − h f ′ (x) + 2
f ′′ (x) − 6
f ′′′ (x) + O(h4 )
h3
f (x + h) − f (x − h) = 2 h f ′(x) + 3
f ′′′ (x) + O(h5 )
f (x+h)−f (x−h) h2
=⇒ 2h
= f ′ (x) + 6
f ′′′ (x) + O(h4 ) (*)
h2 h3
f ′ (x + h) = f ′ (x) + h f ′′ (x) + 2
f ′′′ (x) + 6
f (4) (x) + O(h4 )
h2 h3
f ′ (x − h) = f ′ (x) − h f ′′(x) + 2
f ′′′ (x) − 6
f (4) (x) + O(h4 )
f ′ (x + h) + f ′ (x − h) = 2 f ′ (x) + h2 f ′′′ (x) + O(h4 )
=⇒ f ′ (x + h) − 2 f ′ (x) + f ′ (x − h) = h2 f ′′′ (x) + O(h4)
1 h2
=⇒ 6
[f ′ (x + h) − 2 f ′ (x) + f ′ (x − h)] = 6
f ′′′ (x) + O(h4 )
2
=⇒ f ′ (x)+ 61 [f ′ (x + h) − 2 f ′ (x) + f ′ (x − h)] = f ′ (x)+ h6 f ′′′ (x)+O(h4 ) (**)
f (x+h)−f (x−h)
(*) et (**) =⇒ f ′ (x) + 61 [f ′ (x + h) − 2 f ′(x) + f ′ (x − h)] = 2h
+ O(h4 )
En utilisant δ omme opérateur diéren e du 1er ordre et δ2 elui du 2nd ordre,
ette équation devient :
125
A. FORMULATION À 9 POINTS.
δf (x)
f ′ (x) + 61 δ 2 f ′ (x) = 2h
+ O(h4) et par suite :
δf (x)
f ′ (x) = 4
2 + O(h )
2 h 1 + δ6
Cette dernière équation représente la diéren e entrée impli ite à 3 points au 4ème
ordre pour f ′ (x).
De la même manière nous pouvons aussi montrer que :
δ 2 f (x)
f ′′ (x) = δ2
+ O(h4 )
h 1 + 12
En utilisant ette dernière formule, nous avons :

∂2T δx2 Ti,j ∂2T δy2 Ti,j
= et =

∂x2 ∂y 2 2

2
δx
δy
i,j ∆x2 1+ 12 i,j ∆y 2 1+ 12
et l'équation de Lapla e s'é rit don :
δy2 Ti,j δy2

δx2 Ti,j δx2

2
δx
+ 2
δy
=0 et en multipliant par : 1+ 12
1+ 12
∆x2 1+ 12 ∆y 2 1+ 12
δy2 δy2
h i
δx2 δx2
1+ 12 ∆x2
+ 1+ 12 ∆y 2
Ti,j = 0
δy2 2
h i
δx2 ∆x2 +∆y 2 δx2 δy
=⇒ ∆x2
+ ∆y 2
+ 12 ∆x2 ∆y 2
Ti,j = 0
En développant maintenant ette dernière équation :

Ti−1,j −2 Ti,j +Ti+1,j Ti,j−1 −2 Ti,j +Ti,j+1 ∆x2 +∆y 2 δx2 Ti,j−1 −2 Ti,j +Ti,j+1
∆x2
+ ∆y 2
+ 12 ∆x2 ∆y 2
=0
Le troisième terme se développe de la façon suivante :

∆x2 +∆y 2 Ti−1,j−1 −2 Ti,j−1 +Ti+1,j−1 Ti−1,j −2 Ti,j +Ti+1,j Ti−1,j+1 −2 Ti,j+1 +Ti+1,j+1
C= 12 ∆y 2 ∆x2
−2 ∆x2
+ ∆x2
12 ∆x2 ∆y 2 ∆x
En multipliant par
∆x2 +∆y 2
et en posant en posant β= ∆y
nous aurons :
12 12 β 2
1+β 2
(Ti−1,j − 2 Ti,j + Ti+1,j ) + 1+β 2
(Ti,j−1 − 2 Ti,j + Ti,j+1 ) + ∆x2 C = 0
et en réarrangeant et regroupant les mêmes termes, nous aurons nalement :
- 126 -
2 (5 − β 2 )
Ti−1,j−1 + Ti+1,j−1 + Ti−1,j+1 + Ti+1,j+1 + (Ti−1,j + Ti+1,j ) + . . .
1 + β2
2 (5β 2 − 1)
... (Ti,j−1 + Ti,j+1 ) − 20 Ti,j = 0
1 + β2
- 127 -
A. FORMULATION À 9 POINTS.
- 128 -
Annexe B
Approximation des dérivées aux interfa es
La dis rétisation par la (MVF) de l'équation de diusion impose l'estimation de la
dérivée aux limites du (v.c). Pour ela, nous pouvons utiliser l'approximation en séries
de Taylor de la manière suivante :
δxe 1 dφ δxe 2 1 d2 φ δxe 3 1 d3 φ

+ O(δx4e )

φE = φe + . .( )
2 1! dx e
+ 2
. 2! .( dx2 )e + 2
. 3! .( dx3 )e
δxe 1 dφ δxe 2 1 d2 φ δxe 3 1 d3 φ

+ O(δx4e )

φP = φe − . .( )
2 1! dx e
+ 2
. 2! .( dx2 )e − 2
. 3! .( dx3 )e
En faisant la soustra tion, nous aurons :
dφ dφ φE −φP
+ O(δx3e )

+ O(δx2e )

φE − φP = δxe dx e
=⇒ dx e
= δxe
Cette approximation est souvent appelée approximation entrée et elle est pré ise
au se ond ordre.
De même :
δxw 1 dφ δxw 2 1 d2 φ δxw 3 1 d3 φ
+ O(δx4w )

φW = φw − . .( )
2 1! dx w
+ 2
. 2! .( dx2 )w − 2
. 3! .( dx3 )w
δxw 1 dφ δxw 2 1 d2 φ δxw 3 1 d3 φ

+ O(δx4w )

φP = φw + . .( )
2 1! dx w
+ 2
. 2! .( dx2 )w + 2
. 3! .( dx3 )w
En faisant la soustra tion, nous aurons :
dφ dφ φP −φW
+ O(δx3w )

+ O(δx2w )

φP − φW = δxw dx w
=⇒ dx w
= δxw
La deuxième manière est de s hématiser le prol linéaire puis déduire du graphe la
pente de la droite :
129
B. APPROXIMATION DES DÉRIVÉES AUX INTERFACES.
α β
α β
PSfrag repla ements
W P E
w e x
δx
δxw δxe
Fig. B.1: Approximation de la dérivée aux interfa es.
dφ
φ(x) = a x + b =⇒ dx
= a = tgα
Quel que soit le prol hoisi, nous avons :
φP −φW dφ

tgα = δxw
= dx w
φE −φP dφ

tgβ = δxe
= dx e
- 130 -

Mna 28 12 2022

Transféré par

Droits d'auteur :

Formats disponibles

Mna 28 12 2022

Transféré par

Informations du document

Titre original

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Mna 28 12 2022

Transféré par

Droits d'auteur :

Formats disponibles

Edition 2022 L.

Université Mustafa Benboulaïd Batna 2.

Liste des gures v

1.3 Classi ation des EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.3.1 Cara téristiques des EDP . . . . . . . . . . . . . . . . . . . . . 4

1.4 Diérents types de onditions aux limites (C.L) . . . . . . . . . . . . . 5

2.2 Prin ipe de superposition . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.3 Méthode de séparation des variables . . . . . . . . . . . . . . . . . . . . 8

2.3.1 Equation de diusion ave onditions de Diri hlet . . . . . . . . 9

2.3.2 Equation de Lapla e en oordonnées re tangulaires . . . . . . . 12

2.3.3 Equation d'onde . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Méthodes des diéren es nies 19

3.2 Prin ipe et bases de la MDF . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 Formule de Taylor en oordonnées re tangulaires . . . . . . . . . . . . . 25

3.4 S hémas de dis rétisation . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.5 Notions de stabilité, onsistan e et onvergen e . . . . . . . . . . . . . 25

3.6 Méthodes de dis rétisation . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.6.1 Méthode expli ite . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.6.2 Méthode impli ite . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.6.3 Méthode de Crank-Ni holson . . . . . . . . . . . . . . . . . . . 28

3.7 Appli ation de la méthode expli ite à l'équation de diusion . . . . . . 29

3.8 Appli ation de la méthode impli ite à l'équation de diusion . . . . . . 31

3.9 Formulation matri ielle du s héma expli ite appliqué à l'équation de

3.10 Formulation matri ielle du s héma impli ite appliqué à l'équation de

3.11 Exemple résolu par wxMaxima . . . . . . . . . . . . . . . . . . . . . . 36

3.12.1 Cas de l'approximation entrée . . . . . . . . . . . . . . . . . . 39

3.12.2 Cas de l'approximation dé entrée d'ordre 1 . . . . . . . . . . . . 41

3.12.3 Cas de l'approximation dé entrée d'ordre 2 . . . . . . . . . . . . 43

3.13 Equations non linéaires . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.13.1 Dis rétisation par un s héma expli ite . . . . . . . . . . . . . . 45

3.13.2 Dis rétisation par un s héma impli ite . . . . . . . . . . . . . . 46

3.14 Problèmes multidimensionnels . . . . . . . . . . . . . . . . . . . . . . . 46

3.14.1 Dis rétisation par un s héma expli ite . . . . . . . . . . . . . . 46

3.14.2 Dis rétisation par un s héma impli ite . . . . . . . . . . . . . . 47

3.15 Equation de Lapla e . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.15.1 Dis rétisation par la formule à 5 points (Runge 1908) . . . . . . 48

3.15.2 Conditions limites de Diri hlet variables . . . . . . . . . . . . . 52

3.15.3 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . 53

3.15.4 Conditions aux limites de Neumann . . . . . . . . . . . . . . . . 54

3.15.5 Equation de Poisson . . . . . . . . . . . . . . . . . . . . . . . . 59

3.15.6 Méthode d'ordre élevé . . . . . . . . . . . . . . . . . . . . . . . 60

3.15.7 Cas du domaine non re tangulaire . . . . . . . . . . . . . . . . . 62

3.15.8 Equation non-linéaire et problème tridimensionnel . . . . . . . . 64

3.16 Equation d'onde . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

3.16.1 Résolution de l'équation d'onde . . . . . . . . . . . . . . . . . . 64

3.16.1.1 Dis rétisation . . . . . . . . . . . . . . . . . . . . . . . 65

3.16.1.3 Autres s hémas de dis rétisation . . . . . . . . . . . . 66

3.16.2 Méthodes de résolution . . . . . . . . . . . . . . . . . . . . . . . 67

4 Méthode des Volumes Finis 81

4.2 Equation de diusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.2.1 Problème de diusion unidimenionnel . . . . . . . . . . . . . . . 85

4.2.2 Exemple 1 : Equation de diusion 1D . . . . . . . . . . . . . . . 87

4.2.3 Exemple 2 : Equation de diusion 1D ave terme sour e . . . . 89

4.2.4 Exemple 3 : Equation de diusion 1D ave onve tion à la surfa e

4.2.5 Problème de diusion bidimensionnel . . . . . . . . . . . . . . . 93

4.2.6 Problème tridimensionnel . . . . . . . . . . . . . . . . . . . . . 99

4.2.7 Problèmes symétriques . . . . . . . . . . . . . . . . . . . . . . . 100

4.2.8 Règles de base à satisfaire . . . . . . . . . . . . . . . . . . . . . 100

4.3 Problème de onve tion-diusion . . . . . . . . . . . . . . . . . . . . . 101

Liste des gures v

1.3 Classi ation des EDP . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.4 Diérents types de onditions aux limites (C.L) . . . . . . . . . . . . . 5

2.3.1 Equation de diusion ave onditions de Diri hlet . . . . . . . . 9

3 Méthodes des diéren es nies 19

3.7 Appli ation de la méthode expli ite à l'équation de diusion . . . . . . 29

3.8 Appli ation de la méthode impli ite à l'équation de diusion . . . . . . 31

4.2 Equation de diusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.2.1 Problème de diusion unidimenionnel . . . . . . . . . . . . . . . 85

4.2.2 Exemple 1 : Equation de diusion 1D . . . . . . . . . . . . . . . 87

4.2.3 Exemple 2 : Equation de diusion 1D ave terme sour e . . . . 89

4.2.4 Exemple 3 : Equation de diusion 1D ave onve tion à la surfa e

4.2.5 Problème de diusion bidimensionnel . . . . . . . . . . . . . . . 93

4.3 Problème de onve tion-diusion . . . . . . . . . . . . . . . . . . . . . 101

4.4 Equation de diusion 1D instationnaire . . . . . . . . . . . . . . . . . . 115

2.2 Problème modié. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.3 Diérents types de s hémas de dis rétisation. . . . . . . . . . . . . . . . 25

3.6 Inuen e négative de la simple pré ision. . . . . . . . . . . . . . . . . . 23

3.7 Inuen e positive de la double pré ision. . . . . . . . . . . . . . . . . . . 24

3.8 Inuen e de l'approximation. . . . . . . . . . . . . . . . . . . . . . . . . 24

problèmes de hamp (diusion de la haleur, propagation d'ondes ...et ). Ces te hniques

vent une grande exibilité des paramètres de variation.

modélisent les phénomènes physiques présents (é oulement de uides, transfert de ha-

1.3 Classi ation des EDP

Leur lassi ation se fait de la manière suivante :