Modélisation numérique en génie chimique

GCH2535
Méthode des éléments finis (MEF)

Diapositives adaptées de :
David Vidal
Bruno Blais
François Bertrand
Plan de cours

2
Les tenseurs : objet mathématique utile à la physique
• Ordre 0 : scalaire
→ Température
→ Concentration

• Ordre 1 : vecteur
→ Vitesse
→ Force
→ Déplacement

• Ordre 2 : matrice
→ Contraintes
→ Permeabilité

• Ordre 3 : essayez d’imaginer ça si vous le pouvez …

3
Notation des tenseurs

4
Utilisation du vecteur nabla

5
Utilisation du vecteur nabla

6
Méthode des éléments finis - Objectif

Objectif :
▪ Transformer une équation aux dérivées partielles continue valable sur un
domaine continu en un système linéaire à N équations pour N inconnues,
dont les équations sont associées à un domaine discret appelé maillage.

𝜕𝑢 𝜕2𝑢 𝜕𝑢 𝑘11 𝑘12 ⋯ 𝑘1𝑛 𝑢1 𝑓1

𝐿 𝑢, , 2 , , … + 𝑓 = 0
𝜕𝑥 𝜕𝑥 𝜕𝑡 𝑘21 𝑘22 ⋯ 𝑘2𝑛 𝑢2 𝑓2
⋮ =
+ Conditions Frontières ⋮ ⋮ ⋱ ⋮ ⋮
+ Conditions Initiales (transitoire) 𝑘𝑛1 𝑘𝑛2 ⋯ 𝑘𝑛𝑛 𝑢𝑛 𝑓𝑛

Les inconnues ne sont plus

nécessairement associées à un
maillage cartésien structuré (grille)

Méthode :
▪ Écrire sous forme discrète (i-1, i, i+1 …), non plus l’EDP, mais une forme
« affaiblie » de cette équation par l’entremise de fonctions d’interpolation
exprimées en chaque élément du maillage
2
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

3
Technique d’affaiblissement par la méthode des résidus pondérés

▪ Prenons l’exemple en conduction thermique suivant:

Γ −𝜅𝛻2𝑇 𝒙 = 𝑓(𝒙), ∀𝒙 ∈ Ω
avec
Ω 𝑇 𝒙 = 𝑇0(𝒙), ∀𝒙 ∈ Γ
ou
𝜕𝑇 𝒙
= 𝛼, ∀𝒙 ∈ Γ
𝜕𝑛

Définition : nous appelons résidu (noté Res), l’expression mathématique suivante

reliée au problème étudié :

𝑅𝑒𝑠 𝑇 = 𝑘𝛻2𝑇 𝒙 + 𝑓 𝒙
Ce résidu s’annule quand T(x) est la solution.

4
Méthode des résidus pondérés

1. Pondération du résidu par une fonction-test :

𝜑 𝒙 × 𝜅𝛻2𝑇 𝒙 + 𝑓 𝒙 = 0, ∀𝒙 ∈ Ω et ∀𝜑 𝒙

fonction-test résidu

2. Intégration sur le domaine :

න 𝜑 𝒙 ∙ 𝜅𝛻2𝑇 𝒙 + 𝑓 𝒙 𝑑Ω = 0, ∀𝒙 ∈ Ω et ∀𝜑 𝒙
Ω

−𝜅 න 𝜑 𝒙 ∙ 𝛻2𝑇 𝒙 𝑑Ω = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Ω

Formulation intégrale FORTE

5
Méthode des résidus pondérés

3. Intégration par parties (formule de Green en 2D et 3D) :

−𝜅 න 𝜑 𝒙 ∙ 𝛻2𝑇 𝒙 𝑑Ω = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω , ∀𝒙 ∈ Ω et ∀𝜑 𝒙
Ω Ω

𝜕𝑇 𝒙
𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω − න 𝜑 𝒙 ∙ 𝑑Γ = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Γ 𝜕𝑛 Ω

Formulation intégrale FAIBLE

Motivations :
1. Réduction de l’ordre maximum des dérivées présentes
2. Introduction « naturelle » des conditions frontières

Rappel : intégration par parties en 1D

