Methodes des éléments finis

M´ ethode des ´ el´ ements finis : ´ elasticit´ e plane Yves Debard Universit´ e du Mans Master Mod´ elisation Num´ erique et R´ ealit´ e Virtuelle http://iut.univ-lemans.fr/ydlogi/index.html 24 mars 2006 – 29 mars 2011

Table des mati` eres 1 Rappels 1.1 Hypothèse contraintes planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Hypothèse déformations planes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 3

2 Forme diff´ erentielle

4

3 Forme int´ egrale faible

5

4 Forme discr´ etis´ ee : ´ el´ ements finis 4.1 Discrétisation du domaine : maillage . . . . . . . . . . . . . . . . . 4.2 Représentation élémentaire (ou locale) du champ de déplacements 4.3 Représentation globale du champ de déplacements . . . . . . . . . 4.4 Discrétisation de la forme intégrale faible . . . . . . . . . . . . . . 4.5 Problèmes particuliers . . . . . . . . . . . . . . . . . . . . . . . . . 4.5.1 Problème stationnaire . . . . . . . . . . . . . . . . . . . . . 4.5.2 Modes propres de vibration . . . . . . . . . . . . . . . . . . 4.6 Mise en œuvre pratique : calculs élémentaires et assemblage . . . . 5 Probl` eme ´ elastostatique : ´ energie potentielle 5.1 Calcul des variations . . . . . . . . . . . . . . ´ 5.2 Energie potentielle . . . . . . . . . . . . . . . 5.3 Méthode de Ritz et éléments finis . . . . . . .

et m´ ethode . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . .

7 7 8 8 9 10 10 11 11

de . . . . . .

Ritz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12 12 13 14

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

6 Calculs ´ el´ ementaires : ´ el´ ements isoparam´ etriques ´ ement isoparamétrique : définition . . . . . . . . 6.1 El´ 6.1.1 Représentation de la géométrie . . . . . . . 6.1.2 Maillage conforme . . . . . . . . . . . . . . 6.1.3 Représentation du champ de déplacements 6.2 Bibliothèque d’éléments . . . . . . . . . . . . . . . 6.2.1 Triangle à 3 nœuds . . . . . . . . . . . . . . 6.2.2 Triangle à 6 nœuds . . . . . . . . . . . . . . 6.2.3 Quadrangle à 4 nœuds . . . . . . . . . . . . 6.2.4 Quadrangle à 8 ou 9 nœuds . . . . . . . . . 6.3 Calcul des matrices et des vecteurs élémentaires . . 6.3.1 Transformation des dérivées . . . . . . . . . 6.3.2 Transformation des intégrales . . . . . . . . ´ 6.3.3 Evaluation numérique des intégrales . . . . 6.3.4 Calcul des matrices . . . . . . . . . . . . . . 6.3.5 Calcul des vecteurs . . . . . . . . . . . . . . 6.4 Qualité du jacobien. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

15 15 15 16 16 17 17 18 19 20 21 21 21 22 23 23 25

A Programmes Maple A.1 tri3 int : triangle à 3 nœuds . . . . . . . . . A.2 tri6 int : triangle à 6 nœuds . . . . . . . . . A.3 quad4 int : quadrangle à 4 nœuds . . . . . . A.4 quad8 int : quadrangle à 8 nœuds . . . . . . A.5 quad9 int : quadrangle à 9 nœuds . . . . . . A.6 Qualité du jacobien : quadrangle à 4 nœuds A.7 Qualité du jacobien : triangle à 6 nœuds . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

27 27 27 28 28 29 29 30

R´ ef´ erences

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

32

´ Elasticit´ e plane

1

Introduction Dans ce texte nous présentons la résolution d’un problème d’élasticité plane par la méthode des éléments finis. Nous adopterons les hypothèses suivantes : – Les déplacements et les déformations sont petits. – Le comportement du matériau est élastique et linéaire. – Le matériau est homogène et isotrope. E, ν, α et ρ sont respectivement le module de Young, le coefficient de Poisson, le coefficient de dilatation et la masse volumique du matériau. Le repère {O; x, y, z} est un repère orthonormé. ~ı, ~ et ~k sont les vecteurs unitaires des axes.

1 1.1

Rappels Hypoth` ese contraintes planes

Figure 1 – Plaque sollicitée dans son plan Un solide (figure 1) est en état de contraintes planes par rapport au plan {O; x, y}, s’il existe un repère {O; x, y, z}, tel qu’en tout point M du solide, le tenseur des contraintes soit de la forme :   ~ı ~ composantes sur  ~ k

T~ (M,~ı ) T~ (M, ~ ) T~ (M, ~k)   σxx σxy 0  σyx σyy 0  0 0 0

(1.1)

o` u σxx , σyy et σxy = σyx sont indépendants de z. L’axe ~k est donc, pour tous les points du solide, direction principale et la contrainte principale associée est nulle. Dans la formule (1.1), T~ (M, ~n) est le vecteur contrainte sur la facette ~n en M .

2

Méthode des éléments finis

Le tenseur des déformations se réduit à :   εxx 12 γxy 0 0 [ε(M )] =  12 γxy εyy 0 0 εzz

avec εzz =

−ν (σxx + σyy ) + α ∆T E

(1.2)

o` u ∆T est la variation de température. La loi de comportement s’écrit : {σ} = [D] ({ε} − {εth })

(1.3)

o` u: – {ε} est le vecteur déformation :        

  {ε} =



γxy

∂u ∂x ∂v ∂y

       

 εxx  εyy =      = 2 εxy       ∂v ∂u     + ∂y ∂x

(1.4)

La figure (2) montre la signification des composantes εxx , εyy et γxy du tenseur des déformations.

Figure 2 – Transformation d’un rectangle infiniment petit – {σ} est le vecteur contrainte :

  σxx  {σ} = σyy   σxy

(1.5)

– [D] est la matrice des coefficients élastiques :  [D] =

E 1 − ν2

1 ν ν 1  0 0

 0 0   1−ν 2

– {εth } représente les déformations d’origine thermique :   1 {εth } = α ∆T 1   0 Les déformations et les contraintes ne dépendent que des déplacements suivant x et y : ½ ¾ u(x, y; t) v(x, y; t)

(1.6)

(1.7)

(1.8)


1.2

3

Hypoth` ese d´ eformations planes

Un solide est en état de déformations planes par rapport au plan {O; x, y}, s’il existe un repère {O; x, y, z}, lié au solide, tel qu’en tout point du solide, le champ de déplacements soit de la forme :   u(x, y; t) v(x, y; t)   0

(1.9)

Le tenseur des déformations se réduit à : 

εxx [ε(M )] =  12 γxy 0

 γxy 0 εyy 0 0 0

1 2

(1.10)

Le tenseur des contraintes est alors de la forme :   σxx σxy 0 [σ(M )] = σxy σyy 0  avec σzz = ν (σxx + σyy ) − E α ∆T 0 0 σzz

(1.11)

La loi de comportement s’écrit : {σ} = [D] ({ε} − {εth })

(1.12)

o` u: – {ε} est le vecteur déformation :        

  {ε} =

 γxy

– {σ} est le vecteur contrainte :

∂u ∂x ∂v ∂y

       

 εxx  εyy =      = 2 εxy   ∂u ∂v         + ∂y ∂x

(1.13)

  σxx  {σ} = σyy   σxy

(1.14)

– [D] est la matrice des coefficients élastiques :   λ + 2µ λ 0 λ + 2 µ 0 [D] =  λ 0 0 µ

,

λ=

Eν (1 + ν)(1 − 2 ν)

,

µ=

E 2 (1 + ν)

(1.15)

λ et µ sont les coefficients de Lam´ e du matériau. – {εth } représente les déformations d’origine thermique :   1 {εth } = α ∆T 1   0

(1.16)

4

2


