Institut International de Technologie

Département de Génie Civil

GC2

Cours de Méthode des Eléments Finis

ZOUARI Bassem

[email protected]

2022/2023

1
 Sommaire

Chapitre 1 : Introduction
Chapitre 2 : Formulation d’un problème physique et Principe de la
méthode des éléments finis
Chapitre 3 : Les formulations variationnelles
Chapitre 4 : Les éléments et leurs espaces d’interpolation
Chapitre 5 : La discrétisation
Chapitre 6 : Les éléments 2D : déformations planes et contraintes planes et
Les éléments 2D structuraux : plaques

2
CHAPITRE 1 : INTRODUCTION

L‘activité de conception est un processus complexe de création. Elle consiste à élaborer un produit
ou un système conformément aux exigences d’un client, et dans le respect de certaines règles ou
normes, ce qui revient à borner le domaine de création du produit. Elle se doit en outre de garantir
la rentabilité financière de l’entreprise.

L’élaboration d’un cahier des charges fonctionnel (CdCF) permet d’appréhender la complexité du
projet par une structuration en fonctions et contraintes auxquelles sont associés des critères
d’appréciation, en précisant leur niveau et leur flexibilité. Les contraintes physiques comme le
comportement mécanique, thermique, électrique ou électromagnétique du produit sont l’objet
d’attention particulière de la part de l’ingénieur. Celui-ci doit disposer d’outils et de méthodes pour
concevoir un produit capable de satisfaire aux contraintes physiques.

Si la formulation théorique de la majorité des problèmes physiques a été réalisée dans les
siècles précédents, la résolution de ses problèmes reste toujours d’actualité. Avec le
développement des outils informatiques et des mathématiques appliquées, on trouve aujourd’hui
une large gamme de logiciels capables de fournir des solutions approchées à la majorité des
problèmes physiques. Reste à connaître le degré de précision des solutions fournis. Cette tâche
n’est pas une mince affaire. En effet elle nécessite une connaissance approfondie de la physique
du problème ainsi que les méthodes de résolutions adoptées.

On s’intéresse dans ce cours à une méthode de résolution numérique très utilisée dans le
domaine de la mécanique et du génie civil à savoir la méthode des éléments finis. Le but de ce
cours n’est pas de former des développeurs de la méthode des éléments finis mais plutôt
d’exposer le principe de la méthode en général et d’avoir des utilisateurs avertis. En effet, une
mauvaise utilisation d’un outil de calcul par éléments finis peut aboutir à des résultats erronés.
L’origine de l’erreur sont multiples :
 Une mauvaise modélisation du problème réel : mauvaises conditions aux limites,
mauvaises estimation du chargement, négligence d’un phénomène physique comme la
variation de la température qui induit des contraintes supplémentaires, problème en
grandes transformations,
 Une mauvaise utilisation des options de calcul : intégration réduite, éléments linéaires,
maillage grossier, ….

https://www.youtube.com/watch?v=fmKKFkREu8Q

Exemples de mauvaises utilisations des calculs par éléments finis

Ci-dessous un exemple où une utilisation non avertie d’un outil de calcul par éléments finis peut
fournir des résultats erronés.
Il s’agit d’un problème de flexion d’une poutre console soumise à son propre poids et une
distribution uniforme de force surfacique sur sa surface supérieure.

3
La poutre a pour longueur 1000 mm et largeur 100mm et hauteur 100mm. Le matériau constitutif
présente un module de Young de 30 000 MPa, un coefficient de poisson de 0.3 et une masse
volumique de 3e-6 kg/mm3.
L’accélération de la pesanteur est prise égale à 10 m/s2. La force surfacique présente une valeur
de 0.01 MPa.
PL4
La solution RDM fournie une valeur de la flèche maximale f max 
8EI
La flèche maximale dans ce cas est f max  0.65mm

On a modélisé ce problème très simple en utilisant le logiciel de calcul par éléments finis Abaqus.
Les résultats suivants représentent la déformée de la poutre pour différents maillages, différents
types d’éléments (à interpolation linéaire et à interpolation quadratique) et différents modes
d’intégration (intégration complète et intégration réduite.

Taille de maille 100 mm intégration réduite interpolation linéaire :

4
On remarque que la valeur de la flèche maximale est 100 fois plus grandes que la flèche donnée par
la solution RDM par contre les contraintes longitudinales sont proches de zéro de l’ordre de 10-12.
On ne voit pas l’évolution linéaire de la contrainte avec la hauteur comme dans le cas de la flexion.
Ce résultat est complètement faux.

Taille de maille 100 mm intégration réduite interpolation quadratique :

On remarque que la déformée de la poutre n’est pas réaliste. On observe des modes parasites de
déformations qui sont les modes à énergie nulle dus à l’intégration réduite. La valeur de la flèche
maximale est de l’ordre de 5 fois plus que la flèche donnée par la solution RDM par contre les

5
contraintes longitudinales présentent une évolution linéaire avec la hauteur mais uniquement proche
de l’encastrement. Ce résultat est complètement faux.
Taille de maille 100 mm intégration complète interpolation linéaire :

On remarque que la valeur de la flèche maximale est égale à 0.787 mm qui est proche de la solution
donnée par la théorie de la RDM. Les contraintes longitudinales présentent une évolution linéaire
en fonction de la hauteur
Ce résultat présente une très bonne amélioration par rapport au cas où on a une intégration réduite.
Taille de maille 100 mm intégration complète interpolation quadratique :

6
On remarque que la valeur de la flèche maximale est égale à 0.641 mm qui est très proche de la
solution donnée par la théorie de la RDM. Les contraintes longitudinales présentent une évolution
linéaire en fonction de la hauteur.
Ce résultat présente une très bonne amélioration par rapport au cas où on a une intégration réduite et
par rapport à une interpolation linéaire.

Taille de maille 25 mm intégration réduite interpolation linéaire :

7
On remarque que la valeur de la flèche maximale est égale à 0.69 mm qui est très proche de la
solution donnée par la théorie de la RDM. Les contraintes longitudinales présentent une évolution
linéaire en fonction de la hauteur.
Ce résultat présente une très bonne amélioration par rapport aux maillages grossiers. On remarque
que la structure devient plus flexible.
Taille de maille 25 mm intégration complète interpolation linéaire :

La solution est très proche de la solution analytique.

Remarques :

 Plus le maillage est fin plus on s’approche de la solution théorique.

 L’intégration complète donne des solutions plus proches de la solution analytique que
l’intégration réduite dans la majorité des cas mais avec un coup de calcul plus élevé.
 L’interpolation quadratique est meilleure que l’interpolation linéaire mais avec un coup de
calcul plus élevé (temps CPU)

CPU : en anglais Central Processing Unit (en français unité centrale de traitement, UCT)

Attention : le temps de calcul évolue d’une façon exponentielle avec le nombre de nœuds.
Parfois on est obligé de réduire le nombre d’éléments pour obtenir une solution approchée avec un
temps de calcul raisonnable.

On peut contourner la finesse du maillage par faire un maillage fin uniquement aux endroits où il y
a une forte variation de la contrainte.

8
 Maillage grossier : mauvaise solution Maillage fin : demande beaucoup de CPU

Maillage à taille variable : bon rapport entre qualité de la solution et temps CPU

Un maillage trop fin peut avoir un coût de calcul trop important. Un maillage grossier peut donner
de mauvais résultats. Il peut sous-estimer les contraintes dans certains points de la structure. Une
solution est de faire un maillage avec une taille variable et n’affiner que dans les endroits à fort
gradient de contraintes.

9
CHAPITRE 2 : PRINCIPE DE LA METHODE DES ELEMENTS
FINIS
2.1 Formulation d’un problème physique

La plupart des problèmes de physique peuvent se formuler de la façon suivante :

Trouver un champ (scalaire, vectoriel, ou tensoriel) u(M, t) satisfaisant à un ensemble
d'équations aux dérivées partielles et d'équations ordinaires en tout point M d'un domaine 
et à tout instant t, et respectant des conditions aux limites (éventuellement fonction du temps)
sur la frontière  du domaine.
Remarques :
- Les conditions aux limites sont des relations ou des valeurs imposées à u et/ou à ses
dérivées sur la frontière  .
- Si les conditions sont imposées sur u, on parle de conditions aux limites de type
Dirichlet (exemples : encastrement, appui simple, …)
- Si les conditions sont imposées sur les dérivées du champ u, on parle de conditions aux
limites de type Neumann (exemples : pression, forces de contact, …).
- Il existe aussi des conditions aux limites de type mixte qui font intervenir des relations
entre le champ u et ses dérivées (exemples : appui élastique, convection de chaleur).
- Si le temps t n'apparaît pas comme variable, on dit que le problème est stationnaire,
sinon c'est un problème d'évolution (séisme).

2.1.1 Exemple 1: Le problème thermique stationnaire

Le champ inconnu est un champ scalaire : la température T (M). L'équation différentielle est :
  T M   f M   M  
où f(M) est une fonction donnée (distribution volumique de chaleur) et λ conductivité
thermique du matériau occupant le domaine Ω.

Les conditions aux limites sont par exemple:

- Un champ de température stationnaire imposé sur la partie de frontière 1 :
T (N) = Tf (N)  N  1

- un champ de flux thermique stationnaire  f(N) imposé sur la partie de frontière  2 :

 T  N   n N    f  N   N   2
où n (N) est la normale unitaire extérieure en N à la frontière.

n
Tf

f(M) 1
Ω

Φf
 2

Figure 2.1 : Représentation du problème thermique

10
Remarques :
1)- Pour que la solution soit entièrement déterminée, il faut que 1   2   et
1   2   , c'est à dire qu'il faut imposer des conditions aux limites sur toute la frontière.
2)- Dans un logiciel de calcul par éléments finis spécialisé dans la mécanique, ne pas mettre
de conditions aux limites sur une partie de la frontière est interprété comme une partie sur
laquelle l’effort exercé par le milieu extérieur est nulle. En thermique, la condition similaire
correspond à une paroi adiabatique.

2.1.2 Exemple 2: Le problème élastique linéaire d'évolution

Pour simplifier on se place dans le cas des petites perturbations (petits déplacements et petites
déformations). Le champ inconnu est un champ vectoriel : le champ de déplacement uM , t  en
tout point M de  et à tout instant.
Les équations sont :
2 u
div σ  f   M , t  (Principe fondamental de la dynamique)
t2
où f est une distribution volumique de forces (forces centrifuges, poids, …)

Les conditions aux limites sont :

- Un champ de déplacements u f N , t  (éventuellement nul), fonction du temps, imposé sur la
partie de frontière 1 :
uN , t   u f N , t   N  1 ,
- Un champ de forces surfaciques q f  N , t  (éventuellement nul), fonction du temps, imposés sur
la partie de frontière  2 :
σ N , t   n N   q f  N , t   N   2
où n(N) est le vecteur normal extérieur en N à la frontière.
avec σ, ε respectivement les tenseurs de contraintes et de déformation en tout point de .

La loi de comportement relie le tenseur de contraintes au tenseur de déformation. En faisant

l’hypothèse des petits déplacements et des petites déformations, la loi de comportement
s’écrit : σ   Tr ε I  2 ε . ε  uM , t    T uM , t  (Petites déformations), et  les
1
2
coefficients de Lamé, I le tenseur identité, Tr l’opérateur trace,  l’opérateur gradient, T

l’opérateur de transposition.

n
uf

1
Ω
f(M)
qf
 2

Figure 2.2 : Représentation d’un problème mécanique 3D

11
Remarque :
Compte tenu de la loi de comportement et de la définition de , la dernière condition aux limites se
ramène à des conditions sur les dérivées normales de uM , t  sur la frontière.

2.1.3 Exemple 3 : Barre en traction-compression

Une barre AB de longueur L est soumise à un chargement linéique de traction-compression p(x)
suivant toute sa longueur et une force concentrée F à son extrémité (B). L’autre extrémité (A) est
supposée fixe (figure 2.3). La barre a pour section S. Le module de Young du matériau constituant
la barre est E.

p(x) F

A B

Figure 2.3 : Barre en traction-compression

Le champ inconnu est un champ scalaire qui est le déplacement u(x) en chaque point suivant l’axe
de la barre.
dN ( x)
L’équation d’équilibre local est :  p ( x)  0  x    0, L 
dx
Où N(x) est l’effort normal dans la section située à x.
u (0)  0 conditions aux lim ites type déplacement

Les conditions aux limites sont :  du
 N ( L)  ES dx ( L)  F conditions aux lim ites type effort
du
La loi de comportement est : N ( x)  ES
dx

2.2 Les méthodes approchées : méthode des éléments finis

Si le modèle physique choisi pour modéliser le problème réel nous fournit les équations et les
conditions aux limites, les mathématiques, par contre, sont souvent impuissantes pour en donner
une solution analytique. Les mathématiques prouvent parfois l'existence et l'unicité de la solution.
Il existe plusieurs méthodes de résolution de systèmes d’équations différentielles ou aux dérivées
partielles avec des conditions aux limites, parmi lesquelles on peut citer la méthode des différences
finis, la méthode des volumes finis et enfin la méthode des éléments finis qui fait l’objet de ce
cours. On assiste aujourd’hui à la naissance de nouvelles méthodes numériques comme les
méthodes particulaires (exemple la méthode SPH, Smooth Particule Hydrodynamique, …).

2.3 Principe de la méthode des éléments finis

La plupart des problèmes physiques rencontrés par l’ingénieur n’ont pas de solutions
analytiques vu la complexité de la forme des domaines qu’occupent les pièces (piston, bielle,
pale de turbine, dalle courbe, pont …) et de la nature du chargement et des conditions aux
limites. Pour cela on se dirige vers des solutions approchées numériques.
Soit une pièce occupant un domaine de l’espace . La résolution du problème physique par la
méthode des éléments finis consiste à :
~
Chercher une solution approchée de la solution exacte sous la forme d'un champ FM, t 
défini par morceaux sur des sous domaines de figure.

12
 Sous domaines

Figure 2.4 : Exemple d’une subdivision en sous domaines d’une pièce (maillage)

Les n sous domaines i doivent être tels que :

n
~ ~
 i   et i   j    i  j
i 1

~
où  i désigne l'intérieur de i . Autrement dit, les  i sont une partition de.
~
Les champs fi M, t  , définis sur chaque sous domaines sont des champs choisis parmi une
famille arbitraire de champs (généralement polynomiaux).
La famille de champs locaux est appelée espace des fonctions d'interpolation de l'élément.
~
La famille de champs globaux FM, t  , obtenus par juxtaposition des champs locaux est appelée
espace des fonctions d'interpolation du domaine .
Le champ dans chaque sous domaine  i est déterminé par un nombre fini de valeurs du
champ (ou de valeurs de ses dérivées) en des points choisis arbitrairement dans le sous
domaine, et appelés noeuds. Le champ local est une interpolation entre les valeurs aux noeuds.
Le sous domaine muni de son interpolation est appelé élément.

Chercher une solution par éléments finis consiste donc à déterminer quel champ local on
~
attribue à chaque sous domaine pour que le champ global FM, t  obtenu par juxtaposition de
ces champs locaux soit proche de la solution du problème.

Parmi les contraintes qu'on impose à la solution approchée cherchée, il y a souvent au moins
une continuité simple (C0 ) à la frontière entre les sous domaines.
La figure 2.5 montre une solution approchée discontinue d'un champ scalaire sur un domaine
de dimension 1. La famille de champs locaux est la famille des champs constants par
morceaux.
La figure 2.6 montre une solution approchée continue C0 d'un champ scalaire sur un domaine
 de dimension 1. La famille de champs locaux est la famille des champs polynomiaux de