L
dT ( x ) L
d ( x )
 ( x ) dx = −  T ( x ) dx +  ( x ) T ( x )  0
L

0 dx 0 dx
d 2T ( x ) d ( x ) dT ( x ) dT ( x ) 
L
L L

 ( x )
0 dx 2
dx = − 
0 dx dx
dx +  ( x )
 dx  0

6
Formes forte et faible

Particularité de la méthode des éléments finis (MEF) :

Discrétiser, non pas l’EDP associée au problème (forme forte), mais une forme
« affaiblie » de cette équation

Vocabulaire : cette forme est dénommée diversement :

➢ Forme/formulation faible
➢ Forme/formulation intégrale
➢ Forme/formulation variationnelle
➢ …
Motivation : affaiblir pour réduire certaines contraintes mathématiques
(p.ex. discontinuités) empêchant l'utilisation d'outils classiques
pour sa résolution

Conséquence : la solution d’une forme faible correspond à une solution

approchée ou « faible » en termes de continuité

7
Méthode des résidus pondérés

4. Introduction des conditions frontières :

𝜕𝑇 𝒙
𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω − න 𝜑 𝒙 ∙ 𝑑Γ = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Γ 𝜕𝑛 Ω

▪ Traitement des conditions de Neumann:

𝜕𝑇 𝒙
= 𝛼, ∀𝒙 ∈ Γ
𝜕𝑛

𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω − 𝛼 න 𝜑 𝒙 𝑑Γ = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Γ Ω

▪ Traitement des conditions de Dirichlet, deux possibilités: 𝜕𝑇 𝒙

1. Introduction d’un flux inconnu 𝑞(𝒙) pour ∀𝒙 ∈ Γ: −𝑘𝜑 𝒙 ∙ = 𝜑 𝒙 𝑞(𝒙)
𝜕𝑛

𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω + න 𝜑 𝒙 𝑞(𝒙)𝑑Γ = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Γ Ω

2. Élimination pure et simple du terme et réintroduction manuelle lors de la phase

d’assemblage (voir plus loin) 9
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

10
Bonne préparation de la géométrie

▪ Elle permet de faciliter la réalisation du

maillage et d'abaisser le coût des calculs

=0

symétrie
▪ Quatre points sont à examiner :

𝜕𝑇 𝑥, 𝑦
𝜕𝑥
– le choix de la dimensionnalité du domaine
considéré et/ou du maillage
– l'exploitation des symétries des modèles (p.ex. y
symétrie
axisymétrie?)
– la suppression des détails géométriques 𝜕𝑇 𝑥, 𝑦
=0
superflus 𝜕𝑦

– la délimitation des régions où seront affectés

des modèles de comportements ou
d'environnements différents (p.ex. matériau x
multicouche)

15
Types d’éléments finis (linéaires)
Nomenclature
nœud

1D
Barre
face

2D arête

Triangle Rectangle

3D
Tétraèdre Hexaèdre Pentaèdre Pyramide
ou prisme
16
Types d’éléments finis (quadratiques)
Nomenclature

1D
Barre

nœud d’arête
utilisé pour
l’interpolation
2D quadratique

Triangle Rectangle

3D
Tétraèdre Hexaèdre Pentaèdre Pyramide
ou prisme
17
Autres types d’éléments finis
▪ Éléments finis cubiques : plus riches contenant 2 nœuds par arête (en plus des nœuds
aux coins) – utilisation plus rare

Exemple pour un tétraèdre

▪ Éléments finis quadratiques à arêtes curvilignes pour mieux épouser les géométries
courbes :

Exemple pour un tétraèdre

18
Choix du type et de la richesse des éléments

1. De façon générale, les rectangles et hexaèdres sont plus précis que les triangles,
tétraèdres et prismes à nombre de nœuds égal, car leurs fonctions d’interpolation sont
plus riches et peuvent donc représenter les gradients de façon plus régulière

Interpolation d’un champ u(x,y)

(Gendre, 2013)
par deux éléments triangulaires
linéaires (a) et par un élément
rectangulaire linéaire (b)

2. En revanche, les triangles et tétraèdres permettent de créer des maillages mieux