Forme diff´ erentielle

Le solide V limité par la frontière S est soumis à : ½ ¾ fV x – un champ de forces volumiques : {fV } = fV y

½ ¾ uP – des déplacements imposés sur la frontière Su : {uP } = vP

½ ¾ fσx – des forces surfaciques sur la frontière : Sσ = S − Su : {fσ } = fσy

Figure 3 – Charges et conditions aux limites Résoudre un problème d’élasticité plane consiste à chercher un champ de déplacements : ½ ¾ u(x, y; t) {u(x, y; t)} = v(x, y; t) tel que :

 ∂2u ∂σxx ∂σxy   ρ + + fV x =  2 ∂t ∂x ∂y  ∂2v ∂σyy ∂σxy   ρ + + fV y = 2 ∂t ∂x ∂y

(2.1a)

(2.1b)

en tout point du solide avec : – les relations cinématiques : εxx =

∂u ∂x

,

εyy =

∂v , ∂y

γxy =

∂u ∂v + ∂y ∂x

(2.1c)

– la loi de comportement (ou loi constitutive) : {σ} = [D] ({ε} − {εth }) – les conditions aux limites : ½ ¾ ½ ¾ u uP = v vP

· sur Su

,

σxx σxy σxy σyy

¸ ½

nx ny

(2.1d) ¾ =

½ ¾ fσx fσy

sur Sσ

(2.1e)

S = Su ∪ Sσ , Su ∩ Sσ = ∅ ½ ¾ nx o` u S est la surface du solide et la normale unitaire à S dirigée vers l’extérieur de V . ny Remarque : en pratique, il y a une partition de la surface S pour chaque composante du déplacement.


5

– les conditions initiales à l’instant t = t0 : ¾ ½ ¾ ½ ¾ ½ ¾ ½ u(x, y; t0 ) u0 (x, y) u(x, ˙ y; t0 ) u˙ 0 (x, y) = , = v(x, y; t0 ) v0 (x, y) v(x, ˙ y; t0 ) v˙ 0 (x, y) ½ ¾ r Le vecteur {r} = x dont les composantes sont : ry  ∂ 2 u ∂σxx ∂σxy   − − fV x  rx = ρ 2 − ∂t ∂x ∂y 2    ry = ρ ∂ v − ∂σxy − ∂σyy − fV y ∂t2 ∂x ∂y

(2.1f)

(2.2)

½ ¾ u(x, y; t) est le r´ esidu de l’équation (2.1). Il est nul si le champ de déplacements {u(x, y; t)} = est v(x, y; t) solution de cette équation.

3

Forme int´ egrale faible

Pour résoudre l’équation (2.1) par la méthode des éléments finis, nous utilisons la m´ ethode des r´ esidus pond´ er´ es. Multiplions le résidu ½ ¾ r {r} = x (3.1) ry par un champ de déplacements arbitraire ½ ∗¾ u {u } = v∗ ∗

(3.2)

puis intégrons sur le domaine V : Z Z ∗ ∗ T W(u, u ) = {u } {r} dV = (u∗ rx + v ∗ ry ) dV = 0 ∀ {u∗ } V

(3.3)

V

Intégrons par parties la quantité

Z µ V

∂σxx ∂σxy + ∂x ∂y

¶ u∗ dV

puis utilisons le théorème de la divergence (ou théorème d’Ostrogradski). Il vient : ¶ Z µ ∂σxx ∂σxy + u∗ dV ∂x ∂y V ¶ ¶ Z µ Z µ ∂(u∗ σxx ) ∂(u∗ σxy ) ∂u∗ ∂u∗ = + dV − σxx + σxy dV ∂x ∂y ∂x ∂y V V ¶ Z Z µ ∂u∗ ∂u∗ ∗ = u (σxx nx + σxy ny ) dS − σxx + σxy dV ∂x ∂y S V Imposons la condition u∗ = 0 sur Su . La première intégrale se réduit à : Z u∗ fσx dS

(3.4)

(3.5)

Sσ

De même :

Z V

µ

∂σyy ∂σxy + ∂x ∂y

¶

Z ∗

Z ∗

v dV = Sσ

v fσy dS −

V

µ ¶ ∂v ∗ ∂v ∗ σxy + σyy dV ∂x ∂y

(3.6)

6


En portant les expressions (3.4), (3.5) et (3.6) dans l’équation (3.3), vient : Z Z W(u, u∗ ) = ρ {u∗ }T {¨ u} dV + {ε∗ }T {σ} dV VZ VZ ∗ T − {u } {fV } dV − {u∗ }T {fσ } dS V

o` u l’on a posé :  ∗  εxx  {ε∗ } = ε∗yy  ∗  γxy

,

ε∗xx =

∂u∗ ∂x

,

(3.7)

Sσ

ε∗yy =

∂v ∗ , ∂y

∗ γxy =

∂u∗ ∂v ∗ + ∂y ∂x

,

{¨ u} =

∂2 {u} ∂t2

(3.8)

{ε∗ } est le champ de déformations virtuelles induit par le champ de déplacements virtuels {u∗ }. La forme int´ egrale faible d’un probl` eme d’´ elasticit´ e s’écrit finalement : Trouver le champ de déplacements ½ ¾ u(x, y; t) {u(x, y; t)} = v(x, y; t)

(3.9a)

tel que : Z ∗

W(u, u ) =

Z ∗ T

ρ {u } {¨ u} dV + {ε∗ }T {σ} dV V Z Z ∗ T − {u } {fV } dV − {u∗ }T {fσ } dS = 0 V

V

(3.9b)

Sσ

∀ {u∗ } tel que {u∗ } = {0} sur Su o` u:

  εxx  {ε} = εyy   γxy

,

  σxx  {σ} = σyy   σxy

(3.9c)

avec : – les relations cinématiques : εxx =

∂u ∂x

,

εyy =

∂v ∂y

,

γxy =

∂u ∂v + ∂y ∂x

(3.9d)

– la loi de comportement (ou loi constitutive) : {σ} = [D] ({ε} − {εth })

(3.9e)

{u} = {uP } sur Su

(3.9f)

– les conditions aux limites : – les conditions initiales : {u(x, y; t0 )} = {u0 (x, y)} ,

{u(x, ˙ y; t0 )} = {u˙ 0 (x, y)} ,

{u} ˙ =

∂ {u} ∂t

(3.9g)


7

Remarques : – Les fonctions {u} et {u∗ } doivent être suffisamment régulières pour que les expressions ci-dessus aient un sens. – La fonction {u∗ } est appelée champ de d´ eplacements virtuels. – Le champ de déplacements {u} est dit cin´ ematiquement admissible (CA). – La formulation intégrale (3.9) est l’expression du principe des travaux virtuels. – Dans l’équation (3.3) la fonction {u}] doit être dérivable deux fois et une fois dans l’équation egrale forte et forme int´ egrale (3.9). Ces équations sont dites respectivement forme int´ faible de l’équation différentielle (2.1). – Sous certaines conditions de régularité, les formulations (2.1) et (3.9) sont équivalentes.

4

Forme discr´ etis´ ee : ´ el´ ements finis

La solution analytique de l’équation (3.9) est en général inaccessible. On est donc conduit à chercher une solution approchée par une méthode numérique : la méthode des éléments finis. Cette méthode est un cas particulier de la m´ ethode de Galerkin : le champ de déplacements cherché {u} et les ∗ fonctions test {u } appartiennent au même espace de dimension finie.

4.1

Discr´ etisation du domaine : maillage

Le domaine V est décomposé en sous-domaines V e de forme géométrique simple (les ´ el´ ements) reliés entre eux en des points appelés nœuds (figure 4). Cette opération s’appelle maillage.

Figure 4 – Domaine plan discrétisé en 12 éléments (8 triangles, 4 quadrangles) reliés entre eux par 15 nœuds Le maillage est défini par deux tables : – La table des nœuds contient les coordonnées des nœuds. x1 y1

