Calcul de Strcuture - Introduction

à la méthode des éléments finis-

Master 1 Génie Civil

Dr. Yacouba Simporé

[email protected]

27 novembre 2022
Table des matières

1 Rappels : Élasticité linéaire 4

1.1 Tenseurs fondamentaux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.1 Tenseurs de déplacements . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.2 Tenseur de déformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.1.3 Tenseur des contraintes . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2 Relation entre déplacement et déformation . . . . . . . . . . . . . . . . . . . . . 5
1.3 Relation entre contrainte et déformation . . . . . . . . . . . . . . . . . . . . . . 5

2 Formules intégrales 6
2.1 Méthode des résidus pondérés : Forme intégrale forte . . . . . . . . . . . . . . . 6
2.2 Forme intégrale faible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.1 Intégration par parties (dimension un) . . . . . . . . . . . . . . . . . . . 7
2.2.2 Formule de Green en dimension deux . . . . . . . . . . . . . . . . . . . . 8
2.2.3 Formule de Green en dimension trois . . . . . . . . . . . . . . . . . . . . 8
2.2.4 Forme intégrale faible en élasticité linéaire . . . . . . . . . . . . . . . . . 8

3 Notions générales – Interpolation nodale 12

3.1 Discrétisation du milieu continu en sous domaines . . . . . . . . . . . . . . . . . 12
3.2 Approximation nodale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.2.1 Définition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2.2 Construction nodale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.3 Transformation géométrique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.3.1 Principe des éléments de référence . . . . . . . . . . . . . . . . . . . . . . 13
3.3.2 Interpolation sur l’élément de référence . . . . . . . . . . . . . . . . . . . 14
3.3.3 Fonction forme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.4 Éléments isoparamétriques . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2
 3

4 Matrices de rigidités locales 18

4.1 Découpage de la forme intégrale . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4.2 Matrice Jacobienne [J] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

5 Expansion – Assemblage – Résolution 23

5.1 Représentation du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.1.1 Table de nœud . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.1.2 Tabel de connectivité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
5.2 Matrice de rigidité et vecteurs forces locaux . . . . . . . . . . . . . . . . . . . . 24
5.2.1 Principe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.2.2 Vecteur de localisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
5.2.3 Utilisation des vecteurs de localisation . . . . . . . . . . . . . . . . . . . 25
5.2.4 Utilisation des matrices et vecteurs expansés . . . . . . . . . . . . . . . . 26
5.2.5 Particularités du système finale . . . . . . . . . . . . . . . . . . . . . . . 27

6 Application à l’élasticité linéaire 28

6.1 Application de la forme générale à un élément barre . . . . . . . . . . . . . . . . 28
6.1.1 Calcul de [H] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.1.2 Calcul de la matrice [B] . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
6.2 Élément poutre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
6.2.1 Élément poutre en flexion uniquement . . . . . . . . . . . . . . . . . . . 30
6.3 Élément de poutre en tension et flexion combinées . . . . . . . . . . . . . . . . . 37
6.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
6.3.2 Matrice de rigidité et vecteur force . . . . . . . . . . . . . . . . . . . . . 37
6.4 Extension des formations en dimension 2 et en dimension 3 . . . . . . . . . . . . 38

3
Chapitre 1

Rappels : Élasticité linéaire

1.1 Tenseurs fondamentaux

1.1.1 Tenseurs de déplacements

On considère le tenseur Ui du champ de déplacement qui est un vecteur noté {Ui } et qui a pour
composantes, en trois dimensions  
 u 
{Ui } = v .
w
 

1.1.2 Tenseur de déformation

On considère le tenseur εij du champ des déformations (appelé tenseur de Green-Lagrange

linéarisé) qui est une matrice symétrique notée [ε] admettant, en dimension trois, les composantes
données par  
εx εxy εxz
[ε] = εxy εy εyz  .
εxz εyz εz

Notons qu’on a également :

γxy = 2εxy , γxz = 2εxz , et γyz = 2εyz .

Aussi, on note souvent


⟨ε⟩ = εx εy εz γxy γxz γyz

1.1.3 Tenseur des contraintes

On note le tenseur de champ σij des contraintes (appelé tenseur de Cauchy-Euler ) qui est
une matrice symétrique notée σ définie en dimension trois par
 
σx τxy τxz
[σ] = τxy σy τyz  .
τxz τyz σz

4
1.2. RELATION ENTRE DÉPLACEMENT ET DÉFORMATION 5

On le note souvent

⟨σ⟩ = σx σy σz τxy τxz τyz .

On rappel une surface élémentaire de vecteur normal ⃗n, dans le champ de contrainte [σ] suivie
→
−
des efforts par unité de surface T , tel que

{T } = [σ] · {n} .

1.2 Relation entre déplacement et déformation

On peut déterminer facilement les déformations à partir de la connaissance des déplacements
en utilisant les relations classiques suivantes :
 
∂u ∂u ∂v 1 ∂u ∂v
εx = ∂x
γxy = ∂y
+ ∂x
εxy = 2 ∂y
+ ∂x

∂v ∂u ∂w 1 ∂u ∂w


εy = ∂y
γxz = ∂z
+ ∂x d’où εxz = 2 ∂z
+ ∂x
.

∂w ∂v ∂w
 
εz = ∂z
γyz = ∂z
+ ∂y εyz = 1 ∂v
+ ∂w
2 ∂z ∂y

On peut exprimer cette même relation matricielle en utilisant les relations tensorielles
1
εij = (ui,j + uj,i ) .
2

1.3 Relation entre contrainte et déformation

On peut déterminer facilement les contraintes à partir de la connaissance des déformations en
utilisant la loi de HOOKE
[σ] = [H] · [ε].
La matrice [H] est appelée matrice de HOOKE et vaut, pour un matériau homogène et
isotrope, défini par son module de Young et son coefficient de poisson ν, par
 
1−ν ν ν 0 0 0
 ν
 1−ν ν 0 0 0 
E −  ; où a = 1 − 2ν .
 ν ν 1 ν 0 0 0 
[H] = 
(1 − 2ν)  0
 0 0 a 0 0  2
 0 0 0 0 a 0
0 0 0 0 0 a

La loi de HOOKE s’exprime sous forme tensorielle par

σij = λδij · εkh + 2G · εij ,

où
E
■ G est le module de cisaillement : G = ;
2 (1 + ν)
ν·E
■ λ est le module volumétrique : λ = .
(1 + ν) (1 − 2ν)
5
Chapitre 2

Formules intégrales

On rappelle qu’en élasticité linéaire, le point de départ de l’application de la méthode des

éléments finis (MEF) est un système d’équation aux dérivées partielles (EDP), assorti de
conditions aux limites :
∂τxy
 ∂σ ∂τxz
 ∂x + ∂y + ∂z + Fvx = 0
x

∂τxy
∂x
+ ∂σ ∂y
y
+ ∂τ∂zyz + Fvy = 0
+ ∂τ∂yyz + ∂σ

 ∂τxz
∂x ∂z
z
+ Fvz = 0

Sur Su , on impose le champ de déplacement ; et donc on a

   
 u   ū 
{Ui } = v = v̄ , avec ū, v̄, et w̄ connus.
w w̄
   

→
−
Sur Sf , on applique des efforts de surface fs qui vont influés sur le champ de contraintes à
l’endroit où on applique ses efforts. Cela s’exprime par la relation

 σx · nx + τxy · ny + τxz · nz = fsx

fs = [σ] · {n} ⇐⇒ τxy · nx + σy · ny + τyz · nz = fsy

τxz · nx + τyz · ny + σz · nz = fsz


2.1 Méthode des résidus pondérés : Forme intégrale forte

On rappelle que la solution recherchée en élasticité est un champ de déplacement u à travers
le domaine V . On appelle résidu, la quantité R (→−u ) qui est égale à
∂τxy
 ∂σ ∂τxz
 ∂x + ∂y + ∂z − Fvx
x

R (→
−
u)= ∂τxy
∂x
+ ∂σ ∂y
y
+ ∂τ∂zyz − Fvy
+ ∂τ∂yyz + ∂σ

 ∂τxz
∂x ∂z
z
− Fvz

Bien entendu, si Uex est une solution exacte du problème telle que modélisé, alors on a
 
0
R (Uex ) = 0
0
 

6
2.2. FORME INTÉGRALE FAIBLE 7

Figure 2.1 – Domaine V.

La méthode des résidus pondérés est une méthode qui consiste à exprimer que l’on recherche la
solution sous la forme d’un champ de déplacement U qui annule le résidu et qui par conséquent
annule l’intégrale
y
W = ⟨ψ⟩ {R (→
−
u )} · dv = 0,
V

où, ⟨ψ⟩ est un vecteur de fonction de pondération dont les composantes sont


⟨ψ⟩ = ψu ψv ψw .

