ANALYSE AVANCÉE DES STRUCTURES

ET ÉLÉMENTS FINIS
– CIV6502 –

Prof. Najib Bouaanani, ing., Ph.D.

Département des génies civil, géologique et des mines
École Polytechnique de Montréal

Tous droits réservés

 Najib Bouaanani, École Polytechnique de Montréal, 2004–2011

Cet ouvrage est exclusivement réservé à l’usage du cours

Analyse avancée des structures par éléments finis (CIV6502)

On ne peut reproduire ni diffuser aucune partie

du présent ouvrage, sous quelque forme ou par quelque procédé
que ce soit, sans avoir obtenu au préalable l’autorisation écrite de l’auteur.
Chapitre 6

Problèmes bidimensionnels
Élasticité plane

6.1 Introduction
L’une des premières applications concrètes de la MEF a été la résolution des problèmes d’élasticité en mé-
canique des solides. Dans de tels problèmes, les inconnus sont les déplacements, les déformations et les
contraintes au sein d’un domaine d’étude ayant un comportement élastique et soumis sous à des charges
et/ou des déplacements imposés. En génie civil, un tel domaine d’étude est généralement représenté par une
structure ou un sol.

6.2 Modélisation par éléments finis

Tel que déjà mentionné, résoudre un problème d’élasticité consiste à déterminer les déplacements, les dé-
formations et les contraintes au sein d’un domaine d’étude ayant un comportement élastique, sous l’effet
de charges et/ou de déplacements imposés. Les équations présentées au chapitre 2 permettent de relier les
déplacements, les déformations et les contraintes au sein d’un domaine d’étude vérifiant un état d’élasticité
plane. Ces équations pourraient être utilisées pour résoudre les problèmes d’élasticité planes, mais rares sont
les cas où l’on aboutirait à une solution analytique exacte, notamment à cause de la géométrie généralement
complexe des domaines d’étude.

Dans ce qui suit, nous nous proposons d’utiliser la MEF pour résoudre les problèmes d’élasticité planes.
Les étapes de modélisation sont les mêmes que celles décrites au chapitre 4. Le domaine d’étude sera mo-
délisé par un assemblage d’éléments finis plans, interconnectés par des noeuds. La forme géométrique des
éléments peut être triangulaire ou quadrangulaire. La figure 6.1 illustre un maillage à base d’éléments finis
triangulaires d’un système barrage-fondation. Seule la formulation détaillée des éléments triangulaires sera
présentée, celle-ci étant facilement généralisable à d’autres types d’éléments. La variable d’état primaire
du problème est le déplacement en tout point du domaine d’étude. On définit alors deux DDL par noeuds,
correspondant aux déplacements u et v selon un système d’axes Cartésien x et y. Suite à l’assemblage des
matrices et vecteurs élémentaires, le principe de l’énergie potentielle minimum est utilisé pour construire

1
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-2

un système d’équations linéaires traduisant l’équilibre du domaine étudié sous chargement. La résolution
numérique de ce système fournie les valeurs des déplacements aux noeuds du maillage. Connaissant ces
déplacements, les relations déformations-déplacements et contraintes-déformations, ainsi que les équations
d’équilibre sont utilisées pour déterminer les autres inconnus du problème, à savoir les déformations, les
contraintes et les réactions.

Figure 6.1– Maillage à base d’éléments finis triangulaires d’un système barrage-fondation.

6.3 Éléments finis triangulaires

La figure 6.2 montre un élément finis triangulaire numéroté localement. Chaque noeud i possède deux DDL
matérialisés par les déplacements ui et vi . L’ensemble des DDL de l’élément e sont regroupés dans un
vecteur-colonne u(e)  
u1
 v1 
 
u 
 2
u(e) =   (6.1)
 v2 
 
 u3 
v3

6.3.1 Approximation des déplacements

Nous appliquerons la même démarche étudiée au chapitre 4 et chercherons à établir les déplacements à
l’intérieur de chaque élément en fonction des déplacements de ses noeuds. Faisons l’hypothèse que les dé-

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-3

Figure 6.2– Élément fini triangulaire numéroté

localement.

placements en tout point de l’élément varient linéairement en fonction des coordonnées

u(x, y) = α1 + α2 x + α3 y
(6.2)
v(x, y) = β1 + β2 x + β3 y

où les coefficients αi et βj sont appelés coordonnées généralisées.

L’équation (6.2) s’écrit sous la forme matricielle

 
α1
 
[ ] [ ] α2 
u(x, y) 1 x y 0 0 0 α3 

=   (6.3)
v(x, y) 0 0 0 1 x y  β1 
 
 β2 
β3

L’équation (6.3) s’applique également aux déplacements aux noeuds de l’élément

    
u1 1 x1 y1 0 0 0 α1
 v1  0 0 0 1 x1 y1  α2 
    
u  1 x y 0 0 0  
 2  2 2  α3 
u(e) = =   (6.4)
 v2  0 0 0 1 x2 y2   β1 
    
u3  1 x3 y3 0 0 0   β2 
v3 0 0 0 1 x3 y3 β3

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-4

Cette équation permet d’obtenir le vecteur des coordonnées généralisées αi et βj

   −1  
α1 1 x1 y1 0 0 0 u1
α2  0 0 0 1 x1 y1   v1 
 
   
α  1 x y 0 0 0  
 3  2 2  u2 
 =    (6.5)
 β1  0 0 0 1 x2 y2   v2 
     
 β2  1 x3 y3 0 0 0  u3 
β3 0 0 0 1 x3 y3 v3

d’où en remplaçant dans l’équation (6.3)

 −1  
1 x1 y1 0 0 0 u1
0 0 0 1 x1 y1   v1 
 