x2 y2

x3 y3

x4 y4

x5 y5

... ...

– La table des ´ el´ ements contient, pour chaque élément, le type (triangle à trois nœuds, quadrangle à quatre nœuds, . . . ) et les numéros des nœuds dans le sens trigonométrique. TRI3 1 4 2 —

QUAD4 4 5 3 2

... ... ... ... ...

8


4.2

Repr´ esentation ´ el´ ementaire (ou locale) du champ de d´ eplacements

Le champ déplacements dans chaque élément est défini en fonction des déplacements des nœuds de l’élément : l’approximation est dite nodale. Dans l’élément (e) (figure 5) :  e  u1 (t)        v1e (t)    ¸ ½ ¾ · e   e N1 (x, y) 0 . . . Nne (x, y) 0 u(x, y; t) . . (4.1) = .  0 N1e (x, y) . . . 0 Nnee (x, y)  v(x, y; t)   e   u e (t)      en vne (t) soit sous forme matricielle : {u(x, y; t)} = {N e (x, y)}T {ue (t)}

(4.2)

o` u – – – –

ne est le nombre de nœuds de l’élément. el´ ementaires. les fonctions Nie (x, y) sont les fonctions d’interpolation ´ la matrice [N e (x, y)] est la matrice d’interpolation élémentaire. le vecteur {ue (t)} regroupe les composantes des déplacements des nœuds de l’élément (e).

Figure 5 – Champ de déplacements dans un élément ` a trois nœuds

4.3

Repr´ esentation globale du champ de d´ eplacements

Le champ de déplacements a pour expression sur l’ensemble du domaine V :   u1 (t)        v (t)   ½ ¾ · ¸ 1   u(x, y; t) N1 (x, y) 0 . . . Nn (x, y) 0 . . = . v(x, y; t) 0 N1 (x, y) . . . 0 Nn (x, y)      un (t)       vn (t)

(4.3)

soit sous forme matricielle : {u(x, y; t)} = [N (x, y)] {U (t)} o` u – – – –

n est le nombre de nœuds du maillage. les fonctions Ni (x, y) sont les fonctions d’interpolation (ou fonctions de forme). [N (x, y)] est la matrice d’interpolation. {U (t)} est le vecteur des d´ eplacements nodaux.

(4.4)


4.4

9

Discr´ etisation de la forme int´ egrale faible

De l’expression du champ de déplacements sur le domaine : {u(x, y; t)} = [N (x, y)] {U (t)} on déduit :

(4.5)

¨} {¨ u} = [N ] {U

(4.6a)  ∂N

{ε} = [B] {U } avec [B] =

£

B1 . . . Bi . . . Bn

¤

,

i

 ∂x   Bi =   0   ∂N

i

∂y {u∗ } = [N ] {U ∗ } , ∗

∗

{ε } = [B] {U } ,

 0

  ∂Ni   ∂y   ∂Ni  ∂x

{u∗ }T = {U ∗ }T [N ]T ∗ T

∗ T

T

{ε } = {U } [B]

En portant ces relations dans l’équation (3.9b), il vient : ³ ´ ¨ } + [K] {U } − {F } W({U }, {U ∗ }) = {U ∗ }T [M ] {U o` u:

(4.6b)

(4.6c) (4.6d)

(4.7)

Z [M ] = Z [K] = Z Z {F } = [N ]T {fV } dV + V

ρ [N ]T [N ] dV

(4.8)

[B]T [D] [B] dV

(4.9)

V

V

Sσ

Z [N ]T {fS } dS +

V

[B]T [D] {εth } dV

(4.10)

