Simulation Numérique d’un Écoulement

Monophasique Incompressible dans un Milieu Poreux

El Moutaouakil Hicham Chahid Chadi

1. Dérivation du Modèle
Le modèle repose sur la loi de Darcy pour un fluide incompressible :
Λ
u = − ∇P
µ
Combinée à la conservation de la masse ∇ · u = f , on obtient :
 
ϱ
−∇ · Λ(x) ∇P (x) = f (x)
µ

En posant µ = ρ = 1, on obtient :

−∇ · (Λ(x)∇P (x)) = f (x).

2. Contexte Physique et Applications

Ce type de modèle est largement utilisé dans la simulation de phénomènes physiques
impliquant des écoulements dans des milieux poreux : stockage du CO2 ou de l’hydrogène,
géothermie, pétrolières. Ces modèles permettent de comprendre la distribution de la
pression dans un réservoir souterrain et d’optimiser les stratégies d’exploitation.

3. Discrétisation par Différences Finies

On considère un domaine Ω = [0, 1]2 discrétisé en un maillage cartésien uniforme de pas
h = N1 . En supposant que Λ est constante par cellule, on approxime le laplacien par le
schéma à 5 points :
 
Pi+1,j + Pi−1,j + Pi,j+1 + Pi,j−1 − 4Pi,j
−∇ · (Λ∇P ) ≈ −Λ
h2

Lorsque Λ est variable, on utilise des moyennes harmoniques aux interfaces, par exemple
:   −1
1 1 1
Λi+1/2,j = +
2 Λi,j Λi+1,j

1
4. Schéma Numérique
Le schéma à 5 points avec perméabilité variable s’écrit pour un point intérieur (i, j) :

1 
Λi+1/2,j (Pi+1,j − Pi,j ) − Λi−1/2,j (Pi,j − Pi−1,j ) + Λi,j+1/2 (Pi,j+1 − Pi,j )
h2 
−Λi,j−1/2 (Pi,j − Pi,j−1 ) = fi,j

où les perméabilités aux interfaces sont approchées par des moyennes harmoniques :
  −1   −1
1 1 1 1 1 1
Λi+1/2,j = + , Λi,j+1/2 = +
2 Λi,j Λi+1,j 2 Λi,j Λi,j+1

Les conditions de Dirichlet sont imposées directement dans le système linéaire :

pour tout point (i, j) situé sur le bord du domaine, on remplace l’équation discrète par :

Pi,j = g(xi , yj )
Ce qui correspond, dans l’implémentation numérique, à fixer :

• Ak,k = 1

• bk = g(xi , yj )

où k est l’indice global du point (i, j). Cela garantit l’imposition exacte de la valeur
de pression sur le bord du domaine.

5. Algorithme de Résolution
1. Discrétiser le domaine Ω

2. Construire la matrice A et le vecteur b

3. Appliquer les conditions de Dirichlet

4. Résoudre le système linéaire

5. Évaluer l’erreur

6. Validation par une Solution Exacte

On utilise P (x, y) = xy(1 − x)(1 − y), avec Λ = 1. Le terme source est :

f (x, y) = 2(x(1 − x) + y(1 − y))

Code MATLAB :

2
1 % Parameters
2 Nx = 20; Ny = 20;
3 Lx = 1; Ly = 1;
4 hx = Lx / Nx ; hy = Ly / Ny ;
5 N = Nx * Ny ;
6 Lambda_val = 1;
7

8 % Grid ( cell centers )

9 x = linspace ( hx /2 , Lx - hx /2 , Nx ) ;
10 y = linspace ( hy /2 , Ly - hy /2 , Ny ) ;
11 [X , Y ] = meshgrid (x , y ) ;
12

13 % Exact solution : P (x , y ) = x * y *(1 - x ) *(1 - y )

14 P_exact = X .* Y .* (1 - X ) .* (1 - Y ) ;
15

16 % Source term : f (x , y ) = 2 * Lambda * ( x *(1 - x ) + y *(1 - y ) )

17 f_source = 2 * Lambda_val * ( X .* (1 - X ) + Y .* (1 - Y ) ) ;
18

19 % Assembly
20 A = sparse (N , N ) ;
21 b = zeros (N ,1) ;
22 idx = @ (i , j ) (j -1) * Nx + i ;
23

24 for j = 1: Ny
25 for i = 1: Nx
26 k = idx (i , j ) ;
27 xp = x ( i ) ;
28 yp = y ( j ) ;
29

30 % Dirichlet boundary ( P = 0 on )
31 if i == 1 || i == Nx || j == 1 || j == Ny
32 A (k , k ) = 1;
33 b ( k ) = 0; % since P_exact = 0 on boundary
34 else
35 dx2 = hx ^2; dy2 = hy ^2;
36