degré 1.

13
La figure 2.7 montre une solution approchée continue Cl d'un champ scalaire sur un
domaine  de dimension 1. La famille de champs locaux est la famille des champs
polynomiaux de degré 3.

Figure 2.5 : Solution approchée discontinue

Figure 2.6 : Solution approchée continue C0

Figure 2.7 : Solution approchée continue C1

La qualité de la solution approchée dépend de :

- la division en sous domaines (nombre et dimensions des sous domaines),
- le choix de la famille de champs locaux dans chaque sous domaine,
- les conditions de continuité qu'on impose aux frontières des sous domaines (C0,
Cl ,….).
Une fois ces choix faits, il reste à rechercher, une combinaison de champs locaux qui
satisfait approximativement les équations.

14
Pour résoudre un problème par la méthode des éléments finis, on procède donc par étapes
successives :
1. On se pose un problème physique sous la forme d'une équation différentielle ou aux
dérivées partielles à satisfaire en tout point d'un domaine, avec des conditions aux
limites sur le bord  .
2. On construit une formulation intégrale du système différentiel à résoudre et de ses
conditions aux limites : C'est la formulation variationnelle du problème.
3. On divise  en sous domaines :C'est le maillage. Les sous domaines sont appelés
mailles.
4. On choisit la famille de champs locaux, c'est à dire à la fois la position des noeuds
dans les sous domaines et les polynômes (ou autres fonctions) qui définissent le
champ local en fonction des valeurs aux noeuds (et éventuellement des dérivées),
appelé fonction d’interpolation. La maille complétée par ces informations est alors
appelée élément.
5. On ramène le problème à un problème discret : C'est la discrétisation. En effet, toute
solution approchée est complètement déterminée par les valeurs aux noeuds des
éléments. Il suffit donc de trouver les valeurs à attribuer aux noeuds pour décrire
une solution approchée. Le problème fondamental de la méthode des éléments finis
peut se résumer en deux questions :
(a) Comment choisir le problème discret dont la solution est «proche» de la
solution exacte?
(b) Quelle signification donner au mot «proche»?
6. On résout le problème discret : C'est la résolution.
7. On peut alors construire la solution approchée à partir des valeurs trouvées aux
noeuds et en déduire d'autres grandeurs : C'est le post-traitement (Par exemple, en
élasticité, quand on a obtenu une solution approchée du champ des déplacements, on en
déduis le champ des tenseurs de déformation, puis le champ du tenseur des
contraintes).
8. On visualise et on exploite la solution pour juger de sa qualité numérique et juger
si elle satisfait les critères du cahier des charges : C'est l'exploitation des résultats.

Les étapes 1,2,3,4 et 5 sont souvent rassemblées sous le nom de prétraitement. On

s’intéressera dans ce cours uniquement à ces étapes.
Le travail de ces différentes étapes est assisté par les logiciels. Il reste que pour maîtriser
leur utilisation, il est indispensable de comprendre les fondements de la méthode,
notamment les phases 3 et 4, ne serait-ce que pour comprendre et choisir intelligemment
parmi les options qu'ils proposent.

15
CHAPITRE 3 : LES FORMULATIONS VARIATIONNELLES

Pour résoudre un système différentiel (aux dérivées partielles si  est une surface ou un volume),
modélisant un système physique, il faut le mettre sous une forme intégrale, appelée aussi forme
variationnelle ou encore forme faible. Ces formulations peuvent être déduites du système
différentiel par des raisonnements mathématiques ou par des raisonnements physiques. On traitera à
titre d'illustration des exemples en fin de chapitre.
Pour être suffisamment général, on utilise les notations suivantes :
- On note u(M) le champ inconnu (u est noté en caractères gras par souci de généralité, le
champ inconnu pouvant être scalaire vectoriel ou tensoriel).
- On symbolise le système différentiel par un opérateur différentiel D devant être nul en tout
point d'un domaine .
- On symbolise les conditions aux limites sur le bord par un opérateur C devant être nul sur la
frontière   (l'opérateur C peut être différentiel ou non, selon que les conditions aux limites portent
sur les dérivées de u ou non).
Par exemple, on a :
d 2 f x  df
D: f(x)  2
2  3x
dx dx
 f 0  
C : f(x)  df  
   x 1 
  d x  
On verra d'autres exemples en fin de chapitre.
Tout problème peut donc s'énoncer ainsi :
Trouver le champ u(M) défini sur  tel que

D(u(M)) = 0  M  

C(u(N)) = 0  N   

3.1 Rappel de quelques résultats d'analyse fonctionnelle

Dans des espaces fonctionnels de fonctions définies sur un domaine , avec des conditions sur les
fonctions de cet espace et des conditions sur  qu'on supposera satisfaites, on peut définir un
produit scalaire entre deux fonctions f et g :

f, g   f M  g M  dM


où
f(M) • g(M) = f(M) g(M) : si f et g sont à valeurs scalaires
f(M) • g(M) = f(M)  g(M) : si f et g sont à valeurs vectorielles
f1 g1
f  f2 et g  g2 f, g  

f1 g1  f 2 g 2  f 3 g3dM
f3 g3
f(M) • g(M) = f(M)  g(M) : si f et g sont à valeurs tensorielles du second ordre
16
D'autre part, un produit scalaire a la propriété suivante :
f , g  0  g M   f =0
En utilisant cette propriété, il vient :

 f M  gM  d

 0  gM   f M   0 (3.1)

3.2 Formes variationnelles

En utilisant (3.1), une nouvelle formulation du problème est donc

- Trouver le champ u tel que :

 M 

D(u(M)) dM= 0   M 

- avec :
C(u(N)) = 0  N  
 M  est appelée fonction test.
Cette formulation est une formulation variationnelle «équivalente» au système différentiel initial.
On obtient d'autres formulations variationnelles en transformant les intégrales. En effet, l'opérateur
D fait intervenir des opérateurs tels que gradient, divergence, rotationnel, laplacien, et on donne
dans la section suivante des formules pour transformer les intégrales.
Ces transformations d'intégrales font apparaître des intégrales sur les frontières de , qu'on appelle
intégrales de bord, dans lesquelles on peut tenir compte de tout ou partie des conditions aux
limites : le nombre de conditions aux limites diminue.
D'autre part, après chaque transformation d'intégrale, l'inconnue u apparaît avec un ordre de
dérivation diminué de 1. L'intérêt de cette particularité apparaîtra plus loin.
du dg
C g d l d l   C u d l d l   u g A
B

Finalement, après d'éventuelles transformations d'intégrales, le problème peut se mettre sous la

forme variationnelle générale suivante :
- Trouver le champ u tel que :