[M ] est la matrice de masse (kg). [K] est la matrice de rigidit´ e (N/m). {F } est le vecteur force ´ equivalent aux charges r´ eparties (N). {U } est le vecteur des d´ eplacements nodaux (m). ¨ } est le vecteur des acc´ {U el´ erations nodales (m/s2 ). Remarques : – les matrices [M ] et [K] sont par construction symétriques (car la matrice des coefficients élastiques [D] est symétrique). – dans l’équation (4.7), il convient d’ajouter la contribution de l’amortissement : {U ∗ }T [C] {U˙ } o` u [C] est la matrice d’amortissement (kg/s) et {U˙ } le vecteur des vitesses nodales (m/s). Effectuons une partition des degr´ es de libert´ e en déplacements inconnus (L) et imposés (non nuls : P , nuls : S) ([1], [14], [15]) :      {UL∗ }   {UL } = ?  (4.11a) d’o` u {U ∗ } = {UP∗ } = {0} = {δU } {U } = {UP } 6= {0}   ∗   {US } = {0} {US } = {0}

10


{δU } est une variation quelconque du vecteur {U }. Cette partition induit une partition de [M ], [C], [K] et {F } : 

 [MLL ] [MLP ] [MLS ] [M ] = [MP L ] [MP P ] [MP S ] [MSL ] [MSP ] [MSS ]



,

 [KLL ] [KLP ] [KLS ] [K] = [KP L ] [KP P ] [KP S ] [KSL ] [KSP ] [KSS ]



 [CLL ] [CLP ] [CLS ] [C] = [CP L ] [CP P ] [CP S ] [CSL ] [CSP ] [CSS ]

,

   {FL }  {F } = {FP }   {FS }

(4.11b)

(4.11c)

La forme faible discrétisée s’écrit finalement : Trouver {UL (t)} tel que : W({UL }, {UL∗ })

={UL∗ }T

µ ½ ¾ ½ ¾ ¨L } £ ¤ {U £ ¤ {U˙ L } [MLL ] [MLP ] ¨P } + [CLL ] [CLP ] {U˙ P } {U ½ ¾ ¶ £ ¤ {UL } + [KLL ] [KLP ] − {FL } = 0 ∀ {UL∗ } {UP }

(4.12)

avec les conditions initiales {UL (t0 )} = {UL,0 } , {U˙ L (t0 )} = {U˙ L,0 } Les déplacements nodaux inconnus {UL (t)} sont donc les solutions de l’équation : ¨L } + [CLL ]{U˙ L } + [KLL ]{UL } [MLL ]{U ¨P } − [CLP ]{U˙ P } − [KLP ]{UP } = {FL } − [MLP ]{U

(4.13a)

avec les conditions initiales : {UL (t0 )} = {UL,0 } ,

{U˙ L (t0 )} = {U˙ L,0 }

(4.13b)

Remarque : par construction, les matrices [KLL ] et [MLL ] sont symétriques.

4.5 4.5.1

Probl` emes particuliers Probl` eme stationnaire

Dans un problème stationnaire, l’équation (4.13) se réduit à : [KLL ]{UL } = {FL } − [KLP ]{UP } = {F¯L }

(4.14)

Si le nombre de liaisons est suffisant, la matrice [KLL ] n’est pas singulière (det [KLL ] 6= 0) et les déplacements inconnus sont égaux à : {UL } = [KLL ]−1 {F¯L } Les déplacements étant connus, les actions de liaison sont égales à : · ¸½ ¾ ½ ¾ [KP L ] [KP P ] {UL } {FP } {A} = − [KSL ] [KSP ] {UP } {FS }

(4.15)

(4.16)

´ Elasticit´ e plane 4.5.2

11

Modes propres de vibration

Les modes propres de vibration de la structure sont les solutions de l’équation : ¨L } + [KLL ]{UL } = 0 [MLL ]{U

(4.17)

˜L } sin ω t {UL (t)} = {U

(4.18)

En posant : ˜L } est indépendant du temps, il vient : o` u {U ˜L } = ω 2 [MLL ]{U ˜L } [KLL ]{U

(4.19)

˜L } le vecteur propre associé. o` u ω est une pulsation propre de la structure et {U Les pulsations propres sont les solution de l’équation : ¡ ¢ det [KLL ] − ω 2 [MLL ] = 0

4.6

(4.20)

Mise en œuvre pratique : calculs ´ el´ ementaires et assemblage

Dans la pratique, [M ], [K] et {F } sont construits ´ el´ ement par ´ el´ ement. Cette opération s’appelle assemblage. De l’expression du champ de déplacements dans l’élément (e) : {u(x, y; t)} = [N e (x, y)] {ue (t)}

(4.21)

on déduit : {u∗ } = [N e ] {ue∗ } ,

{u∗ }T = {ue∗ }T [N e ]T

{¨ u} = [N e ] {¨ ue }

 ∂N e 1

{ε} = [B e ] {ue } ,

[B e ] =

£

B1e . . . Bie . . . Bne e

¤

 ∂x   e Bi =   0   ∂N e

,

1

∂y {ε∗ } = [B e ] {ue∗ } ,

 0

  ∂N1e   ∂y   ∂N1e  ∂x

(4.22)

{ε∗ }T = {ue∗ }T [B e ]T

En portant ces expressions dans l’équation (3.9b), il vient : W({U }, {U ∗ }) =

X {ue∗ }T ( [me ] {¨ ue } + [k e ] {ue } − {f e } )

(4.23)

e

o` u:

Z [me ] = V

Z e

[k ] =

(4.24)

[B e ]T [D] [B e ] dV

(4.25)

Ve

Z

Z

e

e T

{f } = V

ρ [N e ]T [N e ] dV e

e

[N ] {fV } dV +

Z e T

Sσe

[N ] {fS } dS +

V

e

[B e ]T [D] {εeth } dV

(4.26)

12


´ ement e Figure 6 – El´ Dans ces formules, V e représente le volume de l’élément (e) et Sσe la partie de Sσ qui appartient à la frontière de l’élément (e) (figure 6). Les matrices et les vecteurs élémentaires sont évaluées numériquement. L’équation (3.9b) s’écrit : W({U }, {U ∗ }) =

X

{U ∗ }T

e

Ã ∗ T

= {U }

³

¨ } + [K e ] {U } − {F e } [M e ] {U

X

¨} + [M ] {U e

X

e

d’o` u: [M ] =

X

[M e ] ,

[K] =

e

e

[K ] {U } −

e

X

[K e ] ,

´

X

!

(4.27)

e

{F }

e

{F } =

X

e

[F e ]

(4.28)

e

Dans les matrices [M e ] et [K e ] et dans le vecteur {F e }, obtenus par expansion respectivement de [me ], [k e ] et {f e }, les seuls termes non nuls sont les termes associés aux degrés de liberté de l’élément (e). Remarques : – La partition des degrés de liberté est effectuée avant la phase d’assemblage. – Dans le logiciel RDM seuls les blocs de matrice (LL) et (LP ) sont assemblés.

5

Probl` eme ´ elastostatique : ´ energie potentielle et m´ ethode de Ritz

Si le problème est indépendant du temps, la forme intégrale faible (3.9) se réduit à : Z Z ∗ ∗ T W(u, u ) = {ε } [D] {ε} dV − {ε∗ }T [D] {εth } dV VZ ZV ∗ T − {u } {fV } dV − {u∗ }T {fσ } dS = 0 V

(5.1)

Sσ

∀ {u∗ } tel que {u∗ } = {0} sur Su

5.1

Calcul des variations

Le problème fondamental du calcul des variations consiste à chercher la fonction u(x) qui rend stationnaire la fonctionnelle (ou fonction de fonctions ) : Z J (u) = a

b

µ ¶ ∂u ∂nu F x, u, , . . . , n dx ∂x ∂x

(5.2)


13

ce qui s’écrit : δJ = 0 ∀ δu

(5.3)

Les principales propriétés de l’opérateur variation δ sont ([2, 9, 17]) :  2 δ (u) = δ(δu) = 0    µ ¶    ∂(δu) ∂u   = δ   ∂x ∂x   µ ¶    ∂u ∂u ∂F ∂F   δF (u, , . . .) = δu + µ ¶ δ + ···   ∂u ∂x ∂u ∂x    ∂   ∂x    δ(F + G) = δF + δG                             

5.2

∂u ∂u , . . .) et G(u, , . . .) sont deux fonctionnelles de u ∂x ∂x δ(F G) = δF G + F δG (règle de Leibniz)

(5.4)

o` u F (u,

δ(F n ) = n F n−1 δF δ(c F ) = c δF o` u c est une constante Z Z δ F dx = δF dx

´ Energie potentielle

Considérons la fonctionnelle : Epot ({u}) = Edef ({u}) − Wext ({u})

(5.5)

o` u: – {u} est un champ de déplacements cinématiquement admissible. – Edef ({u}) est l’énergie de déformation du champ de déplacements {u} : Z Z 1 T Edef ({u}) = {ε} [D] {ε} dV − {ε}T [D] {εth } dV 2 V V – Wext ({u}) est travail des forces appliquées pour le déplacement {u} : Z Z T Wext ({u}) = {u} {fV } dV + {u}T {fσ } dS V

(5.6)

(5.7)

Sσ

– Epot ({u}) est l’énergie potentielle du système pour le déplacement {u}. La condition de stationnarité (5.3) s’écrit : δEpot = δEdef − δWext = 0 ∀ {δu} soit :

Z δEpot =

Z T

{δε} [D] {ε} dV − {δε}T [D] {εth } dV V Z ZV − {δu}T {fV } dV − {δu}T {fσ } dS = 0 ∀ {δu} V

o` u

(5.8)

Sσ

(5.9)

14

Méthode des éléments finis – {δu} est une variation quelconque du champ de déplacements (en particulier : {δu} = {0} sur Su ).   µ ¶   ∂u ∂(δu)          δ            ∂x   ∂x         µ ¶   δεxx    ∂v ∂(δv) = – {δε} = δεyy = δ       ∂y ∂y    δγxy    ¶  µ         ∂u ∂(δu) ∂v ∂(δv)         + + δ    ∂y ∂x ∂y ∂x

Cette équation est identique à (5.1) si on choisit {u∗ } = {δu}. La seconde variation de la fonctionnelle est égale à : Z δ 2 Epot = {δε}T [D] {δε} dV

(5.10)

V

La matrice [D] étant définie positive, on en déduit : δ 2 Epot > 0 ∀ {δu} <> {0}

(5.11)

De plus Epot ({uexact } + {δu}) = Epot ({uexact }) + δEpot |{u}={uexact } + = Epot ({uexact }) +

1 2 δ Epot 2

1 2 δ Epot > Epot ({uexact }) ∀ {δu} <> {0} 2

(5.12)

Le champ de déplacements {u} = {uexact } + {δu} étant cinématiquement admissible (CA), on en déduit : Epot ({uCA }) ≥ Epot ({uexact }) ∀ {uCA } (5.13) On peut donc énoncer le théorème suivant : Parmi l’ensemble des champs de d´ eplacements cin´ ematiquement admissibles, le champ de d´ eplacements exact est celui qui minimise l’´ energie potentielle.

5.3

M´ ethode de Ritz et ´ el´ ements finis

Si on restreint la recherche de la solution aux champs de déplacements définis au paragraphe (4.3), l’énergie potentielle discrétisée est égale à : 1 {U }T [K]{U } − {U }T {F } 2 ½ ¾T · ¸½ ¾ ½ ¾T ½ ¾ 1 {UL } [KLL ] [KLP ] {UL } {UL } {FL } = − [KP L ] [KP P ] {UP } {UP } {FP } 2 {UP }

Epot ({U }) =

et la condition de stationnarité s’écrit : µ ½ ¾ ¶ £ ¤ UL T [KLL ] [KLP ] δEpot ({UL }) = {δUL } − {FL } = 0 ∀ {δUL } UP

(5.14)

(5.15)

d’o` u: [KLL ]{UL } = {FL } − [KLP ]{UP } Cette équation est identique à celle obtenue au paragraphe (4.5.1).

(5.16)


6

15

Calculs ´ el´ ementaires : ´ el´ ements isoparam´ etriques

6.1

´ ement isoparam´ El´ etrique : d´ efinition

` chaque élément réel, on associe un ´ A el´ ement de r´ ef´ erence (figure 7).

Figure 7 – Transformation géométrique 6.1.1

Repr´ esentation de la g´ eom´ etrie

La transformation géométrique (figure 7) qui fait passer de l’élément de référence à l’élément réel possède les propriétés suivantes : – elle est de la forme : x(ξ, η) =

n X

Ni (ξ, η) xi

,

y(ξ, η) =

i=1

n X

Ni (ξ, η) yi

(6.1)

i=1

o` u: – – – – –

n est le nombre de nœuds de l’élément. ξ et η sont les coordonnées d’un point de l’élément de référence. x(ξ, η) et y(ξ, η) sont les coordonnées d’un point de l’élément réel. xi et yi sont les coordonnées du ie nœud de l’élément. les Ni (ξ, η) sont les fonctions d’interpolation ou fonctions de forme

La matrice jacobienne de la  ∂x  ∂ξ [J(ξ, η)] =   ∂x ∂η

transformation est :  P  n ∂N n ∂N P ∂y i i xi yi   ∂ξ  i=1 ∂ξ   = i=1 ∂ξ n n   P P ∂N ∂N ∂y i i  xi yi ∂η i=1 ∂η i=1 ∂η     x 1 y1 ∂N1 ∂Ni ∂Nn  .. ..   · ¸ .   ∂ξ . . . ∂ξ . . . ∂ξ   . J11 J12     xi yi  = = J21 J22 ∂Ni ∂Nn   ∂N1  .. ..  ... ...   . . ∂η ∂η ∂η xn yn

(6.2)

– elle est nodale : un nœud de l’élément de référence devient un nœud de l’élément réel (les deux éléments possèdent donc le même nombre de nœuds) : xi = x(ξi , ηi ) =

n X j=1

Nj (ξi , ηi ) xj

,

yi = y(ξi , ηi ) =

n X j=1

Nj (ξi , ηi ) yj

,

i = 1, . . . , n

(6.3)

16

Méthode des éléments finis o` u (ξi , ηi ) sont les coordonnées du ie nœud de l’élément de référence. On en déduit : ( 0 si i 6= j Nj (ξi , ηi ) = 1 si i = j

(6.4)

– une frontière de l’élément de référence devient une frontière de l’élément réel. – elle est bijective : le d´ eterminant de la matrice jacobienne ne doit pas changer de signe sur l’´ el´ ement. Nous imposerons la condition : det [J(ξ, η)] = J11 J22 − J12 J21 > 0

(6.5)

ce qui implique que l’élément de référence et l’élément réel soient numérotés dans le même sens (en général positif).

6.1.2

Maillage conforme

La transformation géométrique doit assurer la continuité de la géométrie entre les éléments (figure 8).

Figure 8 – Maillage conforme (` a gauche) et non conforme (` a droite] Si deux éléments ont une arête commune : – les deux éléments doivent avoir le même nombre de nœuds sur l’arête. – les coordonnées d’un point de l’arête ne doivent dépendre que des coordonnées des noeuds de l’arête ; elles se réduisent à : x(α) =

na X

Nia (α) xi

,

y(α) =

i=1

na X

Nia (α) yi

,

−1 ≤ α ≤ 1

(6.6)

i=1

o` u na est le nombre de nœuds situés sur l’arête et les Nia (α) sont les fonctions d’interpolation de l’élément à une dimension et na nœuds (figure 15) et (§ 6.2.2 : transformation d’une arête).

6.1.3

Repr´ esentation du champ de d´ eplacements

Les fonctions Ni (ξ, η) qui définissent la transformation géométrique sont les fonctions d’interpolation pour le champ de déplacements (élément isoparamétrique) : u(ξ, η) =

n X

Ni (ξ, η) ui

i=1

,

v(ξ, η) =

n X

Ni (ξ, η) vi

(6.7)

i=1

o` u ui et vi sont les déplacements du nœud i. Critère de complétude : pour que la solution éléments finis converge vers la solution exacte quand la taille des éléments tend vers zéro, l’élément doit pouvoir représenter un champ de déplacements qui


17

correspond à des déformations nulles (mouvement de corps rigide) ou constantes. Considérons donc le champ de déplacements : u(x, y) = ax + bx x + cx y

,

v(x, y) = ay + by x + cy y

(6.8)

d’o` u les valeurs nodales : ui = ax + bx xi + cx yi

,

vi = ay + by xi + cy yi

,

i = 1, . . . , n

(6.9)

Le champ de déplacements sur x s’écrit sous forme paramétrique (équation 6.7) : u(ξ, η) =

n X

Ni (ξ, η) ui =

i=1

= ax

n X

Ni (ξ, η) (ax + bx xi + cx yi )

i=1

n X

Ni (ξ, η) + bx

i=1

n X

Ni (ξ, η) xi + cx

i=1

n X

(6.10) Ni (ξ, η) yi

i=1

En utilisant les relations (6.1), il vient : u(x, y) = ax

n X

Ni (ξ, η) + bx x + cx y

(6.11)

i=1

On retrouve le champ de déplacements (6.8) si : n X

Ni (ξ, η) = 1

(6.12)

i=1

Cette condition est vérifiée par les éléments décrits ci-dessous. Remarque : si le maillage est conforme, le champ de déplacement est continu entre les éléments.

6.2 6.2.1

Biblioth` eque d’´ el´ ements Triangle ` a 3 nœuds

Figure 9 – Triangle ` a 3 nœuds La transformation géométrique est de la forme :   ¤ a b = [P (ξ, η)] {A} x(ξ, η) = a + b ξ + c η = 1 ξ η   c £

[P (ξ, η)] est la base polynomiale de la transformation.

(6.13)

18


La transformation est nodale d’o` u:         1 0 0 x1  x(0, 0) a   x2 = x(1, 0) = 1 1 0 b = [C] {A}       x3 x(0, 1) 1 0 1 c

(6.14)

On en déduit : x(ξ, η) = [P ] [C]−1

    x1  £ ¤ x1  x2 = N1 (ξ, η) N2 (ξ, η) N3 (ξ, η) x2     x3 x3

Il vient pour les fonctions d’interpolation et leurs dérivées (programme tri3 int) :       · ¸T · ¸T −1 −1 1−ξ−η ∂N ∂N  , ξ [N ]T =  = 1  , = 0  ∂ξ ∂η 0 1 η

(6.15)

(6.16)

Figure 10 – Triangle ` a 3 nœuds : fonctions d’interpolation Matrice jacobienne de la transformation  · ¸ x −1 1 0  1 x2 [J(ξ, η)] = −1 0 1 x3

: elle est égale à (équation (6.2)) :  · ¸ y1 x2 − x1 y2 − y1  y2 = x3 − x1 y3 − y1 y3

(6.17)

d’o` u son déterminant : det[J(ξ, η)] = (x2 − x1 ) (y3 − y1 ) − (x3 − x1 ) (y2 − y1 ) = 2 A o` u A est l’aire de l’élément réel.

6.2.2

Triangle ` a 6 nœuds

Figure 11 – Triangle ` a 6 nœuds

(6.18)


19

Base polynomiale :

£ ¤ [P (ξ, η)] = 1 ξ η ξ 2 η 2 ξ η

(6.19)

Fonctions d’interpolation (λ = 1 − ξ − η) (programme tri6 int) :       λ(2 λ − 1) 1 − 4λ 1 − 4λ  4 ξλ        · ¸T 4 (λ − ξ) · ¸T  −4 ξ        ∂N ∂N ξ(2 ξ − 1)  0  −1 + 4 ξ  ,   [N ]T =  , = =  4 ξη   4η   4ξ  ∂ξ ∂η        η(2 η − 1)     −1 + 4 η  0 4 ηλ −4 η 4 (λ − η)

(6.20)

Transformation d’une arˆ ete : considérons l’arête de l’élément de référence passant par les nœuds 3, 4 et 5 ; elle est définie par λ = 1 − ξ − η = 0 d’o` u η = 1 − ξ avec 0 ≤ ξ ≤ 1 ; elle devient dans l’élément réel la courbe d’équation paramétrique :   x1        x2       £ ¤ x3  x(ξ) = 0 0 ξ(2 ξ − 1) 4 ξ(1 − ξ) (1 − ξ)(1 − 2 ξ) 0 x4        x5        x6 = ξ(2 ξ − 1) x3 + 4 ξ(1 − ξ) x4 + (1 − ξ) (1 − 2ξ) x5 y(ξ) = ξ(2 ξ − 1) y3 + 4 ξ(1 − ξ) y4 + (1 − ξ) (1 − 2ξ) y5

,

0≤ξ≤1

d’o` u en posant ξ = (1 − α)/2 : x(α) =

α (α − 1) α (α + 1) ¯1 (α) x3 + N ¯2 (α) x4 + N ¯3 (α) x5 x3 + (1 − α2 ) x4 + x5 = N 2 2

¯1 (α) y3 + N ¯2 (α) y4 + N ¯3 (α) y5 y(α) = N

,

−1 ≤ α ≤ 1

¯i (α) sont les fonctions d’interpolation de l’élément isoparamétrique à une dimension et trois nœuds. Les N Si le nœud 4 est au milieu des nœuds 3 et 5, la transformation se réduit à : x(α) =

1−α 1+α x3 + x5 2 2

,

y(α) =

1−α 1+α y3 + y5 2 2

Remarque : le champ de déplacements des points de l’arête est : ¯1 (α) u3 + N ¯2 (α) u4 + N ¯3 (α) u5 u(α) = N

,

¯1 (α) v3 + N ¯2 (α) v4 + N ¯3 (α) v5 v(α) = N

Il ne dépend que des déplacements des nœuds de l’arête (élément conforme).

6.2.3

Quadrangle ` a 4 nœuds

Figure 12 – Quadrangle ` a 4 nœuds

20


Base polynomiale :

£ ¤ [P (ξ, η)] = 1 ξ η ξ η

(6.21)

Fonctions d’interpolation (programme quad4 int) :   (1 − ξ) (1 − η) 1 (1 + ξ) (1 − η)  [N ]T =  4 (1 + ξ) (1 + η) (1 − ξ) (1 + η) 6.2.4

·

∂N ∂ξ

¸T

  −(1 − η) 1  (1 − η)   =  4  (1 + η)  −(1 + η)

·

∂N ∂η

¸T

  −(1 − ξ) 1 −(1 + ξ)  =  4  (1 + ξ)  (1 − ξ)

(6.22)

Quadrangle ` a 8 ou 9 nœuds

Quadrangle ` a 8 nœuds :

Figure 13 – Quadrangle ` a 8 ou 9 nœuds Base polynomiale :

£ ¤ [P (ξ, η)] = 1 ξ η ξ 2 ξη η 2 ξ 2 η ξη 2

(6.23)

Fonctions d’interpolation (programme quad8 int) :   −(1 − ξ)(1 − η)(1 + ξ + η)   2 (1 − ξ 2 )(1 − η)   −(1 + ξ)(1 − η)(1 − ξ + η)    1 2 (1 + ξ)(1 − η 2 ) T   [N ] =  4 −(1 + ξ)(1 + η)(1 − ξ − η)    2 (1 − ξ 2 )(1 + η)   −(1 − ξ)(1 + η)(1 + ξ − η) 2 (1 − ξ)(1 − η 2 )   −(1 − η)(2 ξ + η)  −4 (1 − η) ξ     (1 − η)(2 ξ − η)    · ¸  ∂N T 1  2 (1 − η 2 )   =  ∂ξ 4  (1 + ξ)(2 ξ + η)    −4 (1 + η) ξ     (1 + η)(2 ξ − η)  2 (1 − η 2 )

 (1 − ξ)(ξ + 2 η)  −2 (1 − ξ 2 )    −(1 + ξ)(ξ − 2 η)   · ¸ ∂N T 1  −4 (1 + ξ) η    =  ∂η 4  (1 + ξ)(ξ + 2 η)   2)   2 (1 − ξ   −(1 − ξ)(ξ − 2 η) −4 (1 − ξ) η

(6.24)



,

(6.25)

Quadrangle ` a 9 nœuds : Base polynomiale :

£ [P (ξ, η)] = . . . ξ i η j

¤ ...

i, j = 0, . . . , 2

(6.26)


21

Fonctions d’interpolation (programme quad9 int) :     DL1 (ξ) L1 (η) L1 (ξ) L1 (η) L2 (ξ) L1 (η) DL2 (ξ) L1 (η)     DL3 (ξ) L1 (η) L3 (ξ) L1 (η)     L3 (ξ) L2 (η) ¸T DL3 (ξ) L2 (η) ·     ∂N  DL3 (ξ) L3 (η) = [N ]T =  L3 (ξ) L3 (η) ,   ∂ξ L2 (ξ) L3 (η) DL2 (ξ) L3 (η)     DL1 (ξ) L3 (η) L1 (ξ) L3 (η)     DL1 (ξ) L2 (η) L1 (ξ) L2 (η)

·

ξ (ξ − 1) [L(ξ)] = 2

6.3 6.3.1

,

DL2 (ξ) L2 (η)

L2 (ξ) L2 (η) o` u:

  L1 (ξ) DL1 (η) L2 (ξ) DL1 (η)   L3 (ξ) DL1 (η)   · ¸T L3 (ξ) DL2 (η)   ∂N L3 (ξ) DL3 (η) =   ∂η L2 (ξ) DL3 (η)   L1 (ξ) DL3 (η)   L1 (ξ) DL2 (η)

1 − ξ2

ξ (ξ + 1) 2

¸

L2 (ξ) DL2 (η) ·

,

2ξ − 1 [DL(ξ)] = 2

−2 ξ

2ξ + 1 2

¸ (6.28)

Calcul des matrices et des vecteurs ´ el´ ementaires Transformation des d´ eriv´ ees

Les dérivées d’une fonction f (x, y) par rapport à ξ et η      ∂f   ∂f ∂x ∂f ∂y  ∂x    +        ∂ξ ∂x ∂ξ ∂y ∂ξ ∂ξ = =  ∂x ∂x ∂f ∂y ∂f ∂f             + ∂η ∂x ∂η ∂y ∂η ∂η

sont :   ∂y  ∂f      ∂ξ   ∂x = [J] ∂y   ∂f      ∂η ∂y

  ∂f       ∂x ∂f       ∂y

On en déduit l’expression des dérivées de f par rapport à x et y :     ∂f  ∂f            ∂x ∂ξ = [J]−1 ∂f  ∂f            ∂y ∂η avec −1

[J] 6.3.2

(6.27)

· ¸ 1 J22 −J12 = det[J] −J21 J11

,

det[J] = J11 J22 − J12 J21

(6.29)

(6.30)

(6.31)

Transformation des int´ egrales

La surface infinitésimale dξ dη au point (ξ,η) de l’élément de référence se transforme en la surface dA au point (x(ξ, η),y(ξ, η)) de l’élément réel.

Figure 14 – Transformation des surfaces On a la relation :

~a ∧ ~b = ~k dA

(6.32)

22


soit

    ∂x  ∂x                 ∂ξ   ∂η  {a} ∧ {b} = ∂y dξ ∧ ∂y dη = {k} det[J(ξ, η)] dξ dη = {k} dA         ∂ξ            ∂η  0 0

(6.33)

dA = det[J(ξ, η)] dξ dη

(6.34)

d’o` u:

L’intégrale sur l’élément réel

Z f (x, y) dV

(6.35)

V

devient donc sur l’élément de référence Z Z t f (x, y) dA = A

t f (x(ξ, η), y(ξ, η)) det [J(ξ, η)] dξ dη

(6.36)

Aref

o` u t est l’épaisseur supposée constante de l’élément. Remarques : – on a la relation

Z det [J(ξ, η)] dξ dη = A

(6.37)

Aref

o` u A est l’aire de l’élément réel. – si l’épaisseur de l’élément est variable et donnée par ses valeurs nodales ti , on a : t(ξ, η) =

n X

Ni (ξ, η) ti

(6.38)

i=1

6.3.3

´ Evaluation num´ erique des int´ egrales

Les intégrales sont évaluées par intégration numérique ([2, 3, 9, 15]) : Z

Z

1

1

f (ξ, η) dξ dη ≈ −1

Z

−1

npi X

1 Z 1−ξ

f (ξ, η) dξ dη ≈ 0

0

wi f (ξi , ηi )

(quadrangle)

(6.39)

i=1 npi X

wi f (ξi , ηi )

i=1

o` u: – npi est le nombre de points d’intégration. – ξi et ηi sont les coordonnées du ie point d’intégration. – wi est le poids du ie point d’intégration.

(triangle)

(6.40)

´ Elasticit´ e plane 6.3.4

23

Calcul des matrices

La matrice de rigidité est égale à : Z Z T [k] = [B] [D] [B] dV = V

≈

t [B(ξ, η)]T [D] [B(ξ, η)] det [J(ξ, η)] dξ dη

Aref

npi X

(6.41)

T

t [B(ξi , ηi )] [D] [B(ξi , ηi )] det [J(ξi , ηi )] wi

i=1

o` u:  ∂N

i

 ∂x   . . . Bn ] avec Bi =   0   ∂N

[B] = [ B1 . . . Bi

i

∂y De même, la matrice de masse s’écrit : Z Z T [m] = ρ [N ] [N ] dV = V

≈

 0

  ∂Ni   ∂y   ∂N 

,

i

    ∂Ni  ∂Ni            ∂x ∂ξ −1 = [J] ∂N  ∂N       i   i  ∂y ∂η

∂x

t ρ [N (ξ, η)]T [N (ξ, η)] det [J(ξ, η)] dξ dη

Aref

npi X

(6.42)

(6.43)

T

t ρ [N (ξi , ηi )] [N (ξi , ηi )] det [J(ξi , ηi )] wi

i=1

o` u:

6.3.5

·

N1 0 . . . Ni 0 . . . Nn 0 [N ] = 0 N1 . . . 0 Ni . . . 0 Nn

¸ (6.44)

Calcul des vecteurs

Force volumique Le vecteur force dˆ u à une force volumique est égal à : Z Z T {f } = [N ] {fV } dV = t [N (ξ, η)]T {fV } det [J(ξ, η)] dξ dη V

≈

npi X i=1

Aref

½ ¾ fV x T t [N (ξi , ηi )] det [J(ξi , ηi )] wi fV y

(6.45)

Gradient thermique Le vecteur force dˆ u à gradient thermique est égal à : Z Z T {fth } = [B] [D] {εth } dV = t [B(ξ, η)]T [D] {εth } det [J(ξ, η)] dξ dη V

≈

npi X i=1

Force surfacique

Aref

T

t [B(ξi , ηi )] [D] {εth } det [J(ξi , ηi )] wi

(6.46)

24


Le champ de déplacements des points situés sur la frontière d’un élément est égal à (élément conforme (§ 6.1.2)) : x(ξ) =

n X

¯i (ξ) xi N

,

y(ξ) =

i=1

u(ξ) =

n X

n X

¯i (ξ) yi N

(6.47a)

¯i (ξ) vi N

(6.47b)

i=1

¯i (ξ) ui N

,

v(ξ) =

i=1

n X i=1

¯i (ξ) sont les fonctions d’interpolation de o` u n est le nombre de nœuds situés sur la frontière et les N l’élément de référence à une dimension et n nœuds (figure 15).

Figure 15 – Transformation des frontières On en déduit : dx =

Ã n X ∂N ¯i i=1

∂ξ

! xi

dξ = Jx dξ

,

dy =

Ã n X ∂N ¯i i=1

ds =

p

∂ξ

! yi

dξ = Jy dξ

(6.48)

q dx2 + dy 2 =

Jx2 + Jy2 dξ = Js dξ

(6.49)

d’o` u l’expression des composantes des vecteurs ~n et ~t (~n est la normale unitaire à Sσ dirigée vers l’extérieur de l’élément ; ~t est le vecteur unitaire tangent à Sσ ) :     dy  dx        ½ ¾ ½ ¾     1 1 Jx Jy ds ds = , {n} = = (6.50) {t} =   dy  Js Jy dx  Js −Jx    −     ds ds Le vecteur force dˆ u à une force surfacique d’intensité f~σ appliquée sur l’une des frontières d’un élément est égal à : Z Z Z 1 T ¯ ¯ (ξ)]T {fσ (ξ)} Js (ξ) dξ {f } = [N ] {fσ } dS = [N ] {fσ } dS = t [N ≈

Sσ npi X

Sσ

−1

(6.51)

¯ (ξi )] {fσ (ξi )} Js (ξi ) wi t [N T

i=1

o` u:

· ¸ ¯1 0 . . . N ¯i 0 . . . N ¯n 0 N ¯ [N ] = ¯1 . . . 0 N ¯i . . . 0 N ¯n 0 N

(6.52)

Pour une force perpendiculaire à la surface d’intensité p : f~σ = p ~n , l’expression ci-dessus se réduit à : ½ ¾ npi X Jy (ξi ) T ¯ {f } ≈ t [N (ξi )] p(ξi ) wi (6.53) −Jx (ξi ) i=1


6.4

25

Qualit´ e du jacobien.

La condition det [J] > 0 impose certaines conditions à la géométrie d’un élément. Exemple 1 : considérons le quadrangle à 4 nœuds représenté sur la figure (16).

Figure 16 – Quadrangle ` a 4 nœuds La matrice jacobienne est égale à (équations (6.2) et (6.22)) (programme § A.6) :

[J] =

1 4

=

1 4

  ¸ d d  −(1 − η) (1 − η) (1 + η) −(1 + η)  L 0   −(1 − ξ) −(1 + ξ) (1 + ξ) (1 − ξ) L L 0 L · ¸ 2 L − d (1 − η) −d (1 − η) −d (1 − ξ) 2 L − d (1 − ξ)

·

d’o` u 8 det [J] = 2 L2 − dL (2 − ξ − η) det [J] est minimal pour ξ = η = −1 : min(det [J]) = L (L − 2 d)/4 La condition det [J] > 0 dans tout l’élément impose donc : d < L/2. Remarque : si d = L/2, les nœuds 1, 2 et 4 sont alignés. Si d > L/2, la transformation qui fait passer de l’élément de référence à l’élément réel génère des points en dehors du quadrangle réel (figure 17).

Figure 17 – Quadrangle ` a 4 nœuds : d > L/2

Exemple 2 : considérons le triangle à 6 nœuds représenté sur la figure (18).

26


Figure 18 – Triangle ` a 6 nœuds ` a bords curvilignes La matrice jacobienne est égale à (équations (6.2) et (6.20)) (programme § A.7) : · ¸ L 4h(1 − 2ξ − η) [J] = 0 L − 4hξ d’o` u det [J] = L ( L − 4 h ξ ) La condition det [J] > 0 dans tout l’élément impose : h < L/4. Remarque : si h = L/4, la parabole qui passe par les nœuds 1, 2 et 3 est tangente à la droite qui passe par les nœuds 3, 4 et 5.

On appelle qualité du jacobien la quantité : qJ =

Aire de l’élément de référence min(det [J(ξ, η)]) Aire de l’élément réel

(6.54)

Remarques : – La qualité du jacobien est comprise entre 0 et 1. La qualité maximale est 1 : dans ce cas, le déterminant du jacobien est constant dans l’élément. – Dans la pratique, on se contente d’évaluer le déterminant du jacobien aux nœuds de l’élément. – D’autres définitions sont possibles. On rencontre souvent celle-ci : qJ =

min(det [J(ξ, η)]) max(det [J(ξ, η)])

(6.55)


A

Programmes Maple

Les programmes suivants sont dans le fichier map elas 2d.txt.

A.1

tri3 int : triangle ` a 3 nœuds

restart:with(linalg): n:=3; # base polynomiale P:=(xi,eta)->[1,xi,eta]; # coordonn´ ees nodales xi_nod:=[0,1,0]: eta_nod:=[0,0,1]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta);

A.2

tri6 int : triangle ` a 6 nœuds

restart:with(linalg): n:=6; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*xi,eta*eta,xi*eta]; # coordonn´ ees nodales xi_nod:=[0,1/2,1,1/2,0,0]: eta_nod:=[0,0,0,1/2,1,1/2]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta);

27

28


A.3

quad4 int : quadrangle ` a 4 nœuds

restart:with(linalg): n:=4; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*eta]; # coordonn´ ees nodales xi_nod:=[-1,1,1,-1]: eta_nod:=[-1,-1,1,1]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta);

A.4


restart:with(linalg): n:=8; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*eta,xi^2,eta^2,xi^2*eta,xi*eta^2]; # coordonn´ ees nodales xi_nod:=[-1,0,1,1,1,0,-1,-1]: eta_nod:=[-1,-1,-1,0,1,1,1,0]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta);