37 A (k , idx ( i +1 , j ) ) = - Lambda_val / dx2 ;

38 A (k , idx (i -1 , j ) ) = - Lambda_val / dx2 ;
39 A (k , idx (i , j +1) ) = - Lambda_val / dy2 ;
40 A (k , idx (i ,j -1) ) = - Lambda_val / dy2 ;
41 A (k , k ) = 2 * Lambda_val * (1/ dx2 + 1/ dy2 ) ;
42

43 b ( k ) = f_source (j , i ) ;
44 end
45 end
46 end
47

48 % Solve
49 P = A \ b;
50 Pnum = reshape (P , Nx , Ny ) ’;
51

3
52 % Error
53 err = abs ( Pnum - P_exact ) ;
54 err_L2 = sqrt ( hx * hy * sum ( err .^2 , ’ all ’) ) ;
55 err_inf = max ( err (:) ) ;
56

57 % Display errors
58 disp ([ ’ L2 ␣ error ␣ ␣ ␣ ␣ ␣ = ␣ ’ , num2str ( err_L2 ) ]) ;
59 disp ([ ’ Linf ␣ error ␣ ␣ ␣ = ␣ ’ , num2str ( err_inf ) ]) ;

Erreurs numériques obtenues :

Nx = Ny h Erreur L2 Erreur L∞
10 0.1 9.6615×10−3 1.1756×10−2
20 0.05 4.8665×10−3 6.0785×10−3
40 0.025 2.4407×10−3 3.084×10−3
80 0.0125 1.2222×10−3 1.5525×10−3
160 0.00625 6.1153×10−4 7.7878×10−4

Table 1: Erreurs entre solution numérique et solution exacte pour différents maillages.

On observe que les erreurs L2 et L∞ diminuent à mesure que le pas de discrétisation

h diminue. Cela confirme la bonne implémentation du schéma numérique et sa capacité
à approcher correctement la solution exacte dans le cadre de la validation. Les résultats
indiquent une convergence régulière et stable du schéma.

7. Simulation sur un Milieu Hétérogène

Code MATLAB :
1 % Param tres
2 Nx = 50; Ny = 50; % nombre de mailles
3 Lx = 1; Ly = 1;
4 hx = Lx / Nx ; hy = Ly / Ny ;
5 N = Nx * Ny ;
6

7 % Grille
8 x = linspace ( hx /2 , Lx - hx /2 , Nx ) ;
9 y = linspace ( hy /2 , Ly - hy /2 , Ny ) ;
10 [X , Y ] = meshgrid (x , y ) ;
11

12 % Matrice de p e r m a b i l i t (x)
13 Lambda = ones ( Ny , Nx ) * 1e -4;
14 Lambda ( Y >= 0.3 & Y <= 0.7 & X >= 0.3 & X <= 0.7) = 1;
15

16 % Assemblage du s y s t m e
17 A = sparse (N , N ) ;
18 b = zeros (N ,1) ;
19

20 idx = @ (i , j ) (j -1) * Nx + i ; % index dans le vecteur

21

4
22 for j = 1: Ny
23 for i = 1: Nx
24 k = idx (i , j ) ;
25

26 % C o o r d o n n e s physiques
27 xp = x ( i ) ;
28 yp = y ( j ) ;
29

30 % Conditions de Dirichlet
31 if i == 1 % x = 0 = > g = 0.5
32 A (k , k ) = 1;
33 b ( k ) = 0.5;
34 elseif j == 1 % y = 0 => g = 1
35 A (k , k ) = 1;
36 b ( k ) = 1;
37 elseif i == Nx || j == Ny
38 A (k , k ) = 1; % autres bords = > g = 0
39 b ( k ) = 0;
40 else
41 % Coefficients aux interfaces
42 Lc = Lambda (j , i ) ;
43 Lr = Lambda (j , i +1) ;
44 Ll = Lambda (j ,i -1) ;
45 Lu = Lambda (j -1 , i ) ;
46 Ld = Lambda ( j +1 , i ) ;
47

48 Ax = hx ; Ay = hy ;
49 dx2 = hx ^2; dy2 = hy ^2;
50

51 % Moyennes harmoniques
52 Lx_p = 2* Lc * Lr / ( Lc + Lr ) ;
53 Lx_m = 2* Lc * Ll / ( Lc + Ll ) ;
54 Ly_p = 2* Lc * Ld / ( Lc + Ld ) ;
55 Ly_m = 2* Lc * Lu / ( Lc + Lu ) ;
56

57 A (k , idx ( i +1 , j ) ) = - Lx_p / dx2 ;