adaptés aux géométries complexes à l’aide d’algorithmes de génération de maillage
automatisée plus performants (p.ex. triangulation de Delaunay)

3. Les éléments quadratiques représentent un bon compromis entre precision et temps de

calcul
− par défaut, COMSOL (version 4) utilise des éléments quadratiques pour la plupart
des types de problèmes physiques (modules)
− deux exceptions: les modules Fluid Flow et Chemical Species Transport qui utilisent
par défaut des éléments linéaires → problèmes dominés par l’advection
19
Raffinement du maillage requis par la physique du problème

134 éléments 1056 éléments 3394 éléments

tCPU < 1 s tCPU  1 s tCPU  3 s

1038 éléments tCPU  1 s

vԦ

▪ Il convient de raffiner le maillage dans les zones de forts gradients du champ d’intérêt (ici
la vitesse du fluide dans une cavité entrainée) pour améliorer la précision à faible coût
20
Raffinement/enrichissement du maillage requis par la géométrie du problème

σ (Pa)

Entrefer
Zoom (120 )
(Comsol, 2013)
՜
𝑉𝑡
Élém. linéaire
Élém. quadratique

Erreur (%)

Éléments par arc de 90

21
Table des connectivités des éléments

▪ Les numéros de nœuds globaux associés à chaque élément sont stockés dans la table
des connectivités (liste de voisins)

Domaine Numérotation des nœuds Numérotation des éléments

No. élément No. des nœuds

1 10 2 34
2 4 33 14
3 1 22 16
… … … …
52 35 19 32
22
Qualité des éléments – éviter les éléments distordus!

 Éléments de bonnes proportions

Meilleur
 Éléments trop aplatis ou étirés

▪ Surviennent souvent lorsque l'on tente de mailler des détails géométriques de

petite taille avec des éléments trop grands

▪ Peuvent nuire à la pertinence du résultat de plusieurs façons :

− fonctions de base d'allure irrégulière → polarisation du champ d’intérêt selon
certaines directions
− possible dégradation du conditionnement du système d'équations résolu

23
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

24
Discrétisation de la forme faible intégrale par un ensemble d’éléments

𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω − 𝛼 න 𝜑 𝒙 𝑑Γ = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Γ Ω

Supposons ici pour

simplifier 𝛼 = 0
Ω
𝜅 න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω = න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω
Ω Ω Γ
𝑁𝑒
Ωj
Ω = ራ Ω𝑗
𝑗=1 Ne = nombre
Ω𝑙 ∩ Ω𝑚 = ∅ d’éléments finis
𝑙≠𝑚

𝑁𝑒 𝑁𝑒

𝜅 ෍ න 𝛻𝜑 𝒙 ∙ 𝛻𝑇 𝒙 𝑑Ω𝑗 = ෍ න 𝜑 𝒙 ∙ 𝑓 𝒙 𝑑Ω𝑗 ∀𝒙 ∈ Ω et ∀𝜑 𝒙
𝑗=1 Ω𝑗 𝑗=1 Ω𝑗
Formulation intégrale faible
On souhaite maintenant se discrétisée
“débarrasser” de ce gradient… 25
Introduction des fonctions d’interpolation en chaque élément

▪ En chaque élément contenant Nn nœuds, nous pouvons interpoler le champ

d’intérêt par :
𝑁𝑛 Température au nœud i 𝑁𝑛
Polynôme 𝒙 − 𝒙𝒋
de Lagrange 𝑇 𝒙 = ෍ 𝑇𝑖 𝜓𝑖 (𝒙) 𝜓𝑖 𝒙 = ෑ
𝒙𝒊 − 𝒙𝒋
𝑖=1 Fonctions de base 𝑗=1
𝑗≠𝑖

Ωk
1) Polynôme de Lagrange de degré 1 - exemple 1D:
𝒙k 𝒙k+1