AMuMB(M   M  (3.2)
- avec C '(u(N)) = 0  N   

ou A et B sont des opérateurs produisant des intégrales sur  et sur   portant sur M  et u(M)
et leurs dérivées, et où C ' est un opérateur produisant les conditions aux limites restantes.

3.3 Transformations d'intégrales

La validité des formules qui suivent est soumise à des conditions de régularité des champs
(scalaires, vectoriels ou tensoriels) définis sur  que nous supposerons satisfaites.

3.3.1 Transformations d'intégrales sur un volume

On note V le volume, S la frontière du volume, et n la normale unitaire extérieure à S.
Dans les formules qui suivent :

17
 - g est un champ scalaire quelconque défini sur V,
- V est un champ vectoriel quelconque défini sur V,
- T est un champ tensoriel du second ordre quelconque défini sur V,
- S est un champ tensoriel du troisième ordre quelconque défini sur V.

Si u est un champ scalaire : grad u est un vecteur, u est un scalaire.

 V  grad u dv    u div V dv   u V  n ds
V V S
(3.3)

 V
g  u dv    grad g  grad u dv 
V  g grad u  n ds
S
(3.4)

si u est un champ vectoriel : grad u est un tenseur du second ordre, div u est un scalaire, rot u est
un vecteur, u est un vecteur.

 V
T  grad u dv    u  div T dv 
V  u  T  n ds
S
(3.5)

 V
g div u dv    u  grad g dv 
V  g u  n ds S
(3.6)
V  rot u dv    u  rot V dv   V  u   n ds
 (3.7)
V V S

 V   u dv    grad V  grad u dv   V  grad u  n ds

V V S
(3.8)

si u est un champ tensoriel du second ordre : grad u est un tenseur du troisième ordre, div u est un
vecteur, u est un tenseur du second ordre.

 V
S  grad u dv    u  div S dv 
V  u  S  n ds S
(3.9)

 V  div u dv
V
   u  grad V dv   V  u  n ds
V S
(3.10)

 T   u dv 
V
  grad T  grad u dv   T  grad u  n ds
V S
(3.11)

3.3.3 Transformations d'intégrales sur une courbe

On note C la courbe, A et B sont les points frontières de la courbe. Dans les formules qui suivent
- g est un champ scalaire défini sur C,
- V est un champ vectoriel tridimensionnel défini sur C.

si u est un champ scalaire :

du dg
 dl   u d l   u g A
B
g
C dl C dl
(3.12)

(c’est la formule classique d’intégration par parties).

si u est un champ vectoriel tridimensionnel :

du dV
 dl   u dl  V uA
B
V
C dl C dl
(3.13)

18
3.4 Exemple 1 : Le problème thermique stationnaire

On suppose que le domaine  est un volume V de frontière S et de normale extérieure n. Le champ

inconnu T est un scalaire.
Le problème s'énonce ainsi :
- Trouver le champ de températures T tel que :
T (M) = f(M) M  V
où f(M) est une source de chaleur volumique connue.
- avec les conditions aux limites :

T(N) =Tf (N)  N  ST

grad T(N)  n(N) = f(N)  N  Sf

(On impose un champ de température stationnaire sur la partie de frontière ST et un champ

de flux thermique stationnaire sur la partie de frontière Sf )
En utilisant (5.1), on obtient une première forme variationnelle :
- Trouver le champ de température T tel que :

   T  f  d v  0  (M) défini sur V

V

- avec les mêmes conditions aux limites, à savoir :

T(N) =Tf (N)  N  ST

grad T(N)  n(N) = f(N)  N  Sf

En utilisant la formule (3.4), on obtient une seconde forme variationnelle :

- Trouver le champ de températures T tel que :

  V
grad   grad T d v    grad T  n d s
S
  V
 f dv  0 

- avec les mêmes conditions aux limites.

Or, l'intégrale de bord peut se décomposer :

 S
 grad T  n d s  ST
 grad T  n d s   Sf
 grad T  n d s

En tenant compte de la condition aux limites de flux sur Sf, il vient :

S
 grad T  n d s   ST
 grad T  n d s  
Sf
 f d s

Finalement, la seconde forme variationnelle du problème est

- Trouver le champ de températures T tel que :
  V
grad   grad T d v  ST
 grad T  n d s   Sf
 f d s   V
f dv 

A(, T) B(
- avec la seule condition aux limites

19
 T (N) = Tf(N)  N  ST

Le problème thermique stationnaire se met bien sous la forme donnée en (3.2).

On peut obtenir une troisième forme variationnelle en utilisant la formule (3.3) :
- Trouver le champ de températures T tel que :

 V
T  d v   S
T grad   n d s   ST
 grad T  n d s 

 Sf
 f d s   V
 f dv  0 

- avec la seule condition aux limites

T (N) = Tf (N)  N  ST

En tenant compte de la condition aux limites de température sur ST, il vient :

 T grad   n d s   T grad   n d s   T grad   n d s
S Sf ST

 Sf
T grad   n d s  ST
Tf grad   n d s
La troisième formulation variationnelle s'écrit donc encore
- Trouver le champ de températures T tel que:

 V
T  d v  Sf
T grad   n d s  ST
 grad T  n d s 

A(, T)
ST
T f grad   n d s  
Sf
 f d s   V
 f dv  0 

B(
- sans conditions supplémentaires.
Bien que pouvant sembler plaisante (il n'y a plus de conditions à respecter sur les frontières), cette
dernière formulation présente des inconvénients.

3.5 Exemple 2 : Statique du solide déformable en petites déformations

Le problème est :
- Trouver le champ de déplacements  en tout point d’un domaine  tel que :

div  + f = 0

D(
- avec les conditions aux limites C() :
- des déplacements imposés (éventuellement nuls)  0 N  sur la partie de frontière Sd :
(N) -  (N) = 0  N  Sd

- des contraintes imposées (éventuellement nulles) q0(N) sur la partie de frontière Sf :

(N)  n(N) – q0(N) = 0  N  Sf

20
 où et  sont respectivement les tenseurs de contraintes et de déformations en tout point
du domaine  vérifiant :  = 2  +  Tr G,   grad   grad T   si le solide déformable est
1
2
élastique isotope linéaire en petites déformations, G est le tenseur identité, Tr est l’opérateur
"Trace" et  sont les coefficients de Lamé.
Cette équation (div  + f = 0) est une équation vectorielle, elle traduit l’équilibre local en tout point
de .

3.5.1 Première forme variationnelle

En utilisant (3.1), on obtient une première formulation équivalente :
- Trouver le champ de déplacements  tel que :

   div   f  d v  0 (3.21)
V

 (M) champ vectoriel dans V.

où  = 2  +  Tr G et 
1
2

grad   grad T  . 
Cette équation est une équation scalaire.
- avec les conditions aux limites :

(N) -  (N) = 0  N  Sd
(N)  n(N) – q0(N) = 0  N  Sf

3.5.2 Seconde forme variationnelle

En utilisant la formule (3.10), on obtient le problème équivalent suivant :
- Trouver le champ de déplacements  tel que :

  σ  grad ψ d v   ψ  σ n ds   ψ  f dv  0  ψ M 
V S V

- avec les conditions aux limites :

(N) -  (N) = 0  N  Sd
(N)  n(N) – q0(N) = 0  N  Sf
L'intégrale de bord peut être décomposée en deux intégrales :

 S
ψ σ n ds   Sd
ψ σ n ds   Sf
ψ σ n ds

On peut alors éliminer la condition aux limites de contrainte imposée en la reportant dans l'intégrale
sur Sf :

 Sf
ψ σ n ds   Sf
ψ  q0 d s

Le problème devient :
- Trouver le champ de déplacements  tel que

21
  σ  grad ψ d v   ψ  f d v   ψ  q0 d s   ψ σ n ds 0  ψ M 
V V Sf Sd

- avec la condition aux limites

(N) -  (N) = 0  N  Sd
Le tenseur des contraintes étant symétrique, on peut écrire

σ  grad ψ  σ  Sym  grad ψ 

Le problème devient :
Trouver le champ de déplacements  tel que

  σ  Sym  grad ψ  d v   ψ  f d v   ψ  q0 d s   ψ σ n ds 0  ψ M 
V V Sf Sd

(3.22)
- avec la condition aux limites

(N) -  (N) = 0  N  Sd
3.5.3 Interprétation mécanique :
Théorème des puissances virtuelles
Appelons «champ de vitesses virtuelles» le champ vectoriel arbitraire .
Remarque : On peut aussi bien l'appeler déplacement virtuel. Dans ce cas, les puissances virtuelles
s'appellent travaux virtuels.
Si  est une vitesse virtuelle, Sym(grad) s'interprète comme un tenseur des vitesses de
déformations virtuel.
L'intégrale   σ  Sym grad ψ  d v s'interprète alors comme une puissance virtuelle des efforts
V

intérieurs dans le champ de vitesses virtuelles .

L'intégrale  ψ  q0 d s s'interprète comme une puissance virtuelle des efforts de surface dans le
Sf

champ de vitesses virtuelles .

L'intégrale  ψ  σ  n d s s'interprète comme une puissance virtuelle des efforts de liaison dans
Sd

le champ de vitesses virtuelles .

L'intégrale  V
ψ  f d v s'interprète comme une puissance virtuelle des efforts de volume dans le

champ de vitesses virtuelles .

L'équation (3.22) traduit le théorème des puissances virtuelles : « la somme des puissances
virtuelles des efforts intérieurs et des efforts extérieurs est nulle dans tout champ de vitesses
virtuelles appliqué à la structure »,
(Si  est appelé «déplacement virtuel», le théorème s'appelle théorème des travaux virtuels). On
pouvait donc déduire cette formulation par des raisonnements mécaniques.
 On réécrit le problème (3.22) en isolant les termes contenant le champ inconnu  :
- Trouver le champ de déplacements  tel que :

  σ  grad ψ d v   ψ  σ  n ds   ψ  f dv  ψ  q0 d s  ψ M 
V Sd V Sf

A(, ) B(

22
où  = 2  +  Tr G et 
1
2
 
grad   grad T  .
avec la condition aux limites :

(N) -  (N) = 0  N  Sd
Le problème du solide élastique se met bien sous la forme donnée en (3.2).

Principe de minimisation de l’énergie potentielle

Si on prend le champ comme étant une variation du champ de déplacement réel det on remplace
le tenseur de contrainte par son expression en fonction de la loi de comportement et du tenseur de
déformation, on obtient :
 C    d v     σ  n d s     f d v     q 0 d s  U    0
V Sd V Sf

En intégrant par rapport à la variable on trouve:

1
U     C    d v     σ  n d s     f d v     q 0 d s
2 V Sd V Sf

Qui représente l’énergie potentielle du système. Le champ de déplacement solution du problème

représente le minimum de l’énergie potentielle de la structure. U    0
Exemple : Cas d’un ressort
Ep
k F

L’équation d’équilibre : kx  F  0 x
1
xkx  F   0 intégration => E p  x   kx 2  Fx
2 x’
Solution qui minimise l’énergie potentielle

E p  x   kx  F  0
x
3.6 Exemple 3 : Barre en traction-compression
Une barre AB de longueur L est soumise à un chargement linéique de traction-compression p(x)
suivant toute sa longueur et une force concentrée F à son extrémité (B). L’autre extrémité (A) est
supposée fixe (figure 1.3). La barre a pour section S. Le module de Young du matériau constituant
la barre est E.

p(x) F

B
A

Figure 4.2 : Barre en traction-compression

Le champ inconnu est un champ scalaire qui est le déplacement u(x) en chaque point suivant l’axe
de la barre.
dN ( x)
L’équation d’équilibre local est :  p ( x)  0  x    0, L 
dx
Où N(x) est l’effort normal dans la section située à x.
23
 u (0)  0 conditions aux lim ites type déplacement

Les conditions aux limites sont :  du
 N ( L)  ES dx ( L)  F conditions aux lim ites type effort
du
La loi de comportement est : N ( x)  ES
dx

1/ Ecrire la première forme variationnelle.

L
 dN ( x) 
0   dx  p( x) dl  0   x    0, L
u (0)  0 conditions aux lim ites type déplacement

Avec  du
 N ( L)  ES dx ( L)  F conditions aux lim ites type effort
2/ Ecrire la seconde forme variationnelle.
d
L

0  N ( x) dx  p( x)dl  N ( x)0  0   x    0, L

L

d
L

  N ( x) dx  p( x)dl   ( L) N ( L)  (0) N (0)  0  x    0, L 

0

d
L

  N ( x) dx  p( x)dl   ( L) F  (0) N (0)  0   x    0, L 

0
La seconde forme variationnelle est :
d
L
du du
0  ES dx ( x) dx  p( x)dl   ( L) F   (0) ES dx (0)  0   x    0, L

Avec u (0)  0

d
L L
du du
 ES ( x) dl  ES (0) (0)   p( x)dx   ( L) F  x    0, L 
0
dx dx dx 0

A(, ) B(
Avec u (0)  0
3/ Ecrire l’expression de l’énergie potentielle de cette structure

24
CHAPITRE 4 : LES ÉLÉMENTS ET LEUR ESPACE DE
FONCTIONS D’INTERPOLATION
~
Les solutions approchées trouvées par la méthode des éléments finis sont une juxtaposition F
~
de champs locaux fi définis dans chaque maille (figure 4.1 : la composante U1 du champ de
déplacement est définie par morceau sur chaque sous domaine).

Figure 4.1 : composante suivant la direction 1 du champ de déplacement.

L'objet de ce chapitre est de montrer comment on définit une interpolation sur une maille de
référence et comment passer ensuite à l’interpolation sur la maille réelle
Pour qu'une maille devienne un élément, il faut :
 Choisir des noeuds dans la maille. En effet, la résolution d'un problème par la méthode des
éléments finis se ramène à calculer les valeurs de la solution approchée aux noeuds du
maillage.
Remarque : Le nombre d'inconnues scalaires par nœud varie selon la nature tensorielle du
champ inconnu et selon la dimension de l'espace physique. Ce nombre est appelé nombre de
degrés de liberté (abréviation courante : ddl). Par exemple, si le champ cherché est un champ
vectoriel dans un espace physique à 3 dimensions, on a 3 ddl par noeud (1 ddl par composante
du champ).
 Choisir une famille de champs locaux, destinés à donner une valeur approchée de la solution
en tout point de la maille en ne connaissant que les valeurs d'une solution approchée aux
nœuds de l'élément. Cette famille de fonctions s'appelle l'espace des fonctions
d'interpolation de la maille.
Remarque : S'il y a plusieurs inconnues par noeud (c'est à dire plusieurs degrés de liberté), on
construit des interpolations pour chacun des degrés de liberté. En général, les interpolations de
chaque degré de liberté appartiennent au même espace de fonctions d'interpolation.
Dans ce chapitre, on va montrer comment on construit les interpolations. On le fera en supposant
qu'il n'y a qu'un seul ddl par nœud. S'il y en a plusieurs on procède de même pour chaque ddl.

4.1 Définition des mailles de référence

Dans un maillage, toutes les mailles ont des formes et des dimensions différentes. On trouve des
mailles linéiques, des mailles surfaciques et des mailles volumiques de toutes formes et de toutes
tailles (figure 4.2). Dans le but d'uniformiser et d'automatiser les calculs, on introduit la notion de
maille de référence. On ne donne ici que les mailles les plus classiques. Par convention
généralement admise (voir figure 4.3) :

25
 Figure 4.2 : maillages : volumique, surfacique et linéique
La maille de référence linéique est le segment x 1   1, 1  .
La maille de référence surfacique triangulaire est le triangle :
x1  0 ; x 2  0 ; x1  x 2  1
La maille de référence surfacique quadrangulaire est le carré :
x 1   1, 1  ; x 2   1, 1  .
La maille de référence volumique tétraédrique est le tétraèdre :
x1  0 ; x 2  0 ; x 3  0 ; x1  x 2  x 3  1 .
La maille de référence volumique pentaédrique est le pentaèdre
(prisme à base triangulaire) : x 1  0 ; x 2  0 ; x 1  x 2  1 ; x 3   1, 1  .
La maille de référence volumique hexaédrique est le cube :
x 1   1, 1  ; x 2   1, 1  ; x 3   1, 1  .
où xl, x2 et x3 sont les coordonnées d'un point courant m de la maille de référence.

Figure 4.3 : Mailles de référence

4.2 Interpolation polynomiale sur les mailles de référence

Les fonctions d'interpolations utilisées dans les logiciels sont pratiquement toujours des
polynômes. Si la fonction à interpoler est vectorielle ou tensorielle, on interpole de la même
manière chacune des composantes. On est donc ramené à des interpolations de champs scalaires.
~
Dans la suite, f r représente donc soit l'interpolation d'un champ scalaire soit l'interpolation d'une
composante de champ vectoriel ou tensoriel.
~
- Pour les mailles linéiques f r est un polynôme de x1,
~
- Pour les mailles surfaciques f r est un polynôme de x1 et x2,

26
 ~
- Pour les mailles volumiques f r est un polynôme de x1, x2 et x3.

4.2.1 Rappels sur les bases de polynômes

L'espace des polynômes de degré d est un espace vectoriel dont la dimension dépend du degré des
polynômes et du nombre de variables. Cet espace possède une base canonique constituée de tous les
monômes de degré non négatif inférieur ou égal à d. Par exemple :
 La base canonique des polynômes à une variable x1 de degré 3 est constituée des

monômes : 1, x 1 , x 12 , x 13 ,
 La base canonique des polynômes à deux variables x1 et x2 de degré 2 est constituée des

monômes : 1, x 1 , x 2 , x 12 , x 22 , x 1 x 2 , 
 La base canonique des polynômes à trois variables x1, x2 et x3 de degré 2 est constituée des
monômes :  1, x 1 , x 2 , x 3 , x 12 , x 22 , x 32 , x 1 x 2 x 2 x 3 x 3 x 1 .
Le tableau qui suit donne les dimensions des espaces de polynôme pour des degrés de 1 à 5 pour les
polynômes à 1, 2 ou 3 variables.

degré 1 variable 2 variables 3 variables

1 2 3 4
2 3 6 10
3 4 10 20
4 5 15 35
5 6 21 56
A partir de la base canonique, on peut engendrer une infinité de bases : Si n est la dimension de
l'espace de polynômes, toute matrice régulière n x n définit une autre base.
Remarque : Certaines bases ou sous bases avec des propriétés ont reçu des noms : Gegenbauer,
Hermite, Jacobi, Laguerre, Tchebitcheff, …
Par exemple les 6 polynômes Pl, P2, P3, P4, P5, P6 sont une base de l'espace des polynômes de degré
2 à 2 variables si la matrice des pij est régulière. :
 P1   p11 p12 p13 p14 p15 p16   1 
P  p p 22 p 23 p 24 p 25 p 26   x 
 2  21  1 
 P3  p 31 p 32 p 33 p 34 p 35 p 36   x2 
      2 
P4  p 41 p 42 p 43 p 44 p 45 p 46   x1 
 P5  p 51 p 52 p 53 p 54 p 55 p 56   x 22 
     
P6  p 61 p 62 p 63 p 64 p 65 p 66   x 1 x 2 
Si on extrait d'une base une sous base de m polynômes de base, on engendre un sous espace
vectoriel de dimension m.

4.2.2 Interpolation polynomiale sur une maille de référence

Soit Er l'espace de référence, et soit m un point courant de Er. Suivant le type de maille (linéique,
surfacique ou volumique), le point m a 1, 2 ou 3 coordonnées. On note dr la dimension de l'espace
de référence (qui est aussi la dimension de la maille réelle).
Pour définir une interpolation, on choisit dans l'espace de référence, n points (n>dr) appelés nœuds
de référence qu'on note m(j), j = 1, ... , n.
~
Une interpolation polynomiale f r sur l'élément de référence Er est définie par :

27
 n
~
f r m    u i  P i  m  (4.1)
i 1
telle que :
   
n
~  j
fr m   u i  P i  m  j  u  j  j 1, ..., n  (4.2)
i 1
où les P(i)(m) sont des polynôme de dr variables, et où les u(j) sont n valeurs données aux nœuds m(j).
~ ~
Autrement dit, une interpolation f r attribue une valeur f r (m) à tout point m, et sa valeur aux n
noeuds m(j) est u(j).
.
Dans une interpolation polynomiale à n noeuds, les P(i) doivent être une n-base de polynôme.
Soit B k m  une n-base de polynômes connue. On peut donc écrire :
n
P i  m    a   B m 
k
i
k (4.3)
k 1

où les a ki  sont le terme général d'une matrice régulière.

La condition (5.2) implique les n2 égalités :
 
P i  m  j   i j (4.4)
soit encore

 a   B m     
n
i j
k k ij (4.5)
k 1

PC  G  (4.5)
où encore sous forme matricielle
 a 11  a k1  a n1   B1 m   
1
 B1 m  j    B1 m n    1  0  0
    
           
 a 1i   a ki   a ni   B k
  
 
m 1  Bk m  j    Bk m n   
  0 
 
1  0

          
a  n   a  n   a  n    B
 1 k n   n  
m 1  Bn m  j    Bn m n   
 0 
 0  1
Ce qui montre que la matrice [P] des coefficients des polynômes P(i) sur la base {Bk} est l'inverse de
la matrice [C] des valeurs des polynômes de la base Bk aux noeuds. La matrice [C] doit donc être
régulière (c.à.d. géométriquement, les k points doivent être distincts, non alignés dans une maille
surfacique et non coplanaires dans une maille volumique) et les coefficients des polynômes de base
de l'interpolation sont donnés par :
P  C 1
Ainsi, les coefficients des polynômes de base de l'interpolation dans l'espace de référence sont
déterminés par les coordonnées des noeuds dans la maille de référence et par le choix de la base
{Bk(m)}. Le calcul de [P] est très simple : il suffit de construire [C] de terme général Ckj = Bk(m(j))
et de l'inverser.
~
En résumé, une interpolation f r sur une maille de référence est caractérisée par :
 la dimension de l’espace de référence (linéique, surfacique, volumique), c'est-à-dire le
nombre de coordonnées dans l'espace de référence,
 le choix de n noeuds par leurs coordonnées dans l'espace de référence,
 le choix de la base {Bk(m)}.
On en déduit la matrice [P] des coefficients des polynômes de base de l'interpolation.
Les polynômes P(i) (polynôme de base de l'interpolation) sont appelés fonctions de forme.

28
Remarque
Dans les logiciels de calcul par éléments finis, on distingue généralement, pour une maille de
référence donnée, deux types d’interpolation. Une interpolation linéaire (ou bilinéaire) et une
interpolation quadratique qui fait intervenir des polynômes de degré 2 en chaque variable.

5.2.3 Exemple 1 (standard)

On prend une maille surfacique (les variables sont donc x1 et x2), et on veut construire une
interpolation à 4 noeuds (voir figure 5.2) dont les coordonnées sont :
      
m 1  x 11 , x 21 ; m 2   x 12  , x 22  ; m 3   x 13 , x 23 ; m 4   x 14  , x 24  
L'interpolation est donc de la forme :
4
~
f r x 1 , x 2    u   P   x , x 
i i
1 2
i 1

Figure 5.2 : Interpolation à 4 noeuds dans un élément de référence surfacique

Les fonctions de forme P(i) sont une base de polynômes à deux variables de dimension 4. Le tableau
page 29 montre qu'on peut travailler dans un sous espace des polynômes de degré 2 (en effet si on
prend des polynômes de degré 1 alors on engendre un espace de dimension 3 <4 ce n’est pas
suffisant. Si on choisit des polynômes de degré 2 alors on obtient un espace de dimension 6>4 on ne
prendra pas les 6 monômes mais on va en extraire 4) . On peut choisir comme sous base canonique
{Bk} les quatre polynômes (B1 = 1, B2 = x1, B3 = x2, B4 = xlx2) de la base canonique.
Remarque : Il s’agit bien d’un choix : quatre polynômes de degré 2 quelconques mais
indépendants peuvent convenir. On s’arrange généralement pour garder des monômes de la base
canonique tels que les rôles des variables aient une certaine «symétrie». L'interpolation obtenue
présente alors une certaine «isotropie». On pouvait prendre aussi une sous base de l'espace des
polynômes de degré 3).
Les conditions (4.2) s'écrivent :

P 1 x   , x     1
1
1
2
1

P 1 x 12  , x 22  0 
P 1 x 13 , x 23  0 
P 1 x 14  , x 24  0
P 2  x   , x     0
1
1
2
1
 1
P 2  x 12  , x 22  
P 2  x 13 , x 23  0   0
P 2  x 14  , x 24 
P 3  x   , x     0
1
1 1
2   0
P 3  x 12  , x 22  
P 3 x 13 , x 23  1   0
P 3  x 14  , x 24 
P 4  x   , x     0
1
1
2
1
  0
P 4  x 12  , x 22  
P 4  x 13 , x 23  0  1
P 4  x 14  , x 24 

En posant : P i   a 1i   a 2i  x 1  a 3i  x 2  a 4i  x 1 x 2 , on construit une 4-base de polynômes P(i) si la
matrice des a ji  est régulière. Les conditions (4.2) s'écrivent alors :

29
  a 11 a 21 a 31 a 41   1 1 1 1  1 0 0 0
 2    x 1
a 1 a 22  a 32  a 42    1 x 12  x 13 x 1 4   0 1 0 0

a 13 a 23  a 33 a 43   x 21 x 22  x 23 x 24   0 0 1 0
 4    1 1   
a 1 a 24  a 34  a 44   x 1 x 2 x 1 2  x 22  x 13  x 23 x 14  x 24   0 0 0 1
On déduit les coefficients des fonctions de forme par inversion :
1
 a 11 a 21 a 31 a 41   1 1 1 1 
 2   
a 1 a 22  a 32  a 42    x 1
1
x 1 2  x 13 x 14  

a 13 a 23  a 33 a 43   x 21 x 22  x 23 x 24  
 4    1 1 2  2  
a 1 a 24  a 34  a 44    x 1 x 2 x 1 x 2 x 13 x 23  x 14  x 24  
Par exemple, les noeuds sont les points (-1, -1), (1, -1), (1, 1), (-1, 1) dans l'espace de référence :
1
 a 11 a 21 a 31 a 41   1 1 1 1 1  1 1 1 
 2   
a 1 a22 
a3 2 
a 42   1 1 1  1 
1 1 1 1  1

 
a 13 a2 3 
a 33 a 43  1  1 1 1 4 1 1 1 1 
 4      
a 1 a 24  a 34  a 44    1  1 1  1 1  1 1  1

Les quatre fonctions de forme de cet élément sont donc

1
1
1  x1  x2  x1 x2   1 1  x1 1  x2 
P 1  x1 , x2  
P  1  1 1 1   1  4 4
 2      x  1 1
P   1 1 1 1 1 P  x1 , x2   1  x1  x2  x1 x2   1  x1 1  x2 
2 
 1  4 4
 P 3   4 1 1 1 1   x2  1 1
 4       P  x1 , x2   1  x1  x2  x1 x2   1  x1 1  x2 
3 
P  1  1 1 1 x1 x 2  4 4
1 1
P  x1 , x2   1  x1  x2  x1 x2   1  x1 1  x2 
4 
4 4

et l'interpolation en fonction des valeurs aux noeuds est :

~
f r x 1 , x 2   u 1 P 1 x 1 , x 2   u 2  P 2  x 1 , x 2   u 3  P 3 x 1 , x 2   u 4  P 4  x 1 , x 2 

4.2.4 Exemple 2 (non standard)

On construit une interpolation à 7 noeuds dans une maille surfacique (voir figure 5.3).

Figure 4.3 : Interpolation à 7 noeuds dans un élément de référence surfacique

L'espace d'interpolation est de dimension 7. Le tableau page 23 montre qu'il faut prendre au
minimum des polynômes de degré 3. On peut prendre, par exemple, comme sous base {Bk} de la
base canonique la suite suivante de monômes :

30
 B 1  1, B 2  x 1 , B 3  x 2 , B 4  x 12 , B5  x 22 , B 6  x 12 x 2 , B 7  x 1 x 22 
La matrice [C] de terme général Ckj = Bk(m(j)) est :

1 1 1 1 1 1 1   5 3 45 27 45 81 27 
 4   
1 2 2 1 1  8 8 8 8 8 8 
3 0 0  1 27 
3 3 3 3   
3 9 27

9

81

 1 2 2 1 1   4 8 8 8 8 8 8 
0 0  dont l'inverse est
 3 3 3 3 3   1 
9

9 9 9 81

27 
1 4 4 1 1   4 8 8 8 8 8 8 
0 0  81 
9

9 9 9 9  P    1 
9

9 9 9

27

1 1 1 1 1   4 8 8 8 8 8 8 
0 0 
 9 9 9 9 9   1 9

3

9 27 27 81
 
0 4 2 1   4 8 8 8 8 8 8 
0 0 0  5 45 3 45 27 27 81 
 27 27 27    
  
2 4 1   4 8 8 8 8 8 8 
0 0 0 0   3 27 27 27 27 27 27 
 27 27 27     
 2 4 4 4 4 4 4 
Ce qui donne les coefficients des fonctions de forme.

4.2.5 Matrice d’interpolation

L’expression (5.1) peut s’écrire sous la forme du produit d’un vecteur ligne contenant les valeurs
aux nœuds par un vecteur colonne contenant les fonctions d’interpolation.
 P 1 m 
 2  
~  n   P m 
n
f r m    u P m   u u ......u    u n N 
i  i  1 2 
i 1   
 P n  m 
u n est le vecteur des degrés de liberté de l’élément.
N  est la matrice d’interpolation.
4.2.6 Exercices
1/ Trouver les fonctions de forme de l’élément linéique à deux nœuds situés à ses deux extrémités.
2/ Trouver les fonctions de forme de l’élément linéique à trois nœuds situés à ses deux extrémités et
à son milieu.
3/ Trouver les fonctions de forme de l’élément surfacique triangulaire à trois nœuds situés à ses
sommets.
4/ Trouver les fonctions de forme de l’élément volumique tétraédrique à quatre nœuds situés à ses
sommets.
5/ Trouver les fonctions de forme de l’élément volumique hexaédrique à huit nœuds situés à ses
sommets. (-1,-1,-1) (1,-1,-1) (1,1,-1) (-1,1,-1) (-1,-1,1) (1,-1,1) (1,1,1) (-1,1,1)

4.3 Continuité C0 inter-éléments : Les éléments conformes

On souhaite généralement que la solution approchée trouvée présente une continuité C0 inter-
élément. Les éléments qui ont cette propriété sont appelés éléments conformes.
Pour les mailles linéiques, la continuité C0 est facile à assurer : Il suffit de mettre des noeuds aux
extrémités de la maille. On est ainsi assuré que les fonctions d'interpolation prennent la même
valeur au nœud commun.
Dans la suite on ne parlera que de la conformité des éléments surfaciques et volumiques.

31
Si on choisit les fonctions d'interpolation telles que : la valeur de l'interpolation sur une frontière
ne dépend que des noeuds de la frontière, on garantit que la solution approchée est continue C0 à
la traversée de la frontière.1

discontinuité
continuité

El1 El2

Figure 4.4 : Conformité partielle d’un élément

Autrement dit, sur une arête d'élément conforme surfacique, l'interpolation doit se réduire à une
interpolation de maille linéique sur les noeuds de l'arête ; et sur une face d'élément conforme
volumique, l'interpolation doit se réduire à une interpolation de maille surfacique sur les noeuds de
la face.

P3 P1
P2
1 1
1

m(1) m(1) m(1)

1 1
1
1 1 m(3)
1 m(3)
m (3) x2 m (2) x2
m(2) x2 m(2) x1
x1
x1
Figure 4.5 : Conformité de l’interpolation : élément triangulaire à trois noeuds
Pour construire des éléments conformes, il faut donc adjoindre les conditions ci-dessus aux
conditions d'interpolation 4.3.
Remarque :
Les principaux éléments de référence conformes qu'on trouve dans les logiciels sont :
- les éléments linéiques dont deux de leurs noeuds sont les extrémités,
- les triangles à 3 noeuds (sommets),
- les triangles à 6 noeuds (sommets + milieux des arêtes),
- les quadrangles à 4 nœuds (sommets),
- les quadrangles à 8 noeuds (sommets + milieux des arêtes),
- les tétraèdres à 4 noeuds (sommets),
- les tétraèdres à 10 noeuds (sommets + milieux des arêtes),
- les tétraèdres à 20 noeuds (sommets + noeuds au tiers des arêtes),

1
La juxtaposition des interpolations des éléments conformes est alors une fonction C0 définie sur , et qui appartient
donc à l'espace de Sobolev H1(). C'est pour des solutions approchées appartenant à cet espace que les théorèmes de
convergence de la méthode sont établis.

32
 - les pentaèdres à 6 noeuds (sommets),
- les pentaèdres à 15 noeuds (sommets + milieux des arêtes),
- les pentaèdres à 24 noeuds (sommets + noeuds au tiers des arêtes),
- les hexaèdres à 8 noeuds (sommets),
- les hexaèdres à 20 noeuds (sommets + milieux des arêtes),
- les hexaèdres à 27 noeuds (sommets + milieux des arêtes + milieux des faces + noeud
central),
- des hexaèdres à 32 noeuds (sommets + noeuds au tiers des arêtes).

4.4 Interpolation dans l'élément réel

Les coordonnées des nœuds de l'élément réel sont différentes de celles de l'élément de référence.

Il convient donc de construire des fonctions d'interpolation sur l'élément réel en utilisant les
~
interpolations f r construites sur les éléments de référence.

4.4.1 Principe de construction

Soient m un point de l'élément de référence, et M un point de l'élément réel.
~ ~
Soient f r m  la fonction d'interpolation connue sur l'élément de référence, et f M  la fonction
d'interpolation cherchée sur l'élément réel.
Considérons une transformation  inversible qui à tout point m de l'élément de référence associe un
point M de l'élément réel :

  : m  élément de référence  M = (m)  élément réel

telle que les noeuds se correspondent :

M(i) = (m(i)) otmail.fr

X3
m M

Elément de référence x1

Elément réel
X1 X2

Figure 5.6 : Transformation de l’élément de référence vers l’élément réel

~
Considérons l'application f définie par :
~ ~ ~
 f : M  élément réel  f (M) = f r m   M = (m)
On a donc :
33
 ~ ~ ~
f M   f m   f r m   m  élément de référence
soit encore :
~ ~ ~ ~
f    f r  f  f r   1
~
L’application f est une interpolation sur l’élément réel car :
n
~ ~
f M   fr m    u i P i  m 
i 1

     
n
~ ~ ~
f M   u i
P i   1 M  et f M i   f r m i   u i 
i 1

~ ~
On doit noter que si f r est une interpolation polynomiale, f n'en est pas une en général, car
 
P i   1 M  n'est pas un polynôme en M en général.

4.4.2 Choix de la transformation  : l'approximation géométrique

On cherche idéalement une transformation  qui fait correspondre les points m de la maille de
référence aux points M de la maille réelle, telle que :
1. la transformation  est inversible
2. la transformation  fait correspondre les frontières
3. la transformation  fait correspondre les noeuds : M(i) = (m(i))
Malgré le nombre de conditions imposées à , il y a une infinité de solutions à ce problème, mais
ces transformations sont compliquées à trouver, surtout lorsque la géométrie de l'élément réel est
compliquée. On va donc chercher des transformations  plus simples.
~
La condition impérative pour que f soit une interpolation est la correspondance des noeuds. On va
choisir des transformations  qui ne respectent que cette condition.
La conséquence, est que si on appelle M' le transformé de m, en général M' n'appartient pas à
l'élément réel, sauf si m est un nœud (figure 5.7).


M
m Approximation de
X3
l’élément réel
M’
Elément de référence x1

Elément réel
X1 X2
Cela signifie que :
- Pour les éléments linéiques, la courbe des points M' =  (m) est différente de la courbe réelle de
l'élément.
- Pour les éléments surfaciques, la surface des points M' =  (m) est différente de la surface réelle
de l'élément, ainsi que la courbe des arêtes.
- Pour les éléments volumiques, les arêtes et les faces transformées sont différentes des arêtes et
des faces réelles.

34
La transformation  transforme l'élément de référence en une approximation géométrique de
l'élément réel. Les noeuds de l'approximation géométrique et de l'élément réel sont confondus, mais
les autres points sont distincts.
~ ~
L'interpolation f = f r o -1 est une interpolation sur une approximation géométrique de la maille
réelle. Autrement dit, dans un calcul par éléments finis, les mailles réelles sont remplacées par leur
approximation géométrique, et cette approximation géométrique est déterminée par le choix de .
On cherche donc une fonction  inversible qui fasse correspondre les noeuds, c'est à dire telle que:
M(i) = (m(i))  le noeud i
Ce problème est un problème d'interpolation sur l'élément de référence : Choisir , c'est choisir
l'interpolation des points M' entre les points M(i). Le problème est donc :
~ ~ ~
pour les mailles linéiques , trouver trois fonctions d'interpolation X1 , X 2 et X 3 telles que :
X 1i   X1 x 1i  
~ ~
X 2i   X 2 x 1i  
~ ~ ~ ~
X 3i   X 3 x 1i   
~ ~ ~
pour les mailles surfaciques , trouver trois fonctions d'interpolation X1 , X 2 et X 3 telles que :
X 1i   X1 x 1i  , x 2i  
~ ~
X 2i   X 2 x 1i  , x 2i  
~ ~ ~ ~
X 3i   X 3 x 1i  , x 2i   
~ ~ ~
pour les mailles volumiques , trouver trois fonctions d'interpolation X1 , X 2 et X 3 telles que :
~ ~

X 1i   X1 x 1i  , x 2i  , x 3i   ~ ~

X 2i   X 2 x 1i  , x 2i  , x 3i   ~ ~

X 3i   X 3 x 1i  , x 2i  , x 3i  
~
où X ji  est la jeme coordonnée dans l'espace physique du nœud réel i , et où x ji  est la jeme
coordonnée dans l'espace de référence du noeud de référence i.
Si on prend des interpolations polynomiales, les fonctions d'interpolation des coordonnées sont de
la forme (On les donne ici pour les éléments volumiques, et sont aisément simplifiées pour les
éléments surfaciques ou linéiques) :
n
~
X 1  X1 x 1 , x 2 , x 3    X 1i  Q i  x 1 , x 2 , x 3  (4.6)
i 1
n
~
X 2  X 2 x 1 , x 2 , x 3    X   Q   x , x
2
i i
1 2 , x3  (4.7)
i 1
n
~
X 3  X 3 x 1 , x 2 , x 3    X   Q   x , x
3
i i
1 2 , x3  (4.8)
i 1

où les Q(i) sont les polynôme de base de l'interpolation géométrique, i  [1, ..., n], et où n est la
dimension de l'espace des fonctions d'interpolation (c'est à dire le nombre de noeuds). Cette
interpolation doit être conforme, pour que les frontières de deux éléments voisins soient les mêmes.
Pour trouver ces interpolations, on emploie les mêmes techniques que précédemment, puisque ce ne
sont que des interpolations sur l'élément de référence.
Quelques définitions :
- Si la méthode d'interpolation des coordonnées est la même que celle des valeurs aux noeuds on dit
que l'élément est isoparamétrique. On a alors [P] = [Q]
- Si l'interpolation des coordonnées est de degré supérieur au degré d'interpolation des valeurs aux
noeuds on dit que l'élément est superparamétrique.
- Si l'interpolation des coordonnées est de degré inférieur au degré d'interpolation des valeurs aux
noeuds on dit que l'élément est subparamétrique.

35
Il reste un dernier problème : la transformation M = (m) doit être inversible. Il faut donc s'assurer
que le déterminant de la matrice jacobienne de la transformation soit non nul en tout point de la
maille de référence :
~ ~ ~
 d X1 d X1 d X1 
 
d~ x1 d x2
~
d x3 
~
d X2 d X2 d X2 
det J   det  0
d x1 d x2 d x3 
 ~ ~ ~ 
d X3 d X3 d X3 
dx d x2 d x3 
 1
Remarque :
Les mécaniciens feront aisément le parallèle avec la transformation d'un milieu continu, d'un
domaine r de référence à un domaine t actuel. Dans ce cas, la matrice jacobienne n'est autre que
la matrice des composantes du gradient de la transformation F dans un système de coordonnées
cartésiennes, et on doit avoir det F  0.
Les équations (4.6 à 4.8) montrent que son calcul se ramène à des dérivées des polynômes de base
Q(i).
Le déterminant est donc une fonction des coordonnées réelles des noeuds X ji  et des variables xi. Il
peut arriver qu'il s'annule pour un certain ensemble de coordonnées nodales ou en certains points
non nodaux (par exemple, si deux noeuds d'un élément réel sont confondus ou trop proches, ou
encore si les trois noeuds d'un triangle réel sont alignés, etc. Il se peut aussi que le jacobien soit nul
en certains points (xr,yr,zr) non nodaux, par exemple si l'approximation géométrique de la maille
linéique, surfacique ou volumique a des points doubles). Tout bon logiciel se doit de faire cette
vérification avant de lancer un calcul.
Lorsque le jacobien est non nul partout, on peut alors trouver -1 en inversant le système (4.4).
On trouve trois fonctions ~x1 , ~
x 2 et ~
x3
x 1 X 1 , X 2 , X 3 
x1  ~ x 2 X 1 , X 2 , X 3  x 3  ~
x2 ~ x 3 X1 , X 2 , X 3  (4.9)
Le calcul de la transformation inverse -1 est difficile, sauf quand l'interpolation des coordonnées est
strictement linéaire (càd. uniquement dans les éléments linéiques à 2 noeuds, les éléments
triangulaires à 3 noeuds et les tétraèdres à 4 noeuds. Dans les autres, l'interpolation ne peut pas être
strictement linéaire. Ce sont les seuls cas ou l'interpolation sur l'élément réel reste polynomiale), car
dans ce dernier cas, l'inversion de  se ramène à des calculs algébriques matriciels. On verra plus
loin que ce calcul n'est pas nécessaire et qu'il suffit de s'assurer que la transformation est inversible.
Un grand nombre des éléments proposés dans les logiciels ont une interpolation linéaire des
coordonnées, mais on trouve aussi des interpolations de coordonnées de degré supérieur, qu'on
appelle souvent éléments courbes. Dans le cas d'interpolation linéaire des coordonnées:
- Les mailles linéiques et les arêtes sont approximées par des segments de droite joignant les
noeuds. S'il y a des noeuds intérieurs, c'est une ligne brisée.
- Les mailles triangulaires et les faces triangulaires à trois noeuds sont approximées par des
triangles plans
- Les mailles triangulaires et les faces triangulaires à plus de trois noeuds sont approximées par des
surfaces polynomiales compliquées à contour polygonal gauche.
- les mailles quadrangulaires et les faces quadrangulaires à quatre noeuds sont approximées par des
paraboloïdes hyperboliques.
- les mailles quadrangulaires et les faces quadrangulaires à plus de quatre noeuds sont approximées
par des surfaces polynomiales compliquées à contour polygonal gauche.

36
On peut remarquer que dans le cas des éléments linéiques à 2 noeuds, des éléments triangulaires à 3
noeuds et des éléments tétraédriques à 4 noeuds, les approximations géométriques des lignes et des
surfaces n'ont jamais de points doubles. Le jacobien ne peut s'annuler que pour des dégénérescences
géométriques (longueurs nulles, surfaces nulles ou volumes nuls). Pour les autres, il convient de
vérifier que le jacobien de  n'est jamais nul non seulement aux noeuds, mais aussi entre les noeuds,
c'est à dire que l'approximation géométrique ne contient pas de point double.
Finalement, l'interpolation sur l'approximation géométrique de l'élément réel est :
n
~
f X1 , X 2 , X 3    u   P   ~x X , X
i 1
i i
1 1 2 , X 3 , ~ x 3 X 1 , X 2 , X 3 
x 2 X1 , X 2 , X 3 , ~ (4.10)

Elle est déterminée quand on connaît les coordonnées des noeuds X ji  et les valeurs aux nœuds u(i).

Quand on a choisi l'interpolation géométrique , on peut calculer à l'avance la matrice [Q] des
~
coefficients des polynômes de base des interpolations des coordonnées X i pour chaque élément de
référence. Par contre, le jacobien dépend des coordonnées réelles, et sa non nullité ne peut être
vérifiée que quand on connaît les coordonnées réelles des noeuds.
On peut toutefois préparer à l'avance la matrice [Q] des coefficients des dérivées des polynômes de
base.

4.4.3 Résumé
Un élément est caractérisé par :
- la forme de sa maille de référence
- la position des noeuds dans la maille de référence (sommets obligatoires pour la conformité)
~
- l'interpolation de référence f r (c'est à dire la matrice [P] des coefficients des polynômes de base
de l'interpolation)
- l'interpolation  des coordonnées des points non nodaux M' (c'est à dire la matrice [Q] des
coefficients des polynôme de base de l'interpolation, et la matrice [Q’] des coefficients des
dérivées de polynôme de base pour le calcul de [J]. L'interpolation des coordonnées détermine
l'approximation géométrique de l'élément.
~
- On en déduit une interpolation (en général non polynomiale) f sur l'élément réel.

4.4.4 Exercices
1/ Trouver l’expression de la transformation  pour un élément linéique à deux nœuds placés à ses
extrémités.
2/ Soit un élément linéique à trois nœuds. Deux sont placés à ses extrémités et un nœud en son
centre.
a) Trouver les fonctions de forme de cet élément
b) Trouver l’expression de la transformation  dans le cas où l’élément est isoparamètrique.

4.5 Dérivées et intégrales dans l'élément réel

Les dérivées et les intégrales dans l'élément réel (les intégrales sont en fait sur l'approximation
géométrique de l'élément) des fonctions d'interpolation apparaissent dans la discrétisation de la
~
formulation variationnelle. On va montrer ici que les dérivées de l'interpolation f dans l'élément
~
réel s'expriment en fonction des dérivées de l'interpolation f r dans l'élément de référence, et qu'on
peut ramener les intégrales sur l'élément réel à des intégrales sur l'élément de référence.

37
4.5.1 Mailles volumiques
Gradient
~ ~ ~
Puisque f = f r o -1 le gradient de l'interpolation f est :
~ ~
grad f = grad f r  (grad -1
~ ~
où grad f et grad f r sont des tenseurs du premier ordre et grad  un tenseur du second ordre (il faut
noter que grad (-1) = (grad )-1 et on verra qu’il n'est donc pas nécessaire de calculer explicitement
les fonctions ~  
x i M j pour calculer la matrice jacobienne de la transformation inverse).
En termes de composantes dans un système de coordonnées cartésiennes :
~ ~
(grad f )i = (grad f r )k (grad -1)ki
~ ~
 
f

 Xi  x k
 fr
J1 k i
avec sommation sur l’indice k, et où on a posé :
~ ~ ~
 d X1 d X1 d X1 
 
d~ x1 d x2
~
d x3 
~
d X2 d X2 d X2 
J   
dx d x2 d x3 
 ~1 ~ ~ 
d X3 d X3 d X3 
dx d x2 d x3 
 1
ou sous forme matricielle
~ ~ ~ 1
 d X1 d X1 d X1 
 
 d x1 d x2 d x3 
~ ~ ~ ~ ~ ~
  f  f  f    f r  f r  f r   d X~ 2 d X~ 2 d X~ 2 
  
  X 1  X 2  X 3    x1  x2  x3   d ~ x1 d x2 d x3 
~ ~ 
d X3 d X3 d X3 
dx d x2 d x3 
 1
~
Le calcul des dérivées de f nécessite donc l'inversion de [J] pour chaque point de l’élément réel.
Ecriture matricielle
En utilisant la notion d’interpolation sur l’élément de référence, l’expression du gradient sur
l’élément de référence est :
 P 1 P 1 P 1 
 x 
 12  x22  x32  
~ ~ ~  P P P 
 f r f r f r  1 2   n   x x2 x3   u 1u  2 .......u  n  B 
   u u .......u  1 
 x1 x2 x3      
    
 P n  P n  P n  
 
 x1 x2 x3 
~ ~ ~
f f f 
  u u .......u B J 
1 2  n  1

  X1  X 2  X 3 
[B] est la matrice d’interpolation des dérivées.

38
Second gradient
~
Le second gradient de l'interpolation f est
~ ~
grad grad f = grad (grad f r  (grad -1)
T ~ ~
 grad    grad grad f r  grad f r  grad grad  
1
 
En termes de composantes dans un système de coordonnées cartésiennes
~
2 f
 
~
 2 fr
~ 1
 
 fr  J k i
 J 1 ki 
 Xi  X j  xk  x j  xk  Xj

avec sommation sur l'indice k.

~
Les dérivées secondes de f r peuvent être préparées d'avance en rangeant dans une matrice [P"] les
~
coefficients des dérivées secondes des polynôme de base de l’interpolation f r . La matrice [J] est
calculée avec la matrice [Q'], mais son inversion doit être faite pour chaque élément réel car elle
dépend des coordonnées réelles des noeuds.

Intégrales
Si on note V' l'approximation géométrique de l'élément réel et si on note Vr l'élément de référence,
l'élément de volume de l'approximation géométrique de l'élément réel est :

dv = det [J] dxldx2dx3

Une intégrale sur l'élément réel se ramène à une intégrale sur l'élément de référence par changement
de variables :
v
~
vr
 ~ ~

 GX1 , X 2 , X 3 d v   G X1 x 1 , x 2 , x 3 , X 2 x 1 , x 2 , x 3 , X 3 x 1 , x 2 , x 3  det Jdx 1dx 2 dx 3
où det [J] est une fonction (en général non polynomiale) de x1, x2, x3.

4.5.3 Mailles linéiques

L'espace de référence est de dimension 1. On note C' l'approximation géométrique de la courbe
~
réelle C engendrée par la transformation . L'interpolation f n'est définie que sur C' et les dérivées
~
de f par rapport aux 3 coordonnées X i' de M n'ont aucun sens.
La courbe C’est naturellement paramétrée par x1 :
~ ~ ~
X 1'  X 1' x 1  X '2  X '2 x 1  X 3'  X 3' x 1 
~
Les trois fonctions d'interpolation X i' sont la définition paramétrique de C'.
La base naturelle tangente de C’est le vecteur tangent non unitaire
~ ~ ~
d M' d X1 d X2 d X3
e1   E1  E2  E3
d x1 d x1 d x1 d x1
Gradient
Le gradient linéique est
~ ~
~ f   fr 
grad f  e1  e1
 x1  x1
La direction de dérivation unitaire u est

39
  e 1 
u  1  e1
e1 ~ 2 ~ 2 ~ 2
 dX 1   dX 2   dX 3 
       
 d x1   d x1   d x1 
La dérivée le long de l'approximation géométrique de la courbe (c'est à dire la dérivée par rapport à
l'abscisse curviligne) est donc :
~
d fr
~
d f ~ ~  d x1
 Du f  grad f  u 
dl ~ 2 ~ 2 ~ 2
 dX 1   dX 2   dX 3 
       
 d x1   d x1   d x1 
Ecriture matricielle
En utilisant la notion d’interpolation sur l’élément de référence, l’expression du gradient sur
l’élément de référence est :
 P 1 
 
 x12  
~  P 
f r  
 u 1u 2 .......u  n   x1   u 1u 2 .......u n  B 
x1   
  
 P n  
 
 x1 
[B] est la matrice d’interpolation des dérivées.

Second gradient

On déduit de ce qui précède

 ~ 
 d fr 
2 ~  
d f d  d x1 

d l2 dl   ~  2  ~ 2  ~  2 
  dX 1    dX 2    dX 3  
  d x1  dx 
 1
dx  
 1 

Remarque : Si l'interpolation des coordonnées est linéaire (l'approximation géométrique de la maille
~
d2Xi
est une ligne brisée) alors les dérivées secondes sont nulles.
d x 12
Intégrales

L'élément de longueur est

~ 2 ~ 2 ~ 2
 dX 1   dX 2   dX 3 
d l       
d x   d x    d x  dx 1
 1   1   1 
On a donc
~ 2 ~ 2 ~ 2
1  dX 1   dX 2   dX 3 
 Gld l  Glx 
l2
      
l1 1
1  d x    d x    d x  dx 1
 1   1   1 

Exercices :
40
1/ Soit un élément tétraédrique à 4 nœuds placés à ses sommets.
a) Trouver les termes de la matrice d’interpolation [N].
b) Trouver les termes de la matrice d’interpolation du gradient [B].
2/ Soit un élément linéique à 2 nœuds placés à ses sommets.
a) Trouver les termes de la matrice d’interpolation [N].
b) Trouver les termes de la matrice d’interpolation du gradient [B].
3/ Soit un élément triangulaire à 3 nœuds placés à ses sommets.
a) Trouver les termes de la matrice d’interpolation [N].
b) Trouver les termes de la matrice d’interpolation du gradient [B].
Exercice 4 : Dérivée dans l’élément réel
Soit un élément linéique a 2 nœuds (aux sommets). Les coordonnées des nœuds de l’élément réel
sont M i   X 1i , X 2i , X 3i  (i=1,2).
1) Trouver l’expression de la transformation  .
2) Trouver le vecteur de la base naturelle e1
3) Les nœuds de l’élément réel sont M 1 0,0,0  et M 2  1,1,1
Calculer la dérivée par rapport à l'abscisse curviligne en fonction des valeurs u i  (i=1, 2) du champ
inconnu aux nœuds M i  .