58 A (k , idx (i -1 , j ) ) = - Lx_m / dx2 ;
59 A (k , idx (i , j +1) ) = - Ly_p / dy2 ;
60 A (k , idx (i ,j -1) ) = - Ly_m / dy2 ;
61 A (k , k ) = ( Lx_p + Lx_m ) / dx2 + ( Ly_p + Ly_m ) / dy2 ;
62 end
63 end
64 end
65

66 % R solution
67 P = A \ b;
68 Pmat = reshape (P , Nx , Ny ) ’;
69

70 % Create the figure and set up subplots

71 figure ;
72

5
 Figure 1: Pression numérique.

73 % Subplot for view (2)

74 subplot (1 , 2 , 1) ;
75 surf (x , y , Pmat ) ;
76 title ( ’ Pression ␣ P (x , y ) ␣ -␣ View ␣ 2 ’) ;
77 xlabel ( ’x ’) ;
78 ylabel ( ’y ’) ;
79 zlabel ( ’P ’) ;
80 shading interp ;
81 view (2) ;
82 colorbar ;
83

84 % Subplot for view (3)

85 subplot (1 , 2 , 2) ;
86 surf (x , y , Pmat ) ;
87 title ( ’ Pression ␣ P (x , y ) ␣ -␣ View ␣ 3 ’) ;
88 xlabel ( ’x ’) ;
89 ylabel ( ’y ’) ;
90 zlabel ( ’P ’) ;
91 shading interp ;
92 view (3) ;
93 colorbar ;

Interprétation : Le profil de la pression numérique obtenue est en accord avec les

attentes physiques. On observe une accumulation de pression dans le coin inférieur gauche
du domaine, à l’interface des conditions P = 1 et P = 0.5. Cette zone agit comme une
région de haute pression, à partir de laquelle la pression décroît dans le reste du domaine.
La décroissance est toutefois fortement influencée par la perméabilité du milieu. En effet,
dans la zone Ω2 = Ω \ Ω1 , la faible perméabilité empêche la propagation efficace du fluide,
ce qui crée des gradients de pression plus importants aux interfaces.
En revanche, dans la région centrale Ω1 , où Λ = 1, le fluide peut circuler plus librement,
ce qui se traduit par une variation plus douce de la pression. Cette zone agit donc comme
un canal préférentiel facilitant le passage de la pression vers le centre du domaine. Les
résultats illustrent bien le comportement typique d’un écoulement en milieu hétérogène,

6
où les zones de haute perméabilité concentrent le transport, tandis que les zones peu
perméables bloquent l’écoulement.
Ces observations confirment que le schéma numérique mis en œuvre est capable de
capturer fidèlement les phénomènes physiques attendus dans un problème de Darcy avec
perméabilité discontinue.

8. Conclusion
Dans ce projet, nous avons mis en œuvre une méthode numérique basée sur les différences
finies pour résoudre une équation de Darcy modélisant un écoulement monophasique
incompressible dans un milieu poreux. Après avoir validé le code dans un cas simple à
perméabilité constante, nous avons simulé un cas plus complexe avec une perméabilité
discontinue, représentant un milieu hétérogène. Les résultats obtenus sont en accord avec
les phénomènes physiques attendus, notamment la canalisation du flux dans les zones plus
perméables et le blocage dans les zones à faible perméabilité.
Cependant, la méthode des différences finies (DFs) présente certaines limitations im-
portantes :

• Sensibilité à la régularité du maillage : La méthode est naturellement adaptée

à des maillages cartésiens réguliers. Cela limite son application à des géométries
complexes, souvent rencontrées dans les problèmes réels d’ingénierie.

• Traitement des discontinuités : En présence de discontinuités dans les coeffi-

cients (comme ici avec la perméabilité), le schéma peut engendrer des erreurs local-
isées si les interfaces ne sont pas bien alignées avec le maillage. Cela peut nécessiter
un raffinement de maillage ou l’utilisation de schémas plus robustes (par exemple
les volumes finis).

• Conservation locale non garantie : Contrairement à d’autres méthodes comme

les volumes finis, les différences finies ne garantissent pas automatiquement la con-
servation locale des flux, propriété essentielle dans les simulations de transport ou
d’écoulement.

• Imposition rigide des conditions aux limites : La méthode impose directement

les conditions de Dirichlet via des substitutions dans la matrice, ce qui peut limiter
la flexibilité, notamment pour des conditions plus complexes (Robin, Neumann,
etc.).

Malgré ces limites, les différences finies constituent une méthode simple, intuitive et
efficace pour les problèmes à géométrie régulière, et permettent une bonne compréhen-
sion des phénomènes physiques en jeu. Elles constituent ainsi une étape précieuse dans
l’initiation à la simulation numérique.