Master 2 IMOI - Méthodes avancées de résolution numérique des EDP 2017-18

Volumes Finis

TP2
Equation de la chaleur

1 Introduction
Le but de ce TP est de résoudre l’équation de la chaleur en 2D par une méthode de Volumes Finis avec
différentes discrétisations en temps. On considère un domaine polygonal Ω ⊂ R2 . Pour T > 0, on cherche
une fonction u = u(x, t) pour x ∈ Ω et t ∈ [0, T ] vérifiant l’équation de la chaleur suivante :

∂u
− ν∆u = f dans Ω×]0, T [ (1)
∂t
u = g sur ∂Ω×]0, T [ (2)
u(·, 0) = u0 dans Ω. (3)

Les fonctions f , g, u0 et le paramètre ν > 0 sont donnés. On associe au domaine Ω un maillage ”Volumes
Finis” dont les volumes de contrôle sont les cellules de Voronoı̈ associées à une triangulation de Delaunay
de Ω. A chaque cellule K, on associe les centres xK qui sont les sommets des triangles (cf. TP1).

xK
e

n K,e S
L
xL

Figure 1 – Volumes de contrôle (Voronoı̈) et triangulation de Delaunay

2 Formulation en Volumes Finis

2.1 Discrétisation en espace
On intègre l’équation (1) sur une cellule K et en utilisant la formule de la divergence, on obtient :
Z Z
d
u(x, t) dx − ν ∇u · nK dΓ = |K|fK
dt K ∂K

1 R
avec fK (t) = f (x, t) dx et nK désigne la normale unitaire dirigée vers l’extérieur de K. On pose
|K| K
1 R
uK (t) = u(x, t) dx et on a, pour t ∈ [0, T ],
|K| K

duK X Z
|K| (t) − ν ∇u · nK,e dΓ = |K|fK (t), (4)
dt e
e∈EK

1
où EK désigne l’ensemble des arêtes de la cellule K et nK,e est la normale unitaire à e dirigée vers
l’extérieur de K.
Z
On approche le flux − ∇u · nK,e dΓ ' FK,e avec
e

(uK − uL )
FK,e = |e| si e 6⊂ ∂Ω. (5)
|xK − xL |
Si K possède une arête sur le bord ∂Ω, on impose
Z
1
uK = ue = g(x, t) dΓ si e ⊂ ∂Ω. (6)
|K| e

Le schéma Volumes Finis de discrétisation en espace s’écrit ainsi, pour t ∈ [0, T ],

duK X
|K| (t) + νFK,e (t) = |K|fK (t). (7)
dt
e∈EK
e6⊂∂Ω

2.2 Discrétisations en temps

Soit n ∈ N∗ et ∆t = T
n> 0 le pas de discrétisation en temps avec tk = k∆t pour k = 0, 1, · · · , n. On note
1 R
l’approximation ukK ' uK (tk ) = u(x, tk ) dx.
|K| K

2.2.1 Schéma explicite

Le schéma explicite s’écrit, pour tout K ∈ T , k = 0, · · · , n,

(uk+1
K − ukK ) X
k k
|K| + νFK,e = |K|fK , (8)
∆t
e∈EK
e6⊂∂Ω

avec
k (ukK − ukL )
FK,e = |e| (e 6⊂ ∂Ω). (9)
|xK − xL |

On écrit les relations (8) pour toute cellule K intérieure i.e. une cellule ne possédant pas d’arête sur le
bord ∂Ω. Pour une cellule K ayant une arête sur ∂Ω, on impose uK = ue donnée par (6). On obtient
ainsi le schéma :
— Pour toute cellule K intérieure,

ν∆t X |e|
uk+1 = ukK − uk − ukL + ∆t fK
k


(10)
K
|K| |xK − xL | K
e∈EK
e=(K|L)

— Pour toute cellule K possèdant une arête e ⊂ ∂Ω, on impose uk+1

K = ue où ue est donnée par (6).

Système linéaire
On prend (pour simplifier) g = 0. Soit N le nombre de cellules de contrôles intérieures. On regroupe les

>
inconnues dans le vecteur uk = (ukK1 , · · · , ukKN )> et on note f k = fK
k
1
k
· · · , fK N
.
Le système linéaire correspondant à (10) s’écrit

uk+1 = IN − ν∆tH −1 A uk + ∆t f k ,


(11)

où IN est la matrice identité d’ordre N et A est la matrice de taille N × N correspondant au Laplacien
(cf. TP1 ”Volumes Finis pour le Laplacien”). La matrice diagonale H de taille N × N est définie par
 
|K1 | 0
 |K2 | 
H= . (12)
 
 . .. 
0 |KN |

2
Le système (11) est explicite au sens où il n’y a pas de système linéaire à résoudre pour déterminer uk+1
à partir de uk . En revanche, le schéma (11) doit vérifier la condition de stabilité
 
 1 X |e| 
 ≤ 1.
λ = ν∆t max  (13)
1≤i≤N  |Ki | |xKi − xL | 
e∈EKi
e=(Ki |L)

2.2.2 Schéma implicite

Le schéma implicite s’écrit, pour tout K ∈ T , k = 0, · · · , n,
k+1
(uK − ukK ) X
k+1 k+1
|K| + νFK,e = |K|fK , (14)
∆t
e∈EK
e6⊂∂Ω

avec
k+1 (uk+1
K − uk+1
L )
FK,e = |e| (e 6⊂ ∂Ω). (15)
|xK − xL |

On écrit les relations (14) pour toute cellule K intérieure i.e. une cellule ne possédant pas d’arête sur le
bord ∂Ω. Pour une cellule K ayant une arête sur ∂Ω, on impose uK = ue donnée par (6). On obtient
ainsi le schéma :
— Pour toute cellule K intérieure,

ν∆t X |e|
uk+1 uk+1 − uk+1 k+1
= ukK + ∆t fK


K + K L (16)
|K| |xK − xL |
e∈EK
e=(K|L)

— Pour toute cellule K possèdant une arête sur le bord ∂Ω, on impose uk+1
K = ue où ue est donnée
par (6).

Système linéaire
On prend g = 0 (pour simplifier). Le système linéaire correspondant à (16) s’écrit

IN + ν∆tH −1 A uk+1 = uk + ∆t f k+1

où A est la matrice de taille N × N correspondant au Laplacien (cf. TP1) et H est définie par (12). On
écrit le système précédent plutôt sous la forme

H + ν∆tA uk+1 = H uk + ∆t f k+1 ,

(17)

la matrice M = H + ν∆tA étant symétrique (ceci peut être avantageux pour la résolution numérique du
système linéaire (17)). On montre également que M est une M -matrice.

Travail demandé.

1. Implémenter les schémas Volumes Finis explicite (10) et implicite (16) en MATLAB. Construire
les systèmes (11) et (17) à partir du TP1 pour la construction de la matrice A du Laplacien.
2. Pour le schéma implicite (17), à chaque itération sur k, on doit résoudre un système linéaire avec
la matrice M = H + ν∆tA. Cette matrice est indépendante de k et on peut donc la factoriser
au début des itérations par une décomposition P MQ = LU où P et Q sont des matrices de
permutation (resp. lignes et colonnes). La résolution du système linéaire Mx = b à partir de cette
décomposition s’écrit en MATLAB :
[L,U,P,Q]=lu(M,’vector’); % factorisation avant les itérations sur n
...
sol = U\(L\b(P));
x=sol(Q);

3
 C’est beaucoup plus rapide d’utiliser cette décomposition (complexité en O(N 2 )) plutôt que de
résoudre à chaque itération le système avec x=M\b (complexité en O(N 3 )).
3. Pour le schéma explicite (11), une fois le maillage construit, calculer la quantité
 
 1 X |e| 
max  .
1≤i≤N  |Ki | |xKi − xL | 
e∈EKi
e=(Ki |L)

Choisir ensuite le pas de temps ∆t pour que la condition de stabilité (13) du schéma explicite soit
satisfaite.
4. Tester votre code avec l’exemple suivant : On choisit Ω =]0, a[×]0, b[, a = 2.5, b = 1, T = 1 et
 2
k22
      
k1 k1 π k2 π
u(x, t) = exp −νπ 2 + t sin x 1 sin x 2 ,
a2 a2 a b

pour x = (x1 , x2 ) ∈ [0, a] × [0, b], t ∈ [0, T ]. On prend k1 = 2, k2 = 1. La fonction u est la solution
de l’équation de la chaleur (1)–(3) avec f ≡ 0 et la donnée initiale
   
k1 π k2 π
u0 (x) = sin x1 sin x2 .
a b

5. Pour le schéma explicite (11), tester la condition de stabilité (13) avec λ > 1.