4.6 Evaluation numérique des intégrales

Le calcul des intégrales nécessaires pour la méthode des éléments finis peut s’effectuer d’une façon
analytique uniquement dans le cas où les termes à intégrer sont des polynômes. Ceci est vrai
uniquement dans le cas où les éléments réels ont une forme simple (pas d’éléments ayant des arêtes
courbes). Dans le cas où les termes à intégrer ne se présentent pas sous forme de polynômes
simples, on utilise des techniques d’intégration numérique.
Les intégrales à calculer par la méthode des éléments finis sont des intégrales ramenées sur le
domaine de référence. Ils s’écrivent sous la forme :
v vr
~ ~ ~

 G X 1, X 2 , X 3  d v   G X1 x1, x2 , x3 , X 2 x1, x2 , x3 , X 3 x1, x2 , x3  det J dx1dx2dx3
4.6.1 Intégration analytique
Les formules suivantes permettent d’intégrer les monômes en xi
1 0 si i impair
 x1i1 
 
1
1 i 1 
Cas d’une dimension :  x dx1     i  1 1   1  2
i
1
1  i  1  1  i  1 si i pair
Cas de deux dimensions :
1 1 0 si i ou j impairs

Pour l’élément de référence carré   x1 x2 dx1dx2  
i j
4
11  i  1 j  1 si i et j pairs