[ ] [ ] 
u(x, y) 1 x y 0 0 0 
1 x2 y2 0 0 0  
 u2 
=     (6.6)
v(x, y) 0 0 0 1 x y 0 0 0 1 x2 y2   v2 
   
1 x3 y3 0 0 0  u3 
0 0 0 1 x3 y 3 v3

On démontre que l’équation (6.6) se simplifie à

 
u1
 v1 
[ ] [ ] 
u 
u(x, y) N1 (x, y) 0 N2 (x, y) 0 N3 (x, y) 0  2
=   (6.7)
v(x, y) 0 N1 (x, y) 0 N2 (x, y) 0 N3 (x, y)  v2 
 
u3 
v3

ou encore sous forme condensée

u = Nu(e) (6.8)
les fonctions de forme Ni étant données par
ai + bi x + ci y
Ni (x, y) = (6.9)
2A
où
a1 = x2 y3 − x3 y2 a2 = x3 y1 − x1 y3 a3 = x1 y2 − x2 y1
b 1 = y2 − y3 b2 = y3 − y1 b3 = y 1 − y 2 (6.10)
c 1 = x3 − x2 c2 = x1 − x3 c3 = x2 − x1
et où A est l’aire de l’élément fini triangulaire

1 x1 y 1
1 1[ ]
A = det 1 x2 y2 = x2 y3 − x3 y2 + x1 (y2 − y3 ) + y1 (x3 − x2 ) (6.11)
2 2
1 x3 y 3

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-5

6.3.2 Déformations

Tel que vu au chapitre 2, les déformations sont reliées aux déplacements par
 
∂
0
   ∂x 
 [ ]
εx   u
   ∂ 
 εy  =  0  (6.12)
 ∂y  v
γxy  
∂ ∂ 
∂y ∂x
soit, en utilisant Eq. (6.7)
   
∂ u1
0
   ∂x   
 [ ]  v1 
εx  ∂  0 N3 0  
    N1 0 N2  u2 
 εy  =  0   
 ∂y  0 N1 0 N2 0 N3  v2 
γxy    
∂ ∂   u3 
∂y ∂x v3
   (6.13)
∂N1 ∂N2 ∂N3 u1
0 0 0
 ∂x ∂x ∂x   v1 
  
 ∂N3   
 ∂N1 ∂N2  u2 
= 0 0 0  
 ∂y ∂y ∂y   v2 
  
 ∂N1 ∂N1 ∂N2 ∂N2 ∂N3 ∂N3  u3 
∂y ∂x ∂y ∂x ∂y ∂x v3
L’équation 6.13, exprimant les déformations en fonction des déplacements aux noeuds de l’élément fini,
s’écrit sous la forme condensée
ε = Bu(e) (6.14)
où, comme pour le cas unidimensionnel, la matrice B est appelée matrice déformations-déplacements. En
utilisant les équations (6.9) à (6.10) explicitant les fonctions de forme, on vérifie aisément l’expression de la
matrice B d’un élément triangulaire en fonction des coordonnées de ses noeuds
 
y2 − y3 0 y3 − y1 0 y1 − y2 0
1 


B=  0 x3 − x2 0 x1 − x3 0 x2 − x1  (6.15)
2A  
x3 − x2 y2 − y3 x1 − x3 y3 − y1 x2 − x1 y1 − y2
On constate que la matrice B est constante pour un élément triangulaire à trois noeuds. Ce résultat était
prévisible puisque les fonctions de forme sont linéaires.

6.3.3 Contraintes

Tel que décrit au chapitre 2, les contraintes et les déformations sont reliées par la matrice constitutive C [Eq. (2.119)]
σ = Cε (6.16)

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-6

Rappelons que cette équation se simplifie aux équations (2.125) et (2.127) pour les états de déformations
planes et de contraintes planes, respectivement. En utilisant Eq. (6.14), on obtient l’expression des contraintes
en fonction des déplacements aux noeuds
σ = CBu(e) (6.17)

6.3.4 Énergie potentielle

Reprenons l’expression générale de l’énergie potentielle

∫ ∫ ∫ ∑
1
Π= σ ε dV −
T
u f dV −
T V
uT f S dS − uTi Pi (6.18)
2 V V S i

où V est le volume du domaine d’étude et S la surface délimitant le domaine.

Dans le cas d’un modèle d’élasticité plane, l’élément de volume dV vaut tdA où t est l’épaisseur constante
du modèle et dA un élément de surface. L’équation (6.18) s’écrit donc
∫ ∫ ∫ ∑
1
Π= σ T ε t dA − uT f V t dA − uT f S t dℓ − uTi Pi (6.19)
2 AT AT LS
| {z } | {z } | {z } i
Π1 Π2 Π3

où AT est l’aire totale du modèle, et où LS est la partie du périmètre du modèle où les forces de surface sont
appliquées. Les termes Π1 , Π2 et Π3 résultent de la contribution de tous les éléments finis du modèle, on
peut donc écrire

∑
m ∫
1
Π1 = σ T ε t dA (6.20)
2 A(e)
e=1
∑m ∫
Π2 = uT f V t dA (6.21)
e=1 A(e)
mS ∫
∑
Π3 = uT f S t dℓ (6.22)
(e)
e=1 LS

où m est le nombre d’éléments finis du modèle, A(e) l’aire de l’élément e, mS le nombre d’éléments finis
(e)
où les forces de surface sont appliquées, et LS la partie du périmètre de l’élément e sur lequel la force de
surface est appliquée.

6.3.5 Matrice de rigidité élémentaire

L’énergie de déformation d’un élément fini e est définie par

∫ ∫
1 1
U (e) = σ T ε t dA = εT σ t dA (6.23)
2 A(e) 2 A(e)

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-7

