Rapport de l’application du MEF sur

Réalisé par : encadré par :

Mr .Oussama BOURIHANE
Kaoutar el habti
Khaoula el alaoui
Ousmane oumarou amani
Fawzi zidouemba
Plan:

I. la définition de DARCY
II. l’utilisation de DARCY
III. Application du MEF sur DARCY
1. la forme variationnelle
2. discrétisation
3. transformation géométrique
4. interpolation
5. contribution élémentaire
6. assemblage
7. Acquisition des conditions aux limites
8. Résolution
IV. programmation sur MATLAB
V. résultats
VI. conclusion
 Introduction
La méthode des éléments finis (MEF) est une technique numérique puissante
utilisée pour résoudre des problèmes complexes dans divers domaines de
l'ingénierie et des sciences appliquées. Développée initialement pour analyser
les structures mécaniques, la MEF s'est progressivement étendue à d'autres
domaines, tels que la mécanique des fluides, la thermodynamique,
l'électromagnétisme et bien d'autres.

Le principe fondamental de la MEF repose sur la division d'un domaine continu

complexe en sous-domaines plus simples appelés "éléments finis". Ces
éléments, connectés entre eux par des points appelés "nœuds", permettent de
représenter mathématiquement des phénomènes physiques en résolvant des
équations différentielles partielles.

Rôle dans la résolution des phénomènes physiques

La MEF joue un rôle crucial dans la modélisation et la simulation de
phénomènes physiques, car elle permet :

1. D'analyser des géométries complexes : Les formes géométriques

irrégulières peuvent être approximées avec précision grâce à des éléments
de différentes formes (triangulaires, quadrilatéraux, etc.).
2. De traiter des matériaux hétérogènes : Elle permet de modéliser des
systèmes composés de matériaux avec des propriétés physiques
différentes.
3. D'obtenir des solutions approximatives : Même dans les cas où des
solutions analytiques exactes sont impossibles, la MEF fournit des
résultats numériques précis.
4. De prédire le comportement des systèmes : Elle est utilisée pour
anticiper les performances structurelles, thermiques ou fluidiques sous
diverses conditions.
1-la définition de DARCY
La loi de Darcy est une relation fondamentale qui décrit le mouvement des
fluides à travers des milieux poreux. Elle est utilisée pour modéliser
l'écoulement de l'eau dans les sols et les roches, ainsi que pour d'autres fluides
Cette équation permet de relier le débit d'un fluide à son gradient de pression,
prenant en compte les propriétés du milieu, telles que la perméabilité.

L’hypothèse de l’écoulement de DARCY :

• Cette loi indique que la vitesse de filtration est proportionnelle
au gradient de pression, et elle dépend des propriétés du fluide
(viscosité) et du milieu poreux (perméabilité).
• Le fluide est incompressible.
• L'écoulement est stationnaire
• Le milieu poreux est saturé
• L'écoulement est laminaire et les vitesses sont faibles
2-L’utilisation de DARCY :
L'équation de Darcy est utilisée dans divers domaines pour modéliser
l'écoulement des fluides à travers des milieux poreux. Son application principale
est dans les domaines de l'hydrogéologie, de l'ingénierie pétrolière, de la
géotechnique et de l'écoulement des fluides dans les matériaux poreux. Pour
mieux comprendre son utilisation, nous allons détailler son application à travers
des exemples et explications dans différents contextes.

1. Application en hydrogéologie :

Dans le domaine de l'hydrogéologie, l'équation de Darcy est utilisée pour

estimer l'écoulement de l'eau souterraine à travers des aquifères ou d'autres
formations géologiques poreuses.
2. Application en génie pétrolier :

Dans l'ingénierie pétrolière, l'équation de Darcy est utilisée pour modéliser

l'écoulement des hydrocarbures dans des réservoirs souterrains, souvent pour
estimer le débit de pétrole ou de gaz naturel extrait.

3. Application en génie civil (géotechnique) :

En géotechnique, l'équation de Darcy est utilisée pour calculer le débit d'eau à

travers les sols dans des situations telles que l'évacuation des eaux de pluie ou
l'analyse de la stabilité des sols.
3-Application de la méthode des éléments finis sur l'équation de
DARCY :
Soit un domaine rectangulaire représenté par les dimensions x et y, où l'écoulement est modélisé par
l'équation de Darcy en régime stationnaire.

h
x

−𝐿 𝐿 −ℎ ℎ
Ω=[ ; ]∗[ ; ]
2 2 2 2

I-Forme variationnelle :

La forme variationnelle est une reformulation de l'équation différentielle d'un

problème en une équation intégrale, rendant le problème compatible avec la
méthode des éléments finis. Elle facilite l'approximation de la solution en
intégrant les conditions aux limites dans un domaine discretisé.

L’équation de DARCY s’écrit sous la forme :

-∇⋅(k∇p) = f

Où :

- k est la perméabilité du milieu (m²)