On appelle cette expression forme intégrale du problème

y 
∂σx ∂τxy ∂τxz

W = ψu + + − Fvx · dv
V ∂x ∂y ∂z
y 
∂τxy ∂σy ∂τyz

+ ψv + + − Fvy · dv
V ∂x ∂y ∂z
y 
∂τxz ∂τyz ∂σz

+ ψw + + − Fvz · dv.
V ∂x ∂y ∂z

2.2 Forme intégrale faible

2.2.1 Intégration par parties (dimension un)

Formule de Green en dimension 1 :

Z x2 Z x2
du dψ
ψ · dx = − u · dx + [ψ · u]xx21 .
x1 dx x1 dx

7
8 CHAPITRE 2. FORMULES INTÉGRALES

2.2.2 Formule de Green en dimension deux

On considère ddeux fonction u(a, y) et ψ(x, y) définies à travers un domaine V du plan (x, y).
On a
x ∂u x ∂ψ I
ψ· dx dy = − · u dx dy + ψ · u · nx · dl
V ∂x V ∂x S
x ∂ψ I
=− · u dx dy + ψ · u · nx dl
V ∂x S
x ∂u x ∂ψ I
ψ· dx dy = − · u dx dy + ψ · u · ny dl
V ∂y V ∂y S

2.2.3 Formule de Green en dimension trois

