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

Volumes Finis

TP3
Volumes Finis pour l’équation de transport

Le but de ce TP est de résoudre l’équation de transport 2D par Volumes Finis. On considère un domaine
Ω ⊂ R2 représentant un récipient rempli par un fluide en mouvement. On suppose connu le champ de
vitesse V(x) en chaque point x ∈ Ω du fluide. La vitesse V est à divergence nulle i.e. div (V) = 0 dans Ω.
A l’instant initial (t = 0), on ajoute un composé (du sel par exemple) dans le fluide avec une concentration
(densité) donnée u0 (x) et on s’intéresse à l’évolution de la concentration u = u(x, t) du composé dans le
fluide. Cette évolution est régie par la conservation de la masse :

∂u
+ div (uV) = 0 dans Ω × R+


∂t

(P ) u(V · n) = 0 sur ∂Ω × R+


u(·, 0) = u0 dans Ω.


Le vecteur n désigne la normale unitaire à ∂Ω, dirigée vers l’extérieur de Ω.

1 Schéma Volumes Finis explicite

Dans un contexte de Volumes Finis, on considère un partitionnement de Ω en cellules de Z contrôle Ki
1
n
auxquelles sont associées des centres ci = xKi . On cherche les approximations ui ' u(x, tn ) dx
|Ki | Ki
où tn = n∆t. Le schéma s’écrit (cf. Cours)
∆t X
(un+1 − uni ) + |eij |Φ uni , unj , nij = 0


i (1)
|Ki |
eij ⊂∂Ki
eij =(Ki |Kj )

où eij désigne l’arête commune aux cellules Ki et Kj ; nij est la normale unitaire à eij dirigée vers
l’extérieur de Ki (cf. Figure 1). Le flux numérique Φ est choisi par décentrement (upwind) :

+
−
Φ (ui , uj , nij ) = ui V · nij + uj V · nij (2)

avec Z

+
− 1
V · nij = max(V · nij , 0), V · nij = min(V · nij , 0), V= V dΓ.
|eij | eij

Sur chaque arête e, on calcule la vitesse moyenne V par :

V(x1 ) + V(x2 )
V= où x1 et x2 sont les deux extrémités de l’arête e. (3)
2

2 Cellules et maillage
Les cellules de contrôle sont les triangles Ki = Ti d’une triangulation de Delaunay de type ”Eléments Fi-
nis”. Les centres ci = xKi des cellules sont les barycentres des triangles. On remarquera qu’en particulier,

Kj

eij cj
Ki nij
ci

Figure 1 – Cellules de contrôle pour les volumes finis

il n’est pas nécessaire que les angles des triangles soient tous plus petits que π/2.

1
Récupérer l’archive suivante :

[Link]
Cette archive contient des scripts MATLAB à compléter ainsi que les mailleurs “Volumes Finis” Delau-
nay/Voronoı̈ contenus dans le répertoire meshTRIVOR/ 1 . Les maillages Delaunay/Voronoı̈ du domaine sont
réalisés avec la fonction mesh2d_trivor.m (voir le TP1 ou bien le fichier ./meshTRIVOR/doc/mesh2d_trivor.pdf
dans l’archive pour les descriptions complètes des structures de données associées aux maillages Delau-
nay/Voronoı̈).

3 Système linéaire
Soit N le nombre de cellules de contrôles (= nombre de triangles de la triangulation). Le schéma (1) avec
(2) s’écrit, pour i = 1, · · · , N ,
 ∆t X +
 ∆t X −
un+1
i = 1− |eij | (V · nij ) uni − |eij | (V · nij ) unj (4)
|Ki | |Ki |
eij ⊂∂Ki eij ⊂∂Ki
eij =(Ki |Kj ) eij =(Ki |Kj )

On pose un = (un1 , · · · , unN )> et le système linéaire correspondant au schéma explicite, s’écrit :
un+1 = (Id − ∆tA)un (5)
où A est la matrice de taille N × N , de la forme

 j1 j2 i j3 