𝜓𝑖(𝒙) 𝑇(𝒙)
𝜓0 𝜓1 𝜓2 𝜓3 𝒙 − 𝒙𝟏
1 𝜓0 𝒙 = 𝑇3
𝒙 𝟎 − 𝒙𝟏
= 𝟏−𝒙
𝑇1
𝒙 − 𝒙0 𝑇2
𝜓1 𝒙 = 𝑇0
0 𝒙1 − 𝒙0
𝒙 = 𝒙 𝒙
0 1 2 3 0 1 2 3
etc…
Ω1 Ω2 Ω3 Ω1 Ω2 Ω3
26
Introduction des fonctions d’interpolation en chaque élément

Ωk
2) Polynôme de Lagrange de degré 2 - exemple 1D:
𝒙k 𝒙k+½ 𝒙k+1
𝜓𝑖(𝒙)
𝜓 0 𝜓½ 𝜓1 𝜓1½ 𝜓2
1
𝒙 − 𝒙½ 𝒙 − 𝒙1
𝜓0 𝒙 =
𝒙𝟎 − 𝒙½ 𝒙𝟎 − 𝒙1
3 1
0 = 2 𝒙𝟐 − 𝒙 +
1 2 2
0 1
2 1 12 2 𝒙
𝒙 − 𝒙𝟎 𝒙 − 𝒙1
Ω1 Ω2 𝜓½ 𝒙 =
𝒙½ − 𝒙𝟎 𝒙½ − 𝒙1
𝑇(𝒙) = −4 𝒙𝟐 − 𝒙
𝑇2
𝑇1½ 𝒙 − 𝒙𝟎 𝒙 − 𝒙½
𝜓1 𝒙 =
𝑇½ 𝒙1 − 𝒙 𝟎𝒙1 − 𝒙½
𝑇1 1
= 2 𝒙𝟐 − 𝒙
𝑇0 2
1
0 1
2 1 12 2 𝒙
etc…
Ω1 Ω2 27
Propriétés de l’interpolation de Lagrange
𝑁𝑛

𝑇 𝒙 = ෍ 𝑇𝑖 𝜓𝑖 (𝒙)
𝑖=1

𝜓𝑖 𝒙𝑗 = 𝜹𝑖𝑗 = ቊ
0 si 𝑗 ≠ 𝑖
1 si 𝑗 = 𝑖 
՜ 𝑇 𝒙𝒊 = 𝑇1 × 0 + ⋯ + 𝑇𝑖 × 1 + ⋯ + 𝑇𝑁𝑛 × 0 = 𝑇𝑖

𝑁𝑛

෍ 𝜓𝑖(𝒙) = 1
𝑖=1

𝑁𝑛

𝛻𝑇 𝒙 = ෍ 𝑇𝑖 𝛻𝜓𝑖 (𝒙)
𝑖=1

Constantes!
28
Introduction des fonctions d’interpolation dans la formulation intégrale faible
discrétisée + remplacement de fonctions-test (Méthode de Galerkine)

𝑁𝑒 𝑁𝑒

𝜅 ෍ න 𝛻𝜑 ∙ 𝛻𝑇 𝑑Ω𝑗 = ෍ න 𝜑 ∙ 𝑓 𝑑Ω𝑗
𝑗=1 Ω𝑗 𝑗=1 Ω𝑗

𝑁𝑛
Pour un élement Ω𝑗 fixé 𝛻𝑇 = ෍ 𝑇𝑖 𝛻𝜓𝑖
𝑖=1

𝑁𝑛

𝜅 න 𝛻𝜑 ∙ ෍ 𝑇𝑖 𝛻𝜓𝑖 𝑑Ω𝑗 = න 𝜑 ∙ 𝑓 𝑑Ω𝑗

Ω𝑗 𝑖=1 Ω𝑗

On choisit 𝜑 tel que… 𝜑 = 𝜓1 , … , 𝜓𝑵𝒏

(méthode de Galerkine) 𝛻𝜑 = 𝛻𝜓1 , … , 𝛻𝜓𝑁𝑛

𝑁𝑛

𝜅 න 𝛻𝜓1 , … , 𝛻𝜓𝑁𝑛 ∙ ෍ 𝑇𝑖 𝛻𝜓𝑖 𝑑Ω𝑗 = න 𝜓1 , … , 𝜓𝑁𝑛 ∙ 𝑓 𝑑Ω𝑗