En tenant compte des équations (6.14) et (6.4.2), Eq. (6.23) devient

∫
1 T
U =
(e)
u(e) BT CBu(e) t dA
2 A(e)
(∫ ) (6.24)
1 (e) T
= u B CB t dA u(e)
T
2 A (e)

ou encore
1 (e) T (e) (e)
U (e) = u k u (6.25)
2
où k(e) est la matrice de rigidité élémentaire de l’élément e
∫
(e)
k = BT CB t dA (6.26)
A(e)

Les matrices B et C étant constantes sur un élément d’élasticité plane triangulaire à trois noeuds, la matrice
de rigidité élémentaire a pour expression
∫
(e) (e) T (e)
k =B CB t dA
A(e) (6.27)
(e) (e) T (e)
= tA B CB

6.3.6 Vecteur des forces de volume

Le terme Π2 de l’équation (6.19) représente la contribution des forces de volume. D’après Eq. (6.8), Π2
devient
∑m ∫
Π2 = uT f V t dA
(e)
e=1 A
∑m ∫ (6.28)
(e) T T V
= u N f t dA
e=1 A(e)

ou encore
∑
m
T (e)
Π2 = u(e) f V (6.29)
e=1
(e)
où f V est le vecteur élémentaire des forces de volume, défini par
∫
V (e)
f = NT f V t dA (6.30)
A(e)

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-8

6.3.7 Vecteur des forces de surface

La contribution des forces de surface est représentée par le terme Π3 de l’équation (6.19). En suivant le
même raisonnement que pour les forces de volume, on peut écrire
∑m ∫
Π3 = uT f S t dℓ
(e)
e=1 L
∑m ∫ (6.31)
(e) T ST S
= u N f t dℓ
(e)
e=1 LS

où NS est la matrice de forme définie à partir de Eq. (6.9), mais évaluée à l’endroit où les forces de surface
sont appliquées. L’équation (6.32) s’écrit également
∑m ∫
(e) T T
Π3 = u NS f S t dℓ
(e)
e=1 LS
(6.32)
∑m
(e) T S (e)
= u f
e=1
(e)
où f S est le vecteur élémentaire des forces de surface défini par
∫
(e) T
fS = NS f S t dℓ (6.33)
(e)
LS

6.3.8 Vecteur des forces concentrées

Comme dans le cas des problèmes unidimensionnels, la contribution des forces concentrées appliquées au
∑
modèle est représentée par le terme i uTi Pi où i est le point d’application de la force concentrée Pi . Seules
les éléments ayant des noeuds confondus avec les points d’application des forces contribuent à ce terme.

6.3.9 Équations de la MEF

En suivant la même démarche que pour les problèmes unidimensionnels, on démontre que l’application du
principe de l’énergie potentielle minimum aboutit à un système d’équations linéaires
KU = FV + FS + FC (6.34)
où U est le vecteur global des déplacements et où
m ∫
∑ ∑
m
K= BT CB t dA = k(e)
e=1 A(e) e=1
m ∫
∑ ∑
m
V T V (e)
F = N f t dA = fV (6.35)
e=1 A(e) e=1
m ∫
∑ ∑
m
S ST S (e)
F = N f t dℓ = fS
(e)
e=1 LS e=1

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-9

Exemple 6.1
Considérer le mur illustré sur la figure 6.1a. La réponse du mur au chargement indiqué est étudiée en uti-
lisant le maillage grossier à base d’éléments finis triangulaires proposé (Figure 6.1b). Déterminer les ma-
trices de rigidité et les vecteurs force élémentaires des éléments 1 et 3 du maillage. Considérer pour l’appli-
cation numérique : ρ = 2400 kg/m3 ; E = 40000 MPa ; L = 1 m ; h = 2 m ; t = 0,20 m ; g = 9.81 m/s2 .

Figure 6.3– Mur en état de contraintes planes : (a) Géométrie 3D ; (b) Maillage grossier 2D à base
d’éléments triangulaires.

Élément 1
– Aire de l’élément
1
A(1) = × 1 × 1 = 0,5
2
– Fonctions de forme
L’application de l’équation Eq. (6.10) aux coordonnées des noeuds de l’élément 1, donne les fonctions
de forme [Eq. (6.9)]
(1) (1) (1)
N1 (x, y) = 1 − y ; N2 (x, y) = x ; N3 (x, y) = −x + y ;

– Matrice déformations-déplacements
D’après Eq. (6.13) ou Eq. (6.15)
 
0 0 1 0 −1 0
 
B(1) =  0 −1 0 0 0 1
−1 0 0 1 1 −1

– Matrice constitutive
L’épaisseur du mur est petite par rapport aux autres dimensions, et le mur est chargé dans son plan. Le
mur peut donc être considéré en état de contraintes planes. La matrice constitutive C est donc donnée

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-10

par    
1 ν 0   4,3956 1,3187 0 
E  ν 1
  
C= 0  
 =  1,3187 4,3956
 × 1010
1 − ν2  
0
   
1−ν
0 0 0 0 1,5385
2
– Matrice de rigidité élémentaire
D’après Eq. (6.26), la matrice de rigidité de l’élément 1 est donnée par
∫
(1)
k = BT CB t dA
A(1)
T
= t A(1) B(1) CB(1)

d’où
 
1,5385 0 0 −1,5385 −1,5385 1,5385
 0 4,3956 −1,3187 0 1,3187 −4,3956 
 
 0 −1,3187 0 −4,3956 1,3187 
 4,3956 
k(1) =  × 103 MN/m
 −1,5385 0 0 1,5385 1,5385 −1,5385 
 
 −1,5385 1,3187 −4,3956 1,5385 5,9341 −2,8571 