− − −
 αij 1
αij 2
1 X
+
αij 3

A = · · · ··· ··· αij ··· · · ·
 ← ligne i

|Ki | |Ki | |Ki | |Ki |
 eij ⊂∂Ki 

± ±
avec αij = |eij | (V · nij ) . Les trois cellules Kj1 , Kj2 , Kj3 sont les trois triangles adjacents au triangle
Ki (cf. Figure ci-dessous).

Kj Kj
1 3
Ki

Kj
2

4 Assemblage
Pour construire la matrice A, on boucle sur les arêtes intérieures des triangles. Pour une arête courante
e 6⊂ ∂Ω, on ajoute les contributions des deux triangles ayant e comme arête commune. On considère ainsi
la matrice élémentaire Aelem suivante :
i j
 
|eij | + |eij | −
 + |K | (V · nij ) +
|Ki |
(V · nij )  i
i
Aelem = 
 

 |eij | − |eij | +
+ (V · nji ) + (V · nji ) j
|Kj | |Kj |
− + + −
Remarque : Puisque nji = −nij , on a les relations (V · nji ) = − (V · nij ) et (V · nji ) = − (V · nij ) .

1. Les mailleurs “Volumes Finis” Delaunay/Voronoı̈ peuvent aussi être récupérés directement à l’adresse [Link]
[Link]/~[Link]/Matlab/[Link]

2
 n ji

Ki

eij n ij

Kj

L’assemblage de la matrice A se fait alors de la façon suivante :

|eij | + |eij | −
A(i, i) = A(i, i) + (V · nij ) , A(i, j) = (V · nij ) ,
|Ki | |Ki |
|eij | + |eij | −
A(j, i) = − (V · nij ) A(j, j) = A(j, j) − (V · nij )
|Kj | |Kj |

5 Condition de stabilité CFL

Pour le schéma ”Volumes Finis” (1) explicite en temps, la condition de stabilité CFL s’écrit :
 1 X 
+
∆t max |eij | (V · nij ) ≤1 (6)
1≤i≤N |Ki |
eij ⊂∂Ki
eij =(Ki |Kj )

Travail demandé.
1. Implémenter le schéma explicite (5) en MATLAB, en complétant le script
vf_transport.m
qui est fourni dans l’archive [Link].
2. Tests numériques. On choisit Ω = [0, 1]2 . Tester votre code avec les données suivantes
 
2
(a) Champ de vitesse constant V(x, y) =
−1
1 3 1 3
Donnée initiale u0 (x, y) = 1D (x, y) où D = [ , ] × [ , ]
8 8 2 4
Que constatez-vous ?
1 1
 
cos(π(x − )) sin(π(y − ))
(b) Champ de vitesse V(x, y) = 
 2 2 
1 1 
− sin(π(x − )) cos(π(y − ))
2 2
1 3 1 3
Donnée initiale u0 (x, y) = 1D (x, y) où D = [ , ] × [ , ]
8 8 2 4
−200y(y − 1)(2y − 1)x2 (x − 1)2
 
(c) Champ de vitesse V(x, y) =
200x(x − 1)(2x − 1)y 2 (y − 1)2
Donnée initiale u0 (x, y) = 1D (x, y) où D est le cercle de centre (0.3, 0.4) et de rayon r0 = 0.15

Vérifier qu’à chaque itération en temps, la masse totale du composé reste constante : le schéma
est globalement conservatif. Tester la stabilité du schéma.
3. Ecrire le schéma Volumes Finis implicite obtenu en prenant le flux numérique

Φ(un+1
i , un+1
j , nij )

dans les relations (1). Ecrire ce schéma sous forme matricielle M un+1 = un . Implémenter ce
schéma en MATLAB. Reprendre les tests numériques précédents et tester la stabilité du schéma
implicite.