Ω𝑗 𝑖=1 Ω𝑗

29
Évaluations des fonctions d’interpolation sur un élément de référence

Motivation : afin d’éviter l’évaluation de chacune des fonctions de base pour

chaque élément, on préfère les évaluer dans un élément de
référence pour lequel ses fonctions sont toujours les mêmes
𝑦 𝒙2(𝑥2, 𝑦2) 𝜁

𝑁𝑛 𝜓෠ 1 𝜉, 𝜁 = 1 − 𝜉 − 𝜁
𝒙 − 𝒙𝒋 Ωj (0,1)
𝜓𝑖 𝒙 = ෑ 𝜓෠ 2 𝜉, 𝜁 = 𝜉
𝒙𝒊 − 𝒙𝒋 𝒙3(𝑥3, 𝑦3)
𝑗=1
෡ 𝜓෠ 3 𝜉, 𝜁 = 𝜁
𝑖≠𝑗 𝒙1(𝑥1, 𝑦1) Ω
𝑥 (0,0) (1,0) 𝜉

Élément réel Élément de référence

Comment ? le changement de base s’obtient de la façon suivante:

𝑥 𝜉, 𝜁 = 𝑥1 𝜓෠ 1 𝜉, 𝜁 + 𝑥2𝜓෠ 2 𝜉, 𝜁 + 𝑥3 𝜓෠ 3 𝜉, 𝜁 = 𝑥1 + 𝑥2 − 𝑥1 𝜉 + 𝑥3 − 𝑥1 𝜁
𝑦 𝜉, 𝜁 = 𝑦1 𝜓෠ 1 𝜉, 𝜁 + 𝑦2𝜓෠ 2 𝜉, 𝜁 + 𝑦3 𝜓෠ 3 𝜉, 𝜁 = 𝑦1 + 𝑦2 − 𝑦1 𝜉 + 𝑦3 − 𝑦1 𝜁

30
Formulation intégrale faible discrétisée avec fonctions d’interpolation exprimées
sur l’élément de référence

𝑁𝑛

𝜅 න 𝛻𝜓1 , … , 𝛻𝜓𝑁𝑛 ∙ ෍ 𝑇𝑖 𝛻𝜓𝑖 𝑑Ω𝑗 = න 𝜓1 , … , 𝜓𝑁𝑛 ∙ 𝑓 𝑑Ω𝑗

Ω𝑗 𝑖=1 Ω𝑗

𝑁𝑛

𝜅 න 𝛻𝜓𝑘 ∙ ෍ 𝑇𝑖 𝛻𝜓𝑖 𝑑Ω𝑗 = න 𝜓𝑘 ∙ 𝑓 𝑑Ω𝑗 ∀𝑘 ∈ 1, … , 𝑁𝑛

Ω𝑗 𝑖=1 Ω𝑗

Ω𝑗 ෡
Ω
𝑑Ω𝑗 = det 𝐉 𝑑Ω ෡
𝛻𝜓𝑙 = 𝛻𝜓෠ 𝑙 𝐉 −1

𝑁𝑛

𝜅 න 𝛻 𝜓෠ 𝑘 𝐉 −1
∙ ෍ 𝑇𝑖𝛻𝜓෠ 𝑖 𝐉 −1 ෡ = න 𝜓෠ 𝑘 ∙ 𝑓 det 𝐉 𝑑Ω
det 𝐉 𝑑Ω ෡ ∀𝑘 ∈ 1, … , 𝑁
𝑛
෡
Ω ෡
Ω
𝑖=1

𝐉 : matrice Jacobienne de transformation de l’élément j

det 𝐉 : déterminant de 𝐉 31
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

32
Assemblage du système matriciel - Prenons un exemple avec Nn=2 et Ne = 4…

0 1 12 13 14 L
x1 x2 x3 x4 x5

𝑁𝑒 𝑁𝑛 𝑁𝑒

𝜅 ෍න 𝛻 𝜓෠ 𝑘 𝐉 −1
∙ ෍ 𝑇𝑖𝛻𝜓෠ 𝑖 𝐉 −1 ෡ = ෍ න 𝜓෠ 𝑘 ∙ 𝑓 det 𝐉 𝑑Ω
det 𝐉 𝑑Ω ෡
෡ ෡
𝑗=1 Ω 𝑖=1 𝑗=1 Ω