1 1 x1
i! j!
  x x dx dx  i  j  2!
i j
Pour l’élément de référence triangulaire 1 2 1 2
0 0
Cas de trois dimensions :
1 1 1 0 si i ou j ou k impairs

Pour l’élément de référence cubique    x x x dx1dx2 dx3  
i j k
1 2 3 8
11 1  i  1 j  1k  1 si i , j et k pairs


41
 1 1 x1 1 x1  x 2
i! j! k!
   x x x dx dx dx  i  j  k  3!
i j k
Pour l’élément de référence tétraédrique 1 2 3 1 2 3
0 0 0

1 1 1 x1 0 si k impair

 0 x x x dx1dx2dx3   2 i! j!
i j k
Pour l’élément de référence prismatique 1 2 3
1 0  k  1i  j  2 ! si k pair

4.6.2 Intégration numérique
Les formules d’intégration numérique permettent d’évaluer les intégrales sous la forme générale
suivante :
   iGr x1i , x2i , x3i 
NPI
~ ~ ~
vr

 r 1 2 3 r 
G x , x , x  d v  G X 
1 1x ,
vr
x 2 , x3 , X 
2 1x , x 2 , x3 , X 
3 1x , x 2 , x3  det J dx1 dx 2 dx3 
i 1
NPI : nombre total de points dits d’intégration numérique sur l’élément de référence
x1i , x2i , x3i  : coordonnées paramétriques du point i
i : coefficient (ou poids) pour le point i
 
