1ère année Master D.F.E.

Département de Physique Université

S2 2021/2022 A.Mira, Béjaïa

Master 1— Dynamique des Fluides et Énergétique —2021/2022

—Examen Final –Analyse numérique– Méthode des Volumes Finis—

Exercice 1. Considérons le problème suivant

((1 + x2 )u′ )′ = x dans Ω =]0, 1[ (1)

u(0) = 0, u′(1) = 1. (2)

Résoudre le problème avec la méthode des volumes finis en subdivisant le domaine en quatre
volumes de contrôle de longueur δ x dont les næuds sont au milieu de chaque volume.

Exercice 2. Soit u une fonction continue sur [0; 3]. On cherche à approcher numériquement la
solution u de l’équation aux dérivées partielles :
 2 2
 ∂ u(x, y) + ∂ u(x, y) = 0,

∀(t, x) ∈ Ω =]0, 3[2 ,
∂ x2 ∂ y2 (3)
 u(x, 0) = u(x, 3) = u(3, y) = 0 et u(0, y) = 1


On suppose u de classe C 2 solution de (12).

1. On pose u(x, y) = φ(x)ψ(y), résoudre analytiquement le problème aux limites 12.

2. On considère une subdivision du domaine [0, 3]2 en neuf (9) volumes de contrôle égaux (9
carré). Par la méthode des volumes finis, on cherche la valeur approchée de u(xi , yj ), que
l’on notera ui,j ≡ u(P ), au nœuds situés au centre de chaque volume de contrôle.

(a) Déterminer une relation que devront vérifier les coefficients ui,j afin de correspondre à
la discrétisation de l’équation 12a.

(b) En parcourant tout les nœuds, trouver le système matriciel discret Au = b qui donne
cette solution approchée du problème 12.
 Corrigé de l’examen
Corrigé de l’exo1

1. Maillage. Nous Subdivisons le domaine de calcul [0, 1] en 4 volumes de contrôle égaux, donc
de longueur δ x = 0.25. Tels que V ci = [xi−1/2 , xi+1/2 ], i = 1, .., 4, avec x1/2 = 0, x3/2 = 0.22,
x5/2 = 0.5, x7/2 = 0.75 et x9/2 = 1.0. Les points noeuds sont au centre des volumes de
contrôle, c’est-à-dire P ∈ {x1 , x2 , x3 , x4 } = {0.125, 0.375, 0.625, 0.875}. (voir la figure 1)
x1/2 x3/2 x5/2 x7/2 x9/2
x1 x2 x3 x4

b=1
a=0 δx = 0.25
W = xi−1 P = xi E = xi+1

w e

Figure 1:

2. Équation discrétisée.

En intégrant l’équation (1) sur un volume de contrôle V c = [w, e], on obtient

 2 e
x
Z Z
2 ′ ′ 2 ′ e

 
(1 + x )u dx = f (x)dx ⇔ (1 + x )u w =
2 w
Vc Vc
 1 2
(1 + e2 )u′ (e) − (1 + w 2 )u′(w) = e − w 2 . (4)


⇔
2

(a) Dans le cas des volumes de contrôle intérieur ou bien les noeuds intérieur, (i = 2, 3),
nous utilisons le schémas centrer pour le calcul des dérivées en w et en e, et elles sont
données par
uP − uW uE − uP
u′ (w) = ; u′ (e) =
δx δx
En substituant dans l’équation (4), l’équation discrétisé aux points noeuds i = 2, 3 est
donnée par
1  1
−(1 + w 2 )uW + (2 + w 2 + e2 )uP − (1 + e2 )uE = − e2 − w 2


(5)
δx 2
 1
Pour δx = 0.25 = l’équation 5et alors
4
1 2
1 2
−(1 + w 2 )uW + (2 + w 2 + e2 )uP − (1 + e2 )uE = δx e − w2 = e − w2


(6)
2 8
D’où en résumé
1 1
i. Pour i = 2 on a w = 0.25 = et e = 0.5 = , donc
4 2
1 1 3
(−17uW + 37uP − 20uE ) = (7)
16 8 16
donc
3
−17uW + 37uP − 20uE = (8)
8
1 3
ii. Pour i = 3 on a w = 0.5 = et e = 0.75 = , donc
2 4
5
−20uW + 61uP − 25uE = (9)
8
(b) Dans le cas du premier volume de contrôle V ci = [x1/2 , x3/2 ] le calcul de la dérivée au
point w = x1/2 = 0 se fait par le schémas en avant comme suite
uP − u(0) 2uP
u′ (w) = =
δx /2 δx
En substituant dans l’équation (4), l’équation discrétisé au point noeud i = 1 est donnée
par
1
49uP − 17uE = . (10)
8
(c) Dans le cas du dernier volume de contrôle V ci = [x7/2 , x9/2 ] (i = 4)
et la valeur de u′ en e = x9/2 = 1 est donnée par la condition au limite u′ (e) = u′ (1) = 1.

En substituant dans l’équation (4), l’équation discrétisé au point noeud i = 4 est donnée
par
9
25(−uW + uP ) = (11)
8

3. Le système Au = b.
    
49 −17 0 0 u1 1
    
 −17 37 −20 0   u2  1  3 
    
AT = b ⇔   =  
 0 −20 61 −25   u3  8  5 
    
    
0 0 −25 25 u4 9
 4. La matrice A est une matrice de bande donc la décomposition de LU donne des matrice de
bande de même longueur. Soit
  
1 0 0 0 49 −17 0 0
 17 1524
  
 
 − 1 0 0 
 0 −11 0
 
A = LU =  49 49
 
 0 245  18341 
− 1 0 
 0
 0 −25 
 381 381 
 9525  220400 
0 0 − 1 0 0 0
18341 18341

On a 


 LY = b

Au = b ⇔ LUu = b ⇔ et


 Uu = Y


La résolution des deux systèmes triangulaire supérieur et inférieur, nous donne

 
t t 91 235 311 11951
u= , , ,
3712 3712 3712 92800

ou bien
t
u = t (0.024, 0.063, 0.083, 0.128)

0.125„ , D’où le graphe Gu = {(xi , ui )} de la solution approché s’écrit

Gu = {(0, 0), (0.125, 0.024), (0.375, 0.063), (0.625, 0.083), (0.875, 0.128)}

Corrigé de l’exo2. Soit

 2 2
 − ∂ u(x, y) − ∂ u(x, y) = 0,

∀(t, x) ∈ Ω =]0, 3[2 ,
∂ x2 ∂ y2 (12)
 u(x, 0) = u(x, 3) = u(3, y) = 0 et u(0, y) = 1


1. On pose u(x, y) = φ(x)ψ(y). L’équation 12a devient




 φ(x)′′ − λ2 φ(x) = 0
φ(x)′′ ψ(y)′′ 
φ(x)′′ ψ(y) + φ(x)ψ(y)′′ = 0 ⇒ =− = λ2 ⇒ et
φ(x) ψ(y) 

 ψ(y)′′ + λ2 ψ(y) = 0


donc 


 φ(x) = K1 cosh(λ x) + K2 sinh(λ x)

et


 ψ(y) = K cos(λ y) + K sin(λ y)

3 4
 avec λ constante. La solution générale est u(x, y) = (K1 cosh(λ x) + K2 sinh(λ x)) (). En
utilisant les conditions aux limites, on obtient

u(x, 0) = 0 ⇒ φ(x)ψ(0) = 0 ⇒ ψ(0) = 0 ⇒ K3 = 0 ⇒ ψ(y) = K4 sin(λ y)

nπ  nπ 
u(x, 3) = 0 ⇒ φ(x)ψ(3) = 0 ⇒ ψ(3) = 0 ⇒ K4 sin(3λ) = 0 ⇒ λ = ⇒ ψ(y) = K4 sin y
3 3

Si on pose φ(x) = K1 cosh(λ (x − 3)) + K2 sinh(λ (x − 3)), alors

nπ
u(3, y) = 0 ⇒ φ(3)ψ(y) = 0 ⇒ ψ(0)φ(3) = 0 ⇒ K1 = 0 ⇒ φ(x) = K2 sinh( (x − 3))
3
n=+∞
nπ nπ
P


Donc u(x, y) = an sinh 3
(x − 3) sin 3
y . Il reste à déterminer les coefficients an .
n=1
On a
n=+∞
X  nπ 
u(0, y) = 1 ⇒ an sinh (−nπ) sin y =1
n=1
3

2 R1 nπ

6 nπ

1 6 nπ



Donc an sinh (−nπ) = sin y dy = cos y = cos − 1 . Finalement
1 0 3 nπ 3 0 nπ 3

n=+∞
6 cos nπ


X
3
−1  nπ   nπ 
u(x, y) = sinh (x − 3) sin y
n=1
nπ sinh (−nπ) 3 3

2. Maillage: Nous Subdivisons alors le domaine de calcul Ω =]0, 3[×]0, 3[ en 3 × 3 = 9

volumes finis (V ci,j )i=1..3;j=1..3, voir figure 2a tels que . V cij =]xi−1/2 , xi+1/2 [×]yj−1/2 , yj+1/2[,
avec x1/2 = 0, x3/2 = 1, x5/2 = 2, y7/2 = 3 et y1/2 = 0, y3/2 = 1, y5/2 = 2, y7/2 = 3 on note
∆xi = xi+1 − xi = 1 et ∆yj = yj+1 − yj = 1
Pour i, j quelconques, on note par P = (xi , yj ) le centre du volume de contrôle V ci,j qui
représente le nœud principale, par W = (xi−1 , yj ) nœ ud du volume V ci−1,j , E = (xi+1 , yj )
nœ ud du volume V ci+1,j , S = (xi , yj−1) nœ ud du volume V ci,j−1 et N = (xi , yj+1) nœ ud du
volume V ci,j+1 (voir figures 2b ).

3. Discrétisation : On note par A = (xi−1/2 , yj−1/2), B = (xi+1/2 , yj−1/2), C = (xi+1/2 , yj+1/2 )

et D = (xi−1/2 , yj+1/2) les sommets du volume de contôle principale dans le sens direct. En
intégrant l’équation (??) sur ce volume de contrôle, on obtient

Z Z Z Z Z
∇ u~nds = ∇ u~nds + ∇ u~nds + ∇ u~nds + ∇ u~nds = 0 (13)
∂Vc [AB] [BC] [CD] [DA]
 N N
u(x, 3) = 0

3 D C
n
P7 P8 P9

w e E
W P
W E
P4 P5 P6 s
u(0, y) = 1 A B
u(3, y) = 0

S
S
P1 P2 P3

0 3
u(x, 0) = 0

Figure 2:

où ~n est le vecteur normale extérieur au V c. On remarque que ~n[AB] = −~j, ~n[BC] = ~i,
~n[CD] = ~j et ~n[DA] = −~i. Alors (13) s’écrit

∂ ∂ ∂ ∂
Z Z Z Z
− uds + uds + uds − uds = 0 (14)
∂y ∂x ∂y ∂x
[AB] [BC] [CD] [DA]

Pour calculer ces intégrales curvilignes, en approxime le gradient de u sur la frontière de

V c par des polynômes d’interpolation. On suppose que le gradient est constant à travers la
frontière ∂ V c. On a, alors
 
∂  ∂  uP − uS
u = u = (15)
∂y [AB] ∂y s δ SP
 
∂  ∂  uE − uP
u = u = (16)
∂x [BC] ∂x e δ PE
 
∂  ∂  uN − uP
u = u = (17)
∂y [CD] ∂y n δ PN
 
∂  ∂  uP − uW
u = u = (18)
∂x [DA] ∂x w δ WP

En substituant (15), (16), (17) et (18) dans (14), on obtient

uP − uS uE − uP uN − uP uP − uW
Z Z Z Z
− ds + ds + ds − ds = 0 (19)
δ SP δ PE δ PN δ WP
[AB] [BC] [CD] [DA]
R R R R R
On note par As = ds, Ae = ds, An = ds, Aw = ds et Su (x, y)dxdy = S̄∆ V ,
[AB] [BC] [CD] [DA] Vc
alors l’équation générale est donnée par
 
As Aw As Aw Ae An Ae An
uS + uW − + + + uP + uE + uN = 0. (20)
δ SP δ WP δ SP δ WP δ PE δ PN δ PE δ PN

Dans notre cas As = Ae = An = Aw = ∆x = 1 . Cette équation peut être réécrite comme

suit
uS + uW − 4uP + uE + uN = 0. (21)

(a) Pour P1 , l’équation est donnée par

−6uP + uE + uN = 1. (22)

(b) Pour P2 , l’équation est donnée par

uW − 5uP + uE + uN = 0. (23)

(c) Pour P3 , l’équation est donnée par

uW − 6uP + uN = 0. (24)

(d) Pour P4 , l’équation est donnée par

uS − 5uP + uE + uN = 1. (25)

(e) Pour P5 , l’équation est donnée par

uS + uW − 4uP + uE + uN = 0. (26)

(f) Pour P6 , l’équation est donnée par

uS + uW − 5uP + uN = 0. (27)

(g) Pour P7 , l’équation est donnée par

uS − 6uP + uE = 1. (28)

(h) Pour P8 , l’équation est donnée par

uS + uW − 5uP + uE = 0. (29)
 (i) Pour P9 , l’équation est donnée par

uS + uW − 6uP = 0. (30)

4. Le système matriciel est donné pa

    
−6 1 0 1 0 0 0 0 0 u1 1
    
 1 −5 1 0 1 0 0 0 0 u2 0 
    
  
    
 0 1 −6 0 0 1 0 0 0 u3 0 
    
  
    
 1 0 0 −5 1 0 1 0 0 u4 1 
    
  
    
Au = b ⇔  0 1 0 1 −4 1 0 1 0 u5 = 0 
    

    
 0 0 1 0 1 −5 0 0 1 u6 0 
    
  
    
u7 1 
    
 0 0 0 1 0 0 −6 1 0   
    
u8 0 
    
 0 0 0 0 1 0 1 −5 1   
    
0 0 0 0 0 1 0 1 −6 u9 0