valable sur tout le domaine Ω ∀𝑘 ∈ 1, … , 𝑁𝑛

Maintenant, pour un det 𝐉 = 𝑥k − 𝑥k+1= ℎ

élément j donné 1
𝐉 −1 =
ℎ

𝑁
1 1 𝑛
𝜅 න 𝛻𝜓෠ 𝑘 ∙ ෍ 𝑇𝑖 𝛻𝜓෠ 𝑖 ℎ𝑑 Ω
෡ = න 𝜓෠ 𝑘 ∙ 𝑓 ℎ 𝑑Ω
෡ ∀𝑘 ∈ 1, … , 𝑁𝑛
෡
Ω ℎ ℎ ෡
Ω
𝑖=1

𝑁𝑛 = 2
2
𝜅
න 𝛻𝜓෠ 𝑘 ෍ 𝑇𝑖 𝛻𝜓෠ 𝑖 𝑑Ω
෡ = ℎ න 𝜓෠ 𝑘 ∙ 𝑓 𝑑Ω
෡ , ∀𝑘 ∈ 1,2
ℎ Ω෡ ෡
Ω 33
𝑖=1
Assemblage du système matriciel (suite)

2
𝜅
න 𝛻𝜓෠ 𝑘 ෍ 𝑇𝑖 𝛻𝜓෠ 𝑖 𝑑Ω
෡ = ℎ න 𝜓෠ 𝑘 ∙ 𝑓 𝑑Ω
෡ , ∀𝑘 ∈ [1,2]
ℎ Ω෡ ෡
Ω
𝑖=1

Nous avons: 𝜓෠ 1 = 1 − 𝜉 𝛻𝜓෠ 1 = −1

𝜓෠ 2 = 𝜉 𝛻 𝜓෠ 2 = 1
Prenons:

𝜅 𝜅 1 −1
෠ ෠ ෡
𝐴𝑗 = 2 න 𝛻 𝜓 𝑘 𝛻 𝜓 𝑖 𝑑 Ω =
ℎ Ω෡ ℎ −1 1

න 1 − 𝜉 ∙ 𝑓(𝜉) 𝑑𝜉 Disons: 𝑓 = 0 c−à−d 𝑘𝛻2𝑇 = 0

𝜉 0
𝑏𝑗 = න 𝜓෠ 𝑘 ∙ 𝑓 𝑑Ω
෡= 𝑏𝑗 =
෡
Ω 0
න 𝜉 ∙ 𝑓(𝜉) 𝑑𝜉
𝜉

𝑇
Donc à l’élément j, nous pouvons écrire: 𝐴𝑗 𝑇 𝑗 = 𝑏 𝑗
𝑗+1
34
Assemblage final – Matrice diffusive

𝜅 𝜅 1 −1
෠ ෠ ෡=
𝐴𝑗 = 2 න 𝛻𝜓𝑘𝛻𝜓𝑖 𝑑Ω
ℎ Ω෡ ℎ −1 1

1 −1 0 0 0
𝜅 −1 1+1 −1 0 0
𝐴 = 0 −1 1+1 −1 0
ℎ 0 −1 1 + 1 −1
0
0 0 0 −1 1

1 −1 0 0 0
𝐴1 𝐴2 𝐴3 𝐴4 𝜅 −1 2 −1 0 0
𝐴 = 0 −1 2 −1 0
ℎ 0 0 −1 2 −1
0 0 0 −1 1

35
Assemblage final – Vecteur élémentaire

𝑏𝑗 , 1 0
𝑏𝑗 = =
𝑏𝑗 , 2 0
𝑏1

𝑏2
𝑏1,1 0
𝑏1,2 + 𝑏2,1 0
𝑏3
𝑏 = 𝑏2,2 + 𝑏3,1 = 0
𝑏3,2 + 𝑏4,1 0
𝑏4 𝑏4,2 0

36
Assemblage final – Système matriciel et imposition des conditions de Dirichlet