Gr x1i , x2i , x3i : valeur numérique de la fonction à intégrer au point i
Choisir un schéma d’intégration consiste à définir les valeurs de i et x1i , x2i , x3i . Les formules les
plus utilisées dans la méthode des éléments finis sont celles de Gauss pour les problèmes à une
dimension, deux dimensions (élément de référence carré), trois dimensions (élément de référence
cubique) et celles de Hammer pour les triangles et les tétraèdres.
Le nombre de points NPI dépond des fonctions à intégrer. Ces formules permettent d’intégrer
exactement des monômes en xi mais fournissent des intégrations approximatives de fonctions
rationnelles.
Par exemple, un schéma de Gauss à n points intègre exactement un monôme en x1 d’ordre 2n-1.
Les tableaux 1 et 2 donnent les poids i et les coordonnées x1i , x2i  pour les schémas d’intégration
de Gauss pour les éléments à une dimension et de Gauss-Hammer pour les triangles.

Pour les éléments de référence linéique, carré et cubique le schéma de Gauss s’écrit :
 
m
Une dimension :  Gr x1 d vr   i Gr x1i où m représente le nombre de points dans la direction
vr
i 1
x1. (NPI=m)
Deux dimensions :  Gr  x1 , x2  d vr   i j Gr x1i , x2j  où m représente le nombre de points dans
m m

vr
j 1 i 1

les directions x1 et x2. (NPI=m2)

42
 25
1  3  7  9 
1  2  3  4  1 81
40 64
2  4  6  8  et 5 
81 81
Trois dimensions :  Gr x1 , x2 , x3  d vr   i jk Gr x1i , x2j , x3k  où m représente le nombre de
m m m

vr
k 1 j 1 i 1

points dans les directions x1, x2 et x3. (NPI=m3)

4.6.3 Intégration réduite

On parlera d’intégration réduite lorsque le nombre de points NPI retenu ne permet pas une
intégration exacte des différents termes.
Un nombre de points suffisant doit être retenu pour éviter la singularité de la matrice de rigidité
après assemblage et prise en compte des conditions aux limites. Cette singularité est associée à la
présence de modes dits parasites (modes à énergie nulle ou Hourglass en anglais (sablier)).
Dans les logiciels de calcul par éléments finis, l’utilisateur a la possibilité de choisir le type
d’intégration numérique.
Il existe trois types qui sont l’intégration complète, l’intégration réduite et l’intégration sélective
(SRI). Dans ce dernier type on intègre d’une façon complète les composantes déviatorique du
tenseur des contraintes et on intègre d’une façon réduite la composante hydrostatique.
Les avantages d’une intégration réduite c’est le gain en temps de calcul de la matrice de rigidité et
l’élimination du blocage numérique dû à l’incompressibilité induite dans les éléments à 4 nœuds
(surfaciques) et à 8 nœuds (volumiques).
43
Exercices :
1
1/ Soit l’intégrale suivante : I   x12  x13dx1
1
a/ Calculer la valeur exacte de I.
b/ Estimer numériquement la valeur de I en utilisant 1 point de Gauss. Calculer l’erreur commise.
c/ Estimer numériquement la valeur de I en utilisant 2 points de Gauss. Calculer l’erreur commise.
d/ Estimer numériquement la valeur de I en utilisant 3 points de Gauss. Calculer l’erreur commise.
1 1
1  x1
2/ Soit l’intégrale suivante : I    dx1dx2
11
2  x2
a/ Calculer la valeur exacte de I.
b/ Estimer numériquement la valeur de I en utilisant 1 point de Gauss. Calculer l’erreur commise.
c/ Estimer numériquement la valeur de I en utilisant 4 points de Gauss. Calculer l’erreur commise.
d/ Estimer numériquement la valeur de I en utilisant 9 points de Gauss. Calculer l’erreur commise.

44
 Tableau 4.3 : Formules d’intégration numériques sur un triangle (Hammer)

4.7 Conclusion
Dans les logiciels, on dispose d'une bibliothèque d'éléments prédéterminés (et on a parfois la
possibilité d'en définir soi-même). Les polynômes de base des interpolations de référence sont
calculés à l'avance et préprogrammés. L'utilisateur ordinaire n'a donc plus à s'en préoccuper.
L'objectif de ce chapitre n'était que de faire comprendre comment on les construit, et surtout de
permettre à l'utilisateur d'un logiciel de les choisir en connaissance de cause.
Il suffit de retenir que choisir un élément dans une bibliothèque c'est à la fois :
- Choisir une forme de maille
- Choisir la place des noeuds dans la maille
- Choisir une interpolation sur l'élément de référence

45
- Choisir une approximation géométrique de la forme de l'élément réel
Les fonctions d'interpolation sur l'approximation géométrique, leurs dérivées, et leurs
intégrales sur l'approximation géométrique sont alors déterminées par les valeurs aux noeuds.
Le choix des éléments est délicat : c'est un compromis entre la qualité de la solution approchée et le
coût du calcul.
Le nombre de noeuds (et donc le nombre d'inconnues égal au nombre de noeuds multiplié par le
nombre d'inconnues par noeud) augmente avec le degré d'interpolation et le nombre de mailles. La
taille du système d'équations à résoudre (et donc le temps de résolution) peut devenir très grande (le
temps de résolution est proportionnel au cube du nombre d'inconnues).
On peut dire que les éléments de faible degré d'interpolation demandent un maillage plus fin que les
éléments de degré plus élevé pour obtenir une «précision équivalente». Il faut donc se préoccuper
de cette question dans la phase de maillage, et l'examen d'une solution approchée peut parfois
amener à recommencer un calcul avec un maillage et/ou des interpolations différents.

46
CHAPITRE 5 : LA DISCRETISATION

Que le problème à résoudre soit formulé sous forme différentielle ou sous forme intégrale, sa
solution exacte est toujours inaccessible. L'objet de ce chapitre est de montrer comment choisir
une solution approchée ; dans la famille F  H1()nddl de fonctions ; construite par les opérations
de maillage et de choix d'éléments. On rappelle que chaque fonction F(M) de F peut s'écrire sous
forme d’une interpolation:
N ddl
FM    u   F   M 
j j

j 1

et que les fonctions de base F  j M  de l'espace F sont des fonctions nulles partout, sauf dans les
éléments qui contiennent le noeud j (On rappelle que ces fonctions de forme ne sont pas
polynomiales en général, sauf si l'interpolation géométrique  est strictement linéaire).

Figure 5.1 : Exemple de fonction de base globale F  j M  sur un maillage avec des éléments triangulaires.

5.1 Principe de la discrétisation

On a vu dans le chapitre précédent que la formulation variationnelle exacte d'un problème peut se
mettre sous la forme :
- Trouver u(M) tel que

A(, u) = B()  (M) défini dans 

- avec les conditions aux limites (éventuellement vides)

C ’ (u) = 0
Si on cherche une solution appartenant à l'espace F, on obtient le problème approché suivant
- Trouver F(M) tel que

A(, F) = B()  (M) défini dans   (5.1)

- avec les conditions aux limites :

C ’(F) = 0

Ce problème n'a pas en général de solution (sauf si par chance la solution exacte appartenait à F).
On va donc remplacer la condition «   (M)» par une condition moins contraignante :

47
On ne cherche à satisfaire l'équation (5.1) que pour certains  (M).
D'autre part, il faut respecter les conditions aux limites C ’(F) = 0 si elles existent

On restreint donc l'espace F aux seules fonctions qui respectent ces conditions. On appelle Fad
l'espace des fonctions admissibles, c'est à dire l'espace des fonctions qui satisfont aux conditions
aux limites C ’(F ) = 0.
Fad est un sous espace de F. En effet, imposer les conditions C ’(F) = 0 revient à imposer NC
relations sur les valeurs des u(i). La dimension de l'espace des fonctions admissibles est donc :
Ninc = Nddl - NC.
Pour déterminer une fonction F(M)  Fad, il faut déterminer un Ninc -uplet de valeurs :

{u(1), u(2), ..., u(i), ..., u(Ninc)}

Si on a choisi correctement Ninc fonctions (i)(M), pour chacune d'elles, l'équation (5.1) va fournir
une équation scalaire :
N inc
A(i, F) - B(i) = 0 où F  u j F  j (5.2)
j 1

et 1≤ i ≤ Ninc
La construction de ce système d'équations s'appelle l'assemblage.
Les fonctions i choisies sont souvent appelées fonctions test.
Pour trouver une solution approchée, on a donc un système de Ninc équations à Ninc inconnues {u(1),
u(2), ..., u(i), ..., u(Ninc) } à résoudre.
Remarques :
1°- Si A (, u) est une application linéaire en u, les termes A (i), F) conduisent à un système
linéaire en u(i). Cette classe de problèmes est appelée problèmes linéaires.
2°- Si A (, u) n'est pas linéaire en u, la solution du système (5.2) est plus complexe. Il existe
cependant des méthodes numériques pour le résoudre qui se ramènent à une succession de
problèmes linéaires.
3°- On peut noter que l'unicité et la convergence de la solution par éléments finis n'est
classiquement établie que pour les problèmes linéaires. Dans presque tous les autres cas, la méthode
des éléments finis reste heuristique.
Tout le problème réside donc dans le choix judicieux des fonctions i(M)
- Une condition impérative est que le choix des i(M) doit être tel que le système
d'équations (5.2) soit régulier.
En particulier, le choix des i(M) et le choix de l'espace des fonctions d'interpolation F(M)
ne sont pas indépendants du choix de la formulation variationnelle : En effet, à chaque
transformation d'intégrale sur , on diminue l'ordre de dérivation de u et on augmente l'ordre de
dérivation de . Il faut donc que les fonctions F et  aient des dérivées d'un ordre suffisant non
nulles partout pour que le système soit régulier (Par exemple, dans la formulation (3.20), la fonction
 apparaît sous la forme de son laplacien. Il faut choisir des fonctions  à laplacien non nul partout.
Une interpolation linéaire ne conviendrait pas).
- On dispose d'algorithmes efficaces pour la résolution de systèmes linéaires symétriques.
On préfèrera donc un choix de i(M) qui conduit à un système symétrique.

48
Il existe plusieurs méthodes pour choisir les fonctions i(M). Elles peuvent conduire à des
systèmes symétriques ou non, et la précision de la solution est plus ou moins sensible au choix des
. On ne développe ici que la méthode la plus utilisée dans les logiciels d'éléments finis : la
méthode de Galerkin.

5.2 La méthode de Galerkin

On choisit comme fonctions i(M) les fonctions de base de l'espace Fad. Ce choix a les
conséquences suivantes :
- Le nombre de fonctions i (et donc le nombre d'équations scalaires) est automatiquement
égal au nombre d'inconnues Ninc.
- On montre en analyse numérique que le système d'équations (5.2) obtenu pour ce choix de i
est nécessairement régulier pour une classe de problèmes appelés problèmes elliptiques : dans ce
cas, A(, u) et B() doivent être linéaires en u et en ,
 N inc

A(i, F) = A  F i  ,  u  j  F  j  
 j 1 
N ddl
et A (F(i),  u   F   ) doit être une forme quadratique définie positive en F(j).
j 1
j j

Pour les autres problèmes, soit le système d'équations n'est pas linéaire en u(j), soit il est
linéaire mais non garanti régulier, et la méthode de Galerkin peut se révéler inadéquate.
- Le fait de choisir des fonctions ) appartenant à Fad peut simplifier certaines intégrales de
bord.
- La méthode de Galerkin impose une contrainte sur le choix de la forme variationnelle : les
fonctions )(M) et les fonctions d'interpolation F(M) appartenant au même espace de fonctions, il
faut choisir une forme variationnelle dans laquelle les dérivées des fonctions  et u soient d'un
ordre tel qu'elles ne s'annulent pas partout.
La recherche de la solution approchée se ramène donc à la résolution du système d'équations
 N inc

A  F i  ,  u  j  F  j   - B(F(i)) = 0
 j 1 
(j)
Les u trouvés déterminent la solution approchée F . Ce système est linéaire et régulier pour les
problèmes elliptiques.
Cette écriture peut être écrite sur chaque sous domaine (élément fini) :

Ecriture sur un élément

On note Pe( i )  X  la fonction de forme liée au nœud i associée à l’interpolation sur l’élément réel (e)
Ne
~
L’interpolation sur l’élément réel e s’écrit f e  X    u  j  Pe j  ( X ) avec Ne le nombre de nœud dans
j 1

l’élément e
 i  N e  j  j )  Nel
 
Nel

 A Pe ,  u Pe    B Pei   0

e 1 

j 1  e 1
Si A est un opérateur linéaire alors :
Ne
   
Nel Nel

  u  j  A Pei , Pe( j )   B Pei   0

e 1 j  1 e 1

 
On note keij  A Pei  , Pe( j ) les termes de la matrice de rigidité élémentaire

49
  
qe(i )  B Pei  la composante du vecteur des forces nodales du nœud (i) associé à l’élément (e)
(un nœud peut appartenir à plusieurs éléments).

Ecriture matricielle
La formulation variationnelle d’un problème physique après transformation des intégrales et
l’utilisation d’une partie des conditions aux limites s’écrit sous la forme :


 D   D dv   E    F  d   f dv   T d
  T

Où D, E et F sont des opérateurs différentiels et f et T des données du problème.

La solution approchée recherchée est définie par morceaux sur les sous domaines générés par le
maillage à savoir     e
e
L’équation précédente devient alors :
Ne Ne

  D   D dv   E    F  d    f dv   T d

e1  e e 1  e
 e  e T

En tenant compte de l’interpolation éléments finis du champ inconnu et de la fonction test (qui sont
de même nature), on peut écrire les quantités suivantes :
D    A p  et D    A p 
Avec [A] la matrice d’interpolation associée à l’opérateur D
   
p 1 .. .. ..  p
T
vecteur des degrés de liberté de l’élément associé au champ inconnu
   
p 1 .. .. ..  p
T
vecteur des degrés de liberté de l’élément associé à la fonction test
L’intégration sur l’élément devient :
 
 D   D dv     A A dv      A Adv      k  
T T T T T
p p p p p e p
e e  e 
Avec ke     A  Adv matrice de rigidité élémentaire (d’un élément).
T

e

On réalise la même démarche pour le terme suivant :

f    p  N  f
T T

L’intégration sur l’élément devient :

 
 f *dv    p  N  fdv   p    N  fdv    p  Fe 
T T T T T

e e  e 
Avec Fe    N  fdv vecteur élémentaire des forces nodales (d’un élément).
T

e

5.3 Application de la méthode de Galerkin au problème thermique stationnaire

On a trouvé trois formes variationnelles pour le problème thermique stationnaire.

5.3.1 Forme l
- Trouver le champ de températures T tel que :

 T  f  d v
V
 0   M  défini sur V

- avec les conditions aux limites

T(N) =Tf (N)  N  ST

