Poly
Poly
Poly
1
Assemblage et conditions aux limites..........................................................................................................54
Exercice ........................................................................................................................................................54
Votre parcours pdagogique .......................................................................................................................56
UTILISATION D'UN LOGICIEL LMENTS FINIS ................................................................................................................ 57
Cration et vrification des donnes:..........................................................................................................57
Excution du calcul: .....................................................................................................................................58
Exploitation des rsultats: ...........................................................................................................................58
Votre parcours pdagogique .......................................................................................................................59
Organigramme d'un logiciel lments finis .................................................................................................60
PROCESSUS DANALYSE & MODLISATION...................................................................................................................... 61
Qu'est-ce qu'un modle ?............................................................................................................................62
Comment estimer les erreurs de discrtisation ?........................................................................................63
Votre parcours pdagogique .......................................................................................................................64
2
Mise en quations des barres 3/64
= D ( )
( ) < Lois de comport em ent > E ( )
= f (u )
Relations Lois de compor tement
gomt riques gnralise
Relat ions
T = n gomtriques
Dans le premier document de cours nous avons tabli la loi de comportement gnralise du modle barre.
H1 : dplacement axial u ( M , t ) = u ( x , t ) xo ==> xx = u, x
H2 : tat de contrainte uni axial xx = E xx
D'o la dfinition de l'effort normal N = ES u, x
Pour terminer la mise en quations des barres, nous pouvons crire une des deux formes du principe de la
mcanique que vous avez vues en MMC :
Le PFD : qui donne un systme d'quations aux drives partielles (formulation locale).
Le PTV : qui est sa forme intgrale ou forme variationnelle et est une forme nergtique globale des
quations du mouvement.
Application du PFD
Nous allons crire les quations de Newton f = ma pour une tranche dpaisseur dx de la barre
Le bilan des efforts extrieurs sur cet lment de matire (figure ci-contre) f
fait apparaitre l'effort normal (torseur des efforts de cohsion)
L'quation de rsultante dynamique dans la direction x donne :
N + dN N + fdx = Sdx u
N N + dN
Soit N, x + f = S u
dx
Compte tenu de la loi de comportement intgre, l'quation locale est : x
x ]0, [ Su ( ESu, x ) = f
,x
3
Mise en quations des barres 4/64
Pour dterminer la rponse dynamique en temps, il faudra se donner les deux conditions initiales:
u ( x, 0) = uo ( x) Dforme et vitesse de dformation
u ( x , 0) =
uo ( x ) initiales de la barre
Application du PTV
Nous allons crire le principe des travaux virtuels u W = A Fo
f
F
pour une barre charge sur sa longueur et ses extrmits. 0 u ( x, t )
Le travail virtuel des quantits dacclration est : A = Su u dx
o
Le travail virtuel des efforts se dcompose en travail virtuel des efforts de cohsion et celui des efforts
extrieurs soit :
Pour les efforts de cohsion Wint = : dV = xx xx dS dx = ES xx xx dx
D 0S 0
xx = u, x
Soit Wint = ESu, x u, x dx
o
Pour les efforts extrieurs Wext = f u dx + Fo uo + F u
o
Le travail des efforts de cohsion Wint peut s'exprimer partir de la variation de l'nergie de dformation
de la barre Wint = Ed
( )
2
avec 2 Ed = : dV = ES u, x dx
D o
4
Mise en quations des barres 5/64
Est quivalente P ( Su ESu, xx f ) dx = 0 P
0
Remarque : si u est une solution approche du problme cette forme intgrale
reprsente le rsidu pondr de lquation locale sur le domaine.
F = N (, t ) = + ES u, x (, t ) N0 0 N
D'o le PTV : P P S u dx + P,x ES u, x dx = P F + Po Fo + P fdx
0 0 0
Le PTV est la forme variationnelle du problme, c'est une formulation globale (notion d'nergie)
u Su u dx = ESu,x u,x dx + f u dx + Fo uo + F u
o o o
Sera utilis pour rechercher les solutions numriques du problme.
Solutions approches
Vous devez tre capable d'crire ces deux formulations pour un problme donn.
Tous les exercices de cours sont corrigs sur le site, mais il faut chercher les rponses aux questions avant
de consulter le corrig.
Exercice 1 : Mise en quations dun barreau en traction
Objectifs : Savoir crire les conditions aux limites pour une barre,
Rsoudre un problme simple en statique,
Pouvoir crire le PTV et savoir passer du PTV au PFD.
1- criture des conditions aux limites.
Donnez les diffrentes conditions aux limites homognes possibles pour une barre.
5
Mise en quations des barres 6/64
Donnez les conditions aux limites correspondantes aux trois figures ci-dessous.
F k x=0
xo xo xo
x= x=0
M
2- Application du PFD.
x
crire le systme d'EDP de ce problme
g Intgrer l'quation diffrentielle en statique
Tracer le diagramme de l'effort normal (analyse type RDM)
3- Application du PTV.
Pour le problme reprsent par la figure ci-dessous
k
x=0 xo x=
Pour un champ de dplacements virtuels cinmatiquement admissible.
Donner lexpression du PTV en ne considrant que la barre.
Retrouver cette expression en considrant la barre et le ressort.
4- quivalence des principes.
Donner lexpression du PFD et passez au PTV (application directe du cours).
Partir du PTV pour retrouver l'quation locale et toutes les CL du problme.
Dmarche inverse celle prsente en cours
Si avec la correction vous n'arrivez pas comprendre la rponse une question, c'est que des lments du
cours ou des pr-requis vous manquent. Revoyez le cours et n'hsitez pas poser la question votre
enseignant, il pourra vous aider rsoudre la difficult.
Pour assimiler le cours il faudra traiter des exercices non corrigs.
6
MEF : tude des treillis 7/64
Le champ de dplacement sur l'lment sera construit sur une approximation polynomiale deux
paramtres de la forme
Approximation linaire du champ de dplacement
a (t )
u ( x, t ) = a1 + a2 x =< 1 , x > 1 Ici les paramtres ai n'ont pas de sens physique
a
2
( t )
1
x N (i ) = 0 N2
N 2 ( x) = vrifie 2
e N2 ( j ) = 1 0 1
x/le
La notion d'approximation nodale est fondamentale dans la mthode des lments finis, elle permet
dutiliser des variables qui ont un sens physique, et sur lesquelles nous pourrons directement imposer les
valeurs donnes par les conditions aux limites de type cinmatique.
7
MEF : tude des treillis 8/64
e
Le travail des quantits d'acclration est : Ae = Su u dx
o
Utilisons lapproximation nodale du champ des dplacements u = < N > {U e }
Le terme u u = u T . u = { U e } < N >T < N > {Ue }
T
e
On peut alors sortir les variables nodales de l'intgrale Ae = { U e } S < N > dx {Ue }
T
<N>
T
o
e
Ae = { U e } [ M e ]{Ue } avec [ M e ] =
T
<N>
T
S < N > dx
o
Nous venons de dfinir la matrice masse lmentaire, le calcul de l'intgrale se fait analytiquement,
on trouve : A titre d'exercice retrouver par le calcul
S e 2 1 les coefficients de cette matrice
[M e ] =
6 1 2
e
Ed = Wint
Le travail des efforts intrieurs est : Wint = ESu, x u, x dx
o
e
ES ( u,x )
2
Pour ce calcul utilisons l'expression de l'nergie de dformation : 2 Ed = dx
o
e
u,2x = uT, x . u, x
2 Ed = {U e } ES < N, x > dx {U e }
T
< N,x >
T
Utilisons l'approximation nodale
o
1 1
[ Ke ]{U e } avec [Ke ] =
ES
Soit pour chaque lment 2 Ed = {U e }
T
e 1 1
f
Le travail des efforts extrieurs est : Wext = f u dx + Fie ui + F je u j Fie F je
(e)
x
o i j
e
2
Ce calcul permet de calculer les charges nodales quivalentes au sens de lapproximation une charge
volumique relle applique la structure
f f e f e
PTV 2 2
i j i j
Charge relle Charge nodale quivalente
8
MEF : tude des treillis 9/64
Exemple Objectif : Dterminer une approximation des premires frquences de rsonnance de la barre avec
un modle lment fini.
Lorsque l'on somme les nergies de chaque lment pour obtenir l'nergie de la structure les matrices
lmentaires s'emboitent les unes avec les autres, en effet
ES 2
Pour l'lment 1 : 2 Ed 1 = (u1 2u1u2 + u22 )
ES 2
Pour l'lment 2 : 2 Ed 2 = (u2 2u2u3 + u32 )
ES 2
Soit pour les deux lments 2 Ed 1+ 2 = (u1 2u1u2 + 2u22 2u2u3 + u32 )
Que l'on peut crire sous forme matricielle 2 Ed1+ 2 = {U }
T
[ K ]{U }
1 1 0
ES
Avec [ K1+ 2 ] = 1 1 + 1 1 sur {U } = < u1 u2
T
u3 >
0 1 1
c'est l'assemblage.
En gnralisant aux n lments on obtient une matrice (n + 1, n + 1) , mais il faut tenir compte de la
condition d'encastrement du premier nud, tous les termes u1 sont nuls, la matrice assemble rduite
est une matrice carre de dimension n
2 1
1 2 1
1 2 1
nES
Matrice raideur assemble rduite ( u1 = 0 ): [K ] = 1 \ \
L
\ \ \
\ 2 1
1 1
De mme pour l'nergie cintique
9
MEF : tude des treillis 10/64
4 1
1 4 1
1 4 1
SL
Matrice masse [M ] = 1 \ \
6n
\ \ \
\ 4 1
1 2
Pour le calcul des pulsations propres (voir fichier MAPLE sur le site)
ES
1 1, 61
SL2 ES
Avec n=2 comparer i = 1, 571 et 4, 712
5, 63 ES SL2
2 SL2
ES
1 1, 589
SL2
ES ES
Pour n=3 2 5,196 anal = 1,571 4, 712 et 7,854
SL 2
SL2
9, 426 ES
3 SL2
La convergence est lente (lments de degr 1)
Avec la matrice modale calcule dans Maple vous pouvez tracer les modes sur la solution analytique, si le premier mode peut
tre assez rapidement approch par des segments, il faudra un maillage trs fin pour approcher la dforme modale des
modes suprieurs.
10
MEF : tude des treillis 11/64
ES 1 1
Reportons ce changement de base dans l'expression de l'nergie de dformation. {U e } {U e }
T
e 1 1
T
ui ui
v T
i C S 0 0 ES 1 1 C S 0 0 vi
2 Ed =
u
j 0 0 C S e 1 1 0 0 C S u j
v j v j
Nous en dduisons lexpression de la matrice raideur lmentaire sur les variables < ui vi uj vj >
[A] [A] C C S
[K e ] = ES
2
avec [A] =
e [ A] [ A] C S S 2
Dans le cas bidimensionnel il est possible de mener les calculs la main.
Ce n'est plus le cas pour les structures tridimensionnelles, c'est pourquoi nous les traiterons exclusivement
du point de vue numrique.
Assemblage et rsolution
Pour chaque lment de la structure nous avons :
Fie
e [ M e ]{Ue } + [ K e ]{U e } = {e } + Les Fie sont les efforts du nud i sur les
F je lments e (effort appliqu l'lment)
e D'o le signe moins pour avoir les efforts
L'assemblage consiste sommer les nergies lmentaires = des lments sur le nud.
D e 0
Pour les efforts nodaux l'quilibre d'un nud quelconque donne Fi Fie = 0
e
Les Fi reprsentent les efforts extrieurs appliqus aux nuds de la structure. Se sont soit des
charges donnes soit des efforts aux appuis (conditions cinmatiques) qui sont des inconnues du
problme.
L'assemblage consiste se donner un ordre de rangement des variables nodales dans le vecteur des
inconnues globales du systme. En pratique ( la main) nous utilisons l'ordre lexicographique pour simplifier
l'criture. La machine (calculateur) utilisera sa propre numrotation pour optimiser la vitesse de traitement
et la taille mmoire utile en fonction des algorithmes de rsolution qu'il utilisera pour traiter les quations,
ces oprations sont transparentes pour l'utilisateur.
11
MEF : tude des treillis 12/64
En statique nous utiliserons une dcomposition du systme matriciel en dplacements inconnus (nuds
ou les charges sont donnes) et dplacements imposs (les charges sont alors inconnues).
[ K11 ] [ K12 ] {U i } {Fd }
=
[ K 21 ] [ K 22 ] {U d } { Fi }
Le premier bloc d'quations nous donne le vecteur des dplacements nodaux inconnus:
{U i } = [ K11 ] {{Fd } [ K12 ]{U d }}
1
Cest le systme rduit
En reportant dans le second nous obtenons le vecteur des efforts de liaison inconnus:
{Fi } = K 22 K 21K111 K12 {U d } + K 21 K111 {Fd }
Dans les exercices trs souvent les dplacements sont imposs nuls, ce qui simplifie les critures et les
{U i } = [ K11 ] {Fd } {Fi } = [ K 21 ]{U i }
1
calculs puis
Post-traitement
Pour effectuer le dimensionnement d'une structure nous avons besoin de calculer l'tat de contrainte dans
la structure, pour un treillis cela revient calculer l'effort normal dans les lments.
Nous utilisons la loi de comportement intgre :
N = ES u, x = ES < N, x > {U e } =
ES
(u j ui ) = Cte
e
L'tat de contrainte est constant dans chaque lment fini
En statique, pour des treillis chargs aux nuds le modle lments finis ne ncessite qu'un lment par
barre du treillis, il donnera la solution analytique exacte. Ce n'est videmment pas le cas ni pour une
colonne charge par son poids propre, ni pour les problmes de dynamique, ou la solution exacte se
dcompose sur une base de fonctions sinusodale (cf chapitre sur les solutions analytiques pour les barres) .
Dans le cas bidimensionnel, ltat de contrainte sur un lment est donn par :
ES ES u j ui
N= (u j ui ) = < C S >
e e v j vi
Exemple
F Analyse
a Nous avons 3 nuds 2 variables par nuds (ui , u j ) les dplacements
a 2 du nud dans le plan. Modle 6 degrs de libert
yo
{U } = {u1 v1 u2 v3 }
T
xo v2 u3
v3 vecteur des dplacements nodaux
u3
3
Les conditions aux limites :
a
u = 0 X
Appui au nud 1 : 1 soit deux efforts inconnus : 1
1
a 2 2 u2
v1 = 0 Y1
Appui glissant au nud 2 : v2 = 0 soit un effort inconnu : Y2
3
F Le travail virtuel des efforts donns et inconnus appliqus la structure
Y1 a Y2 conduit lexpression du vecteur des forces nodales :
{F } = { X 1 0}
T
1 X1 a 2 2
Y1 0 Y2 F
{U } = {0 v3 }
T
pour 0 u2 0 u3
12
MEF : tude des treillis 13/64
xo 1 1 1 1
v3 ES 1 1 1 1
Pour llment 2 (1,3) K2 =
3 u3 2 a 1 1 1 1
a 1 1 1 1
sur {u1 v3 }
v1
= 45 v1 u3
1 u1
1 1 1 1
v3
ES 1 1 1 1
3 u3 Pour llment 3 (2,3) K3 =
2a 1 1 1 1
a
1 1 1 1
v2
= 135
sur {u2 v3 }
2 u2
v2 u3
Lnergie de dformation totale de la structure est la somme des nergies de dformation de chaque
lment, lassemblage des matrices consiste ranger chaque terme dans une matrice globale dfinie sur
le vecteur {U } = {u1 v1 u2 v3 }
T
v2 u3
Do la matrice globale
1 + 2 1 2 0 1 1 les termes de
la matrice K1 sont en bleu
1 1 0 0 1 1 la matrice K2 sont en rouge
ES 1 la matrice K3 sont en vert
[K ] = 2 0 1+ 2 1 1
2a 0 0 1 1 1 1
1 1 1 1 1 + 1 1 1
1 1 1 1 1 1 1 + 1
13
MEF : tude des treillis 14/64
F a
u2 = ES
1 + 2 1 1 u2 0 2
Allure de la dforme
ES F 1
1 2 0 u3 = F u3 = a (1 +
u3
) F
2a ES 2 2
1 0 2 v3 0
F a 1
u2 2
v3 =
ES 2 2
ES 2 F / 2
Y 2 = ( u 2 + u3 v3 )
2a
Post-traitement
Calculons l'effort normal dans les lments Lquilibre de chaque nud
ES est vrifi
N1 = ( u2 ) = F / 2 F
a 2 N2 N3
ES N 2 (2) (3) N3 F / 2
N = ES u, x ==> N 2 = ( u3 + v3 ) = F / 2 F
N1 (1) N1
a 2
F / 2
ES
N3 = ( u3 u2 v3 ) = F / 2
a 2
Remarques
Tous les calculs sont systmatiques et la dmarche suivie sera toujours la mme en statique.
Facilit de programmation de ce type de solution
Seule lanalyse, du problme et des rsultats, reste la charge de lingnieur.
La matrice raideur du systme rduit tait inversible " det( K ) 0 " car les conditions aux limites en
dplacement bloquaient tous les modes rigides de la structure.
Problme statique bien pos
Les efforts calculs aux appuis quilibrent parfaitement le chargement.
Les rsidus d'quilibre sont nuls, car nous travaillons sur la solution analytique de l'quation
matricielle. Dans le cas dune rsolution numrique ces rsidus doivent tendent vers zro (erreur
numrique).
Les contraintes calcules sur les lments quilibrent de faon exacte (aux rsidus prs) les charges
nodales. Ceci est vrai dans ce cas particulier calcul statique dun treillis charg aux nuds car
lapproximation utilise reprsente le champ exact de la solution analytique effort normal constant dans
chaque lment de la structure .
Erreur de discrtisation qui est nulle
En post traitement il est possible disoler un un chaque lment de la structure pour crire lquation
matricielle de lquilibre de llment. Ces calculs permettent de dterminer les efforts internes aux nuds
de la structure, nous en donnons des exemples dans les exercices de cours.
14
MEF : tude des treillis 15/64
Exercices
Les exercices de cours sont corrigs sur le site, il faut chercher les rponses avant de consulter le corrig.
Exercice 8 : Modlisation EF d'une colonne sous son poids propre
Objectifs : Notion d'erreur de discrtisation, et analyse des rsultats EF.
tude de convergence en affinant le maillage.
x 1- tablir l'quation matricielle d'un modle un lment fini
Analyser les rsultats aux nuds (dplacement & efforts)
Tracer les rsultats sur l'lment (dplacement & efforts)
g = 6h
2- Construire le modle utilisant deux lments finis identiques
Analyser les rsultats aux nuds (dplacement & efforts)
Tracer les rsultats sur les lments (dplacement & efforts)
3- Modle 3 lments, pour affiner le maillage dans la zone la plus contrainte
nous utilisons 3 lments de longueur h, 2h, et 3h.
Dduire des calculs prcdents l'quation matricielle du modle.
Tracer les rsultats sur les lments (dplacement & efforts)
Pour amliorer la solution lments finis nous avons augment le nombre d'lments et densifi le
maillage dans la zone la plus charge. Cette mthode dite h convergence demande en gnral un
nombre lev d'lments finis.
La figure suivante prsente les rsultats dun modle lments finis en contraintes planes. Pour quantifier lerreur relative
cette discrtisation, la dmarche est identique celle de cet exercice, elle est base sur lanalyse de la discontinuit du
champ des contraintes entre deux lments adjacents.
en MPa
Discontinuit
145
solution lments finis Lerreur est beaucoup trop importante.
constante par morceau
Ce modle nest pas satisfaisant, il faut
83
affiner le maillage
62
15
MEF : tude des treillis 16/64
16
MEF : tude des treillis 17/64
17
MEF : tude des treillis 18/64
18
Mise en quations des poutres 19/64
Nous avons tabli dans le chapitre prcdent la loi de comportement gnralise du modle poutre.
Application du PFD
Nous allons crire les quations de Newton f = ma pour une tranche dpaisseur dx de la poutre
Le bilan des efforts extrieurs sur cet lment de matire (figure ci-contre) f
fait apparaitre le torseur des efforts de cohsion, l'effort tranchant est associ
T Mf + dMf
aux contraintes de cisaillement qui s'opposent au glissement des sections.
Les quations de rsultante et de moment dynamique sont : Mf x
T + dT T + fdx = Sv dx T + dT
dx
dx dx
(T + dT ) 2 + M f + dM f M f + T 2 0
Soit
On nglige le moment dynamique
x ]0, [ Sv + M f , xx = f de rotation des sections.
T = M f ,x
Compte tenu de la loi de comportement intgre, l'quation locale est : x ]0, [ Sv + EIv, x4 = f
19
Mise en quations des poutres 20/64
Application du PTV
y Fo f
Nous allons crire le principe des travaux virtuels u W = A F
pour une poutre charge sur sa longueur et ses extrmits. Mo
M
0
x
On nglige le moment dynamique
Le travail virtuel des quantits dacclration est : A = Sv v dx de rotation des sections.
o
Le travail virtuel des efforts se dcompose en travail virtuel des efforts de cohsion et celui des efforts
extrieurs soit :
Pour les efforts de cohsion Wint = : dV = xx xx dS dx
D 0S
xx = yv, xx
Soit Wint = EIv, xx v, xx dx xx = E xx
o
Pour les efforts extrieurs Wext = f v dx + Fo vo + F v + M oo + M
o
Le travail des efforts de cohsion Wint peut s'exprimer partir de la variation de l'nergie de dformation
de la barre Wint = Ed
( )
2
avec 2 Ed = : dV = EI v, xx dx
D o
20
Mise en quations des poutres 21/64
Effectuons deux intgrations par partie du terme EIv,xx v,xx dx
o
Fait apparatre les conditions aux
EIv,xx v,xx dx = v, x EIv, x2 0 v, x EI v, x3 dx limites en rotation et moment.
o 0
Fait apparatre les conditions aux
limites en flche et force.
EIv,xx v,xx dx = v, x EIv, x2 0 v EIv,x3 0 + v EI v, x4 dx
o 0
Reportons dans : v Sv v dx = EIv,xx v,xx dx + f v dx + Fo vo + F v + M oo + M
o o o
(
v F EI v
) + ( M + EI v )
v ( Sv + EIv,x )
o , x3 o , x2 o
En regroupant les termes : v
o
f dx =
( ) + ( M EI v )
4
o v F + EI v 3
,x , x3
Le choix de vo 0 et v = 0 sur ]0, ] , nous donne la condition aux limites en force en x=0
(
Fo EI v )
, x3 x =0
=0 Fo = To
Cette condition tient compte de lorientation de
De la mme faon la normale extrieure au domaine
Pour ( v, x ) 0
o (
M o = EI v )
, x 2 x =0
= M f o
Pour v 0 (
F = EI v )
, x3 x =
= T
Pour ( v, x ) 0
(
M = EI v )
, x 2 x =
= M f
Vous devez tre capable de faire la dmonstration dans les deux sens PTV PFD .
Le PTV est la forme variationnelle du problme, c'est une formulation globale (notion d'nergie)
v Sv v dx = EIv,xx v,xx dx + f v dx + Fo vo + F v + M o o + M
o o o
Sera utilis pour rechercher les solutions numriques du problme.
Solutions approches
21
Mise en quations des poutres 22/64
Vous devez tre capable d'crire ces deux formulations pour un problme donn.
Tous les exercices de cours sont corrigs sur le site, mais il faut chercher les rponses aux questions avant
de consulter le corrig.
Exercice 3 : Mise en quations dune poutre en flexion plane
Objectifs : Savoir crire les conditions aux limites pour une poutre,
Rsoudre un problme simple en statique,
Pouvoir crire le PTV et savoir passer du PTV au PFD.
Les hypothses sont celles des poutres longues en petites dformations et petits mouvements. Le
matriau est suppos homogne isotrope lastique
1- criture des conditions aux limites
Exprimer les 4 conditions aux limites homognes suivantes :
g
Dterminer la solution analytique en statique, pour M = 0 .
Calculer la dforme de la poutre
Dterminer le diagramme du moment de flexion
3- Application du PTV.
Pour le problme reprsent par la figure ci-dessous, donner lexpression du PTV correspondant
des champs de dplacements virtuels cinmatiquement admissibles.
yo
g
t
( , E, I, S)
M
xo
Peut-on transformer le PTV pour retrouver l'quation locale et les conditions aux limites.
Les exercices de cours sont corrigs sur le site, il faut chercher les rponses avant de consulter le corrig.
Pour assimiler le cours il faut aussi traiter des exercices non corrigs.
22
Mise en quations des poutres 23/64
23
Mise en quations des poutres 24/64
24
MEF : tude des portiques 25/64
Pour identifier nos quatre variables nodales, nous utilisons une approximation polynomiale cubique de la
forme :
a1 (t )
3 a2 (t ) Approximation de degr 3
v ( x, t ) =< 1 x x
h 2
x > 4 variables
a3 (t )
a4 (t )
Par identification des variables nodales avec lapproximation de la flche et de la rotation aux noeuds, nous
obtenons la relation matricielle suivante :
vi (t ) v ( o, t ) 1
h
0 00 a (t )
(t ) h 1
i ( o, t ) 0 1 0 0 a (t )
v (t ) = h
2
2 3
j v ( , t ) 1 a3 (t )
j (t ) h 1 2 3 2 a4 (t )
(, t ) 0
Inversons cette relation et reportons le rsultat dans l'expression de l'approximation, nous obtenons :
vi (t )
(t )
i
v ( x, t ) = < N >e {U e } = < N1 N 2 N3 N 4 >
h
v (t )
j
j (t )
Avec les fonctions d'interpolation suivantes :
25
MEF : tude des portiques 26/64
N 1 ( s) = 1 3s 2 + 2 s 3 x
o s = 1
N1 ( s ) N3 ( s )
N 3 ( s) = 3s 2 s
2 3
x
N1 et N3 reprsentent la dforme d'une poutre bi - encastre pour laquelle 1
s=
0
on impose un dplacement unit une des deux extrmits
N 2 ( s ) = ( s 2 s 2 + s 3 )
N2 (s)
N 4 ( s ) = ( s 2 + s 3 ) 1 1 s
0 1
N2 et N4 reprsentent la dforme d'une poutre encastre une extrmit. N4 ( s)
Pour laquelle on impose une rotation unit l'autre extrmit.
A = o Sv v dx
L
Partons de u W = A avec Wint = Ed avec 2 Ed = EI ( v, xx ) dx
2
o
L
Wext = f v dx + Fo vo + FL vL + M oo + M L L
o
La poutre pouvant tre modlise par plusieurs lments finis nous calculerons les nergies sur
chaque lment puisque l'approximation nodale est une approximation lmentaire.
o
Utilisons lapproximation nodale du champ des dplacements v, xx = < N, xx > {U e }
En reportant dans l'nergie de dformation, pour chaque lment nous obtenons l'expression
matricielle de l'nergie de dformation lmentaire :
2 Ed = {U e } EI [ N, xx ] dx {U e }
T
[ N,xx ]
T
0
La matrice raideur associe est [ K e ] = [ B ]T EI [ B] dx
0
6 2 6 2
avec [ B ] = < N , xx > = < 2
( 1 + 2 s ) , ( 2 + 3s ) , 2
(1 2s ) , ( 1 + 3s ) >
12 6 12 6
A titre dexercice calculez le terme
EI 6 4 2 6 2
2
Tout calcul fait on trouve : [ K e ] = 3 (1,2) de cette matrice
12 6 12 6
6 2 2 6 4 sur <v , ,v , >
2
i i j j
Cette matrice n'est pas adimensionnelle car v et n'ont pas la mme dimension.
Pour que les coefficients de la matrice soient adimensionnels il faut travailler sur les variables v et
26
MEF : tude des portiques 27/64
Cette expression peut vous
12 6 12 6 permettre de simplifier vos
EI 6 2
calculs numriques.
[ Ke ] = 3z 612 46 12
6
6 2 6 4 sur <v , ,v , >
i i j j
13 11 9 13 420
35 210 70
11 1 13 1
D'o la matrice masse lmentaire est [ M e ] = S 210 105 420 140
9 70 13
420
13
35 11
210
13 420 1140 11 210 1105
sur <vi ,i ,v j , j >
27
MEF : tude des portiques 28/64
Exemple Objectif : Dterminer la rponse statique de la poutre avec un modle lment fini.
y F Modle 1 lment fini
x Ce modle comporte 4 variables : X T =< v1 , 1 , v 2 , 2 >
Les conditions aux limites : (v1 , 1 ) = (0, 0)
Y1
M1
v2 2 dplacements inconnus : X IT =< 0, 0, v2 , 2 >
2
2 efforts inconnus : FI T =< Y1 , M 1 , 0, 0 >
1 2
12 6 12 6 0 Y1
2
EI 6 4 2 6 2 0 M 1
Le PTV appliqu l'lment nous donne l'quation matricielle 3 =
12 6 12 6 v2
F
6 2 2 6 4 0
2
2
EI 12 6 v2 F v2 F 3 1/ 3
Les quations donnant la dforme sont : 2 =
==> =
3 6 4 2 0 2 EI 1/ 2
C'est la solution exacte de la RDM
Les quations donnant les efforts l'encastrement sont :
EI 12 6 v2 Y1 Y1 12 6 1/ 3 F
3 6 2
2 = ==> = F 2 =
2 M1 M1 6 2 1/ 2 F
On vrifie les quations d'quilibre de la structure
Dans cet exemple le modle lment fini donne la solution exacte car celle-ci est un polynme d'ordre 3
comme l'approximation utilise.
Pour calculer ltat de contrainte sur les lments, le diagramme du moment de flexion et celui de l'effort
tranchant, nous utilisons la loi de comportement intgre.
M f = EI v, xx = EI < B, x2 > {U e }
Pour chaque lment nous crirons :
T = EI v, xxx = EI < B, x3 > {U e }
6 2 6 2
Rappel : < B, x2 > = < 2 ( 1 + 2 s ) , ( 2 + 3s ) , 2 (1 2s ) , ( 1 + 3s ) >
Vous notez que le moment de flexion Mf 12 6 12 6
est linaire et que leffort tranchant est < B, x3 > = < 3 , 2 , 3 , 2 >
constant par lment.
28
MEF : tude des portiques 29/64
Modle 1 lment.
Dterminer la matrice raideur, et le vecteur force gnralis associ au poids propre.
crivez le systme rduit des quations, calculez les dplacements nodaux.
Calculer la flche au centre de la poutre, comparer la solution analytique
5 gS 4
v ( / 2) =
384 EI
Calculer les efforts aux appuis, et vrifier l'quilibre global de la structure.
Calculer les efforts sur l'lment, tracer les diagrammes de l'effort tranchant et du moment de
flexion, comparer la solution analytique.
Modle 2 lments.
Dterminer la matrice raideur assemble complte.
Dterminer le vecteur force gnralis associ au poids propre de la structure.
crivez le systme rduit des quations, calculez les dplacements nodaux et comparer la
solution analytique.
Calculer les efforts aux nuds, comparer la solution analytique.
Calculer les efforts sur l'lment et tracer les diagrammes de l'effort tranchant et du moment de
flexion, et comparer la solution analytique.
Comparer la solution analytique..
Rpondez aux mmes questions
Prise en compte de la symtrie
Utiliser la symtrie pour simplifier le modle
Calculer la matrice raideur et retrouver la solution du modle 2 lments.
29
MEF : tude des portiques 30/64
7 F 3 1 F 2 1 F 2
v (C ) = , (C ) = , (B) =
768 EI 128 EI 32 EI
M f ( A) = 3F /16 , M f (C ) = 5 F / 32
30
MEF : tude des portiques 31/64
Il est clair que nous ne manipulerons pas ces matrices manuellement, d'autant que pour effectuer
l'assemblage d'une structure portique il faut effectuer un changement de base pour exprimer toutes les
matrices lmentaires sur une base globale.
Il faut passer aux calculs numriques MEFlab, Cast3M ou Abaqus
ES / S 2
= = lancement de la poutre
EI / 3 I
Pour on tend vers la solution obtenue en ngligeant les dformations de traction
Pour un lment horizontal (orient de i vers j suivant la direction des x) y
La base locale et la base globale correspondent, la matrice raideur est vi vj
i j u
celle donne juste avant sur < ui , vi ,i , u j , v j , j > ui j
j x
i
31
MEF : tude des portiques 32/64
v =0
2 Cette hypothse permet d'crire deux quations de liaison : 2
F u2 u2 u3 = u2
3
2
Le modle ne comporte plus que 2 variables
Calculons directement les matrices lmentaires sur ces 2 variables.
1
EI 12 6
Pour llment 1 (2-1) : [ K1 ] = 6 4 2
3
EI 0 0
Pour llment 2 (2-3) : [ K 2 ] = 0 4 2
3
2 F 3
EI 12 6 u2 F u2 =
==> {U } =
15 EI
D'o le systme rduit des quations :
3 6 8 2 = 0 2
2 = 1 F
2 10 EI
u u
C'est la solution exacte de la RDM
Allure de la dforme
32
MEF : tude des portiques 33/64
2 3 32
< R22 M 22 R32 M 32 >= F < 0,6 0,4 0,6 0,2 >
Ce modle ne nous donne pas toutes les composantes deffort car nous avons nglig les allongements
des lments.
Pour calculer la composante verticale de leffort au
0,6 F
noeud 1, nous pouvons crire les quations d'quilibre F - 0,2 F
de la structure.
33
MEF : tude des portiques 34/64
34
Formulation variationnelle & criture Matricielle 35/64
Formulation Variationnelle
& criture Matricielle
Ce paragraphe consacr la formulation variationnelle dun problme de mcanique, reprend des notions
abordes en Mcanique des Milieux Continus MMC et en Mathmatiques, il est la base de toutes les
simulations en ingnierie mcanique et doit faire partie, ce titre, de votre culture scientifique.
La mthode des rsidus pondrs (annulation d'erreur) utilise comme point de dpart le systme dEDP
(quations diffrentielles) dfini par les quations locales et les conditions aux limites du problme. La
formulation intgrale du problme ainsi obtenue peut tre utilise comme point de dpart c'est La
formulation variationnelle du problme. En mcanique cette forme variationnelle est le PTV (Principe des
Travaux Virtuels).
Problme aux limites
Systme d'EDP formulation forte du problme
Toutes les mthodes d'approximation ont un mme objectif, remplacer un problme mathmatique dfini
sur un milieu continu (quations diffrentielles ou intgrales) par un problme mathmatique discret
(quation matricielle). Problme de dimension finie que l'on sait rsoudre numriquement.
Formulation intgrale
Partons de lquation locale et des conditions aux limites dfinies dans le cadre de la MMC.
M D u div = f PFD appliqu un lment
de matire
M u u = ud
M n = Td
35
Formulation variationnelle & criture Matricielle 36/64
Lquation locale
"EL" P P.( u div( ) - f ) dV = 0 (Formulation forte)
D
Sachant que* : grad s ( P ) = div ( P ) P. div ( )
"EL" P ( P. u + : grad s (P ) - P. f ) dV ( P. n ) dS = 0
D D
En tenant compte des conditions aux limites sur la frontire M n = Td
"EL."
C.L sur
P ( P. u + : grad s (P) - P. f ) dV P. n dS P.Td . dS = 0
D u
La mthode des lments finis utilise cette formulation, avec deux ides fortes
la construction systmatique des fonctions de forme par sous domaine lments finis
utilisation des variables nodales comme paramtres dapproximation ce qui permet dimposer
les conditions aux limites en dplacement du problme.
*
Cette formule utilise la symtrie du tenseur des contraintes et les relations suivantes:
T
: grad ( u ) = div ( u ) u. div ( )
T
: grad (u ) = div ( u ) u. div ( )
La dmonstration de ces relations se fait simplement partir d'une relation indicielle de la forme suivante:
ij ui , j = (ui ij ) , j ui ij , j
36
Formulation variationnelle & criture Matricielle 37/64
L'intrt de ce principe est de fournir directement la forme intgrale sans avoir passer par les quations
aux drives partielles
Pour exprimer le produit : nous utilisons une reprsentation vectorielle des tenseurs
{ } =< xx , yy , zz , xy , xz , yz >
T
La forme matricielle associe aux lois de comportement est alors : { (M) } = [ D( M )]{ ( M ) }
Les relations dformations - dplacements peuvent aussi s'crire sous forme matricielle.
{ (M ) }= {grad u } = [ L] {u
s (M ) (M ) }
[ L] Matrice d'oprateurs diffrentiels correspondant l'expression du gradient
symtrique des dplacements.
*
Cette relation sur le taux de dformation est dmontre en MMC, nous admettons, ici, ce rsultat sans dmonstration.
37
Formulation variationnelle & criture Matricielle 38/64
Cette quation pouvant tre crite quelque soit {q} , nous obtenons l'quation matricielle :
[ M ] {q} + [ K ] {q} = {F }
avec: [M ] = [W ]T [W ] dV
D
[ K ] = [ B] T [ D] [ B ] dV
D
{ F } = [W ]T { f } dV + [W ]T { Td } dS
D
[ ] = 21 22 0
l'intrieur du domaine est voisin de l'tat de
Hypothse contrainte sur les surfaces, donc plan par rapport
0 0 0 la normale la surface
1+
crivons l''inverse de la loi de HOOKE = trace( ) 1 +
E E
Pour dterminer le tenseur des dformations partir du tenseur des contraintes, en ne conservant que les
termes travail non nul.
La dformation 33 qui nest pas prise en compte
11 1 11
0 dans la loi de comportement, peut tre calcule
1
22 = 1 0 22 posteriori par : 33 = ( 11 + 22 )
2 E 0 0 2(1 + ) 12 E
12
11 1 0 11
E
Puis inversons cette relation : 22 = 2
1 0 22 ==> [ D ]
(1 ) (1 ) 212
12 0 0
2
Le modle dformations planes s'applique des pices massives dont les dformations
longitudinales seront suffisamment faibles pour tre ngliges.
11 12 0
Hypothse : [ ] = 21 22 0 soit : { } =< 11
T
22 212 >
0 0 0
38
Formulation variationnelle & criture Matricielle 39/64
= E
(1 + )(1 2 )
crivons la loi de HOOKE : = trace( ) 1 + 2 avec
= 2(1 + )
E
Pour exprimer l'approximation en fonction des dplacements nodaux nous allons identifier les valeurs
nodales des dplacements
u1 1 x1 y1 a1
soit pour le dplacement suivant x ui = u h ( xi , y i ) , aux 3 nuds : u2 = 1 x 2 y 2 a 2
u 1 x y 3 a 3
3 3
39
Formulation variationnelle & criture Matricielle 40/64
a1 23 31 12 u1 A = aire du triangle
1
a 2 = y 23 y 31 y12 u2 avec xij = xi x j et yij = yi y j
a 2 A x x13 x 21 u = x y x y
3 32 3 ij i j j i
40
Formulation variationnelle & criture Matricielle 41/64
W ( x , y ) = ( x 2 a 2 )( y 2 a 2 )
1- Soit : 1 Deux fonctions de forme dfinies sur le domaine
W2 ( x , y ) = W 1 ( x , y ) ( x 2
+ y 2
)
41
Formulation variationnelle & criture Matricielle 42/64
Formulation
Rappeler la forme variationnelle du problme.
Donner la forme matricielle correspondant une approximation de cette quation intgrale.
Calcul de la matrice raideur lmentaire
Les quatre lments tant identiques, on limitera les calculs au premier lment (l'lment 123).
Dfinir les fonctions d'interpolation de cet lment.
Dterminer la matrice raideur et le vecteur force gnralis de cet lment.
Montrer que pour les autres lments les rsultats sont identiques.
Rsolution
Calculer l'approximation du champ de vitesse au centre de la conduite
42
Formulation variationnelle & criture Matricielle 43/64
u ( r , z , t ) = w w ez
z
yo e
1- Exprimer l'oprateur gradient du = grad (u ) dX sur la base b xo
b er
2- En dduire en petites dformations l'oprateur [ L ] ; { } = [ L ]{U }
w T3
j
En dduire l'expression de la matrice [ B ] ; { } = [ B ]{U e }
i
symtrie cylindrique
43
Formulation variationnelle & criture Matricielle 44/64
44
Mthodes numriques 45/64
Mthodes numriques
dans le cadre de la MEF
Les principales tapes de construction d'un modle lments finis sont les suivantes :
Discrtisation du milieu continu en sous domaines.
Construction de l'approximation nodale par sous domaine.
Calcul des matrices lmentaires correspondant la forme intgrale du problme.
Assemblage des matrices lmentaires - Prise en compte des conditions aux limites.
Rsolution du systme d'quations.
Discrtisation du milieu
Cette opration consiste dcomposer le domaine continu en un nombre fini de sous domaines
lments finis .
ne
D = De telle que lim De = D
taille des e 0 e
e =1
Il ne doit y avoir ni recouvrement ni trou entre deux lments ayant une frontire commune. De plus
lorsque la frontire du domaine est complexe, une erreur de discrtisation gomtrique est invitable.
Cette erreur doit tre estime, et ventuellement rduite en modifiant la forme ou en diminuant la taille
des lments concerns.
Modifier la taille des
lments
Pice prsentant des
congs de raccordement
Changer la gomtrie
lments frontire courbe
Erreur de discrtisation
gomtrique
Erreur de discrtisation gomtrique.
Approximation nodale
La mthode des lments finis est base sur la construction systmatique d'une approximation u h du
champ des variables par sous domaine. Cette approximation est construite sur les valeurs approches du
champ aux noeuds de llment. On parle de reprsentation nodale de lapproximation ou plus
simplement dapproximation nodale.
L'approximation par lments finis est une approximation nodale par sous domaines ne faisant intervenir
que les variables nodales du domaine lmentaire
De . M De {u } = [ N
h
(M ) (M ) ] {U n }
{u }
h
Valeur de la fonction approche en tout point M de l'lment
[N] Matrice des fonctions d'interpolation de l'lment
{U n } Variables nodales relatives aux nuds d'interpolation de l'lment
45
Mthodes numriques 46/64
+
lment 2
x
Approximation linaire +
utilisant 3 lments lment 3
46
Mthodes numriques 47/64
0 si i j
Les fonctions d'interpolation satisfont la proprit suivante M i N j ( Mi ) =
1 si i = j
Approximation linaire
base polynomiale utilise est (1, s ) 2 nuds 1 N1 N2
N 1 ( s) = L1 = 1 s 1
0 s
N 2 ( s) = L2 = s 1 2
Approximation quadratique N2 N3
1 N1
Base polynomiale (1, s, s 2 ) 3 nuds
N 1 ( s) = L1 (2 L1 1)
N 2 ( s) = 4 L1 L2 0 1 s
N ( s) = L (2 L 1)
3 2 2 1 2 3
L1
N 1 ( s) = 2 (3 L1 1)(3 L1 2) 1 N1 N2 N3 N4
9
N 2 ( s) = L1 L2 (3 L1 1)
2
9
N 3 ( s) = L1 L2 (3 L2 1) 0
1 s
2
L2 1 2 4
N 4 ( s) = 2 (3 L2 1)(3 L2 2)
3
lments triangulaires t
C'est un triangle rectangle de cot unit. s [ 0,1] et t [ 0,1 s ] (0,1) 3
47
Mthodes numriques 48/64
N1 = 1 s t N2 = s N3 = t
N1 t t N3
N2 3
3 t
3
1
1 1
2
2 2 s
s s
Approximation quadratique t
Base polynomiale (1, s, t , st , s 2 , t 2 ) (0,1) 3
lment 6 nuds, triangle de type "T6".
5 4
On pose : L1 = 1 s t , L2 = s , L3 = t
1 2 s
Les fonctions d'interpolation sont : (0,0) 6 (1,0)
Pour les 3 noeuds sommet: i = 1,2,3 Pour les 3 noeuds d'interface (4,5,6):
N i = Li (2 Li 1) i = 1, 2,3 avec j i, k et k i
3 N i+3 = 4 L j Lk
5 1
1
4 6
1 1 5
6 2
N 4 = 4 L2 L3
N1 = L1 (2 L1 1) 2
4
3
N1 1
= 4 (1 s)(1 t ) N1 s
1
3
N2 = 4 (1 + s)(1 t ) 1
2
1
1
N3 = 4 (1 + s)(1 + t ) (-1,-1) (1,-1)
1 2
N 1
= 4 (1 s)(1 + t )
4
Approximation quadratique "Q8" t
7 6 5
Pour viter davoir des nuds internes, on utilise des bases polynomiales
incompltes mais symtriques contenant tous les monmes dun mme degr. s
8
4
2 2 2 2
Base polynomiale (1, s, t , st , s , t , s t , t s )
1 2 3
Nous avons donn les fonctions dinterpolation des lments les plus simples pour lesquels il est possible
de faire un certain nombre de calculs analytiquement. Signalons que les expressions des fonctions
d'interpolation de nombreux autres lments de rfrence sont donnes dans le livre de Dhatt et Touzot
"Une prsentation de la mthode des lments finis". Vous y trouverez aussi des exemples de programmes
pour la construction systmatique de ces fonctions dinterpolation.
48
Mthodes numriques 49/64
linaire e
quadratique
cubique
Transformations gomtriques d'lments une dimension
Ces transformations gomtriques utilisent les fonctions d'interpolation linaire,
quadratique et cubique dfinies plus haut.
lments deux dimensions.
Pour ces lments les transformations gomtriques conduisent respectivement des frontires linaires,
quadratiques.
49
Mthodes numriques 50/64
s s
(0,0) (1,0) (-1,-1) (1,-1)
Linaire (3) Quadratique (6) Linaire (4) Quadratique (8)
x y z
s s s s x x
y z = [ J ]
t =
x
t t t y y
x y z
u u u
u z
z
x =< N g ( s , t , u ) > { x n }
Compte tenu de lexpression de la transformation gomtrique y =< N g ( s , t , u ) > { y n }
z =< N ( s , t , u ) > z
g { n}
< N g >
s
s
[
[ J ] = t < x y z >= < N g > {x n } { y n } {zn }
t
]
u < N g >
u
[ J ] est le produit d'une matrice (3, n) par une matrice (n, 3) connues.
Pour chaque lment, la matrice Jacobienne s'exprime en fonction des drives des fonctions
de la transformation gomtrique et des coordonnes des nuds gomtriques de l'lment
rel.
50
Mthodes numriques 51/64
La relation inverse permet alors de calculer les drives premires par rapport aux coordonnes relles des
fonctions d'interpolation.
x s
1 1
y = [ J ] t La transformation devant tre une bijection [ J ] existe
z u
Une singularit de J peut apparatre lorsque l'lment rel est trop "distordu" par rapport l'lment de
rfrence lment dit dgnr . De faon gnrale on vite lors du maillage d'utiliser des lments
trop disproportionns, car ils nuisent la prcision numrique du modle
De plus en plus de logiciels de pr-traitement proposent des outils de contrle de la qualit
du maillage bas sur la taille, les proportions et le calcul du Jacobien des lments utiliss.
Pour le calcul des drives secondes 2 et 2 des fonctions d'interpolation (problmes de flexion), la
x 2 xy
dmarche est identique mais les calculs sont plus complexes, reportez vous au livre de D&T (pages 55-57).
Cette dernire intgrale ne peut tre value analytiquement que dans des cas extrmement simples. En
gnral, la fonction intgrer est une fraction rationnelle polynomiale complique. Le calcul de l'intgrale
sur l'lment de rfrence est donc effectu numriquement.
Les formules d'intgration numrique permettent d'valuer l'intgrale sous la forme gnrale suivante :
Npi
f dv f (i )i
Dref i =1
*
Pour un problme de dynamique lintgration rduite peut introduire des modes parasites (phnomnes dhourglass)
51
Mthodes numriques 52/64
Les schmas d'intgration les plus utiliss pour les lments 2D sont :
points Coordonnes i Poids i
1 1
4 s= t= =1
3 3
1 1
1 s= t= = 1/ 6
3 3
1/ 6 1/ 6
s = 2/3 t = 1/ 6 = 1/ 6
3
3
1 2 1/ 6 2/3
Vous trouverez dans le D&T (pages 280 300) les tableaux et les figures donnant la
position et les poids d'intgration pour diffrents schmas d'intgration.
Organisation des calculs numriques
pour chaque lment e
52
Mthodes numriques 53/64
x x y2 y1 x21 y21
[ J ] = x2 x1 =
y3 y1 x31 y31
Cette matrice est
une constante
3 1
A est laire de l'lment rel.
1 1 y31 y21
Son inverse : [ J ]
J =2A est le jacobien de la transformation.
=
2 A x31 x21
x 1 y31 y21 s
Sachant que : =
y 2 A x31 x21 t
Nous en dduisons :
N1, x 1 1 1 y31 + y21 1 y23
= J = =
N1, y 1 2 A x31 x21 2 A x32
N 2, x 1 1 1 y31
= J =
N 2, y 0 2 A x31
u
N 3, x 1 0 1 y21
xx
,x
= J = { } = = v = [ B ] {U e }
N
3, y 1 2 A x21 yy , y
2 xy u, y + v, x
D'o lexpression de la matrice [ B ( x, y ) ]
N 0 N 2, x 0 N3, x 0
1, x
[ B ( x, y ) ] = 0 N1, y 0 N 2, y 0 N 3, y
N1, y N1, x N 2, y N 2, x N3, y N3, x
y23 0 y31 0 y12 0
1
[ B ( s, t ) ] = 0 x32 0 x13 0 x21
2A
x32 y23 x13 y31 x21 y12 e paisseur suppose
uniforme de llment
Ces rsultats sont ensuite utiliss pour le calcul des matrices lmentaires :
1 1 s
Matrice raideur: [K e ] = [ B( M )]T [ D( M )] [ B( M )] dv = e [ B ( s, t )]T [ D] [ B( s , t )] 2 A dsdt
De 0 0
Les termes de cette matrice sont des constantes.
1 1 s
Matrice masse: [M e ] = [N ( x, y )]T [N ( x, y )] dv = e [N (s, t )] [N ( s, t )] 2 A dsdt
T
De 0 0
1 1 s
Dans le cas dune charge de volume f d {Fde } = e [ N ( s, t )]T { f d } 2 A dsdt
0 0
Nayant que des polynmes intgrer le calcul analytique ou numrique de ces matrices ne
pose pas de problme, et il conduira des rsultats exacts qui pourront tre rutiliss.
53
Mthodes numriques 54/64
e =1
Cette opration traduit simplement que lnergie associe au domaine tudi est la somme des nergies
lmentaires des sous domaines. Cela consiste ranger dans une matrice globale les termes des matrices
lmentaires. Le programme dfinira lordre des variables globales pour optimiser la place mmoire
(disque) de la matrice globale, mais aussi le temps de calcul en fonction des algorithmes de rsolution
utiliss*. Deux mthodes classiques matrices bandes et matrices ligne de ciel sont prsentes dans
le livre de D&T.
Aprs assemblage, nous obtenons la forme matricielle du principe des travaux virtuels :
[ M ] {U} + [ K ] {U } = { Fd } + { Fi }
Soit N quations pour N+P inconnues. Pour rsoudre, il faut tenir compte des P conditions
aux limites cinmatiques associes aux P composantes inconnues du vecteur { Fi } .
Dans le cas dun calcul statique, La mthode directe de rsolution par blocs, peut se prsenter sous la
forme suivante en regroupant les termes des matrices.
Dans le cas particulier ou {U d } = {0} seul les termes de [ K11 ] et [ K 21 ] sont utiles
{U i } = [ K11 ] { Fd 1}
1
En effet
{ Fi } = [ K 21 ]{U i } {Fd 2 }
Nous nabordons pas ici les mthodes numriques de rsolution de ces quations matricielles. Ces
mthodes sont vues en analyse Numrique. Pour les problmes de statique vous trouverez dans les codes
EF deux types de mthodes
Les mthodes directes : limination de Gauss, dcomposition LDU, ou Cholesky
Les mthodes itratives de type Gauss-Seidel.
Exercice
Les exercices de cours sont corrigs sur le site, il faut chercher les rponses avant de consulter le corrig.
Exercice 1 : EF-T3 pour llasticit plane
Objectifs : Assimiler les techniques de calcul au niveau lmentaire mise en uvre dans un
modle lment finis pour un problme de mcanique.
*
Le choix de lalgorithme dassemblage est un problme numrique & informatique. Les algorithmes dpendent de la taille du
systme matriciel, de la nature du problme (dynamique, statique, linaire, non linaire, etc ), de la machine utilise (place
mmoire disponible, espace de stockage, paralllisme, . ).
54
Mthodes numriques 55/64
55
Mthodes numriques 56/64
4 0 2b
b
1 e
2
xo
En dduire lexpression de [ J ]
1
Calculer la drive premire par rapport aux coordonnes relles des fonctions d'interpolation.
En dduire lexpression de la matrice [ B ] en fonction de s et t
Est il possible de calculer analytiquement la matrice raideur dun lment rectangulaire ?
Calculs numriques
Analyser le script Q4_ke qui utilise lintgration numrique
Le diaporama Q4 propos sur le site vous aidera faire le lien avec le cours.
Utiliser MEFLAB pour raliser un modle en contrainte plane d'une poutre console.
F=200Kg
Poutre en acier
h=20cm
yo L=3m section
e=1cm
rectangulaire
xo
56
Utilisation d'un logiciel lments finis 57/64
57
Utilisation d'un logiciel lments finis 58/64
Excution du calcul:
Ce bloc, le plus coteux en temps machine est souvent excut en tche de fond. Un fichier de rsultats
permet gnralement de vrifier que les diffrentes phases de calculs se sont correctement droules :
- Interprtation des donnes, vrification des paramtres manquants
- Construction des matrices (espace utile pour les gros problmes)
- Singularit de la matrice raideur (problme de Conditions aux limites ou de dfinition des
lments)
- Convergence, nombre d'itrations, etc. ...
Ce fichier peut ventuellement contenir les rsultats du calcul (dplacements, rsidus, contraintes,...) ce
qui lui confre dans ce cas un volume gnralement trs important.
Il peut arriver que le calcul choue. Les principales sources d'erreurs, gnralement observes ce niveau,
sont les suivantes:
Diffrentes vrifications doivent tre effectues pour valider les rsultats. Ces vrifications entranent dans
la plupart des cas remettre en cause le modle pour en crer un nouveau, dont on espre quil
amliorera la solution prcdente.
Pour valider une solution, il faut procder dans lordre
dans un premier temps, estimer la prcision du modle.
Puis procder sa validation. Vrification (et remise en cause) des hypothses du modle.
58
Utilisation d'un logiciel lments finis 59/64
Les indicateurs sur la prcision du modle sont gnralement locaux, ils concernent des informations
lmentaires calcules aux nuds ou aux points dintgration. Or ces informations sont souvent
extrapoles ou lisse pour tre reprsent en valeur moyenne sur llment ou en courbe diso valeurs
sur le domaine.
Les indicateurs locaux sur la prcision dun modle mcanique peuvent tre :
Discontinuit des contraintes entre des lments adjacents. Le plus simple, pour un matriau
isotrope, est de visualiser la contrainte quivalente de Von Mises, cela permet davoir une ide des
zones fortement charges et ayant un fort gradient de contrainte. Ces zones seront lobjet de toute
notre attention.
Valeur du tenseur des contraintes sur les bords libres ou chargs (certaines valeurs sont alors
connues). En pratique il faudra estimer ces valeurs partir des valeurs obtenues aux points
dintgration.
Densit dnergie interne de dformation sur chaque lment, lidal tant davoir un cart le plus
faible possible.
Ayant quantifi la qualit de la solution numrique (prcision), diffrents contrles vous permettrons de
valider votre modle :
Ordre de grandeur des rsultats obtenus
Vrification des hypothses du modle
Par exemple en lasticit linaire il faut vrifier que lamplitude des dplacements reste
faible par rapport aux dimensions de la structure, que les dformations et les contraintes
observes respectent les hypothses de linarits de la loi de comportement.
Que les choix de dparts sont justifis.
La comparaison et lanalyse des rsultats des diffrentes modlisations que vous aurez ralises, vous
permettra d'amliorer puis de valider un modle "final" fiable. Ltude mene vous amnera conclure sur
ladquation entre la structure et le cahier des charges.
La synthse de ces calculs prliminaires est indispensable car elle vous permet de justifier et de dfinir
les limites (du ou) des modles retenus.
Avant de passer la pratique, prcisons comment se droule une tude base sur l'utilisation d'un logiciel
lments finis.
59
Utilisation d'un logiciel lments finis 60/64
60
Processus danalyse & modlisation 61/64
Il lui faudra :
Formaliser les non dits et les rflexions qui justifient les choix explicites ou implicites de son analyse
du problme (dfinition de son modle).
valuer la confiance qu'il accorde aux rsultats produits.
Analyser les consquences de ces rsultats par rapport aux objectifs viss.
Ne perdez jamais de vue que l'analyse des rsultats ncessite une bonne comprhension des diffrentes
tapes mathmatiques utilises lors de l'approximation, pour pouvoir estimer l'erreur du modle
numrique par rapport la solution exacte du problme mathmatique. N'oubliez pas non plus que le
modle numrique ne peut fournir que des rsultats relatifs aux informations contenues dans le modle
mathmatique qui dcoule des hypothses de modlisation.
De faon gnrale, les diffrentes tapes d'analyse d'un problme physique s'organisent suivant le
processus schmatis par la figure suivante.
Processus danalyse
Problme physique
Hypothses de modlisation
Modle mathmatique
P Hypothses de discrtisation
r
o Modle numrique
c volution du
modle
d Erreur de discrtisation (valuation) numrique
u prcision sur les grandeurs dintrt
r volution du
e modle
EF Analyse des rsultats mathmatique
Vrification des hyp. de modlisation
61
Processus danalyse & modlisation 62/64
62
Processus danalyse & modlisation 63/64
En rsum, les questions essentielles auxquelles l'ingnieur devra rpondre s'il veut effectuer une analyse
par un modle numrique dans de bonnes conditions, sont :
Quel modle mathmatique utiliser ?
Quel modle numrique faut-il lui associer ?
Quelle est l'erreur d'approximation commise ? (Cette question est aborde ci-dessous)
Peut-on amliorer le modle numrique ?
Faut-il changer le modle mathmatique ?
63
Processus danalyse & modlisation 64/64
Il faut savoir qu'une erreur locale importante n'entrane pas que les rsultats sur le reste de la structure
sont faux, il peut y avoir convergence globale avec des erreurs locales importantes. (il faut estimer la zone
concerne par l'erreur locale et voir s'il est raisonnable de considrer qu'elle n'aura pas d'influence sur les
autres rsultats, l'exprience vous y entrainera)
Dans tous les cas il faut comprendre l'origine de l'erreur locale, elle peut simplement tre lie au maillage,
il est alors simple de la rduire si besoin. Ou bien venir d'une singularit du modle en quel cas un maillage
grossier la fera disparaitre, et un maillage fin la rendra encore plus importante et pourra cacher les autres
rsultats sur la pice tudie.
Un modle lment finis classique ne peut pas donner de rsultats satisfaisants au voisinage d'une
singularit. Il faut si c'est ncessaire changer le modle pour remplacer la singularit par un modle
rgulier, il sera alors possible de calculer la grandeur avec la prcision souhaite.
A vous de passer la pratique sur les diffrents projets proposs sur le site, nous vous conseillons de
dbuter par les projets d'initiation qui sont plus simple modliser. Les projets industriels ont deux
principaux objectifs :
L'analyse du problme conduisant proposer des modles simplifis
L'analyse des rsultats pour valider (ou non) ces modles.
64