1 −1 0 0 0 𝑏1,1 0
𝜅 −1 2 −1 0 0 𝑏1,2 + 𝑏2,1 0
𝐴 = 0 −1 2 −1 0 𝑏 = 𝑏2,2 + 𝑏3,1 = 0
ℎ 0 0 −1 2 −1 𝑏3,2 + 𝑏4,1 0
0 0 0 −1 1 𝑏4,2 0

1 −1 0 0 0 𝑇1 0 Supposons:
−1 2 −1 0 0 𝑇2 0 • T1 = Ta
0 −1 2 −1 0 𝑇3 = 0 • T5 = Tb
0 0 −1 2 −1 𝑇4 0
0 0 0 −1 1 𝑇5 0

1 0 0 0 0 𝑇1 Ta
−1 2 −1 0 0 𝑇2 0
0 −1 2 −1 0 𝑇3 = 0
Pour les conditions de 0 0 −1 2 −1 𝑇4 0
Neumann, cf. diapo 9 0 0 0 0 1 𝑇5 Tb 37
Assemblage final – Système matriciel et imposition des conditions de Dirichlet

1 0 0 0 0 𝑇1 Ta
−1 2 −1 0 0 𝑇2 0
0 −1 2 −1 0 𝑇3 = 0 Forme non-compacte
0 0 −1 2 −1 𝑇4 0
0 0 0 0 1 𝑇5 Tb

2 −1 0 𝑇2 Ta
−1 2 −1 𝑇3 = 0 Forme compacte
0 −1 2 𝑇4 Tb

Rappel:

d 2T  2 − 1 0  T2  Ta + (h  ) f 2 
2

−  2 = f dans 
dx MDF − 1 2 − 1 T  =  (h2 ) f

T ( x = a) = Ta   3    3

 0 − 1 2  T4  Tb + (h 2  ) f 4 
T ( x = b) = Tb

38
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

39
Quel algorithme doit-on choisir pour résoudre un système linéaire?

Réponse:
▪ On choisit l’algorithme en fonction de la taille N du système et du type de matrice (p.ex.
diagonale, creuse,…) à résoudre.

▪ Méthodes disponibles:
1. Inversion de matrice et formule de Cramer (calculs de déterminant)
2. Méthodes directes “avancées”(p.ex. factorisation LU, algorithme de Thomas,…)
3. Méthodes itératives (p.ex. méthodes de gradient conjugué, Gauss-Seidel, GMRES,
BiCGStab…)

N= 5 10 25 105 106
Formule de Cramer 1300 flops 4108 flops 1015 ans  

Factorisation LU 100 flops 900 flops  O(msec)  O(heure)  O(mois)

Méthode itérative - - -  O(min)  O(heure)

40
Comparaison MDF vs. MEF

Différences finies (rappel) Éléments finis

− EDP + C.F. & C.I.. − = Forme FORTE

− Obtention forme FAIBLE intégrale
− Création de la grille cartésien − Création du maillage du domaine
➢ Nœuds équidistants ➢ Nœuds
➢ Éléments (connectivité)
− Obtention de l’équation discrète − Discrétisation de la forme faible intégrale sur
➢ Schémas « prêts à l’emploi » chaque élément + fonction d’interpolation
➢ Idem pour les C.F.

− Construction du système − = Assemblage (introduction des CF)

− Résolution du système − Résolution du système
− Post-traitement − Post-traitement

41
Toujours verifier la convergence de son schéma numérique!!!

Erreur

10 1 0.1 0.01 0.001

h (m)

▪ Utiliser 3 maillages de tailles de maille différentes (grossier, moyen & fin)

▪ Pour calculer l’erreur, si vous n’avez pas de solution analytique, comparez

alors les résultats obtenus par rapport au maillage le plus fin ou à
l’extrapolation de Richardson

▪ Un ordre 2 est souhaitable 42

Précision des éléments finis

1
Erreur L2 = ‫׬‬ 𝑇𝑖 − 𝑇(𝑥𝑖 ) 2 𝑑Ω = 𝑇𝑖 − 𝑇(𝑥𝑖) 2
Ω Ω