1,5385 −4,3956 1,3187 −1,5385 −2,8571 5,9341

– Vecteur des forces de volume

L’intensité de la force de volume est f V = ρg, agissant vers le bas (sens contraire de l’axe y), d’où
[ ]
0
fV =
−ρg

Le vecteur des forces de volume de l’élément 1 est donné par l’équation (6.30)
∫
(1)
fV = NT f V t dA
A(1)
 
0
 ∫ 1 (∫ 1 ) 
 
 −t ρg(1 − y) dy dx 
 
 0 x 
 
 
 0 

= ∫ ( ∫ ) 
1 1 
 −t ρgx dy dx 
 
 0 x 
 
 0 
 ∫ (∫ ) 
 1 1 
−t ρg(−x + y) dy dx
0 x

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-11

soit, tous calculs faits    

0 0
1 1
   
(1) tρg 
0
 0
 
fV =−   = −0,785   (kN)
6 1 1
   
0 0
1 1
– Vecteur des forces de surface
Un des bords de l’élément 1 est soumis à un chargement variant linéairement en fonction de y
1
f S (y) = y (kN/m2 ) pour 0 ≤ y ≤ 1
2
d’où  
1
y
fS = 2 
0
La matrice de forme NS est obtenue en évaluant N à l’endroit où les forces de surface sont appliquées,
soit en posant dans le cas de l’élément 1, x = 0
[ ]
1 − y 0 0 0 y 0
NS =
0 1−y 0 0 0 y

Le vecteur élémentaire des forces de surface est obtenu à partir de l’équation (6.33)
∫
S (1) T
f = NS f S t dℓ
(1)
L
 S∫ 1 
1
t (1 − y)y dy 
 0 2 
 
 0 
 
 
 
 0 
=



 0 
 ∫ 1 
 
 1 2 
 t y dy 
 2 
 0 
0

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-12

d’où  t   
0,0167
 12   
   
0  0 
   
   
   
0  0 
=    ((kN))
(1)
fS  = 
0  0 
   
 t   
  0,0333
   
6  
0 0

Élément 3
On vérifie aisément que l’élément 3 a les mêmes fonctions de forme, la même matrice déformations-
déplacements et donc la même matrice de rigidité que l’élément 1. Les deux éléments ont également le
même vecteur force de volume. Leur vecteurs force de surface sont cependant différents.
– Vecteur des forces de surface
Deux bords de l’élément 3 sont soumis à des forces de surface : l’une variant linéairement, et l’autre
constante. Le vecteur force appliqué associé à la force variant linéairement s’écrit
 
1
(y + 1)
faS =  2  kN/m2
0

La matrice de forme NS est obtenue dans ce cas en évaluant N en x = 0

[ ]
1 − y 0 0 0 y 0
NS =
0 1−y 0 0 0 y

D’où le vecteur élémentaire associé à la force variant linéairement

∫
S (3) T
fa = NS f S t dℓ
(3)
L
 S∫ 1 
1
t (1 − y ) dy 
2
 0 2 
 
 0 
 
 
 
 0 

= 

 0 
 ∫ 
 1 
 1 
t (1 + y)y dy 
 0 2 
 
0

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-13

soit, tous calculs faits  t   

0,0667
3  
   
0  0 
   
   
   
0  0 
faS
(3)

= =  

0  0 
   
 5t   
  0,0833
   
 12   
0 0

Par ailleurs, le vecteur force appliqué associé à la force constante a pour expression
[ ]
0
fbS = kN/m2
−0,80

La matrice de forme NS est obtenue cette fois-ci en évaluant N en y = 1

[ ]
0 0 x 0 −x + 1 0
NS =
0 0 0 x 0 −x + 1

D’où le vecteur élémentaire associé à la force constante

∫
S (3) T
fb = NS f S t dℓ
(3)
LS
 
 0 
 
 
 0 
 
 
 0 
 ∫ 

= 1 

 −t 0,80x dx 
 
 0 
 
 0 
 ∫ 1 
 
t (0,80x − 0,80) dx
0

soit    
0 0
 0   0 
   
 0   
(3)    0 
fbS = =  (kN)
−0,40t −0,0800
   
 0   0 
−0,40t −0,0800

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-14

Finalement, le vecteur force de surface associé à l’élément 3 est la somme des contributions des forces
linéaire et constantes
(3) (3) (3)
fS = faS + fbS
 
0,0667
 0 
 
 
 0 
=  (kN)
 −0,0800 
 
 0,0833 
−0,0800

L’élément fini triangulaire étudié jusqu’à présent (figure 6.2) est dit linéaire puisque la variation des dépla-
cements à l’intérieur d’un tel élément est linéaire [Eq. (6.2)]. C’est historiquement le premier élément fini
développé pour résoudre des problèmes d’élasticité plane.

6.3.10 Comportement de l’élément triangulaire linéaire – CST

En utilisant Eq. (6.4.2), on peut calculer les contraintes au sein d’un élément fini triangulaire linéaire. La
matrice B étant constante [Eq. (6.15)], les contraintes le sont également au sein de chaque élément fini de ce
type. L’élément est connu dans la littérature sous le nom élément fini triangulaire à déformations constantes
ou élément CST 1 . Le calcul des contraintes au sein des éléments CST est illustré à l’aide de l’exemple 6.2.

Exemple 6.2
La résolution du système d’équations linéaires obtenus pour le mur de refend de l’exemple 6.1 aboutit au
vecteur des déplacements globaux
 
0,000
 0,000 
 
 
 0,000 
 
 0,000 
 
 0,377 
 
 −7,880 
 
U=  × 10−7 (m)
 1,963 
 
 −9,029 
 
 
 3,407 
 
 −10,021 
 
 3,806 