y ∂u y ∂ψ {
ψ· dx dy dz = − · u dx dy dz + ψ · u · nx ds
V ∂x V ∂x S
y ∂u y ∂ψ {
ψ· dx dy dz = − · u dx dy dz + ψ · u · ny ds
V ∂y V ∂y S
y ∂u y ∂ψ {
ψ· dx dy dz = − · u dx dy dz + ψ · u · nz ds
V ∂z V ∂z S

2.2.4 Forme intégrale faible en élasticité linéaire

y ∂σx ∂τxy ∂τxz

 y 
∂τxy ∂σy ∂τyz

W = ψu + + dx dy dz + ψv + + dx dy dz
V ∂x ∂y ∂z V ∂x ∂y ∂z
y 
∂τxz ∂τyz ∂σz
 y
+ ψw + + dx dy dz − (ψu Fvx + ψv Fvy + ψw Fvz ) dx dy dz,
V ∂x ∂y ∂z V

avec
y 
∂σx ∂τxy ∂τxz
 y 
∂ψu ∂ψu ∂ψu

A= ψu + + dx dy dz = − σx + τxy + τxz dx dy dz
V ∂x ∂y ∂z V ∂x ∂y ∂z
{
+ (ψu · σx · nx + ψu · τxy · ny + ψu · τxz · nz ) ds,
S

y 
∂τxy ∂σy ∂τyz
 y 
∂ψv ∂ψv ∂ψv

B= ψv + + dx dy dz = − · τxy + · σy + · τyz dx dy dz
V ∂x ∂y ∂z V ∂x ∂y ∂z
{
+ (ψv · τxy · nx + ψv · σy · ny + ψv · τyz · nz ) ds
S

et
y 
∂τxz ∂τyz ∂σz
 y 
∂ψw ∂ψw ∂ψw

C= ψw + + dx dy dz = − · τxz + · τyz + · σz dx dy dz
V ∂x ∂y ∂z V ∂x ∂y ∂z
{
+ (ψw · τxz · nx + ψw · τyz · ny + ψw · σz · nz ) ds.
S

8
2.2. FORME INTÉGRALE FAIBLE 9

d’où
y ∂ψu

∂ψv ∂ψw

W =− · σx + · σy + · σz dx dy dz
V ∂x ∂y ∂z
y   ∂ψ ∂ψw
 
∂ψu ∂ψv
 
∂ψv ∂ψw

u
− τxz + + τxy + + τyz + dx dy dz
V ∂z ∂x ∂y ∂x ∂z ∂y
{ {
+ ψu (σx · nx + τxy · ny + τxz · nz ) ds + ψv (τxy · nx + σy · ny + τyz · nz ) ds
{S yS
+ ψw (τxz · nx + τyz · ny + σz · nz ) ds − (ψu Fvx + ψv Fvy + ψw Fvz ) dx dy dz
S V

Pour ce qui est de l’intégrale sur S (intégrale de surface), on rappelle que la frontière du domaine
est divisée en deux parties : Su (déplacements imposés) et Sf (effort de surface).
Ainsi
{ {
ψu (σx · nx + τxy · ny + τxz · nz ) ds + ψv (τxy · nx + σy · ny + τyz · nz ) ds
S S
{
+ ψw (τxz · nx + τyz · ny + σz · nz ) ds
{S
= ⟨ψ⟩ {fs } ds : conditions aux limites.
S

Finalement y {
W = A1 − ⟨ψ⟩ {fv } dx dy dz + ⟨ψ⟩ {Fs } ds.
V S

Hypothèse de Galerkine

Le choix des fonctions de pondération dans le cadre de l’hypothèse de Galerkine consiste à

prendre les variations des composantes du champ de solution recherchée :
⟨ψ⟩ = ⟨ψu , ψv , ψw ⟩ = ⟨δu , δv , δw⟩
où δ est l’opérateur de variation.

Propriétés de l’opérateur δ


δ(u + v) = δu + δv ; δ(cu) = cδu ; δ(u · v) = uδv + vδu ; δ(P ) = δP1 δP2 δP3 ;
 
∂u ∂(δu) y y 
δ = ; δu dx dy dz = δ u dx dy dz ;
∂x ∂x V V


avec u et v deux champs scalaires, et ⟨P ⟩ = P1 P2 P3 un champ vectoriel.

Forme intégrale faible avec l’hypothèse de Galerkine


On a ⟨ψ⟩ = δu δv δw , et donc
y  ∂(δu) ∂(δv) ∂(δw)

W =− σx + σy + σz dx dy dz
V ∂x ∂y ∂z
y   ∂(δu) ∂(δv)  
∂(δu) ∂(δw)
 
∂(δv) ∂(δw)

− τxy + + τxz + + τyz + dx dy dz
V ∂y ∂x ∂z ∂x ∂z ∂y
y 
  {

− δu δv δw τyz dx dy dz + δu δv δw {fs } ds
V S

9
10 CHAPITRE 2. FORMULES INTÉGRALES

en appliquant la dernière propriété de l’opérateur δ, on a :

y
W =− [δ(εx ) σx + δ(εy ) σy + δ(εz ) σz + δ(γxy ) τxy + δ(γxz ) τxz + δ(γyz ) τyz ] dx dy dz
y V {
− ⟨δu⟩ {Fv } dx dy dz + ⟨δu⟩ {fs } ds
V S
or


⟨ε⟩ = εx εy εz γxy γxz γyz =⇒ ⟨δε⟩ = δεx δεy δεz δγxy δγxz δγyz
et

⟨σ⟩ = σx σy σz σxy σxz σyz
Donc on a :
y y {
W =− ⟨δε⟩ {σ} dx dy dz − ⟨δu⟩ {Fv } dx dy dz + ⟨δu⟩ {fs } ds
V V S


avec ⟨δu⟩ = δu δv δw
La forme intégrale faible devient
y y {
W =− ⟨δε⟩ {σ} dx dy dz − ⟨δu⟩ {Fv } dx dy dz + ⟨δu⟩ {fs } ds
V V Sf

=0
où y
■ ⟨δε⟩ {σ} dx dy dz est le travail virtuel interne ;
yV −
→
■ ⟨δu⟩ {Fv } dx dy dz est le travail virtuel des efforts de volume Fv ;
{ V

■ ⟨δu⟩ {fs } ds est le travail virtuel des efforts de surface.

Sf

Théorème 2.1
La forme inétgrale faible avec l’hypothèse de Galerkine pour la plupart des problèmes
d’élasticité, correspond en fait aux théorèmes des travaux virtuels dans sa formulation de
déplacement :
y y {
− ⟨δε⟩ {σ} dx dy dz − ⟨δu⟩ {Fv } dx dy dz + ⟨δu⟩ {fs } ds = 0.
V V Sf

Prise en compte des forces ponctuelles

L’expression déterminée au paragraphe précédent s’applique que dans le cas où il n’y a pas de
forces ponctuelles qui sont appliquées.
Si tel est le cas, on peut adapter facilement cette expression en introduisant dans l’expression
le travail virtuel m
X
δUpi · Fpi ,
i=1
des m forces ponctuelles, ce qui donne :
y y { m
X
− ⟨δε⟩ {σ} dx dy dz − ⟨δu⟩ {Fv } dx dy dz + ⟨δu⟩ {fs } ds + δUpi · Fpi = 0
V V Sf
i=1

10
2.2. FORME INTÉGRALE FAIBLE 11

Exercice 2.1. On considère ci-dessous une zone dans l’espace réel (x, y).
On considère un problème physique décrit sur cette zone pour l’équation aux dérivées partielles
suivantes (équation de la chaleur ou équation de Poisson avex k = 1 et q = 1)

∂ 2T ∂ 2T
+ + 1 = 0.
∂x2 ∂y 2
T = 0 imposée sur le bord 2-3 de la zone
∂T
= 0 (flux nul) sur le bord 1-2, 1-3 de V .
∂n
La solution recherchée est la champ scalaire T (x, y) à travers V .
1. Établir la forme intégrale forte du problème.
2. Établir la forme intégrale faible du problème avec l’hypothèse de Galerkine.

Exercice 2.2. On considère une barre de longueur L = 1 et de module de Young constant E,

défini le long de l’axe x, entre x = 0 et x = 1.
Cette barre ne travaille qu’en tension ou en compression le long de l’axe x. La section de la
barre A est variable le long de l’axe et sa variation est linéaire entre A(0) = e et A(1) = 1.
On applique sur cette barre une force horizontale F en x = 1 et un blocage du déplacement u(x)
en x = 0.
La solution recherchée est le champ de déplacement u(x) entre x = 0 et x = 1.
Formaliser ce problème par une équation différentielle et des conditions de frontière.
Établir la forme intégrale forte du problème.
Établir la forme intégrale faible du problème avec l’hypothèse de Galerkine.

11
Chapitre 3

Notions générales – Interpolation nodale

Introduction
Les principales étapes de construction d’un modèle Éléments finis sont les suivantes :
■ Discrétisation du milieu continu en sous domaine ;
■ Construction de l’approximation nodale ;
■ Calcul des matrices élémentaires correspondant à la forme intégrale du problème ;
■ Assemblage des matrices élémentaires et prise en compte des conditions aux limites ;
■ Résolution du système d’équation linéaire.

3.1 Discrétisation du milieu continu en sous domaines

Cette opération consiste à décomposer le domaine continu en un nombre fini de sous domaines,
“éléments finis”. Le domaine D en question sera :
ne
X
D= De : lim (∪De ) = D.
e→0
e=1

Il faut donc pouvoir représenter au mieux la géométrie du domaine étudié, souvent complexe
par des éléments finis de forme géométrique simple. Il ne doit y avoir ni recouvrement, ni
trou entre deux éléments et un frontière commune. De plus, lorsque la frontière du domaine
est complexe, une erreur de discrétisation géométrique est inévitable. Cette erreur doit être
exprimer et éventuellement réduite en modifiant la forme ou en diminuant la taille des éléments
concernés.

3.2 Approximation nodale

La méthode des éléments finis est basée sur la construction systématique d’une approximation
U ∗ (U ∗ est une fonction qui approxime U ) du champ des variables U par sous domaines.
Cette approximation est construite sur les valeurs approchées du champ au nœud de l’élément.
On parle de représentation nodale ou simplement d’approximation nodale.

12
3.3. TRANSFORMATION GÉOMÉTRIQUE 13

3.2.1 Définition

L’approximation par éléments finis est une approximation nodale par sous domaines ne faisant
intervenir que les variables nodales du domaine.
∀M ∈ De , {U ∗ (M )} = [N [M ]] {Un } ,
■ {U ∗ } est la valeur de la fonction approchée en tout point M de l’élément De ;
■ [N ] matrice des fonctions d’interpolation de l’élément ;
■ {Un } variables nodales relatives aux nœuds d’interpolation de l’élément.

3.2.2 Construction nodale

L’approximation nodale est construite à partir d’une approximation générale :

∀M : U ∗ (M ) = ⟨Φ(M )⟩ {a} ,
■ ⟨Φ(M )⟩ est une base de fonctions connues indépendantes ;
■ {a} : est le vecteur de n paramètres de l’approximation.

Les n nœuds Mn sont des points de l’élément pour lesquels l’approximation U ∗ est identifiée à
la valeur du champ de variable.
   
Φ(M1 ) Φ(M1 )
{Un } = {U ∗ (Mn )} =  ...  {a} ; avec [T ] =  ..
   
. 
Φ(Mn ) Φ(Mn ) .
Si [T ] est inversible, on a
{a} = [T ]−1 {Un } .
En reportant ce résultat dans l’approximation, nous obtenons la matrice des fonctions d’interpolation
U ∗ (M ) = ⟨Φ(M )⟩ [T ]−1 {Un }
= [N (M )] {Mn } ; avec [N (M )] = ⟨Φ(M )⟩ [T ]−1 .

Base polynomiale


1D ■ Linéaire : 1 x ;


■ Quadratique : 1 x x2 ;


2D ■ Linéaire : 1 x y ;


■ .Quadratique : 1 x y x2 xy y 2 ;


3D ■ Linéaire : 1 x y z ;


■ Quadratique : 1 x y z x2 xy y 2 xz yz z 2 .

3.3 Transformation géométrique

3.3.1 Principe des éléments de référence

Le principe consiste à associer à chaque type d’éléments réels un élément unique de forme
constante, défini dans un espace de reférence (s, t, ν) ayant le même nombre de nœud.

13
14 CHAPITRE 3. NOTIONS GÉNÉRALES – INTERPOLATION NODALE

(x3 , y3 )

E. réel
(s3 , t3 )
(x1 , y1 ) (x2 , y2 )

Elément réference

(s1 , t1 ) (s2 , t2 )

Figure 3.1 – Elément réel-élément de réference.

3.3.2 Interpolation sur l’élément de référence

On rappelle qu’en dimension trois, l’interpolation sur l’élément réel d’une fonction u∗ (x, y, z) à
partir de n valeurs au nœud Ui s’établit par l’expression

NnR (x, y, z)⟩ u1 ... un

n o
u∗ (x, y, z) = N1R (x, y, z) · · ·

Toujours en 3D, l’interpolation nodale sur l’élément de référence dans l’espace de coordonnées
(s, t, ν) est
 
 u1 
 
u (s, t, ν) = ⟨N1 (s, t, ν) · · · Nn (s, t, ν)⟩ ...
∗

u  
n

où les Ni (s, t, ν) sont les facteurs d’interpolation sur l’élément de référence.

3.3.3 Fonction forme

L’utilisation d’élément de référence demande la connaissance de la relation entre élément réel

et élément de référence. Cette relation est représentée en 3D par trois fonctions

x(s, t, ν) , y(s, t, ν) , et z(s, t, ν) .

Les fonctions de forme, notées N i (s, t, ν), permettent d’établir facilement cette relation entre
les deux espaces de coordonnées et ceci pour tout type d’éléments finis.

14
3.3. TRANSFORMATION GÉOMÉTRIQUE 15

Pour une transformation en trois dimensions, on a toujours pour un élément à n nœud :

  

  x1 
  . 
..


x(s, t, ν) = N 1 (s, t, ν) · · · N n (s, t, ν)





 
 

x



 n





  
 x1 



  . 
..

y(s, t, ν) = N 1 (s, t, ν) · · · N n (s, t, ν)


 
 



 x n





  
 x1 




  . 
..

z(s, t, ν) = N 1 (s, t, ν) · · · N n (s, t, ν)






 
 

 x n

3.3.4 Éléments isoparamétriques

Pour la plupart des éléments finis utilisés dans la pratique, les fonctions d’interpolation sur
l’élément de référence sont égales aux fonctions formes.
Ces éléments sont appelés isoparamétriques. Nous utiliserons essentiellement les éléments
isoparamétriques dans ce cours. Dans ce cas, en dimension trois, on a les résultats suivants :

  
 x1 

 
..

 
x(s, t, ν) = ⟨N (s, t, ν) · · · N (s, t, ν)⟩



 1 n
 .
   
xn









  
 x1 


 
..
 
y(s, t, ν) = ⟨N1 (s, t, ν) · · · Nn (s, t, ν)⟩ .

 
  


 xn





  
 x1 



 
..
 
z(s, t, ν) = ⟨N1 (s, t, ν) · · · Nn (s, t, ν)⟩ .





 
  
 x n

avec les Ni les fonctions d’interpolation sur l’élément de référence.

Exemple 3.1. Pour l’élément de référence, les fonctions d’interpolation sont

N1 (s, t) = 1 − s − t ; N2 (s, t) = s ; N3 (s, t) = t .

15
16 CHAPITRE 3. NOTIONS GÉNÉRALES – INTERPOLATION NODALE

(0, 1)

Elément réference

(0, 0) (1, 0)

Figure 3.2 – Elément réel

Les relations qui lient l’élément de référence aux fonctions formes sont :
  


 D x1 
E 
x(s, t) = − −



 1 s t s t x 2
  

x3

 


  
y1 




 D E 



 y(s, t) = 1−s−t s t y2
  

y3
 

En résumé, pour exprimer l’interpolation sur l’élément de référence en dimension 2 et 3, on

procède de la manière suivante :
■ On exprime l’interpolation sur l’élément de référence (Fiche) ;
■ On exprime les relations entre les coordonnées cartésiennes (x, y, z) et (s, t, ν) à l’aide des
fonctions Ni , pour obtenir les relations

x(s, t, ν) ; y(s, t, ν) ; et z(s, t, ν) ;

■ On remplace les relations s(x, y, z), t(x, y, z) et ν(x, y, z) dans l’expression de l’interpolation
sur l’élément de référence pour obtenir finalement l’expression de l’interpolation sur
l’élément réel.

Exercice 3.1. Soit un élément triangulaire défini en dimension deux dans le plan (x, y). Les
coordonnées des nœuds de l’élément sont les suivantes :

nœud 1 : (0, 0) ; nœud 2 : (1, 0) ; nœud 3 : (1, 1).

16
3.4. CONVERGENCE 17

On interpole ensuite un champ de température T (x, y) à l’intérieur de cet élément. Les températures
connues sont respectivement T1 , T2 et T3 . En passant par l’élément de référence associé à cet
élément réel, déterminer l’expression de l’interpolation T (x, y).

Exercice 3.2. On considère une barre étudiée en transfert de chaleur le long de l’axe x. La
barre est définie entre x = 0 et x = 1. On mesure la température de cette barre en trois points :
■ nœud 1, pour x = 0, on a : T1 = 20 ;

■ nœud 2, pour x = 0.5, on a : T2 = 25 ;

■ nœud 3, pour x = 1, on a : T3 = 22.

On cherche à interpoler à partir de ces mesures, un champ de température T (x) pour x ∈ [0, 1].
Déterminer l’interpolation T (x) en utilisant un maillage avec un seul élément (1D quadratique) ;
la fonction d’interpolation est :
⟨N (s)⟩ = − 21 s(1 − s) 1 − s2 12 s(1 + s) .



3.4 Convergence
On peut garantir dans certaines conditions que la méthode converge vers la solution exacte du
problème (tel que modélisée) en déplacement, déformation ou contrainte.

Définition 3.1
On dira que l’on a convergence sur la solution, c’est-à-dire que u∗ (x, y, z) tend vers la solution
exacte uex (x, y, z), lorsque l’on fait tendre la taille des éléments vers 0, si on a

lim ∥u∗ (x, y, z) − uex (x, y, z)∥ = 0,

taille(De )→0

où l’opérateur ∥ · ∥ est une norme qui reste à définir.

Définition 3.2
On dira que l’on a la convergence sur les dérivées premières de la solution si on a :
∗
∂u ∂u ex
lim (x 1 , x2 , x3 ) − (x 1 , x2 , x3 ) = 0, pour i ∈ {1, 2, 3} .
taille(De )→0 ∂xi ∂xi

Théorème 3.1
■ Avec les tétraèdres linéaires en élasticité, on a convergence des déplacements, des
déformations et des contraintes.
■ Avec les tétraèdres linéaires en élasticité, on a continuité entre élément voisin sur les
déplacements, mais pas sur les déformations et les contraintes.
■ Pour avoir théoriquement convergence et continuité sur les déplacements, les
déformations, les contraintes en élasticité 3D, il faut au minimum utiliser des tétraèdres
quadratiques en 10 nœuds.

17
Chapitre 4

Matrices de rigidités locales et vecteurs

forces locaux

4.1 Découpage de la forme intégrale

On procède à la discrétisation de la forme intégrale faible en procédant de la façon suivante.
C’est-à-dire le découpage en intégrales sur les éléments et l’introduction de l’interpolation nodale
qui a pour objet de faire apparaître dans la formulation, les valeurs du champ de solution
recherchée au nœud du maillage.
Ainsi, on procède à la discrétisation de la forme intégrale, en écrivant que les intégrales sur le
domaine V et la surface Sf sont la sommation d’intégrale sur chaque élément.
y y {
W =− ⟨δε⟩ · {σ} dv − ⟨δu⟩ · {Fv } dv + ⟨δu⟩ · {fs } ds
V V Sf
N
X
= W(i) = 0
i=1

avec y y {
W(i) = − ⟨δε⟩ · {σ} dv − ⟨δu⟩ · {Fv } dv + ⟨δu⟩ · {fs } ds.
V (i) V (i) Sf

Comme {σ} = [H] · {ε}, donc

y y {
Wi = ⟨δε⟩ [H] {ε} dv − ⟨δu⟩ {Fv } dv + ⟨δu⟩ {fs } ds.
V (i) V (i) Sf (i)

On introduit maintenant l’interpolation nodale du champ de solution recherchée.

 
u
Soit n le nombre de nœuds de l’élément (i) {U } = v ;
w
 
     
u 1 v   w1 
  .1 
  

.

 ..

u = N1 · · · Nn .. ; v = N1 · · · Nn .. ; w = N1 · · · Nn .

u   
v   
w  
n n n


On introduit le vecteur Un(i) :


 
Un(i) = u1 v1 w1 · · · un vn wn ; {U } = [N ] Un(i) ,

18
4.1. DÉCOUPAGE DE LA FORME INTÉGRALE 19

avec  
N1 0 0 N2 0 0 N3 ··· Nn 0 0
 0 N1 0 0 N2 0 0 ··· 0 Nn 0 
[N ] = 
 ··· ··· ··· ··· ··· ··· ···

··· ··· ··· ··· 
0 0 N1 0 0 N2 0 ··· 0 0 Nn
De cette expression, on tire

{δU } = [N ] · δUn(i)

⟨δU ⟩ t {U } = t [N ] · δUn(i) = δUn(i) t[N ]







δUn(i) = δu1 δv1 δw1 · · · δun δvn δwn

Matrice B

L’interpolation nodale conduit également à la définition de la matrice B telle que


{ε} = [B] · Un(i) .

La matrice B est une matrice de taille 6 × 3n (avec n le nombre de nœud) qui est globalement
formé avec les termes de la matrice [J] (à définir). Ainsi

⟨ε⟩ = t {ε} = t [B] · Un(i) = Un(i) t[B],





de même
⟨δε⟩ = δUn(i) t [B].



Matrice de rigidité locale et vecteur force local

On a donc
y t y
δUn(i) t [N ] {Fv } dv



W(i) = − δUn(i) [B][H][B] Un(i) dv −
V (i) V (i)
{
t
+ δUn(i) [N ] {fs } ds
Sf (i)
 y y {
 
t t t


= δUn(i) − [B][H][B] Un(i) dv − [N ] {Fv } dv + [N ] {fs } ds
V (i) V (i) Sf (i)
 y
y {
 
t t t

 
= δUn(i) − [B][H][B] dv Un(i) − [N ] {Fv } dv + [N ] {fs } ds
V (i) V (i) Sf (i)

Posons
y y {
t t t
[Ki ] = [B][H][B] dv , et {Fi } = − [N ] {Fv } dv + [N ] {fs } ds ;
V (i) V (i) Sf (i)

alors Wi devient

  
Wi = δUn(i) −[Ki ] Un(i) + {Fi } .

19
20 CHAPITRE 4. MATRICES DE RIGIDITÉS LOCALES

Définition 4.1
La matrice [Ki ] est appelée matrice de rigidité locale pour l’élément (i) et vaut
y
t
[Ki ] = [B][H][B] dv.
V (i)
Le vecteur {Fi } est appelé vecteur force local pour l’élément (i) et vaut
y {
t t
{Fi } = − [N ] {Fv } dv + [N ] {fs } ds.
V (i) Sf (i)

Remarque 4.1. 1. La matrice [Ki ] est toujours symétrique en elasticité.

2. Le vecteur {Fi } est constitué de deux termes, le premier terme représente la contribution
des forces de volume au vecteur force et le second terme représente la contribution des
forces de surface au vecteur force.

4.2 Matrice Jacobienne [J]

Les relations de la transformation géométrique en 3D est égale à
  
 x1 

 E 
.
 D 
Nn · ..

x (s, t, ν) = N1 ···




 
  
xn









  
 y1 


 E 
.
 D 
y (s, t, ν) = N1 ··· Nn · ..

  
 


 yn





  
 z1 



 E 
..
 D 
z (s, t, ν) = N1 ··· Nn ·




  .
   
 zn

∂x
On en déduit que l’on peut exprimer par
∂s
 
x 
∂x
∂N1 ∂Nn
  .1 
= ∂s · · · ∂s
· .. .
∂s   
xn

Cela amène à une nouvelle définition de la matrice Jacobienne qui est particulière à l’élément
fini : ∂  ∂x ∂y ∂z 
 ∂s 
  ∂s ∂s ∂s

 
  

  ∂x ∂y ∂z 

∂
[J] = ∂t · x y z =   ∂t ∂t ∂t 


 
  

 
∂  ∂x ∂y ∂z
∂ν ∂ν ∂ν ∂ν

20
4.2. MATRICE JACOBIENNE [J] 21

 ∂N1 ∂Nn

∂s
··· ∂s  
  x1 y1 z1
∂Nn   .. .. .. 
 ∂N 
 ∂t · · ·
[J] =  1
·
∂t   . . .
  x y z
n n n
∂N1 ∂Nn
∂ν
··· ∂ν
La relation inverse permet alors de calculer les dérivées premières par rapport aux coordonnées
réelles des fonctions d’interpolation.
∂ ∂
 ∂x
  
  ∂s 
 

  
   
  
∂ −1 ∂
∂y = [J ] · ∂t

 
 
 


   
∂  
∂ 
∂z ∂ν

Connaissant les dérivées premières

∂Ni ∂Ni ∂Ni
, et
∂s ∂t ∂ν
ces relations participent à la construction de la matrice [B(s, t, ν)] à partir de la matrice
[B(x, y, z)], qui s’exprime en fonction des dérivées premières
∂NiR ∂NiR ∂NiR
, et .
∂x ∂y ∂z
Pour le calcul des dérivées secondes
∂2 ∂2
, , ···
∂x2 ∂y 2
des fonctions d’interpolation (problème de flexion), la démarche est identique, mais les calculs
sont plus complexes.

Calcul numérique des intégrales

y y
f (x, y, z) dx dy dz = f (s, t, ν) . |det [J]| ds dt dν
V (i) Dr(i)

Exercice 4.1. On considère une barre définie le long ed l’axe x, entre x1 et x2 .

Les paramètres caractérisant le matériau et la géométrie de la barre sont les suivants :
■ ρ est la masse volumique du matériau ;

■ E est le module de Young du matériau ;

■ A est la section de la barre ;

■ L est la longueur de la barre.

On considère deux types de chargement à la barre (F1 et F2 ) et ensuite une force par unité de
longueur q(x) répartie le long de la barre. L’équation d’équilibre est donnée par :

d
AE dU



 dx dx
+ q(x) = 0





dU
dx
(x1 ) = −FAE
1
= εx





 dU (x ) = F2 = ε

dx 2 AE x

21
22 CHAPITRE 4. MATRICES DE RIGIDITÉS LOCALES

1. Déterminer les formulations forte et faible par hypothèse de Galerkine pour le problème.
2. Déterminer la matrice de rigidité locale et le vecteur force local de l’élément barre.

22
Chapitre 5

Expansion – Assemblage – Résolution

5.1 Représentation du maillage

Définition 5.1
■ Le nombre d’élément du maillage est noté N .
■ Le nombre de nœud du maillage est noté Nd .
■ Le nombre de nœud par élément est noté n.

5.1.1 Table de nœud

Indice du nœud Coordonnées x Coordonnées y Coordonnées z
1 x1 y1 z1
2 x2 y2 z2
.. .. .. ..
. . . .
Nd xNd yNd zNd

La table des nœuds contient tout simplement les coordonnées cartésiennes des nœuds sous la
forme ci-dessus. Le numéro de nœud de cette table (allant de 1 à Nd ) est un numéro global
du nœud au sein du maillage, qu’il ne faut pas confondre avec l’indice d’un nœud au sein d’un
élément.

5.1.2 Tabel de connectivité

Cette table représente essentiellement la manière dont les nœuds de la liste précédente sont
reliés les uns aux autres pour former des éléments.
On considère ici, pour simplifier, que tous les éléments sont du même type et donc qu’ils ont
tous le même nombre de nœuds.
Cette table se présente sous la forme d’un tableau à deux entrées CONNEC (i, j) ; où i (indice
de ligne du tableau) correspond à l’indice de l’élément considéré (et donc i varie entre 1 et N ),
et j (indice de colonne) correspond à l’indice du nœud au sein de l’élément (et donc varie entre
1 et n).

23
24 CHAPITRE 5. EXPANSION – ASSEMBLAGE – RÉSOLUTION

Ensuite, le contenu du tableau est la liste des numéros des nœuds qui forment l’indice de
l’élément j.

5.2 Matrice de rigidité et vecteurs forces locaux

5.2.1 Principe

  
On a : W (i) = δUn(i) [−K(i)] Un(i) + {Fi } .
De manière générale, la taille des matrices de rigidité locale est

(n × N DDL) × (n × N DDL)

où
■ n est le nombre de nœud par élément ;
■ N DDL est le nombre de degré de liberté par nœud.

Définition 5.2 (Les degrés de liberté au nœud )

Pour les éléments volumiques classiques (tétraèdres et hexaèdres, par exemple), on a
N DDL = 3. Ce qui correspond aux composantes de déplacement dans les trois directions
de l’espace.
Certains éléments peuvent également avoir (c’est le cas des éléments finis poutres et des
coques, par exemple) des degrés de liberté connu en rotation.

Règle 5.1.
1. La taille des matrices de rigidité locales est : (n × N DDL) × (n × N DDL).
2. La taille des vecteurs forces locaux est : n × N DDL.
3. La taille des matrices de rigidité globales et des matrices de rigidité locales expansées est :
(Nd × N DDL) × (Nd × N DDL).
4. La taille des vecteurs forces globales et des vecteurs forces locaux expansés est : Nd ×
N DDL.

5.2.2 Vecteur de localisation


On définit pour chaque élément du maillage, un vecteur de localisation noté : LOC(i) pour
l’élément (i).
Cet outil permet facilement d’envoyer les termes des matrices locales et vecteurs locaux au bon
endroit dans les matrices locales expansées et vecteurs locaux expansés. La taille de ces vecteurs
de localisation est : n × N DDL.
La Ke composante du vecteur de localisation de l’élément (i) est définit grâce à la formule

LOC(i) (k) = (numéro nœud global − 1) N DDL + indice ddl,

où indice ddl est l’indice du degré du nœud considéré, il varie donc entre 1 et n × N DDl.

24
5.2. MATRICE DE RIGIDITÉ ET VECTEURS FORCES LOCAUX 25

On forme les valeurs des composantes du vecteur de localisation de l’élément (i) en bouclant
sur les nœuds de l’élément (i), et sur les degrés de liberté de chaque nœud selon l’algorithme.
Dans cet algorithme, k correspond à l’indice du composant du vecteur de localisation et sera
incrémenté de 1 à n × N DDL.

Algorithme 5.1 (Formation du vecteur de localisation pour l’élément (i))

k=1
Pour j allant de 1 à n
Pour l allant de 1 à N DDL
LOC(i) (k) = (CONNEC(i, j) − 1) × N DDL + l
k =k+1
Fin Pour
Fin Pour

5.2.3 Utilisation des vecteurs de localisation

Comme mentionné ci-dessous, on utilise les vecteurs de localisation pour former les matrices et
les vecteurs locaux expansés.

Définition 5.3
On utilise, dans la suite du texte, les notations suivantes :
■ [Ki ] est la matrice de rigidité locale de l’élément (i) ;
■ {Fi } est le vecteur force local pour l’élément (i) ;
■ [Kie ] est la matrice de rigidité locale expansée pour l’élément (i) ;
■ {Fie } est le vecteur force local expansé pour l’élément (i).

On commence par initialiser toutes les composantes des matrices locales expansées et des
vecteurs locaux expansés à 0.

Algorithme 5.2 (Formation du vecteur local expansé pour l’élément (i))

Pour k allant de 1 à n ×

N DDL
e
Fi LOC(i) (k) = F(i) (k)
Fin Pour

Exemple 5.1.

25
26 CHAPITRE 5. EXPANSION – ASSEMBLAGE – RÉSOLUTION

Algorithme 5.3 (Formation de la matrice local expansé pour l’élément (i))

Pour k allant de 1 à n × N DDL

Pour l allant de 1 à n × N DDL

Kie LOC(i) (k) , LOC(i) (l) = K(i) (k, l)
Fin Pour
Fin Pour

 
Le terme à l’intersection de la ligne k et la colonne
 e  l de K (i) est envoyé à l’intersection de la
ligne LOC(i) (k) et de la colonne LOC(i) (l) de K(i) .

Exemple 5.2.

5.2.4 Utilisation des matrices et vecteurs expansés

De cette manière, la contribution de chaque (i) à la forme intégrale faible peut maintenant
s’écrire   e   e 
W(i) = ⟨δUNd ⟩ − K(i) {UNd } + F(i) ,
où le vecteur {UNd } est tel que


⟨UNd ⟩ = u1 v1 w1 · · · un vn wn .

Ce vecteur regroupe tout simplement toutes les composantes des déplacements au nœud du
maillage. On a globalement


⟨δUNd ⟩ = δu1 δv1 δw1 · · · δun δvn δwn .

Matrice de rigidité globale et vecteur force global

On a :
N
X N
X  e 
{UNd } + {Fie } = 0


W = W(i) = 0 =⇒ ⟨UNd ⟩ − K(i)
i=1 i=1
" N N
#
X  e  X
⇐⇒ ⟨δUNd ⟩ − K(i) {UNd } + {Fie } = 0.
i=1 i=1

Et finalement, on introduit [K] et {F }, ce qui donne

N
X
W(i) = ⟨δUNd ⟩ [−[K] {UNd } + {F }] .
i=1

 e globale [K] et le vecteur force global {F } àe partir

Règle 5.2. On calcule la matrice de rigidité
des matrices de rigidité locales expanées K(i) et les vecteurs forces locaux expansés {Fi } de la
manière suivante :
X N XN
 e   e
[K] = K(i) et {F } = F(i) .
i=1 i=1

26
5.2. MATRICE DE RIGIDITÉ ET VECTEURS FORCES LOCAUX 27

5.2.5 Particularités du système finale

Ici, on a un système linéaire particulier puisque :

−[K] {UNd } + {F } = 0 ⇐⇒ [K] {UNd } = {F } .

■ La matrice [K] est complétement connue.

■ Le vecteur {UNd } est majoritairement inconnu, car on impose des déplacements et donc
certaines composantes de {UNd } sont connues.
■ Le vecteur {F } est majoritairement connu, car les seuls composantes de {F } qui sont
inconnues sont les composantes correspondant au déplacement imposé.
En général, les déplacements imposés correspondent aux appuis.

27
Chapitre 6

Application à l’élasticité linéaire

6.1 Application de la forme générale à un élément barre

On considère l’exercice 1 sur l’élément barre.
Nous avons vu que la forme générale de la matrice de rigidité locale et le vecteur force local est
Z Z I
t t t
  
K(i) = [B][H][B] dv et F(i) = − [N ] {Fv } dv + [N ] {fs } ds.
V (i) V (i) Sf (i)

Nous allons maintenant appliqué la forme générale de ces équations à l’élément barre.

6.1.1 Calcul de [H]

On a


⟨ε⟩ = εx εy εz 0 0 0 et ⟨σ⟩ = σx 0 0 0 0 0
donc
⟨δε⟩ {σ} = δεx · σx .
Ainsi, pour simplifier, on considère ici que les vecteurs lignes ⟨ε⟩ et ⟨σ⟩ se réduisent à des
scalaires :
⟨ε⟩ = εx et ⟨σ⟩ = σx .
Ainsi [H] = E.

6.1.2 Calcul de la matrice [B]

On rappelle que la matrice [B] est telle que


{ε} = [B] Un(i) ,




avec le vecteur UN (i) qui est ici égal à Un(i) = U1 U2 .
On introduit ici la matrice J −1 telle que
 

 
dU
{ε} = J −1 {εs }
 
avec {εs } = .
ds

28
6.1. APPLICATION DE LA FORME GÉNÉRALE À UN ÉLÉMENT BARRE 29

L 2
donc J −1 = .
 
Ici [J] =
2 L
On introduit ensuite la matrice [Bs ] telle que :
   
 du 
 u1
⟨εs ⟩ = [Bs ] · Un(i) et {εs } = ; avec Un(i) = N1 (s) N2 (s)
ds u2
 

dN dN  u1
=⇒ {εs } = ds1 ds2
u2

donc
dN dN2

[Bs ] = ds
1
ds
.
Comme
1−s 1+s
[Bs ] = − 12 1


N1 (s) = et N2 (s) = , il vient 2
.
2 2
Or
 
u
{ε} = [B] 1 = J −1 {εs }
 
u2
 
u
[Bs ] 1 .
 −1 
= J
u2

Finalement,
[B] = J −1 [Bs ]
 

et dans le cadre de l’application

2
1 1
 1

[B] = −2 2
= −1 1 .
L L

⟨N (x)⟩ = N1 2x−xL2 −x1 N2 2x−xL2 −x1


L+2x−x2 −x1
= L−2x+x 2 +x1


2L 2L

Ici, la force repartie q(x) est considérée comme une force de volume qui est égale à

q(x)
{Fv } =
A
I
t
Pour ce qui est du deuxième terme : [N ] {fs } ds.
Sf (i)

Il n’y a pas d’effort de pression de surface appliquée, et donc {fs } = 0. Cela conduit donc à
notre force {Fi } qui est égale à
Z I
t t
{FI } = [N ] {Fv } dv + [N ] {fs } ds + {· · · }
V (i) Sf (i)
Z  
t q(x) F1
= [N ] dv +
V (i) A F2
1 
Z (L − 2x + x 2 + x 1 )
1 2L
=   q(x) dv .
A V (i) 1
2L
(L + 2x − x2 − x1 )

29
30 CHAPITRE 6. APPLICATION À L’ÉLASTICITÉ LINÉAIRE

6.2 Élément poutre

6.2.1 Élément poutre en flexion uniquement

Introduction

On considère maintenant une poutre définie le long de l’axe x, entre x1 et x2 , ne travaillant

qu’en flexion ; c’est-à-dire que l’on n’a pas de chargement sur cette poutre, qu iinduit un effort
normal, mais uniquement une effort normal et uniquement un effort tranchant et un moment
fléchissant.
Les paramètres caractérisant le matériau et la géométrie de la barre sont les suivants :
■ I est le moment d’inertie de la section de la poutre ;
■ E est le module de Young du matériau ;
■ A est la section de la poutre ;
■ L est la longueur de la poutre.
On considère deux (02) types de chargement appliqués à la poutre tel que illustrés. D’abord les
forces ponctuelles F1 et F2 , verticales (positives vers le haut) appliquées aux deux extrémités
et ensuite les moments ponctuels M1 et M2 appliquées aux deux extrémités (positifs dans le
sens trgonométrique anti-horaire).

Équations fondamentales

La solution recherchée est N (x) (pour x ∈ [x1 , x2 ]), le champ de déplacement vertical (positif
vers le haut), le long de l’axe x. On considère également dans la résolution, la dérivée de v(x)
qui correspond physiquement, lorsque la poutre fléchit, à la rotation de la section de la poutre
au point de coordonnées x. Cette rotation, anti-horaire, est classiquement notée

dv
θ(x) = .
dx
Ce phénomène physique, avec les hypothèses mentionnées ci-dessus, est décrit par l’équation
différentielle
d2 d2 v
 
EI 2 = 0,
dx2 dx
et, pour E et I constantes, on a
d4 v
EI= 0.
dx4
Les conditions aux limites, imposées par l’application des forces F1 et F2 , ainsi que les moments
M1 et M2 , se traduisent de la manière suivante :
on rappelle que l’expression du moment fléchissant le long de l’axe x est

d2 v
M (x) = EI
dx2
et que l’expression de l’effort tranchant le long de l’axe x est

dM d3 v
V (x) = − , i .e. V (x) = −EI .
dx dx3
30
6.2. ÉLÉMENT POUTRE 31

Sur la surface, on a
V = −F1 et M = −M1 ,
c’est-à-dire
dM d3 v
(x1 ) = F1 ⇐⇒ EI 3 (x1 ) = F1
dx dx
d2 x d2 v M1
=⇒ M (x1 ) = EI 2 (x1 ) = −M1 =⇒ 2
(x1 ) = − .
dx dx EI
De même, sur la face droite, on a :

d3 v F2 d2 v M2
3
(x2 ) = − et 2
(x2 ) = .
dx EI dx EI

La forme intégrale forte est donnée par

Z x2  4 
dv
EI ψ (x) dx = 0, ∀ψ : fonction de pondération.
x1 dx4

La forme intégrale faible

x2 x2 x2 !
d4 v dψ d3 v
Z   Z  3
d v
W = EI ψ (x) dx = EI − dx + ψ 3 (x)
x1 dx4 x1 dx dx3 dx x1

x2 x2 !
dψ d3 v
Z  3
d v
= −EI dx − ψ 3 (x)
x1 dx dx3 dx x1

x2 x2  3 x2 !
d2 ψ d2 v dψ d2 v
Z 
d v
= −EI − dx + (x) − ψ 3 (x)
x1 dx2 dx2 dx dx2 x1 dx x1

x2 x2  3 x2
d2 ψ d2 v dψ d2 v
Z 
d v
= EI dx − (x) + ψ 3 (x)
x1 dx2 dx2 dx dx2 x1 dx x1

x2
d2 ψ d2 v d2 v d2 v
Z  
dψ dψ
= EI dx − EI
(x2 ) 2 (x2 ) − (x1 ) 2 (x1 )
x dx2 dx2 dx dx dx dx
1 3 3

d v d v
+ EI ψ (x2 ) 3 (x2 ) − ψ (x1 ) 3 (x1 )
dx dx

x2
d2 ψ d2 v
Z  
dψ dψ
= EI dx − M 2 (x 2 ) + M 1 (x 1 ) + (−F2 ψ (x2 ) − F1 ψ (x1 )) .
x1 dx2 dx2 dx dx

Ainsi donc, la forme intégrale faible est

Z x2 2 2  
d ψd v dψ dψ
W = EI 2 2
dx − M2 (x2 ) + M1 (x1 ) − (F2 ψ (x2 ) + F1 ψ (x1 )) = 0.
x1 dx dx dx dx

Forme intégrale faible avec hypothèse de Galerkine

 
∂v ∂
ψ = δv, δ = (δv)
∂x ∂x

31
32 CHAPITRE 6. APPLICATION À L’ÉLASTICITÉ LINÉAIRE

x2
d2 d2 v
Z    
dv dv
W = EI (δv) dx − M2 δ (x2 ) − M1 δ (x1 ) − F2 δv (x2 ) − F1 δv (x1 ) .
x1 dx2 dx2 dx dx

dv dv
Posons (x2 ) = θ (x2 ) et (x1 ) = θ (x1 ).
dx dx
La forme intégrale faible avec hypothèse de Galerkine est donc
x2
d2 (δu)
Z    

 F1
 M1
W = EI dx − δv (x1 ) δv (x2 ) − δθ (x1 ) δθ (x2 ) =0
x1 dx2 F2 M2

Maillage et passage de l’élément de reférence

La matrice Jacobienne dans ce cas est identique à la matrice Jacobienne de l’élément barre

1 L L  −1  2
[J] = (x2 − x1 ) = ; det[J] = et J = .
2 2 2 L
Ce qui donne pour la transformation des dérivées

dv  −1  dv ∂  ∂
= J −1

= J sachant que
dx  ds ∂x ∂s
2
 
dv d dv d dv
= J −1 [J]−1
 
=⇒ 2
=
dx dx dx ds ds
2 2
dv 4 dv
donc 2
= 2 2.
dx L ds
L
En appliquant la transformation de l’intégrale avec dx =
ds, dx = det[J] ds, il vient
2
Z 1  2  2     
8 dv dv
 M1
 F1
W(i) = 3 EI δ − δθ(−1) δθ(1) − δv(−1) δv(1)
L −1 ds2 ds2 M2 F2

avec
d2 (δv) d2 v 4 d2 v 4 d2 v L
 
dx = δ × × ds.
dx2 dx2 L2 ds2 L2 L2 2

Introduction à l’interpolation nodale

Ensuite, la grande particularité c’est que cet élément n’est pas isoparamétrique car on introduit
une interprétation qui n’est pas de degré 1.
On introduit donc une interpolation particulière v(s) qui est égale à

v(s) = a0 + a1 s + a2 s2 + a3 s3 .

On introduit ensuite dans l’interpolation les grandeurs nodales (aux nœuds) v1 , θ1 , v2 , θ2 .

Pour que l’interpolation soit cohérente, on doit avoir

dv dv
v(−1) = v1 , v(1) = v2 , (x1 ) = θ (x1 ) = θ1 , et (x2 ) = θ (x2 ) = θ2 .
dx dx
32
6.2. ÉLÉMENT POUTRE 33

dv 2 dv
En remarquant que (x1 ) = , cela donne
dx L ds
dv 2 dv dv L
(x1 ) = (−1) = θ1 =⇒ (−1) = θ1
dx L ds ds 2
dv 2 dv dv L
(x2 ) = (1) = θ2 =⇒ (1) = θ2 .
dx L ds ds 2
Ces conditions sont introduites dans l’interpolation, et donnent alors


 v(−1) = v1 = a0 − a1 + a2 − a3 
1 −1 1 −1
  
a0 v1


 v(1) = v2 = a0 + a1 + a2 + a3 0 1 −2 3  a1   L θ1 
⇐⇒ 
1 1
  = 2 
dv
(−1) = L
θ = a1 − 2a2 + 3a3 1 1  a2   v2 
2 1

 dsdv
 L
0 1 2 3 a3 θ

L 2 2
ds
(1) = θ
2 2
= a1 + 2a2 + 3a3

notant  
1 −1 1 −1
0 1 −2 3 
D=
1 1

1 1
0 1 2 3
il s’ensuit      
a0 v1 
 a0 
 a1  L   
  = D−1  2 θ1  , a1

 
 a2   v2  et V (s) = 1 s s2 s3

 a2 
L  
a3 θ a3

2 2

 

 v1  

2 3
 −1  L θ1 
i .e. V (s) = 1 s s s D 2
 v2 

L 
θ

2 2
 
1 L 3 L
= (v1 + v2 ) − (θ2 − θ1 ) + (v2 − v1 ) − (θ1 + θ2 ) s
2 8 4 8
 
L 1 L
+ (θ2 − θ1 ) s + − (v2 − v1 ) + (θ1 + θ2 ) s3
2
8 L 8

On a finalement cette interpolation

 
v1 
 

 θ1 
V (s) = N1 (s) N2 (s) N3 (s) N4 (s)

 v2 
 
θ2


avec 
1


 N1 (s) = 2
− 34 s + 14 s3





L
 N2 (s) = − L8 s − L8 s2 + L8 s3


8



N3 (s) = 21 + 34 s − 14 s3









 N (s) = − L − L s + L s2 + L s3

4 8 8 8 8

33
34 CHAPITRE 6. APPLICATION À L’ÉLASTICITÉ LINÉAIRE

On peut également remarquer que l’on a

 
v 1
dv
dN1 dN2 dN3 dN4
  
= ds ds ds ds
θ v
1 2
ds
θ2
 

Ce qui donne ensuite

 
 v1 
d2 v
d2 N1

 
θ1

d2 N2 d2 N3 d2 N4

= ds2 ds2 ds2 ds2
ds2 
 v2 
 
θ2


avec
 
dN1 d2 N1


 ds
(s) = − 34 + 43 s2 

 ds2
(s) = 3
2
s

 


 

d2 N2
 
dN2
(s) = − L8 − L4 s − 38 Ls2 (s) = − L4 + 34 Ls

 

ds ds2

 

et
d2 N3
 dN3 3 3 2

(s) = − s (s) = − 32 s

 



 ds 4 4 

 ds2

 


 

dN4 d2 N4
= − L8 + L4 s + 38 Ls2 L
+ 34 Ls
 

ds
(s) 
ds2
(s) = 4

 2 
d N1 
ds2 



 


 

 d 2N 
 2
d2 v
   

  ds2 
δ = δv1 δθ1 δv2 δθ2
ds2 
 d2 N3 

ds2 

 

 


 

 d2 N4 
 
ds2

De même
 
F 1 
 

 0 
δV1 .F1 = δv1 δθ1 δv2 δθ2

 0
 
0

 
0
 

 0 
δv2 .F2 = δv1 δθ1 δv2 δθ2

 F2 
 
0

 
 0 
 

 M1 
δθ1 M1 = δv1 δθ1 δv2 δθ2

 0 
 
0

 
 0 
 

 0 
δθ2 M2 = δv1 δθ1 δv2 δθ2

 0 
 
M2


34
6.2. ÉLÉMENT POUTRE 35

 2 
d N1
ds2 

 

 

   
v1 

 d2 N 

 2
Z 1  
  ds2  D E
 
8EI
d2 N1 d2 N2 d2 N3 d2 N4 θ1

=⇒ W(i) = δv1 δθ1 δv2 δθ2 ds2 ds2 ds2 ds2
ds
L3 −1 
 d2 N3 
 
v2 
 
ds2  θ2

  

 


 

 2 

d N4 
ds 2
 
 F1 
 

 M1 
− δv1 δθ1 δv2 δθ2
 F2 
 
M2
 
  2  
d N1
 ds2 

 

 
     
 2 
v1   F1 
  
d N
 
 
 2 

2 D      
 8EI ds θ1 M1 
 E 
d2 N1 d2 N2 d2 N3 d2 N4

= δv1 δθ1 δv2 δθ2  3 ds2 ds2 ds2 ds2

−
 L  v F2 

2 2
 ddsN23 
  
  
     
θ2 M2 
  
 
 
 

 
 
 
 d2 N4 
 
ds2

Ainsi

3
   
2s  F1 

 8EI
Z 1
− L + 3 Ls M1 

3
− L4 + 34 Ls − 23 s L
+ 34 Ls ds


K(i) = 3  4 34  et F(i) =
L  − s  2s 4 F2 
−1 2 

 
L 3
4 + 4 Ls
M2


Finalement
 
12 6L −12 6L
  EI  6L 4L2 −6L 2L2 
K(i) = 3  
L −12 −6L 12 −6L
6L 2L2 −6L 4L2

Ajout d’une force repartie

On ajoute maintenant au modèle précédent une force par unité de longueur q(x), repartie le
long de la poutre (ce chargement q(x) est homogène au niveau des unités en des N/m).
Ce chargement reparti q(x) permet par exemple de modeliser le poids propre de la poutre dans
le champ de gravité. De manière plus générale, q(x) peut avoir une variation quelconque.
Introduisant une force repartie, le modèle devient

d4 v d4 v
EI = q(x) ⇐⇒ EI − q(x) = 0.
dx4 dx4
35
36 CHAPITRE 6. APPLICATION À L’ÉLASTICITÉ LINÉAIRE

Notons {G} l’apport de force de volume.

Z

G(i) = ⟨δU ⟩ {Fv } dv ; ici Fv = q
V (i)
 
Z 
 N1 (x)


 N2 (x)
= δv1 δθ1 δv2 δθ2 q(x) dx
V (i) 
 N3 (x)

N4 (x)
 
   
Z 
 N1 (x)
 
Z 1
 N1 (x)


N2 (x)
 L N2 (x)

q(x) dx = q(x(s)) ds
N3 (x)
V (i)   2 −1 
 N3 (x)

N4 (x) N4 (x)
   
 
N1 (x)
L 1 N2 (x)
Z  
= q(x(s)) ds
2 −1 
N3 (x)

N4 (x)
 

Si q est constante, alors  


Z 1 N1 (x)

 L 
N2 (x)

G(i) = q ds .
2 −1  N3 (x)

N4 (x)
 

Finalement  

 F1 
 
M1
 
{Fi } = + G(i) .

 F2 
 
M2


Exemple 6.1. On considère un cas pratique (figure ci-dessous) de poutre encastré, de longueur
1 m soumise à une charge repartie constante q = 120 N/m, ainsi qu’à un moment C = 50 N.M
à son extrémité.
On considère également que EI = 1 000 N.m2 .

Résultat important

Si on considère un élément poutre ne travaillant qu’en flexion, de longueur L, le moment d’inertie

I de section A et de module de Young E, chargée par deux forces ponctuelles F1 et F2 et deux
moments ponctuels M1 et M2 .
La matrice de rigidité locale est donnée par
 
12 6L −12 6L
  EI  6L 4L2 −6L 2L2 
K(i) = 3  
L −12 −6L 12 −6L
6L 2L2 −L 4L2
et le vecteur force local est
   1 3 1 3

 F1  2
− 4
s + 4
s
 M1  L 1  L − L s − L s2 + L s3 
  Z
F(i) = +  8 18 3 8 1 3 8 
 F2 
  2 −1  2
+ 4s − 4s 
L L L 2 L 3
M2 −8 − 8s + 8s + 8s
 

36
6.3. ÉLÉMENT DE POUTRE EN TENSION ET FLEXION COMBINÉES 37

Pour q constante
   
 F
  1  1
 M1  qL  L 
F(i) = +  6 

 F 2  2  1 
 
M2 − L6


6.3 Élément de poutre en tension et flexion combinées

6.3.1 Introduction

On considère maintenant un epoutre définie le long de l’axe x, entre x1 et x2 , travaillant à la

fois en tension et en flexion, c’est-à-dire que l’on a sur cette poutre un chargement qui induit à
la fois un effort normal, un effort tranchant et un moment fléchissant.
Les paramètres caractérisant le matériau et la géométrie de la barre sont les suivants : I, A, E
et L.
On considère deux types de chargement appliquées à la poutre tels que illsutré sur la figure.
D’abord, les forces ponctuelles Fx1 et Fx2 appliquées aux deux extrémités, induisant un effort
normal. Ensuite, on a les forces ponctuelles Fy1 et Fy2 appliquées également aux extrémités,
induisant un effort tranchant et un moment fléchissant. On considère également les moments
M1 et M2 appliqués aux deux extrémités.
On a trois composantes d’effort par nœud (Fx1 , Fy1 et M1 , par exemple pour le nœud 1). Cela
entraîne que l’on a trois degrés de liberté
■ u1 , v1 , θ1 pour le nœud 1 ;
■ u2 , v2 , θ2 pour le nœud 2.

6.3.2 Matrice de rigidité et vecteur force

Pour un élément de barre, en négligeant les forces reparties, on a vu que l’on a la réalisation
suivante
    
EA 1 −1 u1 Fx1
=
L −1 1 u2 Fx2

Pour un élément de poutre en flexion uniquement, toujours en négligeant les forces reparties,
on a vu que l’on a la relation
    
12 6L −12 6L  v1  Fy 
EI  2 2  
   1

6L 4L −6L 2L θ M1

1
  = .
L3 −12 −6L 12 −6L   v2 
 
 Fy2 

6L 2L2 −8L 4L2 θ2 M2
   

On en déduit facilement la relation globale reliant les efforts au degré de liberté pour l’élément

37
38 CHAPITRE 6. APPLICATION À L’ÉLASTICITÉ LINÉAIRE

de poutre superposant ces deux comportements

EA
− EA
    
L
0 0 L
0 0 u1 Fx1 
 
    




12EI 6EI
− 12EI 6EI  
     
 0 0 v

F

L2  1  y
 
 L3 L2 L3 
 1 

     

    
 

6EI 4EI
 0 0 − 6EI 2EI  
θ M
 
1 1
 
 L2 L L2 L    
   =
 EA EA
   
−
u2  
0 0 0 0    Fx2 
 L L   
    




     
 0 − 12EI − 6EI 0 12EI 6EI  
− L2  v2   F
 
y2 
 
 L3 L2 L3 
 

     

6EI 2EI  

0 0 − 6EI 4EI
θ M

L2 L L2 L 2 2
| {z }| {z } | {z }
A U F

Résultat important

Si on considère un élément de poutre travaillant en flexion et en tension, de longueur L, du

moment d’inertie de section I, de section A et de matériau de module de Young E, chargé par
des forces ponctuelles Fx1 , Fy1 , Fx2 et Fy2 , et deux moments ponctuels M1 et M2 .
Le comportement de ce type d’élément est représenté par l’égalité AU = F .

6.4 Extension des formations en dimension 2 et en dimension

3

Dans les paragraphes précédents, nous avons vu les modèles de barre et de poutre confinés le
long de l’axe x, entre les coordonnées x1 et x2 .
Dans une structure réelle, les éléments sont repartis à travers le plan (pour un treilli, pour une
structure en deux dimensions) ou à travers l’espace (pour un treilli, pour une structre en trois
dimensions).
Prenons l’exercice dans le plan, pour un élément de barre tel que illustré.
On définit deux repères de coordonnées, le repère (x, y) et le repère global qui représente le
plan dans lequel la structure est définie.
Le repère (x, y) quant à lui est un repère qui est local à l’élément barre.

Passage en deimension deux

Soit u1 , v1 , u2 et v2 les composantes de déplacement aux nœuds 1 et 2 exprimés dans le repère

global. Soit Fx1 , Fy1 , Fx2 et Fy2 les composantes des forces aux nœuds 1 et 2 exprimés dans le
repère global.

38
6.4. EXTENSION DES FORMATIONS EN DIMENSION 2 ET EN DIMENSION 3 39

Posons
   
1 0 −1 0 cos ϕ sin ϕ 0 0
 ′  EA  0 0 0 0 − sin ϕ cos ϕ 0 0 
K(i) =   et [P ] =  
L −1 0 1 0  0 0 cos ϕ sin ϕ 
0 0 0 0 0 0 − sin ϕ cos ϕ
 ′ 
et soit Kixy la matrice de rigidité locale de l’élément de barre dans le repère global (X, Y ).
On montre que
K(i)xy = [P ]−1 K(i)
 ′   ′ 
[P ].
Finalement
   
 ′  [A] −[A] C 2 ϕ CϕSϕ
K(i)xy = ; avec [A] = , où Cϕ = cos ϕ et Sϕ = sin ϕ
−[A] [A] CϕSϕ S 2 ϕ

Passage en trois dimensions

Si on note (X1 , Y1 , Z1 ) et (X2 , Y2 , Z2 ) les coordonnées des nœuds 1 et 2 dans le repère global
(X, Y, Z), on a :
L2 = (X2 − X1 )2 + (Y2 − Y1 )2 + (Z2 − Z1 )2
et la matrice de rigidité local de l’élément de barre dans le repère global est donnée par
 2
αγ −α2 −αβ −αγ

α αβ
 
 2 2

 αβ
 β βγ −αβ −β −βγ  
 
 2 2

 ′  
 αγ βγ γ −αγ −βγ −γ 

K(i)XY Z =   

 −α2 −αβ −αγ α2 αβ αγ 
 
 
 
−αβ −β 2 −βγ αβ β 2
βγ 
 
 
2 2
−αγ −βγ −γ αγ βγ γ

où
X 2 − X1 Y2 − Y1 Z2 − Z1
α= , β= et γ = .
L L L

39
Bibliographie

[1] Jean-Christophe Cuillère

Introduction à la méthode des éléments finis

40