𝑁𝑒𝑙𝑒
ℎele ℎele
1
= ෍ 𝑇𝑖 − 𝑇(𝑥𝑖 ) 2 ≤ 𝐶ℎ𝑘+1
𝑁𝑒𝑙𝑒
𝑖=1

avec ℎ = max ℎ𝑒𝑙𝑒 (taille de maille) et k = ordre des éléments utilisés

𝑒𝑙𝑒

Exemple:
Erreur

k +1  2

1 / hD
10 1 0.1 0.01 0.001
h (m) 43
 Hypothèse
Extrapolation de Richardson généralisée principale :
(lorsque la solution analytique n’est pas connue) Les solutions obtenues
avec les maillages fin
i=1
p1 estimé et grossier tendent de
façon asymptotique
Tdx −Tdx
fin grossier vers une solution
TRich. =Tdx + p
fin dxgrossier i
−1
dxfin
Autrement, la
non → pi = pi+1
Procédure procédure itérative ne
→ i=i+1
itérative converge pas !
pi+1 = pi ?
oui
(%)

ϵ = 𝒪 dxpi+1
Tdx −TRich.

1. Ordre de convergence “p”

TRich.

2. Estimé TRich.de la solution

grossier numérique pour dx → 0

fin
ϵ=

44
Remarques finales sur la méthode des éléments finis
• Permet de résoudre une variété d’EDP:
• elliptique (probl. stationnaire de diffusion): 𝛼𝛻 2 𝑈 + 𝑓 = 0
𝜕𝑈
• parabolique (probl. instationnaire de diffusion): − 𝛼𝛻 2 𝑈 + 𝑓 = 0
𝜕𝑡
• hyperbolique (probl. instationaire d’advection): 𝜕𝑈
− 𝑣 ∙ 𝛻𝑈 + 𝑓 = 0
𝜕𝑡
• mixte transport (probl. instationnaire 𝜕𝑈
d’advection-diffusion): − 𝑣 ∙ 𝛻𝑈 + 𝛼𝛻 2 𝑈 + 𝑓 = 0
𝜕𝑡

• Avantages:
1. Flexibilité → maillages non-structurés
2. Prise en compte “naturelle” des conditions frontières
3. Contrôle de l’erreur

• Inconvénients:
1. Mise en place et programmation plus ardues que MDF
2. Parallélisation par décomposition de domaines ardue → système matriciel
3. Génération de maillage pour des géométries 3D peut devenir compliquée

45
Méthode des Volumes Finis (MVF)
n
Théorème de flux-divergence (ou de Green-Ostrogradski):

ම 𝜵 ∙ 𝑭 𝑑Ω = ඾ 𝑭 ∙ 𝑑𝒏
Ω Γ
La divergence d’un champ vectoriel F sur un volume Ω
est égale au flux de ce champ à travers la frontière Γ de
ce volume (intégrale de surface).

Prenons maintenant: 𝑭 = 𝑘𝜵𝑇 ՜ ම 𝜵 ∙ 𝑘𝜵𝑇𝑑Ω = ඾ 𝑘𝜵𝑇 ∙ 𝑑𝒏

Ω Γ 𝑘=constante
Volume i
( ici cas “1D”) ම 𝑘𝜵2 𝑇𝑑Ω = ඾ 𝑘𝜵𝑇 ∙ 𝑑𝒏
Ω Γ Discrétisation
𝑇𝑖 − 𝑇𝑗 du flux par une
Centroïde
ම 𝑘𝜵2 𝑇𝑑Ω =඾𝑘 ∙ 𝑑𝒏 différence finie
× × Ω Γ ℎ
du volume i j

h La méthode consiste donc à sommer tous les flux entrant

par toutes les frontières d’un volume et de discrétiser ce
𝑇𝑖 − 𝑇𝑗 flux au moyen d’un schéma de différences finies
𝜵𝑇 =
ℎ
46
Lectures suggérées

Chapitre 8
47
Autres lectures pour les plus passionnés

https://onlinelibrary.wiley.c https://www.amazon.com
om/doi/book/10.1002/978 /The-Finite-Element-
0470510858 Method-
Engineering/dp/04864118
18
48