−13,204
1. De l’anglais Constant-Strain Triangle.

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-15

Ces déplacements sont donnés selon la numérotation globale montrée sur la figure 6.1. Calculer les
contraintes au sein du mur de refend.

La numérotation locale illustrée sur la figure 6.4 est adoptée. Les déplacements u(e) aux noeuds de chaque
élément e du modèle sont
   
0,000 0,000
 0,000   0,000 
   
 1,963   0,000 
   
u(1) =   × 10−7 (m); u(2) =   × 10−7 (m)
 −9,029   0,000 
   
 0,377   1,963 
−7,880 −9,029
   
0,377 0,377
 −7,880   −7,880 
   
   1,963 
 3,806  −7  
(3)
u =  × 10 (m); (4)
u =  × 10−7 (m)
 −13,204   −9,029 
   
 3,407   3,806 
−10,021 −13,204

Figure 6.4– Numérotation locale des

noeuds de chaque élément du modèle.

Les matrices B(e) reliant les déformations au sein de chaque élément e à ses déplacements aux noeuds
sont calculées à l’aide de Eq. (6.15)
 
0 0 1 0 −1 0
 
B(1) = B(3) =  0 −1 0 0 0 1
−1 0 0 1 1 −1
 
−1 0 1 0 0 0
 
B(2) = B(4) = 0 0 0 −1 0 1
0 −1 −1 1 1 0

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-16

En considérant un état de contraintes planes, la matrice constitutive C est donnée par

   
1 ν 0   4,3956 1,3187 0 
E  ν 1
  
0  = 0 
 × 10
10
C=    1,3187 4,3956
1−ν  2
  
1−ν
0 0 0 0 1,5385
2
Les contraintes au sein de chaque élément sont ensuite calculées en utilisant Eq. (6.4.2)
 
−3,42
 
σ (1) = CB(1) u(1) =  −32,55  (kPa)
−1,19
 
−11,91
 
σ (2) = CB(2) u(2) =  −39,69  (kPa)
3,02
 
−1,07
 
σ (3) = CB(3) u(3) =  −8,88  (kPa)
−0,24
 
1,47
 
σ (4) = CB(4) u(4) =  −16,26  (kPa)
1,07

La figure 6.5 montre la variation des contraintes σy dans le mur de refend.

Figure 6.5– Variation des contraintes σy dans le mur.

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.3. ÉLÉMENTS FINIS TRIANGULAIRES 6-17

Les résultats de l’exemple 6.2 montre que l’élément CST n’est pas approprié pour modéliser des zones à
grandes variations de contraintes et de déformations. L’élément est également sensible au verrouillage en
cisaillement et doit être utilisé avec précaution. Ce comportement défectueux est illustré ci-après.

Considérons un assemblage de deux éléments CST, soumis à un moment de flexion M tel que montré que
la figure 6.6a. Supposons que sous l’effet du moment de flexion, le noeud 3 du maillage se déplace des
quantités u3 et v3 selon les axes x et y, respectivement (figure 6.6b). Les noeuds 1 et 2 sont fixes. Étudions
le comportement de l’élément 1. La configuration déformée de cet élément est montrée sur la figure 6.6b.
L’interpolation linéaire des déplacement au sein de l’élément donne [Eq. (6.2)]

u(x, y) = α1 + α2 x + α3 y
v(x, y) = β1 + β2 x + β3 y

Figure 6.6– Comportement d’un élément fini CST.