50
 grad T(N)  n(N) = f(N)  N  Sf

La fonction inconnue T(M) apparaît par son Laplacien. Il faut donc choisir un espace d'interpolation
F tel que le laplacien de ses fonctions de base soit non nul, ce qui interdit les interpolations
linéaires. D'autre part, la construction de l'espace Fad respectant la condition de flux est difficile : Il
faut choisir des éléments qui contrôlent la dérivée normale au bord, ce qui est difficile
Par contre, l'espace des fonctions test (qui est aussi Fad) paraît "trop luxueux", car  n'est pas dérivé
dans cette formulation.
La forme A (, T ) est bien linéaire en  et T . Le système en T(i) est donc linéaire et régulier. Bien
que la méthode de Galerkin appliquée à cette formulation paraisse «déséquilibrée» du point de vue
des dérivées de  et T , elle est praticable si on dispose d'éléments adéquats.

5.3.2 Forme 2
- Trouver le champ de températures T tel que
  grad   grad T d v    grad T  n d s     f d s   f dv  
V ST Sf V

A(, T) B(
- avec la seule condition aux limites
T (N) = Tf (N)  N  ST
Les fonctions  et T n'apparaissent que sous la forme de leur gradient. On peut donc choisir un
espace d'interpolation F plus simple : il suffit que les fonctions de base de F aient un gradient non
nul. On peut donc prendre des éléments plus simples : des interpolations linéaires sont suffisantes
(on peu prendre des interpolations de degré supérieur pour améliorer la qualité de la solution
approchée).
L'espace Fad est plus facile à construire : il suffit d'imposer des valeurs aux T(j) sur la frontière ST.
D'autre part, puisque   Fad, on a

ST
 grad T  n d s   ST
Tf grad T  n d s

Cette intégrale est plus facile à calculer.

La forme A (, T ) est toujours linéaire en  et T . Le système algébrique

N inc N inc
  grad Fi  
V
 T  jgrad F j d v 
j1
 ST
Tf  T  grad F   n d s 
j1
j j

  F i   f d s   F i  f d v (5.3)
Sf V

est linéaire en T(i). Cette formulation variationnelle est bien adaptée à la méthode de Galerkin. Les
intégrales sur V se ramènent à une somme d'intégrales sur les mailles. Par exemple :
Nm

V
F i  f d v  
m 1
m
F i  f d v

où Nm est le nombre de mailles. Il est à noter que les intégrales de volume et de surface de
l'équation i ne portent finalement que sur les fonctions de base de l'espace d'interpolation (On
rappelle aussi que, si on utilise des éléments de référence, les intégrales sur un élément réel se

51
ramènent à une intégrale sur l’élément de référence, et qu’elles peuvent être partiellement préparée
à l’avance, voir section 4.5).
Ces intégrales sur les mailles ne sont non nulles que sur les mailles qui contiennent le noeud i. On
ne les calcule donc que sur les mailles qui contiennent ce noeud. Les inconnues qui interviennent
dans l'équation i sont donc uniquement les T(k) tels que le noeud k appartient à une maille contenant
le noeud i. Il en est de même pour les intégrales de surface. Il reste donc :

V
F i  f d v  
mI
m
F i  f d v

où I est l'ensemble des mailles contenant le noeud i.

5.4 Application de la méthode de Galerkin au problème de statique de solide

déformable

Pour les mêmes raisons que dans l'exemple précèdent, on choisit la formulation variationnelle :
- Trouver le champ de déplacements  tel que

    grad  d v       n d s      f d v     q 0 d s   M 
V Sd V Sf