- p estLa pression hydraulique (en [Pa])
- f est le terme source (m/s)
- ∇p : Le gradient de pression (en [Pa/m]

on suppose que :
-K est cst

-Le fluide est incompressible

-L’écoulement est laminaire

∫ 𝛿𝑃 𝛻(𝑘𝛻𝑃)𝑑𝛺 + ∫ 𝛿𝑃𝑓 ⅆ𝛺 = 0
𝛺 𝛺

La formule de GREEN :

⃗ ) = 𝜵(𝝀)𝒗
𝜵(𝝀𝒗 ⃗⃗⃗⃗
⃗ + 𝝀𝜵(𝒗)

On considère :

δp=λ et ⃗⃗⃗⃗𝑣 =k∇p

∇(δp.k∇p)=∇(δp).k∇p + δp ∇(k∇p)

Donc ,

∫ 𝛅𝐏𝛁(𝐤𝛁𝐏) ⅆ𝛀 = ∫ 𝛁(𝛅𝐏𝐤𝛁𝐏) ⅆ𝛀 − ∫ 𝛁𝛅𝐏(𝐤𝛁𝐏) ⅆ𝛀

𝛀 𝛀 𝛀

L'équation générale devient :

⃗ ⅆ𝛛𝛀 + ∫ 𝛅𝐏 𝐟 ⅆ𝛀
∫ 𝛁𝛅𝐏(𝐤𝛁𝐏) ⅆ𝛀 = ∫ 𝛅𝐏(𝐤𝛁𝐏) 𝐧
𝛀 𝛀 𝛀

II. Discrétisation
La discrétisation transforme le domaine continu en un ensemble fini de points (nœuds), ce qui
permet de résoudre l’équation de manière numérique. La section détaille également l’interpolation
de Lagrange, utilisée pour approximer les valeurs entre les nœuds du maillage.
Sens géométrique : Pour simplifier les calculs, on applique une transformation géométrique τ aux
éléments, en changeant de variable. Cette transformation facilite l'application de la MEF en
réduisant la complexité du domaine en un cadre géométrique plus simple. on doit faire un
changement de variable :

y 𝜂

x 𝜉
𝑝(𝑥 ) 𝑝(𝜉 )

Pour modéliser la déformation d’un élément et permettre les calculs entre les nœuds,
l’interpolation de Lagrange est utilisée. Cela permet de déterminer des valeurs
approximatives entre les points connus.

x = a0, + ξa1, + ηa2, + ξηa3,

{ (*)
y = b0, + ξb1, + ηb2, + ξηb3,

Sens mathématique

Pour le premier terme : ∫𝛺 𝛻𝑢 𝛻𝑝 𝑑Ω = ∑ ∫𝛺𝑒 𝛻𝑢 𝛻𝑝 𝑑Ω

Pour le deuxième terme : ∫𝛺 𝑢 𝑄𝑔 𝑑Ω = ∑ ∫𝛺𝑒 𝑢 𝑄𝑔 𝑑Ω 𝑒

On utilise l’élément de référence et on travaille sur un rectangle

D’après l’interpolation de Lagrange on a :

𝑝(𝑥 ) = ∑𝑁𝑖 𝑝𝑖 = 0 𝑎𝑣𝑒𝑐 𝑁𝑖 = 𝑁𝑖 (𝑥 ) = 𝑁𝑖 (𝑥 , 𝑦 )

𝑝(𝜉 ) = ∑𝑁𝑖 𝑝𝑖 = 0 𝑎𝑣𝑒𝑐 𝑁𝑖 = 𝑁𝑖(𝜉 ) = 𝑁𝑖(𝜉 , 𝜂 )
D’après la transformation géométrique :

𝜏(−1, −1) = (𝑥1, 𝑦1)

𝜏 (1, −1) = (𝑥2, 𝑦2 )
{
𝜏 ( 1 , 1 ) = (𝑥3, 𝑦3 )
𝜏 (−1, 1) = (𝑥4, 𝑦4 )

III. Interpolation
On utilise une base polynomiale en lien avec le nombre de nœuds, ici quatre pour chaque
élément quadrilatère. Cette base est déterminée en fonction du triangle de Pascal,
permettant une approximation adéquate de la fonction de potentiel.

On prend la base < b > = < 1 x y xy >

On aura :
1 −1 −1 1 1 1 1 1
1 1 −1 −1 −1 1 1 −1
V = [< b (ξ⃗⃗𝑖 ) >] = V = [ ] 𝑒𝑡 𝑉 −1 = [ ]
1 1 1 1 −1 −1 1 1
1 −1 1 −1 1 −1 1 −1

1
N1 = (1 − ξ − η + ξη)
4
1
N2 = (1 + ξ − η − ξη)
4
1
N3 = (1 + ξ + η + ξη)
4
1
{ N4 = (1 − ξ + η − ξη)
4

Ces fonctions sont fondamentales en MEF, car elles définissent comment chaque point dans un
élément contribue à la solution globale.

𝑝( 𝜉) = ∑𝑁𝑖𝑝𝑖 =< 𝑁 > {𝑝}

∂p
=< N, ξ > {p}
∂ξ
→ ⃗⃗⃗⃗
𝛻𝜉 𝑝 = [ 𝑝, 𝜉 𝑝, 𝜂 ] → ⃗⃗⃗⃗
𝛻𝜉 𝑝 = [ 𝑁, 𝜉 𝑁, 𝜂 ]{𝑝}
𝜕𝑝
=< 𝑁, 𝜂 > {𝑝}
{𝜕𝜂
 1
a0, = 4
(x3, + x4, + x1, + x2, )
x1, = a0, − a1, − a2, + a3, 1
x2, = a0, + a1, − a2, + a3, a1, = 4
(x3, − x4, − x1, + x2, )
{ 1
x3, = a0, + a1, + a2, + a3, a2, = (x3, + x4, − x1, − x2, )
x4, = a0, − a1, + a2, + a3, 4
1
{ a3 = 4
(x3, − x4, + x1, − x2, )

Remplaçons dans (*) on obtient :

Pour le vecteur de transformation 𝑝 (𝜉 ) = 𝑥 , on a :

{ 𝑥 =< 𝑁(𝜉) > {𝑥𝑖 } 𝑦 =< 𝑁(𝜉) > {𝑦𝑖 } Avec i = (1, 2, 3, 4)

𝑥 = [ < 𝑁(𝜉 ) > < 𝑁(𝜉 ) > ]{𝑋𝑖 }

IV. Contribution élémentaire

Chaque élément a sa propre contribution, qui dépend de la transformation géométrique et des
fonctions de forme. Ces contributions sont ensuite assemblées pour former la matrice globale.

On utilise l’interpolation de Lagrange on a :

Et on a : 𝑝(𝑥 ) = 𝑝(𝜉 ) =< 𝑁(𝜉 ) > {𝑝}𝑒

𝑢(𝑥 ) = 𝑢(𝜉 ) = < 𝑁(𝜉 ) > {𝑢}𝑒 = < 𝑢 >𝑒 {𝑁(𝜉 )}

Donc :

∫𝛺𝑒 𝑢 𝑄𝑔 𝑑Ω = ∫𝛺𝑒 < 𝑢 >𝑒 {𝑁(𝜉 )} 𝑄𝑔 𝑑Ω

Avec 𝑑Ω = (𝐽) 𝑑𝜉 𝑑𝜂

∫𝛺𝑒 𝑢 𝑄𝑔 𝑑Ω = ∫𝛺𝑒 < 𝑢 >𝑒 {𝑁(𝜉 )} 𝑄𝑔 𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

Transformation d’intégrale : On utilise la transformation géométrique pour déterminer 𝛻𝜉 p et [𝐽] la

matrice jacobien :

𝜕𝑝 𝜕𝑝 𝜕𝑥 𝜕𝑇 𝜕𝑦 𝜕𝑝 𝜕𝑝 𝜕𝑝 𝜕𝑝
= +
𝜕𝜉 𝜕𝑥 𝜕𝜉 𝜕𝑦 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝜉 𝜕𝑥
→ = 𝜕𝑝
𝜕𝑝 𝜕𝑝 𝜕𝑥 𝜕𝑇 𝜕𝑦 𝜕𝑝 𝜕𝑝 𝜕𝑝
= +
{𝜕𝜂 𝜕𝑥 𝜕𝜂 𝜕𝑦 𝜕𝜂 {𝜕𝜂} [𝜕𝜂 𝜕𝜂] {𝜕𝑦}

𝜕𝑝 𝜕𝑝
𝜕𝜉 𝜕𝜉
→ 𝛻𝜉𝑝 = [𝐽]𝛻𝑝 Avec [J] = [𝜕𝑝 𝜕𝑝
]
𝜕𝜂 𝜕𝜂
 ∂p
=< N, ξ > {xi }
∂ξ
𝜕𝑝
=< 𝑁, 𝜂 > {𝑥𝑖 }
{𝜕𝜂

Donc : [𝐽] = [ < 𝑁, 𝜉 > < 𝑁, 𝜂 > ]{𝑥𝑖 , 𝑦𝑖 } = [𝐺]{𝑋}

Partie tangente :

On peut l’écrire en fonction de G :

⃗⃗⃗⃗
𝛻𝜉 𝑝 = [𝐽] 𝛻𝑝

Or: ⃗⃗⃗⃗
𝛻𝜉 𝑝 = [ < 𝑁, 𝜉 > < 𝑁, 𝜂 > ]{𝑝}𝑒 → ⃗⃗⃗⃗
𝛻𝜉 𝑝 = [𝐺] {𝑝}𝑒

Donc : 𝛻𝑝 = [𝐽]−1 [𝐺] {𝑝}𝑒

Ainsi on aura :

𝑘 ∫Ω 𝛻𝑢 𝛻𝑝 𝑑Ω = 𝑘∑ ∫Ω𝑒 < 𝑢 > 𝑒 [𝐺]𝑇 [𝐽]−𝑇 [𝐽]−1 [𝐺] {𝑝}𝑒 𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

On dérive les fonctions de forme et on trouve :

1 −1 + 𝜂 1− 𝜂 1 + 𝜂 −1 − 𝜂
[G] = [ ]
4 −1 + 𝜉 −1 − 𝜉 1 + 𝜉 −1 − 𝜉

−1 + 𝜂 −1 + 𝜉
1 1 − 𝜂 −1 + 𝜉
[G]T = [ ]
4 1 + 𝜂 1 + 𝜉
−1 − 𝜂 1 − 𝜉

𝑥1 𝑦1
1 −1 + 𝜂 1− 𝜂 1 + 𝜂 −1 − 𝜂 𝑥2 𝑦2
[𝐽] = [𝐺]{𝑋} = [ ]∗[ ]
4 −1 + 𝜉 −1 − 𝜉 1 + 𝜉 −1 − 𝜉 𝑥3 𝑦3
𝑥4 𝑦4

1 (−1 + 𝜂)𝑥1 − (−1 + 𝜂)𝑥2 + (1 + 𝜂)𝑥3 − (1 + 𝜂)𝑥4 (−1 + 𝜂)𝑦1 − (−1 + 𝜂)𝑦2 + (1 + 𝜂)𝑦3 − (1 + 𝜂)𝑦4
= [ ]
4 (−1 + 𝜉)𝑥1 − (1 + 𝜉)𝑥2 + (1 + 𝜉)𝑥3 − (−1 + 𝜉)𝑥4 (−1 + 𝜉)𝑦1 − (1 + 𝜉)𝑦2 + (1 + 𝜉)𝑦3 − (−1 + 𝜉)𝑦4

𝐽11 𝐽12
[J] = [ ]
𝐽21 𝐽22

1 𝐽22 −𝐽12 𝑗11 𝑗12

[J−1 ] = [ ] = [𝑗] = [ ]
𝐽11 𝐽22 − 𝐽12 𝐽21 −𝐽21 𝐽11 𝑗21 𝑗22
On pose :

𝑗11 𝑗12 1 −1 + 𝜂 1− 𝜂 1 + 𝜂 −1 − 𝜂
𝑀 = [𝑗][𝐺] = [ ]∗ [ ]
𝑗21 𝑗22 4 −1 + 𝜉 −1 − 𝜉 1 + 𝜉 −1 − 𝜉

1 (−1 + 𝜂)𝑗11 − (1 − 𝜉)𝑗12 (1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12 (1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12 −(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12
= [ ]
4 (−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22 (1 − 𝜂)𝑗21 − (1 + 𝜉)𝑗22 (1 + 𝜂)𝑗21 + (1 + 𝜉)𝑗22 −(1 + 𝜂)𝑗21 + (1 − 𝜉)𝑗22

(−1 + 𝜂)𝑗11 − (1 − 𝜉)𝑗12 (−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22

1 (1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12 (1 − 𝜂)𝑗21 − (1 + 𝜉)𝑗22
[M]T = [ ]
4 (1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12 (1 + 𝜂)𝑗21 + (1 + 𝜉)𝑗22
−(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12 −(1 + 𝜂)𝑗21 + (1 − 𝜉)𝑗22

𝑚11 𝑚12 𝑚13 𝑚14

1 𝑚21 𝑚22 𝑚23 𝑚24
[M]T [M] = [ ]
16 𝑚31 𝑚32 𝑚33 𝑚34
𝑚41 𝑚42 𝑚43 𝑚44

Avec :

𝑚11 = ((−1 + 𝜂)𝑗11 − (1 − 𝜉)𝑗12)2 + ((−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22) 2

𝑚12 = 𝑚21 = [(1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12] ∗ [(−1 + 𝜂)𝑗11 − ((1 − 𝜉)𝑗12] + [(1
− 𝜂)𝑗21 − (1 + 𝜉)𝑗22] ∗ [(−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22)]

𝑚13 = 𝑚31 = [(1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12] ∗ [(−1 + 𝜂)𝑗11 − (1 − 𝜉)𝑗12] + [(1

+ 𝜂)𝑗21 + (1 + 𝜉)𝑗22] ∗ [(−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22]

𝑚14 = 𝑚41 = [−(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12] ∗ [(−1 + 𝜂)𝑗11 − (1 − 𝜉)𝑗12] + [−(1

+ 𝜂)𝑗21 + (1 − 𝜉)𝑗22] ∗ [(−1 + 𝜂)𝑗21 − (1 − 𝜉)𝑗22]

𝑚22 = ((1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12) 2 + ((1 − 𝜂)𝑗21 − (1 + 𝜉)𝑗22)2

𝑚23 = 𝑚32 = [(1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12] ∗ [(1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12] + [(1 − 𝜂)𝑗21
− (1 + 𝜉)𝑗22 ] ∗ [(1 + 𝜂)𝑗21 + (1 + 𝜉)𝑗22]

𝑚24 = 𝑚42 = [(1 − 𝜂)𝑗11 − (1 + 𝜉)𝑗12] ∗ [−(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12] + [(1

− 𝜂)𝑗21 − (1 + 𝜉)𝑗22] ∗ [−(1 + 𝜂)𝑗21 + (1 − 𝜉)𝑗22]

𝑚33 = ((1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12)2 + ((1 + 𝜂)𝑗21 + (1 + 𝜉)𝑗22) 2

𝑚34 = 𝑚43 = [(1 + 𝜂)𝑗11 + (1 + 𝜉)𝑗12] ∗ [−(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12] + [(1

+ 𝜂)𝑗21 + (1 + 𝜉)𝑗22] ∗ [−(1 + 𝜂)𝑗21 + (1 − 𝜉)𝑗22]

𝑚44 = (−(1 + 𝜂)𝑗11 + (1 − 𝜉)𝑗12)2 + (−(1 + 𝜂)𝑗21 + (1 − 𝜉)𝑗22)2

 𝑚11 𝑚12 𝑚13 𝑚14
𝑘 𝑚21 𝑚22 𝑚23 𝑚24
𝐾𝑒 = 𝑘 ∫𝛺𝑒 [𝑀]𝑇 [𝑀] ∗ 𝑑𝑒𝑡(𝐽) = ∫𝛺𝑒 [ ]
16 (𝐽11 𝐽22 − 𝐽12 𝐽21) 𝑚31 𝑚32 𝑚33 𝑚34
𝑚41 𝑚42 𝑚43 𝑚44

Cas particulier (élément rectangulaire) :

Maintenant on va traiter le cas où on a un rectangle avec des points non définie

𝑎 𝑏
1 −1 + 𝜂 1− 𝜂 1 + 𝜂 −1 − 𝜂 𝑎+1 𝑏
[𝐽] = [𝐺]{𝑋} = [ ]∗[ ]
4 −1 + 𝜉 −1 − 𝜉 1 + 𝜉 −1 − 𝜉 𝑎+1 𝑏+ℎ
𝑎 𝑏+ℎ
1 𝐿 0
→ [𝐽] = [ ]
2 0 ℎ
2 0 𝐿+ℎ
Avec : → [𝐽]−1 = [ 𝐿 2] et → 𝑑𝑒𝑡(𝐽) = 4
0 ℎ

Donc on peut calculer :

−1 + 𝜂 1− 𝜂 1 + 𝜂 −1 − 𝜂
1
𝑀 = [𝑗][𝐺] = [ 𝐿 𝐿 𝐿 𝐿 ]
2 −1 + 𝜉 −1 − 𝜉 1 + 𝜉 −1 − 𝜉
ℎ ℎ ℎ ℎ

Et sa transposée égale :

−1 + 𝜂 −1 + 𝜉
𝐿 ℎ
1 − 𝜂 −1 + 𝜉
1 𝐿 ℎ
[M]T =
4 1 + 𝜂 1 + 𝜉
𝐿 ℎ
−1 − 𝜂 1 − 𝜉
[ 𝐿 ℎ ]

Ce qui implique : 𝑀𝑇 𝑀 = [𝐽−1 ]𝑇 [𝐺]𝑇 [𝐽−1 ] [𝐺]

 1 1 1 1 1 1 1 1
(𝜂 − 1)2 + 2 (𝜉 − 1)2 − (1 − 𝜂)2 − 2 (𝜉² − 1) (𝜂² − 1) + (𝜉² − 1) − (𝜂² − 1) − 2 (𝜉 − 1)2
𝐿2 ℎ 𝐿2 ℎ 𝐿2 ℎ2 𝐿2 ℎ
1 2
1 1 1 1 1 1 1
1 − 𝐿2 (𝜂 − 1) − ℎ2 (𝜉² − 1) (1 − 𝜂)2 + 2 (𝜉 + 1)2 (1 − 𝜂²) − 2 (𝜉 + 1)2 − 2 (1 − 𝜂²) − 2 (1 − 𝜉²)
𝑀𝑇 𝑀 = 𝐿2 ℎ 𝐿2 ℎ 𝐿 ℎ
4 1 1 1 1 1 1 1 1
(𝜂² − 1) + 2 (𝜉² − 1) (1 − 𝜂²) − 2 (𝜉 + 1)2 (𝜂 + 1)2 + 2 (𝜉 + 1)2 − 2 (𝜂 + 1)2 + 2 (1 − 𝜉²)
𝐿2 ℎ 𝐿2 ℎ 𝐿2 ℎ 𝐿 ℎ
1 1 2
1 1 1 1 1 1
[− 𝐿2 (𝜂² − 1) − ℎ2 (𝜉 − 1) − 2 (1 − 𝜂²) + 2 (𝜉² − 1) − 2 (𝜂 + 1)2 − 2 (𝜉² − 1)2 (𝜂 + 1)2 + 2 (1 − 𝜉)2 ]
𝐿 ℎ 𝐿 ℎ 𝐿2 ℎ

𝑘𝑒 = 𝑘 ∫𝛺𝑒 [𝐺]𝑇 [𝐽]−𝑇 [𝐽]−1 [𝐺] 𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

𝑘𝑒 = 𝑘 𝑑𝑒𝑡(𝐽) ∫𝛺𝑒 𝑀𝑇 𝑀𝑑𝜉 𝑑𝜂

𝑘𝑒
𝑙ℎ
= 𝑘
4
1 1 1 1 1 1 1 1
(𝜂 − 1)2 + 2 (𝜉 − 1)2 − (1 − 𝜂)2 − 2 (𝜉² − 1) (𝜂² − 1) + (𝜉² − 1) − (𝜂² − 1) − 2 (𝜉 − 1)2
𝐿2 ℎ 𝐿2 ℎ 𝐿2 ℎ2 𝐿2 ℎ
1 1 1 1 1 1 1 1
1 − 2 (𝜂 − 1)2 − 2 (𝜉² − 1) (1 − 𝜂)2 + 2 (𝜉 + 1)2 (1 − 𝜂²) − 2 (𝜉 + 1)2 − 2 (1 − 𝜂²) − 2 (1 − 𝜉²)
+ ∫ 𝐿 ℎ 𝐿2 ℎ 𝐿2 ℎ 𝐿 ℎ 𝑑𝜉 𝑑𝜂
4 1 1 1 1 1 1 1 1
(𝜂² − 1) + 2 (𝜉² − 1) (1 − 𝜂²) − 2 (𝜉 + 1)2 (𝜂 + 1)2 + 2 (𝜉 + 1)2 − 2 (𝜂 + 1)2 + 2 (1 − 𝜉²)
𝐿2 ℎ 𝐿2 ℎ 𝐿2 ℎ 𝐿 ℎ
1 1 2
1 1 1 1 1 1
[− 𝐿2 (𝜂² − 1) − ℎ2 (𝜉 − 1) − 2 (1 − 𝜂²) + 2 (𝜉² − 1) − 2 (𝜂 + 1)2 − 2 (𝜉² − 1)2 (𝜂 + 1)2 + 2 (1 − 𝜉)2 ]
𝐿 ℎ 𝐿 ℎ 𝐿2 ℎ

D’où la matrice du terme de conduction est :

2 2 −2 1 −1 1 1 2
+ + + −
𝐿2 ℎ2 𝐿2 ℎ2 𝐿2 ℎ2 𝐿2 ℎ2
1 1 2 2 1 2 −1 1
𝑙ℎ 𝐿2 + ℎ2 𝐿2
+ 2
ℎ 𝐿2
− 2
ℎ
−
𝐿2 ℎ2
𝑘𝑒 = 𝑘
6 −1 1 1 2 2 2 −2 1
− − 2 + 2 +
𝐿2 ℎ2 𝐿2 ℎ 𝐿2 ℎ 𝐿2 ℎ2
1 2 −1 1 −2 2 2 2
[ 𝐿2 − ℎ2 −
𝐿2 ℎ2
+
𝐿2 ℎ2
+
𝐿2 ℎ2 ]

𝜕𝑝 𝑙ℎ
𝑘∫ 𝑢 𝑑𝑆 = 𝑘 ℎ(𝑝𝑒𝑥𝑡 − 𝑝𝑖𝑛𝑡) < 𝑢 >𝑒
𝜕𝑛 4

∫𝛺𝑒 𝑄𝑔 𝑢 𝑑Ω = ∫𝛺𝑒 𝑄𝑔 < 𝑢 >𝑒 {𝑁(𝜉)}𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

= 𝑄𝑔 < 𝑢 >𝑒 ∫ {𝑁(𝜉)}𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

(1 − 𝜉 − 𝜂 + 𝜉𝜂)
𝐿ℎ (1 + 𝜉 − 𝜂 − 𝜉𝜂)
∫𝛺𝑒 𝑄𝑔 𝑢 𝑑Ω = 𝑄𝑔 < 𝑢 >𝑒 ∫𝛺𝑒 ( )
4 (1 + 𝜉 + 𝜂 + 𝜉𝜂)
(1 − 𝜉 + 𝜂 − 𝜉𝜂)
 1
𝐿ℎ 1 𝑒
∫𝛺𝑒 𝑄𝑔 𝑢 𝑑Ω = 𝑄𝑔 < 𝑢 > ( )
4 1
1

V. Assemblage
L’assemblage consiste à combiner les contributions de tous les éléments pour former une matrice de
rigidité globale qui représente le système complet. Les termes de conduction et de convection sont
additionnés dans cette matrice. La matrice globale est essentielle pour les calculs finaux car elle
intègre toutes les interactions entre éléments et prend en compte la structure entière du domaine.

−𝑘 ∫𝛺 𝛻𝑢 𝛻𝑝 𝑑Ω + ∫𝛺 𝑄𝑔 𝑢 𝑑Ω = 0

Cette équation contient deux termes, commençons par le premier terme exprimant la rigidité,

P𝑜𝑢𝑟 𝑙’é𝑙é𝑚𝑒𝑛𝑡 𝒆𝟏, nous avons ∶

𝑘11 𝑘12 𝑘13 𝑘14 𝑝1

𝑘21 𝑘22 𝑘23 𝑘24 𝑝2
𝐾𝑒1 = < 𝑢1 𝑢2 𝑢5 𝑢4 > [ ]{ }
𝑘31 𝑘32 𝑘33 𝑘34 𝑝3
𝑘41 𝑘42 𝑘43 𝑘44 𝑝4

𝐾𝑒1 = < 𝑢 > [𝑘1 ] {𝑝}

𝑘11 𝑘12 0 𝑘14 𝑘13 0 𝑝1

𝑘21 𝑘22 0 𝑘24 𝑘23 0 𝑝2
0 0 0 0 0 0 𝑝3
𝐾𝑒1 = < 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢6 >
𝑘41 𝑘42 0 𝑘44 𝑘43 0 𝑝4
𝑘31 𝑘32 0 𝑘34 𝑘33 0 𝑝5
[ 0 0 0 0 0 0] {𝑝6}

Pour l’élément 𝐞𝟐, nous avons ∶

𝑘11 𝑘12 𝑘13 𝑘14 𝑝2

𝑘21 𝑘22 𝑘23 𝑘24 𝑝3
𝐾𝑒2 = < 𝑢2 𝑢3 𝑢6 𝑢5 > [ ]{ }
𝑘31 𝑘32 𝑘33 𝑘34 𝑝6
𝑘41 𝑘42 𝑘43 𝑘44 𝑝5

0 0 0 0 0 0 𝑝1
0 𝑘11 𝑘12 0 𝑘14 𝑘13 𝑝2
0 𝑘21 𝑘22 0 𝑘24 𝑘23 𝑝3
𝐾𝑒2 = < 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢6 >
0 0 0 0 0 0 𝑝4
0 𝑘41 𝑘42 0 𝑘44 𝑘43 𝑝5
[0 𝑘31 𝑘32 0 𝑘34 𝑘33] {𝑝6}
Finalement, on assemble e1 et e2 :

𝐾𝑒1 + 𝐾𝑒2 = < 𝑢1 𝑢2 𝑢2 𝑢4 𝑢5 𝑢6

𝑘11 𝑘12 0 𝑘14 𝑘13 0 𝑝1
𝑘21 𝑘11 + 𝑘22 𝑘12 𝑘24 𝑘14 + 𝑘23 𝑘13 𝑝2
0 𝑘21 𝑘22 0 𝑘24 𝑘23 𝑝3
>∗
𝑘41 𝑘42 0 𝑘44 𝑘43 0 𝑝4
𝑘31 𝑘41 + 𝑘32 𝑘42 𝑘34 𝑘44 + 𝑘33 𝑘43 𝑝5
[ 0 𝑘31 𝑘32 0 𝑘34 𝑘33] {𝑝6}

Pour le deuxième terme exprimant la force,

∫𝛺 𝑄𝑔 𝑢 𝑑Ω = ∫𝛺𝑒 𝑄𝑔 < 𝑢 >𝑒 {𝑁(𝜉)} 𝑑𝑒𝑡(𝐽) 𝑑𝜉 𝑑𝜂

Pour l’élément 𝐞1, nous avons ∶

1
𝐿ℎ 1
𝑄𝑒 1 = < 𝑢1 𝑢2 𝑢5 𝑢4 > ( ) 𝑄𝑔
4 1
1
1
1
1
𝐿ℎ 0
𝑄𝑒 = < 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢6 > 𝑄𝑔
4 1
1
( 0)

Pour l’élément 𝐞𝟐, nous avons ∶

1
2
𝐿ℎ 1
𝑄𝑒 = < 𝑢2 𝑢3 𝑢6 𝑢5 > ( ) 𝑄𝑔
4 1
1
0
1
𝐿ℎ 1
𝑄𝑒 2 = < 𝑢2 𝑢3 𝑢6 𝑢5 > 𝑄𝑔
4 0
1
( 1)

Finalement, on assemble
 1
2
𝑙ℎ 1
𝑄𝑒 1 + 𝑄𝑒 2 = < 𝑢1 𝑢2 𝑢3 𝑢4 𝑢5 𝑢6 > 𝑄𝑔
4 1
2
( 1)

VI. Acquisition des conditions aux limites

Les conditions aux limites définissent ce qui se passe aux bords du domaine. Elles
représentent les contraintes physiques, comme la hauteur ou le flux d’eau.

Type :

1. Dirichlet :

Imposent une valeur fixe pour p(x,y).

Exemple : Une hauteur d’eau constante de 10 cm sur un bord.

𝑝(𝑥, 𝑦 = 0) = 𝑝𝑒𝑥𝑡
{
𝑝(𝑥 = 0, 𝑦) = 𝑝𝑒𝑥𝑡

2. Neumann :
𝝏𝒑
Imposent une valeur pour le flux ( ).
𝝏𝒏

Exemple : Un débit d’eau constant traversant une frontière.

𝑄(𝑥, 𝑦 = 0) = 0
{
𝑄(𝑥 = 0, 𝑦) = 0

VII. Résolution
La résolution consiste à résoudre un système d’équations de la forme :

𝐾𝑥 = 𝐹

• K : Matrice globale, issue des contributions élémentaires. Elle représente les interactions
entre les nœuds.
• x : Vecteur des inconnues. Il contient les valeurs du potentiel hydraulique p aux nœuds.
• F : Vecteur global de charges, représentant les sources externes et les conditions aux limites.
 programmation sur MATLAB :

%======================== Parameters =========================%

• Lx = 2; Ly = 2;

• nElx = 20; % Number of elements according to x

• nEly = 20; % Number of elements according to y

• s = @(x,y) 0; % Source term (e.g., fluid injection/production rate)

• nPg = 4; % Number of Gauss Points

• K = 1e-12; % Permeability of the medium (m^2)

• % function Mesh

• [coor, conn, nN, nDofEl, nDofT, nNEl, nEl] = generateMesh(Lx, Ly, nElx, nEly);

• % conditions aux limites generateMesh : Fonction qui crée

le maillage.
• [B1, B2, B3, B4] = nddlpn(Lx, Ly, coor, nElx, nEly); nddlpn : Fonction qui identifie les
nœuds situés sur les bords .
• bcdof = [B1; B2; B3]; • bcdof : Regroupe les indices des
nœuds soumis à des conditions aux
• bcval = [10*ones(size(B1)); 0*B2; 0*B3]; limites.
• bcval : Valeurs imposées
• % calcule des matrices et vecteurs

• [Xi, Eta, W] = gaussPoints(nPg); •

• F = zeros(nDofT, 1);

• kT = zeros(nDofT);

• for iel = 1:nEl

• % Elementary Mat Vec gaussPoints : Retourne les

coordonnées et poids des points de
• xEl = coor(conn(iel,:), 1); Gauss.
shapFunctions : Calcule les
• yEl = coor(conn(iel,:), 2); fonctions de forme (N) et leur gradient
(G).
• ke = zeros(nDofEl);

• Fe = zeros(nDofEl, 1);

• for pg = 1:nPg
• xi = Xi(pg); eta = Eta(pg); w = W(pg);

• [N, G] = shapFunctions(xi, eta);

• J = G * [xEl yEl];

• x = N * xEl; y = N * yEl;

• j = inv(J);

• ke = ke + (K * (G' * j' * j * G)) * det(J) * w;

• Fe = Fe + N' * s(x, y) * det(J) * w;

Assemblage : Ajoute les contributions
• end locales dans les matrices globales :
• nods : Nœuds de l'élément
• % assemblage global courant.
• index : Correspondance entre
• nods = conn(iel, :); indices locaux et globaux.
• assVecMat : Fonction
• index = globalIndex(nods); d'assemblage.
• [F, kT] = assVecMat(F, kT, Fe, ke, index);

• end

• % imposition des conditions aux limites

• for idof = 1:length(bcdof)

Pour chaque degré de liberté
• dof = bcdof(idof);
soumis à une condition aux limites
• F(dof) = bcval(idof); :

• kT(dof, :) = 0; • Impose la valeur dans le

vecteur F.
• kT(dof, dof) = 1; • Modifie la matrice kT pour
refléter une condition de
• End Dirichlet.

• % partie solution du système linéaire

• p = kT \ F;

• % post traitement (visualisation des résultats)

• X = reshape(coor(:, 1), [nElx+1, nEly+1]);

• Y = reshape(coor(:, 2), [nElx+1, nEly+1]);

• P = reshape(p, [nElx+1, nEly+1]);

 • Reshape : Réorganise les
• figure; données pour correspondre au
maillage 2D.
• contourf(X, Y, P); • Contourf : Trace la distribution
des pressions en utilisant des lignes
• colorbar; de contour.
• title('Pressure Distribution');
• Colorbar : Ajoute une légende
pour les valeurs de pression.
• xlabel('x'); ylabel('y’);

• % Used Functions

• function [coor, conn, Nnoeud, Nddlpelt, Nddlt, Nnpelt, nelt] = generateMesh(Lx, Ly,
nelx, nely)

• % Nodes Quadrilateral Element Mesh

• Nx = nelx + 1; Ny = nely + 1;

• hx = Lx / nelx; hy = Ly / nely;

• nelt = nelx * nely;

• Nnoeud = Nx * Ny; % Total number of nodes

• Nnpelt = 4; % Number of nodes per element

• Nddlpn = 1; % Number of dof per node

• Nddlpelt = Nddlpn * Nnpelt; % Number of dof per element

• Nddlt = Nddlpn * Nnoeud; % Total number of dof

• k = 1;

• for j = 1:Ny

• for i = 1:Nx

• coor(k, 1) = (i - 1) * hx;

• coor(k, 2) = (j - 1) * hy;

• coor(k, 3) = 0.0;

• k = k + 1;

• end
• end

• start = 1;

• for j = 1:nely

• for i = 1:nelx

• conn(start, 1) = (j - 1) * Nx + i;

• conn(start, 2) = (j - 1) * Nx + i + 1;

• conn(start, 3) = j * Nx + i + 1;

• conn(start, 4) = j * Nx + i;

• start = start + 1;

• end

• end

• end

• function [B1, B2, B3, B4] = nddlpn(Lx, Ly, coor, nelx, nely)

• % Bords Definition

• % Left edge : x = 0

• B1 = find(abs(coor(:,1)) < Lx / (10 * nelx));

• % Bottom edge y = 0

• B2 = find(abs(coor(:,2)) < Ly / (10 * nely));

• % Right edge x = Lx

• B3 = find(abs(coor(:,1) - Lx) < Lx / (10 * nelx));

• % Top edge y = Ly

• B4 = find(abs(coor(:,2) - Ly) < Ly / (10 * nely));

• end

• function index = globalIndex(nods)

• % Global indicies

• n = length(nods); % = 4
• index = zeros(n, 1);

• for nod = 1:n

• index(nod) = nods(nod);

• end

• end

• function [F, kT] = assVecMat(F, kT, Fe, ke, index)

• % Assembly

• F(index) = F(index) + Fe;

• kT(index, index) = kT(index, index) + ke;

• end

• function [xi, eta, w] = gaussPoints(npg)

• if npg == 1

• xi = 0; eta = 0; w = 4;

• elseif npg == 4

• x1 = -1/sqrt(3); x2 = 1/sqrt(3);

• w1 = 1; w2 = 1;

• xi(1) = x1; eta(1) = x1; w(1) = w1 * w1;

• xi(2) = x1; eta(2) = x2; w(2) = w1 * w2;

• xi(3) = x2; eta(3) = x1; w(3) = w2 * w1;

• xi(4) = x2; eta(4) = x2; w(4) = w2 * w2;

• end

• end

• function [N, G] = shapFunctions(xi, eta)

• p = @(xi, eta) [1 xi eta xi * eta];

• pn = [p(-1, -1); p(1, -1); p(1, 1); p(-1, 1)];

• invPn = inv(pn);

• N = p(xi, eta) * invPn;

• pxi = @(xi, eta) [0 1 0 eta];

• gxi = pxi(xi, eta) * invPn;

• peta = @(xi, eta) [0 0 1 xi];

• geta = peta(xi, eta) * invPn;

• G = [gxi; geta];

• end
4-Resultats :
En considerant que :

- B1=bord gauche
- B2 =bord du bas
- B3= bord de droite
- B4= bord du haut

1er cas: B1=10MPa ; B2=B3=0MPa ;

Nous avons appliqué une préssion de 10Mpa au bord droit : C’est là ou le fluide
a été injecté ; et une pression de 0Mpa à gauche et en bas: c’est la sortie du
fluide . On trouve que les lignes de courant sont exponentielles dans la figure ci-
dessous :

Distribution de la pression
2
9
1.8
8
1.6
7
1.4
6
1.2
y (m)

5
1

4
0.8

0.6 3

0.4 2

0.2 1

0
0 0.5 1 1.5 2
x (m)

2eme cas : B1=0; B3=0

- F=0

Ici on a appliqué une préssion de 10 en B1 et 0 en B3 ;

Les bords supérieurs et inférieurs sont libres(pas de conditions aux limites) :

c’est comme ci ils ont été isolés ;

Nous trouvons alors des lignes de courants linéaires suivant x et verticales

sur la figure suivante :

Distribution de la pression
2
9
1.8
8
1.6
7
1.4
6
1.2
y (m)

5
1

4
0.8

0.6 3

0.4 2

0.2 1

0
0 0.5 1 1.5 2
x (m)

- F=1
En ajoutant une force F=1 on trouve que la pression est plus elévé au
centre
 11
Distribution de la pression x 10
2
4.5
1.8
4
1.6
3.5
1.4
3
1.2
y (m)

2.5
1

2
0.8

0.6 1.5

0.4 1

0.2 0.5

0
0 0.5 1 1.5 2
x (m)

3eme cas : B2=B4=B3=0

- F=0
Dans ce cas les bords droit ,bas,haut sont de pressions nulles et l’entré du
fluide est au bord gauche ; Donc ici le fluide peut sortir dans tous les 3
bors de pression nulle ;
On trouve des lignes de courant circulaire
 Distribution de la pression
2
9
1.8
8
1.6
7
1.4
6
1.2

5
1

4
0.8

0.6 3

0.4 2

0.2 1

0
0 0.5 1 1.5 2

- F=1

En appliquant une force constante F=1, la pression est plus forte au centre .
 11
Distribution de la pression x 10
2

1.8 2.5

1.6

1.4 2

1.2
1.5
y (m)

0.8
1
0.6

0.4
0.5
0.2

0 0
0 0.5 1 1.5 2
x (m)

- 4eme cas : B1=10 ; B3=0

En prenant une perméabilité hétérogène avec

K = (x, y) 1e-12 * (1 + 0.5 * sin (2 * pi * x) .* cos(2 * pi * y));

Ici nous avons pris une permeabilité variable pour pour illustrer le cas du sol car
ce dernier est hétérogène ; il est composé du sable , le gravier, l’argile… etc ; Et
tous ont leur permeabilité différente

Sur la figure de droite nous observons la distribution de la permeabilté ; nous

avons pris une fonction trigonometrique pour illustrer les vides dans le
milieu ;donc dans la partie blanche on a une permeabilité élevée dans le fluide
passe plus facilement que dans la partie noire ou la permeabilité est faible

Sur la figure de gauche nous observons la distribution de pression on constate

que les lignes isobares deviennent curviligne.
 -13
Distribution de pression Distribution de la permeabilite x 10
2 2
14
1.8 1.8
13
1.6 1.6
12
1.4 1.4
11
1.2 1.2
y (m)

y (m)
10
1 1

0.8 0.8 9

0.6 0.6 8

0.4 0.4 7

0.2 0.2
6

0 0
0 0.5 1 1.5 2 0 1 2
x (m) x (m)

Conclusion :
L’étude de la loi de Darcy, combinée à des outils tels que la méthode des
éléments finis, est essentielle en conception mécanique pour comprendre et
maîtriser les interactions entre fluides et structures poreuses. Que ce soit dans
le développement de systèmes de refroidissement, de filtres industriels ou de
matériaux innovants comme les mousses métalliques ou les composites, la
capacité à modéliser les écoulements dans les milieux poreux est cruciale.

Cette approche favorise également l’innovation en permettant la conception de

solutions sur mesure, optimisées pour répondre à des défis techniques
spécifiques. Par exemple, dans les domaines de l’aéronautique et de
l’automobile, la gestion des flux à travers des matériaux poreux contribue à
l’amélioration de la performance et à la réduction des coûts énergétiques.

Enfin, la maîtrise de ces phénomènes ouvre la voie à des matériaux intelligents

et durables, renforçant ainsi l’impact de la mécanique dans la résolution de
problématiques complexes et dans la transition vers des technologies plus
respectueuses de l’environnement. L’intégration de ces concepts dans les
projets d’ingénierie est un véritable moteur d’innovation et de progrès
technique