A.5


restart:with(linalg): n:=9; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*eta,xi^2,eta^2,xi^2*eta,xi*eta^2,xi^2*eta^2]; # coordonn´ ees nodales xi_nod:=[-1,0,1,1,1,0,-1,-1,0]: eta_nod:=[-1,-1,-1,0,1,1,1,0,0]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta);

A.6

Qualit´ e du jacobien : quadrangle ` a 4 nœuds

restart:with(linalg): n:=4; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*eta]; # coordonn´ ees nodales xi_nod:=[-1,1,1,-1]: eta_nod:=[-1,-1,1,1]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta); dN:=matrix([dNxi,dNeta]); coord:=matrix([[d,d],[L,0],[L,L],[0,L]]); J:=multiply(dN,coord):J:=simplify(%); detJ:=det(J);

29

30


A.7

Qualit´ e du jacobien : triangle ` a 6 nœuds

restart:with(linalg): n:=6; # base polynomiale P:=(xi,eta)->[1,xi,eta,xi*xi,eta*eta,xi*eta]; # coordonn´ ees nodales xi_nod:=[0,1/2,1,1/2,0,0]: eta_nod:=[0,0,0,1/2,1,1/2]: # fonctions d’interpolation C:=matrix([seq(P(xi_nod[i],eta_nod[i]),i=1..n)]); N:=multiply(P(xi,eta),inverse(C)); N:=map(factor,N); dNxi:=map(diff,N,xi); dNeta:=map(diff,N,eta); dN:=matrix([dNxi,dNeta]); coord:=matrix([[0,0],[L/2,h],[L,0],[L/2,L/2],[0,L],[0,L/2]]); J:=multiply(dN,coord):J:=simplify(%); detJ:=det(J);