A(,) B(
 = 2  +  Tr G et   grad   grad T  .
1
où
2
- avec la condition aux limites :

(N) -  (N) = 0  N  Sd
C ’()
Dans cette formulation, l’inconnue  et les fonctions  apparaissent tous les deux sous la forme de
leur gradient.
D’autre part,  appartenant à Fad, on a :

Sd
    n ds  Sd
0    n d s

Le problème discret est donc :

Trouver les Ninc valeurs (i) tel que :

~  sym grad Fi  d v  i   

 0 ~  n d s    F  f d v   F  q 0 d s  F base de Fad (6.4)
i  i  i 
 
V Sd V Sf
N inc
~    jSym grad F  j .
où ~  2  ~   Tr ~ G
 et  j 1

Là encore, les intégrales de volume et de surface de l’équation i se réduisent à des intégrales sur les
seules mailles qui contiennent le nœud i.

5.5 L’assemblage dans la méthode de Galerkin

L’assemblage est l’opération qui consiste à construire le système d’équations à résoudre.

Le fait de choisir comme fonctions  les fonctions de base de l’espace Fad simplifie beaucoup le
calcul des intégrales. Comme on l’a vu dans les exemples précédents, l’équation relative à la

52
fonction de base F(i) ne fait intervenir en fait que les intégrales sur les mailles qui contiennent le
nœud i. On en déduit une méthode pour construire la matrice [K] du système à résoudre :

 u 1   b 1 
K        
u  Ninc   b  Ninc  
   

Pour chaque maille, on calcule les intégrales de la formulation. Pour chaque F(i) dans la maille, on
porte le coefficient de u(j) en ligne i colonne j de [K] en l’ajoutant à ce qui s’y trouve déjà, et on
porte en ligne i de [b] le résultat des intégrales du second membre toujours en l’ajoutant à ce qui s’y
trouve déjà.
L’exemple suivant illustre l’opération d’assemblage pour des matrices élémentaires 6x6
N
 Wi   q em K em q em  t
t t

m 1

k111 1
k12 k131 k141 1
k15 k161   u1  k112 k122 k132 k142 k152 k162   u 2 
 1 1 1 1 1     
 k 22 k 23 k 24 k 25 k 26   v1   k 222 k 232 k 242 2
k 25 k 262   v 2 
 1
k 33 1
k 34 1
k 35 1 
k 36     k 332 k 342 k 352 2 
k 36  2 
u1 v1 1 u 2 v 2  2  1 1 1 
1
  u 2 v 2  2 u 3 v3  3      .. 
 k 44 k 45 k 46  2 u  k 442 2
k 45 k 462   u 3 
 1
k 55 1 
k 56 v 2   k 552 2 
k 56 v3 
     

1
k 66   2   k 662  3 
k11i k12i k13i k14ik16i   u i 
k15i k11N k12N k13N k14N k16N   u n 1 
k15N
 i     

i
k 22 i
k 23 i
k 24k 26   vi 
i
k 25  k 22N k 23N k 24N k 25N
k 26N   vn 1 
 k 33 k 34 k 35 k 36  i 
i i i i   k 33N N N N
k 34 k 35 k 36  n 1 

 u i vi  i u j v j  j   i i i 
u
  ...  u n 1 vn 1  n 1 u n vn  n   
 k 44 k 45 k 46   j   k 44N k 45N k 46N   u n 
 k 55i k 56i   v j   k 55N k 56N   v n 
     
  j  k 66N    n 
i
 k 66 
u 1 v1 1 u 2 v2  2 u 3 v3  3 u i vi  i . u j v j  j .u n 1 vn 1  n 1 u n vn  n 
k 1
11 k 1
12 k 1
13 k 1
14
1
k
15 k 1
16 
 1 1 1 1 1 
 k 22 k 23 k 24 k 25 k 26 
 k 1
33 k 1
34 k 1
35 k 1
36

 

1
k 44  k112 1
k 45  k122 1
k 46  k132 k142 k152 k162 
 k k1 2
k k
1 2
k 2
k 2
k 262 
 55 22 56 23 24 25

 k k
1
66
2
33 k 2
34 k 2
35 k 362 
 k 2
k 2
k 462 
 44 45 
 k 2
55 k 562 
 
 k 662 
 
 
 k11i k12i k13i k14i k15i k16i 
 i
k 22 i
k 23 i
k 24 i
k 25 i
k 26 
 i i

 k 33 k 34 k 35i i
k 36 
 
 
 
 
 
 
 
 
 i
k 44 i
k 45 i
k 46 
 
 k 55i i
k 56 
 i
k 66 
 
 
 k11N k12N k13N k14N k15N k16N 

 k 22N k 23N k 24N k 25N k 26N 
 
 k 33N k 34N k 35N k 36N 
 k 44N k 45N N 
k 46
 
 k 55N k 56N 
 N
k 66 


u1 v1 1 u2 v2 2 u3 v3 3 .... u i vi i .... ... u j vj j u n 1 vn 1  n 1 un vn n 

53
5.6 Résolution
Après assemblage de la matrice de rigidité et du vecteur des forces nodales on trouve un système à
résoudre sous la forme : K q   F  avec q  le vecteur des ddl de la structure entière

D’une façon générale, on décompose d’une part le vecteur déplacements en deux vecteurs
qI connu ( qI  q0  ) et q II inconnu et d’autre part le vecteur forces en deux vecteurs
F I inconnu (les réactions) et F II connu :
q F  K  K 12   qI   F I 
q   I  et F    I  . Le système à résoudre devient :  11   
q II  F II  K 21 K 22  q II  F II 
K  q   K 12 q II  F I
Equivalent à  11 I (S)
K 
 21 I q   K 22 q II  F II

Puisque le vecteur q I  q0  connu alors le deuxième système devient : K 22 q II  F II  K 21 q0 
(le second membre est connu)
Pour retrouver les réactions, on utilise la première équation de (S)
F I  K 11 q0   K 12 q II

Remarque :

Les conditions aux limites sont souvent plus complexes que celles qui ont été données dans les
exemples proposés précédemment. On peut par exemple imposer des relations entre certains degrés
de libertés (par exemple pour simuler des nœuds liés à un solide indéformables). Le sous espace Fad
est alors plus difficile à construire. En fait on ne le construit pas et on rajoute des équations pour
forcer les conditions aux limites. Le nombre d’inconnus est donc Nddl et non Ninc. Il existe en
analyse numérique des méthodes pour résoudre un système algébrique avec des contraintes sur la
solution : on peut en citer (entre autres) la méthode des multiplicateurs de Lagrange et la méthode
de pénalisation. Certains logiciels offrent ces choix en option.

54
CHAPITRE 6 : Eléments 2D : Elasticité plane
L’objectif de ce chapitre est la mise en œuvre des étapes de la M.E.F pour la résolution des
problèmes d’élasticité plane (2D) en mécanique des solides déformables en utilisant les éléments
surfaciques de Lagrange de continuité C0.

6.1 Les équations des problèmes plans :

On distingue deux types de problèmes plans :

Les contraintes planes où le tenseur des contraintes présente la forme suivante :
  11  x1 , x2   12  x1 , x2  0 
 
  12  x1 , x2   22 x1 , x2  0 
 0 0 0 

Ce type de problème est rencontré lorsqu’une dimension de la pièce est très faible devant les deux
autres dimensions (plaques minces) et le vecteur chargement en tout point est perpendiculaire à
cette direction (figure 6.1-a).
Les déformations planes où le tenseur des déformations se présente sous la forme suivante :
  11  12 0 
 
   12  22 0 
 0 0 0 

Ce type de problème est rencontré lorsqu’une dimension de la pièce est très grande devant les deux
autres (barrage) et le vecteur chargement en tout point est perpendiculaire à cette direction (figure
6.1-b).

Figure 6.1 : Exemples de problèmes plans : (a) contraintes planes ; (b) déformations planes

Soit un solide occupant un domaine D soumis à des forces de volume {B} en tout point et des
forces surfaciques {F} sur une partie de la frontière SF . L’autre partie de la frontière SU est soumise
à un déplacement imposé (figure 6.2).

55
 Figure 6.2 : description du problème plan
 
L’équation d’équilibre du milieu continu : div  B  0 sur V  ij , j  Bi  0 sur V
Les conditions aux limites :

Type efforts   n  F  sur SF  ij n j  Fi sur SF
 
Type déplacement u  d sur SU ui  di sur SU
Lois de comportement sous forme matricielle:

 xx 
 
Les composantes du tenseur contrainte écrites sous forme vectorielle σ   yy 
 
 xy 
Relation déformation-déplacement :
  
Le vecteur déplacement présente deux composantes : u  ui  vj
1  
 
Le tenseur des déformations est   u T u . Il s’écrit sous forme vectorielle
2
  xx   u x 
  
ε    yy    v y 
  2  u y  v x 
 xy xy   
Pour un problème en contraintes planes E*  E et  *  
Pour un problème en déformations planes E*  E 1  2 et  *   1   2    
Formulation variationnelle :
Le champ virtuel pris dans ce cas un champ de déplacement   u pour retrouver le principe
des travaux virtuels :

La première forme variationnelle est :   ij , j  Bi .ui dV  0 
V

Après transformation de l’intégrale :    ij .ui , j dV    ij .n j .ui dS   Bi .ui dV  0

V S V
L’introduction des conditions aux limites, sachant que S=SF+SU :
 
  ij .ui , j dV   Fi .ui dS    ij .n j .ui dS   Bi .ui dV  0
V SF SU V

Vu la symétrie du tenseur contrainte :   ij .ui , j dV    ij . ij dV

V V
Ecriture sous forme matricielle en adoptant un champ virtuel cinématiquement admissible à

zéro u  0 sur SU:
On a :  ij . ij   xx . xx   yy . yy   xy . xy   yx . yx   xx . xx   yy . yy  2 xy . xy  ε σ
T

    dV   u F dS   u BdV  0 avec      yy  xy 

T T T T
xx
V SF V

Cette équation représente l’équation intégrale forme faible où l’équation Principe des Travaux Virtuels
(PTV).

56
Puisque la loi de comportement n’est pas encore prise en considération, cette équation est valable pour
de nombreux problèmes 2D.
Introduisons la loi de comportement, on obtient l’équation suivante :
T
 
   E    dV   u F dS   u BdV  0
th T T

V SF V

6.2 Elément fini triangulaire à trois nœuds T3

On va appliquer la méthode des éléments finis à la résolution d’un problème de plaque sollicitée
en cisaillement (figure 6.3). Il s’agit d’un problème en contraintes planes.
Deux maillages seront étudiés : un maillage grossier contenant deux éléments triangulaires
(figure 6.4-a) (pour faire des calculs à la main) et un maillage plus fin avec24 éléments (figure
6.4-b).

Figure 6.3 : plaque mince sollicitée en cisaillement

Figure 6.4 : maillages de la plaque : (a) maillage grossier ; (b) maillage fin

L’élément triangulaire à 3 nœuds présente six degrés de liberté qui sont les deux translations
suivant les direction x et y en chaque nœuds (figure 6.5).

Figure 6.5 : Degrés de liberté de l’élément triangulaire T3.

6.2.1. Interpolation des déplacements

L’élément de référence est le triangle rectangle isocèle de côté 1 (figure 6.6).

Figure 6.6 : élément de référence triangulaire

Les fonctions d’interpolation sur l’élément de référence sont : (voir TD chapitre 3)
N1  1  x1  x2 N 2  x1 N 3  x2

57
L’interpolation est linéaire (polynôme de degré 1)
Les deux composantes du vecteur déplacement sont interpolées sur l’élément de référence
comme suit :
3
u  x1 , x2    ui N i  x1 , x2   u1 N1 x1 , x2   u 2 N 2  x1 , x2   u3 N 3 x1 , x2 
i 1
3
v x1 , x2    vi N i x1 , x2   v1 N1 x1 , x2   v2 N 2  x1 , x2   v3 N 3  x1 , x2 
i 1

Si on met le vecteur des déplacements nodaux sous la forme suivante (appelé aussi vecteur des ddl):

u1   u1 
v  v 
 1  1
u  u  x1 , x2   N1 0 N2 0 N3 0  u2 
u   2  alors le vecteur déplacement s’écrit :    
v2  v x1 , x2   0 N1 0 N2 0 N 3  v2 
u3  u3 
   
 v3   v3 
u  x1 , x2 
   N u N : matrice d’interpolation des déplacements sur l’élément de référence
v x1 , x2 

Pour construire l’interpolation sur l’élément réel, on doit définir une transformation entre
l’élément de référence et l’élément réel. L’élément T3 est un élément isoparamètrique, on
choisit la transformation géométrique qui utilise les mêmes fonctions d’interpolation pour
approcher la géométrie et le champ de déplacement :
3
x  x x1 , x2    x i  N i   x1 , x2   x1 1  x1  x2   x 2 x1  x 3 x2 
i 1

 : mx1 , x2   M x, y   x  x1 x  x1   x2 x 3  x1   x1  x1 x 21  x2 x 31
1 2

   
3
y  y x1 , x2    y i  N i   x1 , x2   y1  x1 y 2  y 1  x2 y 3  y 1  y 1  x1 y 21  x2 y 31
i 1

6.2.2. Interpolation des déformations

Les dérivées sur l’élément de référence sont :
 u1 
v 
 u   1
u ,  
 1   x1   1,1N 0 N 0 N 0  u2   1 0 1 0 0 0
    u   
2 ,1 3,1
    B u     0 0 1 0
u
u,2     N1, 2 0 N 2, 2 0 N 3, 2 0 v2   1 0
 x2  u3 
 
 v3 
u1 
v 
 v   1
 
 v,1   x1  0 N1,1 0 N 2,1 0 N 3,1  u2  0  1 0 1 0 0
    v        B' u     0 0 0 1
u
v,2    0 N1, 2 0 N 2, 2 0 N 3, 2  v2  0  1
 x2  u3 
 
 v3 

58
Avec B  matrice d’interpolation des dérivées sur l’élément de référence.
Les dérivées sur l’élément réel sont, d’après le chapitre 4:
1
dx dx 
~ ~ ~ ~ ~ ~ 1
  f  f    f r  f r   d x1 d x2    f r  f r   x 21 x 31 
       21 
  x  y    x1  x2   d y d y    x1  x2   y y 31 
 d x1 d x2 
~ ~ ~ ~
  f r  f r  1  y 31  x 31    f r  f r  T
   21   j 
  x1  x2  2 A  y x 21    x1  x2 
1
 
Avec A  x 21 y 31  x31 y 21 est l’aire de l’élément réel (obtenu par produit vectoriel)
2
 u 
u , x   x  u ,1  1 y
23
0 y 31 0 y12 0
   
u ,

 u
  j 
u ,
  j B  
 u  B  u  
2 A
 32 13 21 u
   
y  2   x 0 x 0 x 0 
 y 
 v 
v, x   x   v,1  1 0 y
23
0 y 31 0 y 12 
  
v ,

 v
  j  
v ,
  j B '  u    B 'u  
2 A
 32 13 21 
u
   y   2   0 x 0 x 0 x 
 y 
La matrice d’interpolation des déformations [D] est définie par:
  xx   u , x   y 23 0 y 31 0 y12 0 
ε    yy    v, y   Du  1  0 x32 0 x13 0 x 21 u
2  u ,  v,  2 A 32
x y 23 x13 y 31 x 21 y12 
 xy   y x  
 y 23 0 y 31 0 y12 0 
D  1  0 x 32 0 x13 0 x 21 
2 A 32
x y 23 x13 y 31 x 21 y12 

6.2.3. Matrice de rigidité élémentaire
Introduisons les expressions des déformations et des déplacements issus de l’interpolation éléments
finis dans la formulation faible, on trouve :
T
 th T
 
   E    dV   u F dS   u BdV  0
T

V SF V

Or le domaine V est divisé en plusieurs sous domaines Ve tel que V

e 1.. N
e V et  Ve  
e 1.. N

 
     E     dV   u BdV   u F dS   0
T th T T

e  Ve Ve Se SF 
u  x1 , x2 
On choisit, d’après la méthode de Galerkin, ε  Du et    N u
vx1 , x2 

  u D EDudV    u N BdV    u D σ dV   uNl FdS
e Ve
T T

e Ve
T T

e Ve
T T th

ef S e  S F
T

A partir de cette équation on définit les termes suivant :

ke    DT E DdV matrice de rigidité élémentaire
Ve

59
f    N  BdV vecteur élémentaire des forces volumiques
e
0 T

Ve

f    D  dV vecteur élémentaire des forces volumiques

e
th T th
induites par le champ de
Ve

température.
 
f e   Nl  F dS vecteur des forces nodales surfaciques.
T

Se S F

Le principe des travaux virtuels s’écrit alors :

   
 uT keu   uT f e0   uT f eth  u f e  
e e e ef

s11  s12  s13 

k e    D E DdV   D E D det( J )dVr  s12  s22  s23 
T T

Ve Vr
s13  s 23  s33 
  1  *   * 1  * 
 b
  i jb  a a
i j  
 b a
i j  ai b j 
 
*
Et   2   2 
sij 

4 A 1   *
*2
 1  *
  1  *
 
 b j ai  a j bi   ai a j  bi b j  
 2   2  
Avec :
ai   x j  xk et bi  y j  yk (a1=x3-x2=x32 b1=y2-y3=y23=-y32)

Lk est la longueur du côté contenant les nœuds Ni et Nj (voir figure 6.7)

 f10     f1th     f1   
  0  0 th
   
th
 
 th 
   
T  
f e   N BdV   f 2  f e   D   dV   f 2  f e   Nl  F dS   f 2 
T
 
Ve  f0 
 3 
Ve
   f th 
 3   
Se F  f 
 3   

Figure 6.7 : chargement surfacique et volumique sur un élément

Pour trouver L’équation de rigidité globale pour tout le domaine on peut utiliser la méthode
d’assemblage qui sera illustrée à travers l’exemple suivant :

60
Connectivité des éléments :
Elément 1 : nœud 1, nœud 2, nœud 3 (il faut parcourir l’élément dans le sens trigonométrique)
Elément 2 : nœud 1, nœud 3, nœud 4
Elément1 Elément 2

Matrice globale :

Introduction des conditions aux limites u1=v1=u4=0

Réaction selon x au nœud 1 ?

Réaction selon y au nœud 1 ?

Réaction selon x au nœud 4 ?

Devant chaque ddl connu on a une composante de réaction inconnue

D’une façon générale, on décompose d’une part le vecteur déplacements en deux vecteurs
qI connu ( qI  q0  ) et q II inconnu et d’autre part le vecteur forces en deux vecteurs
F I inconnu (les réactions) et F II connu :

61
 q  F   K  K 12   q I   F I 
q   I  et F    I  . Le système à résoudre devient :  11   
q II  F II  K 21 K 22  q II  F II 
K  q   K 12 q II  F I
Equivalent à  11 I (S)
K 21 q I  K 22 q II  F II
Puisque le vecteur q I  q0  connu alors le deuxième système devient : K 22 q II  F II  K 21 q0 
(le second membre est connu)
qII  K 221 F II  K 21q0 
Pour retrouver les réactions, on utilise la première équation de (S)
F I  K 11q0   K 12 q II  K 11q0   K 12 K 221 F II  K 21q0 
Résolution :
Ecrire le système à résoudre en éliminant les lignes et les colonnes 1, 2 et 7. Résoudre et obtenir les
déplacements u2, v2, u3, v3 et v4.

Détermination des contraintes :

Elément 1 :

Elément 2 :

Application :
On considère une plaque chargée dans son plan par 2 forces ponctuelles comme montrée sur la figure ci-
dessous :

La plaque est maillée par deux éléments T3

62
 on remarque une forte discontinuité de la
composante  yy du tenseur des contraintes

63
CHAPITRE 7 : Eléments finis plaques
Une plaque est un solide défini par une surface de référence, ou surface moyenne, plane (plan x y
par exemple) et une épaisseur h(x,y) . Cette épaisseur doit être petite par rapport aux autres
dimensions (extensions, rayons de courbure) de la structure à modéliser. La figure 7.1 illustre ces
propos. Cette plaque peut être constituée d’un matériau homogène ou être obtenue par l’empilement
de différentes couches de matériaux orthotropes.

Figure 7.1 : Géométrie d’une plaque.

La théorie des plaques adoptée dans la suite est une théorie dite de premier ordre basée sur les
hypothèses suivantes:
 Les sections droites ou planes : les points situés sur une normale à la surface moyenne
restent sur une droite dans la configuration déformée.
 La déformation transversale  zz nulle (pas de variation d’épaisseur).

7.1 Cinématique des plaques de Reissner- Mindlin

Les plaques de Reissner Mindlin sont basées sur la cinématique suivante :

Un point de la surface moyenne présente cinq variables cinématiques indépendantes qui sont :
 Les trois déplacements de membrane u et v les déplacements dans le plan moyen et w le
déplacement transversal
 
 Les deux rotations  x et  y de la fibre normale à la surface moyenne autour des axes i et j
(figure 7.2).

64
 Figure 7.2 : Cinématique de la déformation d’un élément plaque.
Soit un point P de la surface moyenne. Ce point va se transformer au point P’ tel que :
u  x, y 

PP'  u  x, y   v x, y 
w x, y 

Un point M de la plaque défini par PM  zk va se transformer après déformation au point M’
 x 0 z y

tel que P' M '  zk   y  0   z x
0 z z
Le vecteur déplacement est alors :

z y u  x , y   z y
 
MM '  MP  PP'  P' M '   zk  u  x, y    z x  v x, y   z x
z w
u  x , y   z y

On note U  v x, y   z x le vecteur déplacement d’un point quelconque de la plaque.
w

7.3 Champs de déformations et de contraintes

Le champ de déformation dans une plaque est défini par dans l’hypothèse des petittes
perturbations (HPP) :

 u , x  z y , x
1
u, y  z y , y v, x  z x , x  1  y  w, x 
2 2
1  t  1 
2
 
  U  U   u, y  z y , y  v, x  z x , x  v , y  z x , y
1
w, y  x  
2 2 


1
 y  w, x 
1
 w, y  x  0 

2 2

 u, x
1
u, y v, x  0   y , x 1
 y , y  x , x  0  0 0
1
 y  w, x 
2 2 2
1  1   
  u , y v, x  0  z   y , y  x , x  w, y  x  
1
v, y x ,y 0   0 0
2  2   2 
 0 0 0   0 0 0     w, 
1 1
w, y  x  0 
     2 y x
2 

Ce tenseur des déformations est divisé en déformations membranaires, déformations de flexion

et déformations de cisaillement transverse.
En écriture vectorielle :
   xx  yy 2 xy  e z 
e  u, x v, y u, y  v, x : déformation membranaire
z   z  y ,x   x , y  y , y  x , x : déformation de flexion
   y  w, x w, y  x : déformation de cisaillement transversal

Le champ de contrainte est déterminé à partir de la loi de comportement qui peut être isotrope
ou orthotrope.

65
  xx  xy  xz 
 
   xy  yy  yz 
 xz  yz 0 
 
Sous forme vectorielle :    xx  yy  xy et    xz  yz respectivement les contraintes
dans le plan et les contraintes de cisaillement transverse.

7.4 Deuxième forme variationnelle : principe des travaux virtuels

L’application de la deuxième forme variationnelle ou du principe des travaux virtuels en

considérant un champ de déplacement virtuel (fonction test) qui a la même forme que la
cinématique de déformation de la plaque, on obtient :
Wint  Wext  0
Avec Wint     : dV         dV
V V

Wext   U  f v dV   U Fs ds   u  f    mdS   u  f s    ms dS

V A Af
Sf

Où fv sont les forces de volumes, Fs les forces surfaciques.

h f vx  x, y, z  h
2
 f 
 f x, y     f vy x, y, z dz ,
2
mx, y    z  vy dz sont des efforts par unité de surface
h 
2 f x, y , z 
 h  f vx 
 vz  2

moyenne.
h  Fsx 
 f s    Fsy dz sont les forces surfaciques ramenées au contour de la surface moyenne où les
2

h  
2 F
 sz 
efforts surfaciques sont appliqués.
Le travail des efforts intérieurs s’écrit en tenant compte de la forme particulière du champ de
déplacement virtuel qui a la même forme que la cinématique des plaques, on obtient :
h2 h 2

Wint     : dV           dV     e    z       dV

V A A
h 2 h 2

   e N   z  M    T dV
A

 N x  h2  x   M x  h 2  xx  h
2 
Avec : N    N y     y dz , M    M y    yy  zdz et T        xz dz
        Tx 
  Ty  h 2  yz 
 N   h 2   M   h 2  
 xy   xy   xy   xy 
Nx, Ny, Nxy : efforts résultants de membrane (en N/m)
Mx, My, Mxy : efforts ou moments résultants de flexion (en mN/m)
Tx, Ty : efforts résultants de cisaillement ou efforts tranchants (en N/m)

La loi de comportement pour un matériau orthotrope s’écrit :

   H     0  et    H    0 
 H11 H 12 H 13 
 H 44 H 45 
Avec H    H 22 H 23  et H     sont les matrices de comportement relatifs aux
  H 55 
 H 33 
contraintes dans le plan et au cisaillement transverse.

 0  et  0  sont les contraintes d’origine thermique.

66
 Les relations entre efforts et déformations sont données par (en négligeant les contraintes
d’origine thermique) :

 N x  h2  x  h
   
 
2
N    N y      y dz   H e  H  zdz H m e  H mf  
 N  h 2   h
 xy   xy  2

 M x  h 2  x  h
   
   
2
M    M y      y  zdz   H ez  H  z 2 dz  H mf e  H f  
M   h 2   h
 xy   xy  2

h h
T  2
  2
T    x     xz dz   H  dz H c  
Ty  h 2  yz  h
2

7.5 Plaques de Kirchhoff : éléments DKT et DKQ

La théorie de Kirchhoff pour les plaques minces est basée sur l’hypothèse dite de conservation des
normales : les points matériels situés sur une normale à la surface xy avant déformation restent sur
une normale à la surface moyenne déformée.
On admet ainsi que la rigidité de cisaillement est très grande par rapport à la rigidité de flexion.
Cela implique que les déformations de cisaillement transversal sont négligeables par rapport aux
autres composantes :
   y  w, x w, y  x  0 0 et   0 0
 y   w, x et  x  w, y
u  x, y   zw, x

Le champ de déplacement est alors U  v x, y   zw, y
w

Les éléments DKT (Discrete Kirchhoff Triangle) et DKQ (Discrete Kirchhoff quadrilateral) ont 3
ou 4 nœuds et 3 degrés de liberté par nœud :
u n  wi  xi  yi i  1, n n=3 pour DKT et n=4 pour DKQ qui sont le déplacement hors plan w et
Les deux rotations suivant la normale à l’arrête et à l’élément  s et dans la direction de l’arrête  n
(figure 7.3-a).
Le déplacement hors plan est interpolé par rapport aux nœuds situés aux sommets de l’élément.
Les rotations sont interpolées par rapport à des nœuds situés aux milieux des arrêtes.
Le tableau de la figure 7.3-b montre les fonctions de forme.

67
Figure 7.3 : (a) ddl de rotation des éléments DKT et DKQ ; (b) fonctions de forme

68
69
Références Bibliographiques

Modélisation des structures par éléments finis

Jean-Louis Batoz, Gouri Dhatt
Hermès, Paris, 1995

La méthode des éléments finis

Jean Garrigues
Janvier 2002

Finite Element Methods and Their Applications

Zhangxin Chen
Springer, New York, 2005

Méthodes numériques appliquées à la conception par éléments finis

David Dureisseix
Université Montpellier 2
Septembre 2003

70