Les coordonnées généralisées αi et βj sont à déterminer à partir des déplacements aux noeuds en écrivant
{ { {
u1 = α1 = 0 ; u2 = α1 + α3 h = 0 ; u3 = α1 + α2 b
v1 = β1 = 0 ; v2 = β1 + β3 h = 0 ; v3 = β1 + β2 b

d’où
u3
α1 = 0 ; α2 = ; α3 = 0
b
v3
β1 = 0 ; β2 = ; β3 = 0
b

Les déplacements u et v au sein de l’élément 1 deviennent donc

u3
u(x, y) = x
b
v3
v(x, y) = x
b

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-18

d’où les déformations

∂u u3 ∂v ∂u ∂v v3
εx = = ; εy = = 0; γxy = + =
∂x b ∂y ∂y ∂x b

La formulation de l’élément CST implique donc la présence d’une déformation de cisaillement γxy qui
n’existe pas dans la réalité car l’assemblage de la figure 6.6 est soumis à un moment de flexion pure. Ce
cisaillement parasite absorbe de l’énergie de déformation et rend l’élément CST artificiellement plus rigide.
Soulignons que le verrouillage en cisaillement conduit à des résultats non-conservateurs, d’où la nécessité
de bien évaluer son importance dans un modèle donné.

Pour minimiser les effets des défauts associés à l’élément CST, il faut s’assurer que la densité du maillage
est suffisante pour capter adéquatement les variations des contraintes et des déformations. L’élément CST
est facile à programmer et les maillages à base de triangles peuvent être construits automatiquement grâce
notamment à la technique de triangulation de Delaunay. Ces deux points constituent les principaux attraits
faisant en sorte que l’élément CST est encore très utilisé dans la pratique en dépit de ses défauts.

6.4 Éléments finis rectangulaires

La figure 6.7 montre un élément finis rectangulaire à quatre noeuds numérotés localement. Les huit DDL de
l’élément sont regroupés dans un vecteur-colonne u(e)
 
u1
v 
 1
 
 u2 
 
 v2 
u =
(e)
u 
 (6.36)
 3
 
 v3 
 
 u4 
v4

La démarche suivie au paragraphe 6.3 s’applique. Nous faisons l’hypothèse que les déplacements à l’intérieur
de chaque élément varient linéairement en fonction des coordonnées

u(x, y) = α1 + α2 x + α3 y + α4 xy
(6.37)
v(x, y) = β1 + β2 x + β3 y + β4 xy

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-19

Figure 6.7– Élément fini rectangulaire numéroté

localement.

ou sous forme matricielle

 
α1
α 
 2
 
[ ] [ ]α 
 3
u(x, y) 1 x y xy 0 0 0 0 α4 


= (6.38)
v(x, y) 0 0 0 0 1 x y xy  
 β1 
 
 β2 
 
 β3 
β4
Appliquée aux noeuds de l’élément, Eq. (6.38) donne
    
u1 1 a b ab 0 0 0 0 α1
v   0  α 
 1  0 0 0 1 a b ab   2
    
u2   1 −a b −ab 0 0 0 0  α3 
    
  
v2   0 0 0 0 1 −a b −ab  α4 
u(e) = u  =  1 −a −b
 
 β  (6.39)
 3  ab 0 0 0 0   1
    
 v3   0 0 0 0 1 −a −b ab   β2 
    
 u4   1 a −b −ab 0 0 0 0   β3 
v4 0 0 0 0 1 a −b −ab β4
d’où le vecteur des coordonnées généralisées αi et βj
   −1  
α1 1 a b ab 0 0 0 0 u1
α   0  v 
 2  0 0 0 1 a b ab   1
     
α3   1 −a b −ab 0 0 0 0  u2 
     
α4   0 0 0 0 1 −a b −ab   v2 
 =    (6.40)
 β   1 −a −b  u 
 1  ab 0 0 0 0   3
     
 β2   0 0 0 0 1 −a −b ab   v3 
     
 β3   1 a −b −ab 0 0 0 0  u4 
β4 0 0 0 0 1 a −b −ab v4

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-20

En remplaçant dans l’équation (6.38)

 −1  
1 a b ab 0 0 0 0 u1
  v 
 0 0 0 0 1 a b ab   1
   
[ ] [ ]

1 −a b −ab 0 0 0 0 

 u2 
 
u(x, y) 1 x y xy 0 0 0 0 
 0 0 0 0 1 −a b −ab 

 v2 
 
=
v(x, y) 0 0 0 0 1 x y xy 
 1 −a −b ab 0 0 0 0 

u 
 3
   
 0 0 0 0 1 −a −b ab   v3 
   
 1 a −b −ab 0 0 0 0   u4 
0 0 0 0 1 a −b −ab v4
(6.41)

L’équation (6.41) se simplifie à

 
u1
v 
 1
 
[ ] [ ] u 
 2
u(x, y) N1 0 N2 0 N3 0 N4 0  
 v2 
= (6.42)
v(x, y) 0 N1 0 N2 0 N3 0 N 4  u3 


 
 v3 
 
 u4 
v4

ou encore sous forme condensée

u = Nu(e) (6.43)
les fonctions de forme Ni étant données par

(a + x)(b + y) (a − x)(b + y)
N1 (x, y) = ; N2 (x, y) =
4ab 4ab
(6.44)
(a − x)(b − y) (a + x)(b − y)
N3 (x, y) = ; N4 (x, y) =
4ab 4ab

6.4.1 Déformations

En suivant la même procédure qu’au paragraphe 6.3.2, on démontre que

  
∂N1 ∂N2 ∂N3 ∂N4 u1
0 0 0 0
   ∂x ∂x ∂x ∂x   v1 
  
εx  ∂N4   
   ∂N1 ∂N2 ∂N3   u2 
 εy  =  0 0 0 0   (6.45)
 ∂y ∂y ∂y ∂y   v2 
γxy   
 ∂N1 ∂N1 ∂N2 ∂N2 ∂N3 ∂N3 ∂N4 ∂N4  u3 
∂y ∂x ∂y ∂x ∂y ∂x ∂y ∂x v3

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-21

ou sous forme matricielle

ε = Bu(e) (6.46)
où la matrice B est donnée par
 
(b + y) 0 −(b + y) 0 −(b − y) 0 (b − y) 0
1  


B=  0 (a + x) 0 (a − x) 0 −(a − x) 0 −(a + x)
4ab  
(a + x) (b + y) (a − x) −(b + y) −(a − x) −(b − y) −(a + x) (b − y)
(6.47)
Les éléments de la matrice B sont donc des fonctions linéaires de x et y.

6.4.2 Contraintes, matrice de rigidité et vecteur des forces élémentaires

Le vecteur des contraintes σ et la matrice de rigidité k(e) de l’élément fini rectangulaire sont données par les
mêmes équations (6.4.2) et (6.26)
σ = CBu(e)
et ∫ ∫ ∫
b a
(e) T
k = B CB t dA = BT CB t dxdy (6.48)
A(e) −b −a

où A(e) est l’aire de l’élément e et t sont épaisseur. Les termes de forces sont établies de la même façon
que pour l’élément triangulaire. L’assemblage des matrices et des vecteurs élémentaires pour constituer les
équations du modèle se fait également en suivant la même procédure.

6.4.3 Comportement de l’élément rectangulaire linéaire – Q4

L’élément fini rectangulaire que l’on vient de décrire est qualifié de linéaire en vertu de la variation linéaire
des déplacements s’y développant. L’élément est connu dans la littérature sous le nom d’élément Q4.

Considérons un tronçon d’une poutre soumis à un moment de flexion pure (figure 6.8a). En supposant un
état de contraintes planes, les déformations au sein du tronçon valent

θ θ
εx = − y; εy = ν y; γxy = 0 (6.49)
2a 2a
Le tronçon est modélisé par un élément fini Q4 dont la configuration déformée est illustrée sur la figure 6.8b.
En utilisant les équations (6.46) et (6.47), on démontre que les déformations dans l’élément Q4 sont

θ̃ θ̃
εx = − y; εy = 0 ; γxy = − x (6.50)
2a 2a

où l’angle θ̃ est défini par la configuration déformée de l’élément fini (figure 6.8b). Comme dans le cas
de l’élémet CST, la formulation de l’élément Q4 entraîne l’apparition d’une déformation de cisaillement
parasite. On vérifie par ailleurs que sous l’effet d’un même moment appliqué M = M f, le rapport entre

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-22

l’énergie de déformation Ue de l’élément Q4 et l’énergie de déformation U du tronçon de poutre est donné

par
Ue θ̃ 2(1 − ν 2 )
= = (6.51)
U θ 1−ν
2+ 2
r
où r est le ratio géométrique r = b/a. Ce résultat signifie que plus l’élément Q4 est étiré (r → 0), plus il
devient verrouillé en cisaillement en absorbant plus d’énergie de déformation.

Figure 6.8– Vérouillage en cisaillement de l’élément fini Q4 : (a) Tronçon de poutre soumis à un moment de flexion
pure ; (b) Modélisation avec un élément fini Q4.

6.4.4 Élément rectangulaire avec modes incompatibles – Q6

Considérons un élément fini rectangulaire Q4, et réécrivons l’équation Eq. (6.43) en ajoutant deux nouvelles
fonctions de forme N5 et N6
 
[ ] β1
β 
N 5 0 N 6 0  2
u = Nu(e) +   (6.52)
0 N5 0 N6 β3 
β4

où
x2 y2
N5 (x, y) = 1 − ; N6 (x, y) = 1 − (6.53)
a2 b2

et où les paramètres β1 à β4 sont des degrés de liberté artificiels car ils ne sont associés à aucun noeud de
l’élément fini. Les fonctions de forme quadratiques N5 et N6 sont montrées sur la figure 6.9. L’élément fini
est connu sous le nom d’élément fini Q6.

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.4. ÉLÉMENTS FINIS RECTANGULAIRES 6-23

Figure 6.9– Fonctions de forme quadratiques d’un élément fini rectangulaire avec modes incompatibles : (a) fonction
de forme N5 ; (b) fonction de forme N6 .

La figure 6.9 illustre la configuration déformée de cet élément sous l’effet d’un moment de flexion pure.
Il apparaît alors que cet élément permet de mieux représenter un comportement flexionnel. Notons que ce
mode de déformation est possible grâce à l’ajout des fonctions de forme quadratique. Nous verrons plus tard
que ce mode met en défaut la compatibilité des déformations au sein d’un maillage, d’où le nom de mode
incompatible.

Figure 6.10– Déformation d’un élément fini Q6

sous l’effet d’un moment de flexion pure.

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-24

Le champs des déformations résultant de l’équation (6.52) est donné par

4 (
∑ )
∂Ni 2x
εx = ui − β1
∂x a2
i=1
∑4 ( )
∂Ni 2y
εy = vi − 2 β4 (6.54)
∂y b
i=1
∑4 ( )
∂Ni ∂Ni 2x 2y
γxy = ui + vi − 2 β2 − 2 β3
∂y ∂x a b
i=1

On voit donc que, contrairement à l’élément Q4, l’ajout des fonctions de forme quadratiques permet d’an-
nuler le cisaillement parasite γxy et d’obtenir la déformation due à l’effet de Poisson εy = −νεx . L’élément
Q6 est donc recommandé pour des calculs où le comportement flexionnel est prédominant. Nous discuterons
plus tard de la performance de cet élément en terme de convergence.

6.5 Éléments finis d’ordre supérieur

6.5.1 Formulation isoparamétrique

Tel que mentionné précédemment, la formulation isoparamétrique est devenu le standard utilisé dans les
logiciels d’éléments finis. En effet, elle se généralise facilement au cas bidimensionnel et la figure 6.11
montre les axes des coordonnées naturelles d’éléments finis triangulaires ou quadrilatères.

Figure 6.11– Axes des coordonnées naturelles d’éléments finis courants : (a) Élément triangulaire ; (b) Élément qua-
drilatère.

Dans le cas plus général d’un élément fini e à n noeuds, et 2n degrés de libertés ui et vi , les équations (4.133)
et (4.134) se généralisent à

∑
n ∑
n
u(r, s, t) = Ni (r, s, t)ui ; v(r, s, t) = Ni (r, s, t)vi (6.55)
i=1 i=1
∑n ∑n
X(r, s, t) = Ni (r, s, t)Xi ; Y (r, s, t) = Ni (r, s, t)Yi (6.56)
i=1 i=1

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-25

où u, v sont les déplacements selon les axes X et Y et où r et s sont les coordonnées naturelles, avec
−1 ≤ r ≤ 1 et −1 ≤ s ≤ 1.

Les fonctions de forme des éléments unidimensionnels de 2 à 4 noeuds ont été présentés au tableau 4.2.
L’équation (4.147) se généralise dans le cas bidimensionnel à

∂X ∑ ∂Ni
n
∂Y ∑ ∂Ni n
= Xi ; = Yi
∂r ∂r ∂r ∂r
i=1 i=1
(6.57)
∂X ∑
n
∂Ni ∂Y ∑
n
∂Ni
= Xi ; = Yi
∂s ∂s ∂s ∂s
i=1 i=1

où les dérivées des fonctions de forme dans les systèmes de coordonnées naturel et global sont reliées par
   
∂Ni ∂Ni
   
 ∂r   
  = J  ∂X  (6.58)
 ∂Ni   ∂Ni 
∂s ∂Y
où J est la matrice Jacobienne définie par
 
∂X ∂Y
 
 ∂r 
J =  ∂r  (6.59)
 ∂X ∂Y 
∂s ∂s

Notons J le déterminant de la matrice Jacobienne

J = det J (6.60)

et dA(XY ) , dA(xy) et dA(rs) les éléments de volume exprimés respectivement dans les systèmes d’axes
global, local et naturel

dA(XY ) = dXdY ; dA(xy) = dxdy ; dA(rs) = drds (6.61)

On démontre que
dA(XY ) = dA(xy) = J dA(rs) (6.62)

Compte tenu de l’équation (6.62), l’intégration de toute fonction F sur un élément de surface dans le système
d’axes global ou local peut s’effectuer en intégrant sur un élément de surface dans le système d’axes naturel
∫ ∫ ∫ 1 ∫ 1
F dXdY = F dxdy = F J drds (6.63)
A(e) A(e) −1 −1

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-26

En particulier, la matrice de rigidité élémentaire a pour expression

∫ ∫ 1 ∫ 1
(e) T
k = B CB t dA = BT C B J t drds (6.64)
A(e) −1 −1

On démontre également la relation suivante pour les éléments de ligne

dℓ = J ℓ dr (6.65)

où √( )2 ( )2
ℓ ∂x ∂y
J = + (6.66)
∂r ∂r

Une intégrale de ligne se calcule alors en écrivant

∫ ∫ ∫ 1
F dX = F dx = F J ℓ dr (6.67)
ℓ(e) ℓ(e) −1

6.5.2 Intégration numérique

Tel que discuté précédemment, les matrices et vecteurs élémentaires de la méthode des éléments finis ne
peuvent pas toujours être calculés de façon explicite exacte comme c’était le cas à l’exemple 6.1. L’évalua-
tion d’une intégrale surfacique peut s’effectuer en appliquant l’équation (4.153) à chacune des différentes
directions. On démontre que dans ce cas
∫ 1 ∫ 1 ∑
nG ∑
mG
I= F(r, s) drds ≈ wi wj F(ri , sj ) (6.68)
−1 −1 i=1 j=1

où nG et mG sont les nombres de points d’intégration de Gauss utilisés dans chaque direction. On parle alors
d’intégration d’ordre nG × mG .

Exemple 6.3
Évaluer l’intégrale suivante en utilisant la quadrature de Gauss
∫ 1∫ 1
I= rs4 (r + s) drds
−1 −1

La fonction à intégrer est une fonction polynôme en r et s. Le nombre de points d’intégration mini-
mum pour que l’intégration soit exacte dans la direction r est nG = 2. Le nombre de points d’intégration
minimum pour que l’intégration soit exacte dans la direction s est mG = 3. L’intégration numérique sera

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-27

donc effectuée à l’ordre 2 × 3. Posons F(r, s) = rs4 (r + s). En se référant au tableau 4.3, on peut donc
écrire
∑
nG ∑
mG
I≈ wi wj F(ri , sj )
i=1 j=1

Pour mieux visualiser les différents termes, la figure 6.12 montre les positions des points de Gauss en deux
dimensions. Le domaine d’intégration est matérialisé par un élément fini quadrilatère.

Figure 6.12– Intégration numérique en deux

dimensions.

Notons F(r, s) = rs4 (r + s). Les différents termes intervenant dans l’évaluation numérique de l’inté-
grale sont résumées au tableau 6.1.

Tableau 6.1– Points d’intégrations de Gauss et coefficients de pondération.

Point de Gauss ri s F(ri , sj ) wi wj wi wj F(ri , sj )
√ √j
P1 −1/√3 0,6 −0,0410 1,0000 5/9 −0,0228
P2 −1/√3 0 0,0000 1,0000 8/9 0,0000
√
P3 −1/√3 − 0,6 0,2810 1,0000 5/9 0,1561
√
P4 1/√3 0,6 0,2810 1,0000 5/9 0,1561
P5 1/√3 0 0,0000 1,0000 8/9 0,0000
√
P6 1/ 3 − 0,6 −0,0410 1,0000 5/9 −0,0228

d’où

I = −0,0228 + 0,0000 + 0,1561 + 0,1561 + 0,0000 − 0,0228

= 0,2667

Soulignons que cette valeur se compare très bien au résultat de l’intégration exacte
∫ 1∫ 1
4
I= rs4 (r + s) drds =
−1 −1 15

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-28

En suivant la même démarche décrite à l’exemple 6.3, on peut déterminer le nombre de points d’intégration
minimum pour que l’intégration numérique des matrices et des vecteurs élémentaires d’un élément fini donné
donne des résultats exacts, ou très proches des valeurs exactes.

Notons que l’évaluation d’une intégrale surfacique sur un triangle ne s’effectue pas selon deux axes comme
c’est le cas pour des éléments quadrangulaires. L’intégration numérique s’écrit dans ce cas
∫ 1 ∫ 1−r ∑
nG
I= F(r, s) drds ≈ wi F(ri , si ) (6.69)
0 0 i=1

où nG sont les nombres de points d’intégration utilisés. On parle alors d’intégration d’ordre nG . Plusieurs
chercheurs ont proposé des techniques pour déterminer les positions optimales ri et si des points d’intégra-
tion, ainsi que les coefficients de pondération wi . On distingue notamment les techniques de Hammer, de
Stroud ou de Cowper. La figure 6.13 illustre les points d’intégration rapportés par Cowper. Le nombre des
points d’intégration minimum pour une intégration exacte est noté n∗G .

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani
6.5. ÉLÉMENTS FINIS D’ORDRE SUPÉRIEUR 6-29

Figure 6.13– Intégration numérique sur un élément fini triangulaire.

Analyse avancée des structures et éléments finis - Hiver 2011  Najib Bouaanani