31

R´ ef´ erences [1] J. H. Argyris et H.-P. Mlejnek – Die methode der finiten elemente, Band I. Verschiebungsmethode in der statik, Vieweg, 1986. [2] K.-J. Bathe – Finite element procedures in engineering analysis, Prentice Hall, 1996. [3] J.-L. Batoz et G. Dhatt – Modélisation des structures par éléments finis, Volume 1. Solides élastiques, Hermès, 1990. [4] M. Bonnet et A. Frangi – Analyse des solides déformables par la méthode des éléments finis, ´ ´ Editions de l’Ecole polytechnique, 2007. [5] L. Chevalier – Mécanique des systèmes et des milieux déformables. Cours, exercices et problèmes corrigés, Ellipses, 2004. [6] R. D. Cook, D. S. Malkus et M. E. Plesha – Concepts and applications of finite element analysis, 3 éd., Wiley, 1989. [7] M. A. Crisfield – Finite elements and solution procedures for structural analysis, Pineridge Press, 1986. [8] G. Dhatt et G. Touzot – Une présentation de la méthode des éléments finis, Maloine, 1984. [9] G. Dhatt, G. Touzot et E. Lefran¸ cois – Méthode des éléments finis, Hermès, 2005. [10] D. Euvrard – Résolution des équations aux dérivées partielles de la physique, de la mécanique et des sciences de l’ingénieur. Différences finies, éléments finis, problèmes en domaines non bornés, 3 éd., Masson, 1994. [11] F. Frey et J. Jirousek – Traité du génie civil, Volume 6. Méthode des éléments finis, Presses Polytechniques et Universitaires Romandes, 2001. [12] R. H. Gallagher – Introduction aux éléments finis, Pluralis, 1976. [13] L. Gallimard et J.-P. Pelle (éds.) – Estimateurs d’erreur pour les analyses éléments finis (revue européenne des éléments finis vol.12 no 6/2003), Hermès, 2003. [14] T. J. Hughes – The finite element method. Linear static and dynamic finite element analysis, Dover, 2000. [15] J.-F. Imbert – Analyse des structures par éléments finis, 3 éd., Cépaduès, 1995. `ze et J.-P. Pelle – La maˆıtrise du calcul en mécanique linéaire et non linéaire, [16] P. Ladeve Hermès, 2001. [17] A. Le Pourhiet – Résolution numérique des équations aux dérivées partielles. Une première approche, Cépaduès, 1988. [18] R. H. MacNeal – Finite elements. Their design and performance, Dekker, 1994. [19] N. Ottosen et H. Petersson – Introduction to the finite element method, Prentice Hall, 1992. [20] J.-P. Pelle, P. Beckers et L. Gallimard – Estimations des erreurs de discrétisation et analyses adaptatives. Application ` a l’automatisation des calculs éléments finis, Cours IPSI, 1996. [21] A. Portela et A. Charafi – Finite elements using Maple. A Symbolic Programming Approach, Springer, 2002. ´ et I. Babuˇ [22] B. Szabo ska – Finite element analysis, Wiley, 1991. ´ ements finis pour l’ingénieur. Grands principes et petites recettes, Tec & Doc [23] P. Thomas – El´ (Collection EDF R&D), 2006. [24] P. Trompette – Mécanique des structures par la méthode des éléments finis, Masson, 1992. [25] C. Wielgoz – Cours et exercices de résistance des matériaux : élasticité, plasticité, éléments finis, Ellipses, 1999.

32


[26] O. C. Zienkiewicz – La méthode des éléments finis appliquée ` a l’art de l’ingénieur, Ediscience, 1973. [27] — , The finite element method, McGraw-Hill, 1977. [28] O. C. Zienkiewicz et R. L. Taylor – La méthode des éléments finis. Formulation de base et problèmes linéaires, AFNOR, 1989. [29] — , The finite element method, Volume 1. The basis, Butterworth-Heinemann, 2000. [30] — , The finite element method, Volume 2. Solid mechanics, Butterworth-Heinemann, 2000.

Methodes des éléments finis

Recommend Documents