Poly [PDF]

Zitiervorschau

MEF Méthode des Éléments Finis

Polycopié de cours Apprentis ECN Responsable du cours : Hervé Oudin

Objectifs : Aborder les notions fondamentales de la méthode des éléments finis à partir de l'étude des structures treillis et portiques. Puis généraliser à différents problèmes de la physique pour comprendre et savoir utiliser les modèles numériques de type éléments finis dans le cadre de problèmes plus complexes utilisant des modèles de l'ingénieur. Ce cours arbore la notion de modélisation indispensable pour dimensionner et valider un produit avec un logiciel de calcul industriel. Pré requis : En Mathématiques : Algèbre linéaire, géométrie, intégration et équations différentielles. En Physique : Mécanique des solides et mécanique des milieux continus. Des connaissances en vibration des structures sont un plus. Supports du cours : Tous les éléments de ce polycopié peuvent être consultés sur le site https://meefi.pedagogie.ec-nantes.fr/MEF/MEF.htm Vous y trouverez entre autre les corrigés des exercices de cours mais aussi d'autres supports pédagogiques qui peuvent vous servir pour compléter votre formation. Les notions abordées : Modèle barre application aux treillis Modèle poutre application aux portiques Méthodes variationnelles et méthodes numériques Notion de modélisation application à l'utilisation d'un code éléments finis Méthodes et outils d'analyse des résultats d'un modèle éléments finis Déroulement : Séances de 4h S1 Cours /TD 1 : Étude des treillis S2 Cours /TD 2 : Étude des portiques S3 TP1 MEFlab en salle informatique S4 Cours /TD 3 : FV & méthode numérique S5 TP2 MEFlab en salle informatique S6 TP Abaqus en salle informatique S7 TP Abaqus en salle informatique S8 Cours /DS

TA1 : exercice à rédiger sur les treillis TA2 : exercice sur les portiques TA3 : rédiger le CR de ce TP TA4 : rédiger le CR de ce TP TA5 : rédiger le CR de ce TP

Pour les TD et TP MEFlab : Chaque thème est abordé en TD de façon analytique pour montrer ce qu'il est possible de traiter à la main et établir des solutions de référence. Les TP sont basés sur des exercices simples à réaliser avec MATLAB. Ces TP sont l'occasion d'utiliser des outils numériques pour voir comment les calculs sont abordés pour des structures plus complexes. Les TA "Travail en Autonomie" L'objectif est pédagogique, ils vous permettent d'avoir une évaluation pédagogique de votre travail et de résoudre vos difficultés au fur et à mesure. Ils doivent être rendus lors de la séance suivante dans le déroulement du cours. Le TP ABAQUS : L'objectif est de savoir analyser les résultats obtenus pour différentes modélisations d'un même problème. Être capable d'en faire une synthèse et de la présenter dans un rapport. Évaluation : deux notes avec un même poids Les Travaux en Autonomie seront notés : (note collective le travail peut être rédigé par binôme) Le DS donnera lieu à une note individuelle.

Table des matières MISE EN ÉQUATIONS DES BARRES........................................................................................................................................... 7

Application du PFD.........................................................................................................................................7 Application du PTV.........................................................................................................................................8 Équivalence des principes..............................................................................................................................8 Bilan & exercice..............................................................................................................................................9 MODÈLE ÉLÉMENTS FINIS POUR L'ÉTUDE DES TREILLIS.............................................................................................. 11

L’élément fini barre......................................................................................................................................11 Modèle éléments finis d'un treillis ..............................................................................................................14 Exercices.......................................................................................................................................................19 MISE EN ÉQUATIONS DES POUTRES EN FLEXION PLANE ............................................................................................... 23

Application du PFD.......................................................................................................................................23 Application du PTV.......................................................................................................................................24 Équivalence des principes............................................................................................................................24 Bilan & exercice............................................................................................................................................25 MODÈLE ÉLÉMENTS FINIS POUR L'ÉTUDE DES PORTIQUES 2D ................................................................................... 29

L’élément fini de flexion plane ....................................................................................................................29 Application aux portiques............................................................................................................................34 Statique des portiques plans simples ..........................................................................................................35 FORMULATION VARIATIONNELLE & ÉCRITURE MATRICIELLE .............................................................................. 39

Formulation intégrale ..................................................................................................................................39 Écriture matricielle du PTV ..........................................................................................................................41 Application aux modèles de l’ingénieur ......................................................................................................42 Exercices : Formulations Variationnelles en physique ................................................................................44 MÉTHODES NUMÉRIQUES DANS LE CADRE DE LA MEF.................................................................................................. 47

Discrétisation du milieu ...............................................................................................................................47 Approximation nodale .................................................................................................................................47 Calcul des matrices élémentaires ................................................................................................................51 Assemblage et conditions aux limites..........................................................................................................56 Exercice ........................................................................................................................................................56 UTILISATION D'UN LOGICIEL ÉLÉMENTS FINIS ................................................................................................................ 59

Création et vérification des données:..........................................................................................................59 Exécution du calcul: .....................................................................................................................................60 Exploitation des résultats: ...........................................................................................................................60 Votre parcours pédagogique .......................................................................................................................61 Organigramme d'un logiciel éléments finis .................................................................................................62 PROCESSUS D’ANALYSE & MODÉLISATION...................................................................................................................... 63

Qu'est-ce qu'un modèle ?............................................................................................................................64 Comment estimer les erreurs de discrétisation ?........................................................................................65 Votre parcours pédagogique .......................................................................................................................66 TEXTE DES TD.............................................................................................................................................................................. 67

TD1 Modèle éléments finis d'un treillis .......................................................................................................67

TD2 Étude d'un portique par les éléments finis ..........................................................................................67 TD3 Modèle EF bidimensionnel en thermique ............................................................................................68 TEXTE DES TP AVEC MEFLAB ................................................................................................................................................ 71

TP1 Prisse en main de MELlab et application aux treillis et portiques........................................................71 TP2: Illustration des méthodes numériques avec MEFlab ..........................................................................72 PRÉSENTATION DE MEFLAB / MEFTAVE.............................................................................................................................. 75

Analyse des scripts éléments finis ...............................................................................................................75 Description des scripts de données .............................................................................................................79 PREMIERS CALCULS AVEC ABAQUS....................................................................................................................................... 83

Prise en main du logiciel ..............................................................................................................................83 Analyse des résultats ...................................................................................................................................85 TP : ÉTUDE DE CAS...................................................................................................................................................................... 85

Pour aller un peu plus loin ...........................................................................................................................87

Mise en équations des barres

7/87

Mise en équations des barres Les quatre champs inconnus de la « MMC » sont : Champs vectoriels u : Déplacements f : Forces

Champs tensoriels

ε : Déformations σ : Contraintes

Les relations entre ces champs peuvent être représentées par la figure suivante

Σ (σ ) Relations géomét riques

σ = D (ε ) < Lois de comport em ent >

E (ε )

Lois de compor tement généralisée

ε = f (u ) Relat ions géométriques

T =σ n F ( f ) < Pr incipe de la dynamique > U (u ) Figure 1 : Relatio ns entre les ch amps de la MMC .

Dans le premier document de cours nous avons établi la loi de comportement généralisée du modèle barre. H1 : déplacement axial u ( M , t ) = u ( x , t ) xo ==> ε xx = u, x H2 : état de contrainte uni axial D'où la définition de l'effort normal

σ xx = Eε xx N = ES u, x

Pour terminer la mise en équations des barres, nous pouvons écrire une des deux formes du principe de la mécanique que vous avez vues en MMC : Le PFD : qui donne un système d'équations aux dérivées partielles (formulation locale). Le PTV : qui est sa forme intégrale ou forme variationnelle et est une forme énergétique globale des équations du mouvement.

Application du PFD Nous allons écrire les équations de Newton f = ma pour une tranche d’épaisseur dx de la barre Le bilan des efforts extérieurs sur cet élément de matière (figure ci-contre) fait apparaitre l'effort normal (torseur des efforts de cohésion) L'équation de résultante dynamique dans la direction x donne : N + dN − N + fdx = ρ Sdx uɺɺ Soit N, x + f = ρ S uɺɺ Compte tenu de la loi de comportement intégrée, l'équation locale est : ∀x ∈ ]0, ℓ[ ρ Suɺɺ − ( ESu, x ) = f

f

N + dN

N dx

x

,x

Les conditions aux limites aux extrémités de la barre peuvent être, u = ud (t ) , en déplacement imposé : ou en force imposée : ESu, x = N d (t ) . Ces 2 conditions permettent de fixer les deux constantes d'intégration en x 7

Mise en équations des barres

8/87

Pour déterminer la réponse dynamique en temps, il faudra se donner les deux conditions initiales: u ( x, 0) = uo ( x) Déformée et vitesse de déformation  initiales de la barre ɺ ɺ u ( x , 0) = u ( x ) o 

Application du PTV Nous allons écrire le principe des travaux virtuels ∀δ u δ W = δ A pour une barre chargée sur sa longueur et à ses extrémités.

f

Fo 0

Fℓ

u ( x, t )

ℓ

ℓ

Le travail virtuel des quantités d’accélération est : δ A = ∫ ρ Suɺɺ δ u dx o

Le travail virtuel des efforts se décompose en travail virtuel des efforts de cohésion et celui des efforts extérieurs soit : ℓ

ℓ

Pour les efforts de cohésion δ Wint = − ∫ σ : δ ε dV = − ∫ ∫ σ xx δε xx dS dx = − ∫ ES ε xx δε xx dx 0S

D

0

ℓ

ε xx = u, x

δ Wint = − ∫ ESu, x δ u, x dx

Soit

o ℓ

δ Wext = ∫ f δ u dx + Fo δ uo + Fℓ δ uℓ

Pour les efforts extérieurs

o

Le PTV conduit à l’équation intégrale suivante : ℓ

∀δ u

ℓ

ℓ

∫ ρ Suɺɺ δ u dx = − ∫ ESu,x δ u,x dx + ∫ f δ u dx + Foδ uo + Fℓδ uℓ o

o

o

C’est la forme variationnelle du problème. Les deux derniers termes correspondent au travail virtuel des efforts appliqués aux extrémités du barreau. Dans le cas ou la condition aux limites porte sur le déplacement, l’effort à l’extrémité est une inconnue du problème. Pour déterminer l'équation du mouvement il faudra tenir compte des conditions aux limites en déplacements. Restreindre le choix des déplacements virtuels à des champs virtuels admissibles, permet d'éliminer l'effort de liaison inconnu de la forme variationnelle. Si u = ud (t ) respectée alors δ u = 0 et le F δ u est éliminé de la Formulation Le travail des efforts de cohésion δ Wint peut s'exprimer à partir de la variation de l'énergie de déformation de la barre δ Wint = −δ Ed ℓ

avec

( )

2 Ed = ∫ σ : ε dV = ∫ ES u, x D

2

dx

o

Équivalence des principes Partons du PFD pour retrouver le PTV. La démarche, présentée de façon générale dans le cours de MMC est utilisée ici dans le cas particulier du modèle barre, sur des équations mono dimensionnelles.

8

Mise en équations des barres L'équation locale

9/87

ρ Suɺɺ − ESu, xx − f = 0

∀x ∈ ]0, ℓ[

ℓ

Est équivalente à

∫ P ( ρ Suɺɺ − ESu, xx − f ) dx = 0

∀P

0

Remarque : si u est une solution approchée du problème cette forme intégrale représente le résidu pondéré de l’équation locale sur le domaine.

Effectuons une intégration par partie du terme en u,xx ℓ

ℓ

ℓ

∫ P ES u, xx dx =  P ESu, x  0 −

∫ P, x ES u,x dx

0

0

Nous obtenons ℓ

∀P

ℓ

ℓ

ℓ

∫ P ρ S uɺɺ dx + ∫ P,x ES u, x dx =  P ESu,x  0 + ∫ P fdx 0

0

0

Introduisons les conditions aux limites en force aux extrémités de la barre Fo Fo = − N ( o, t ) = − ES u, x ( o, t ) Fℓ = + N ( ℓ, t ) = ES u, x ( ℓ, t ) ℓ

∀P

D'où le PTV :

0

N0

ℓ

Fℓ

ℓ Nℓ

ℓ

∫ P ρ S uɺɺ dx + ∫ P,x ES u, x dx = Pℓ Fℓ + Po Fo + ∫ P fdx 0

0

0

Bilan & exercice Le PFD donne un système d’équations aux dérivées partielles (EDP), c'est une formulation locale

équation locale : ∀x ∈ ]0, ℓ[ ρ Suɺɺ − ESu, xx = f PFD  2 conditions aux limites en x = 0 et en x = ℓ Il est utilisé pour rechercher la solution analytique d'un problème. Les conditions initiales ne servent que pour résoudre l'équation différentielle en temps Réponse dynamique d'une structure. Le PTV est la forme variationnelle du problème, 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 numériques du problème. Solutions approchées Vous devez être capable d'écrire ces deux formulations pour un problème donné.

9

Mise en équations des barres

10/87

Tous les exercices de cours sont corrigés sur le site, mais il faut chercher les réponses aux questions avant de consulter le corrigé.

Exercice 1 : Mise en équations d’un barreau en traction Objectifs : Savoir écrire les conditions aux limites pour une barre, Résoudre un problème simple en statique, Pouvoir écrire le PTV et savoir passer du PTV au PFD. 1- Écriture des conditions aux limites. Donnez les différentes conditions aux limites homogènes possibles pour une barre. Donnez les conditions aux limites correspondantes aux trois figures ci-dessous.

F x=ℓ

x=0

k xo

xo

x=0

xo

M

2- Application du PFD. x

Écrire le système d'EDP de ce problème g

ℓ

Intégrer l'équation différentielle en statique Tracer le diagramme de l'effort normal (analyse type RDM)

3- Application du PTV. Pour le problème représenté par la figure ci-dessous k xo x=0 x=ℓ Pour un champ de déplacements virtuels cinématiquement admissible. Donner l’expression du PTV en ne considérant que la barre. Retrouver cette expression en considérant la barre et le ressort.

4- Équivalence des principes. Donner l’expression du PFD et passez au PTV (application directe du cours). Partir du PTV pour retrouver l'équation locale et toutes les CL du problème. Démarche inverse à celle présentée en cours

Si avec la correction vous n'arrivez pas à comprendre la réponse à une question, c'est que des éléments du cours ou des pré-requis vous manquent. Revoyez le cours et n'hésitez pas à poser la question à votre enseignant, il pourra vous aider à résoudre la difficulté. Pour assimiler le cours il faudra traiter des exercices non corrigés.

10

Modèle Éléments finis pour l'étude des Treillis

11/87

Modèle éléments finis pour l'étude des treillis Un treillis est constitué d'éléments barres qui ne travaillent qu'en traction compression. Nous allons utiliser la méthode des éléments finis pour modéliser ces structures. Nous débutons par la présentation de l'élément fini barre, en détaillant le calcul des matrices élémentaires permettant d'exprimer le principe des travaux virtuel sous forme matricielle. Puis nous verrons comment utiliser ces résultats pour modéliser des treillis bidimensionnels

L’élément fini barre uj

ui

Approximation : Considérons un élément de longueur ℓ e

(e) i Le repère local de l'élément est orienté du nœud i vers le nœud j. Les deux variables nodales sont les déplacements notés (ui , u j ) dans la direction x

j

x

Le champ de déplacement sur l'élément sera construit sur une approximation polynomiale à deux paramètres de la forme Approximation linéaire du champ de déplacement  a (t )  Ici les paramètres ai n'ont pas de sens physique u ( x, t ) = a1 + a2 x =< 1 , x >  1  a2 (t )  L'approximation nodale sera construite en identifiant les déplacements nodaux à la valeur de l'approximation soit : u (0, t ) = ui (t ) = a1 en x = 0 a1 = ui  u j − ui en x = ℓ e u (ℓ e , t ) = u j (t ) = a1 + a2ℓ e nous en déduisons  a = 2  ℓe D’où l'approximation nodale u ( x , t ) =< 1 −

x x  ui (t )  , >  ℓ e ℓ e u j (t ) 

Les paramètres ui , u j ont ici un sens physique

Cette approximation sera notée : u = < N > {U e } Les fonctions d’interpolation de l’approximation nodale sont : x N1 ( x) = 1 − ℓe

 N (i ) = 1 vérifie  1  N1 ( j ) = 0

1

x N 2 ( x) = ℓe

 N (i ) = 0 vérifie  2  N2 ( j ) = 1

1

N1

0

1

x/le

N2 0

1

x/le

La notion d'approximation nodale est fondamentale dans la méthode des éléments finis, elle permet d’utiliser des variables qui ont un sens physique, et sur lesquelles nous pourrons directement imposer les valeurs données par les conditions aux limites de type cinématique.

Expression matricielle des énergies élémentaires Nous devons calculer, le travail des quantités d'accélération ainsi que le travail des efforts intérieurs et celui des efforts extérieurs associé à notre élément, en utilisant l'approximation nodale. 11

Modèle Éléments finis pour l'étude des Treillis

12/87

Le travail des quantités d'accélération est : δ Ae =

ℓe

∫ ρ Suɺɺ δ u dx o

Utilisons l’approximation nodale du champ des déplacements u = < N > {U e } Le terme uɺɺ δ u = δ u T . uɺɺ =

{δ U e }

T

< N >T < N > {Uɺɺe }

On peut alors sortir les variables nodales de l'intégrale δ Ae = {δ U e }

T

ℓe

ρ S < N > dx {Uɺɺe }

∫

T

o

δ Ae = {δ U e } [ M e ]{Uɺɺe } avec [ M e ] = T

ℓe

∫

T

ρ S < N > dx

o

Nous venons de définir la matrice masse élémentaire, le calcul de l'intégrale se fait analytiquement, on trouve : A titre d'exercice retrouver par le calcul les coefficients de cette matrice ρ S ℓ e 2 1  [M e ] =   6 1 2  ℓe

δ Ed = −δ Wint

Le travail des efforts intérieurs est : δ Wint = − ∫ ESu, x δ u, x dx o

Pour ce calcul utilisons l'expression de l'énergie de déformation : 2 Ed =

ℓe

∫ ES ( u,x )

2

dx

o

2 Ed = {U e }

T

Utilisons l'approximation nodale

ℓe

∫ < N,x >

T

u,2x = uT, x . u, x

ES < N, x > dx {U e }

o

Soit pour chaque élément 2 Ed = {U e }

T

[ Ke ]{U e } avec [Ke ] =

ES ℓe

 1 − 1 − 1 1   

ℓ

Le travail des efforts extérieurs est : δ Wext = ∫ f δ u dx + Fie δ ui + F je δ u j

f

Fie

(e)

i

o

F je j

Pour la densité de charge f appliquée sur l'élément nous devons tenir compte de l'approximation nodale pour exprimer le travail virtuel.

δ We =

ℓe

∫

f ( x ) δ u dx

Compte tenu de l’approximation δ We = {δ U e }

T

0

ℓe

∫ < N ( x) >

T

f ( x ) dx

0

Nous pouvons effectuer le calcul de l'intégrale si la répartition de charge sur l'élément est connue. ℓe    T Pour une charge f = Cte . On trouve δ We = {δ U e } {φe } , avec {φe } = f  2  ℓe   2  Ce calcul permet de calculer les charges nodales équivalentes au sens de l’approximation à une charge volumique réelle appliquée à la structure f i

PTV

f ℓe

j

i

Charge réelle

f ℓe

2 j

Charge nodale équivalente

12

2

x

Modèle Éléments finis pour l'étude des Treillis

13/87

Objectif : Déterminer une approximation des premières fréquences de résonnance de la barre avec un modèle élément fini.

Exemple

Tester la méthode avec un modèle à 1 élément

ℓ

xo

Généraliser à n éléments (maillage régulier) Modèle à 1 élément fini 2

ES  1 −1 , et [ K ] = sur < u1 u2 >   L  −1 1  6 1 2  Les vibrations de la barre sont modélisées par un système à 1 DDL

uɺɺ2 +

ES u2 = 0 L

1

[M ] =

xo

La condition : u1 = 0

ρ SL 3

ρ SL  2 1 

ω1 = 3

ES

ρ SL2

ES

≅ 1, 732

ρ SL2

À comparer à la solution analytique ωi =

π

ES

2

ρ SL

2

≅ 1,571

ES

ρ SL2

L’erreur d’approximation 10% sur la première pulsation propre est importante car on utilise une approximation linéaire pour une fonction sinusoïdale.

Modèle à N éléments finis 1

2

3

n

n+1

Pour tout élément [ K e ] =

ES  1 L / n  −1

−1   1

et [ M e ] =

ρ SL  2 6n  1

1

 2

Lorsque l'on somme les énergies de chaque élément pour obtenir l'énergie de la structure les matrices élémentaires s'emboitent les unes avec les autres, en effet ES 2 Pour l'élément 1 : 2 Ed 1 = (u1 − 2u1u2 + u22 ) ℓ ES 2 Pour l'élément 2 : 2 Ed 2 = (u2 − 2u2u3 + u32 ) ℓ ES 2 Soit pour les deux éléments 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  T Avec [ K1+ 2 ] = −1 1 + 1 −1 sur {U } = < u1 u2  ℓ  0 −1 1 

u3 >

c'est l'assemblage. En généralisant aux n éléments on obtient une matrice (n + 1, n + 1) , mais il faut tenir compte de la condition d'encastrement du premier nœud, tous les termes u1 sont nuls, la matrice assemblée réduite est une matrice carrée de dimension n  2 −1   −1 2 −1      −1 2 −1 nES   Matrice raideur assemblée réduite ( u1 = 0 ): −1 \ \ [K ] =   L   \ \ \   \ 2 −1   −1 1   De même pour l'énergie cinétique 13

Modèle Éléments finis pour l'étude des Treillis

14/87

4 1  1 4 1     1 4 1  ρ SL   Matrice masse 1 \ \ [M ] =   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  ρ SL2 ω ≅ 5, 63 ES  2 ρ SL2   ES ω1 ≅ 1, 589 ρ SL2   ES ES  Pour n=3 ωanal = 1,571 4, 712 et 7,854 ω2 ≅ 5,196 2 ρ SL ρ SL2   ω ≅ 9, 426 ES  3 ρ SL2 La convergence est lente (éléments de degré 1) Avec la matrice modale calculée 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 très fin pour approcher la déformée modale des modes supérieurs.

Modèle éléments finis d'un treillis La démarche générale pour traiter un problème par une modélisation éléments finis est la suivante : Analyse du problème choix de discrétisation (définition des inconnues) Boucle sur les éléments Calcul des matrices élémentaires et charges généralisées Assemblage & C. Limites équation matricielle à résoudre Résolution déformée de la structure et efforts aux appuis 14

Modèle Éléments finis pour l'étude des Treillis Post-traitement

15/87

contraintes dans les barres et efforts aux nœuds.

Pour un treillis bidimensionnels

j

(e)

α

yo i

Soit une barre formant un angle α avec l’axe xo du repère global.

uj

vj

Pour effectuer l'assemblage nous devons exprimer le déplacement axial u en fonction de ses composantes sur la base globale ( u , v ). u  u  u =< cos α sin α >   =< Cα Sα >   v  v 

uj

x

xo

Ce n'est que la première ligne d'un changement de base classique.

Appliquons ce changement de base aux nœuds de l’élément  ui  Cα  = u j   0

Sα 0

0 Cα

 ui    0   vi    Sα  u j   v j 

u   Cα  = v   − Sα

Reportons ce changement de base dans l'expression de l'énergie de déformation. {U e }

T

T

 ui  v   i  C 2 Ed =    α u  j  0 v j   

Sα

0

0 Cα

T

0  ES  1 −1 Cα Sα  ℓ e  −1 1   0

Sα

0

0 Cα

[A] − [A]  − [ A] [ A] 

ℓe

avec

 Cα2 Cα S α

[A] = 

ES  1 −1   {U e } ℓ e  −1 1 

 ui    0   vi    Sα  u j  v j   

Nous en déduisons l’expression de la matrice raideur élémentaire sur les variables < ui

[K e ] = ES 

Sα  u    Cα   v 

vi

uj

vj >

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 numérique.

Assemblage et résolution Pour chaque élément de la structure nous avons :  Fie  ∀e [ M e ]{Uɺɺe } + [ K e ]{U e } = {φe } +    F je 

Les

Fie sont les efforts du nœud i sur les

éléments e (effort appliqué à l'élément)

L'assemblage consiste à sommer les énergies élémentaires

ℓe

∫

=∑∫

D

e

D'où le signe moins pour avoir les efforts des éléments sur le nœud.

0

Pour les efforts nodaux l'équilibre d'un nœud quelconque donne Fi − ∑ Fie = 0 e

Les Fi représentent les efforts extérieurs appliqués aux nœuds de la structure. Se sont soit des charges données soit des efforts aux appuis (conditions cinématiques) qui sont des inconnues du problème. L'assemblage consiste à se donner un ordre de rangement des variables nodales dans le vecteur des inconnues globales du système. En pratique (à la main) nous utilisons l'ordre lexicographique pour simplifier l'écriture. La machine (calculateur) utilisera sa propre numérotation pour optimiser la vitesse de traitement 15

Modèle Éléments finis pour l'étude des Treillis

16/87

et la taille mémoire utile en fonction des algorithmes de résolution qu'il utilisera pour traiter les équations, ces opérations sont transparentes pour l'utilisateur. En statique nous utiliserons une décomposition du système matriciel en déplacements inconnus (nœuds ou les charges sont données) et déplacements imposés (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 déplacements nodaux inconnus:

{U i } = [ K11 ] {{Fd } − [ K12 ]{U d }} −1

C’est le système réduit

En reportant dans le second nous obtenons le vecteur des efforts de liaison inconnus: {Fi } =  K 22 − K 21K11−1 K12  {U d } +  K 21 K11−1  {Fd } Dans les exercices très souvent les déplacements sont imposés nuls, ce qui simplifie les écritures et les

{U i } = [ K11 ] {Fd } −1

calculs

puis

{Fi } = [ K 21 ]{U i }

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 éléments. Nous utilisons la loi de comportement intégrée : ES N = ES u, x = ES < N, x > {U e } = (u j − ui ) = Cte ℓe L'état de contrainte est constant dans chaque élément fini

En statique, pour des treillis chargés aux nœuds le modèle éléments finis ne nécessite qu'un élément par barre du treillis, il donnera la solution analytique exacte. Ce n'est évidemment pas le cas ni pour une colonne chargée par son poids propre, ni pour les problèmes de dynamique, ou la solution exacte se décompose sur une base de fonctions sinusoïdale (cf chapitre sur les solutions analytiques pour les barres) . Dans le cas bidimensionnel, l’état de contrainte sur un élément 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 nœuds à 2 variables par nœuds (ui , u j ) les déplacements du nœud dans le plan.

a 2

yo

{U }

T

xo

Modèle à 6 degrés de liberté = {u1 v1 u2

u3

a 1

a 2

u3

v3 }

vecteur des déplacements nodaux

v3 3

v2

2

u2

Les conditions aux limites : u = 0 Appui au nœud 1 :  1  v1 = 0 Appui glissant au nœud 2 : v2 = 0

X soit deux efforts inconnus :  1  Y1 soit un effort inconnu : Y2

Le travail virtuel des efforts donnés et inconnus appliqués à la structure conduit à l’expression du vecteur des forces nodales :

16

Modèle Éléments finis pour l'étude des Treillis 3

Y2

X1 a 2

1

(2)

(3) (1) a 2

1

2

xo v3 3

0

v3 }

u3

Calculons la matrice raideur [K] de cette structure. u1 2 1 Pour l’élément 1 (1,2) a 2 ES  1 −1 K1 =   sur {u1 u2 } a 2  −1 1  1 1  ES  1 1 K2 = Pour l’élément 2 (1,3) 2 a  −1 −1   −1 −1

u1

Pour l’élément 3 (2,3)

v3

u3

3

0 u3

Les équations 1, 2 et 4 nous donnerons les efforts aux appuis en fonction de ces déplacements.

α = 45°

1

u2

Les équations 3,5, et 6 nous permettent de déterminer le champ de déplacement de la structure (sa déformation).

a

v1

0}

Nous avons donc 6 inconnues pour 6 équations

3

yo

pour

F

Y1 0 Y2

2

 0   X1   0  Y     1 u   0  [ K ]  02  =  Y     2 u3   F       v3   0 

a

{F } = { X 1 T {U } = {0 T

F

a

Y1

17/87

a

v2

α = 135°

2

u2

−1 −1 −1 −1 1 1  1 1 sur {u1

v1

u3

v3 }

 1 −1 −1 1    ES  −1 1 1 −1 K3 = 2a  −1 1 1 −1    1 −1 −1 1  sur {u2

v2

u3

v3 }

u2

L’énergie de déformation totale de la structure est la somme des énergies de déformation de chaque élément, l’assemblage des matrices consiste à ranger chaque terme dans une matrice globale définie sur le vecteur {U } = {u1 v1 u2 T

D’où la matrice globale 1 + 2   1 ES  [K ] =  − 2 2a  0   −1  −1

1 1

v2

u3

− 2 0

0 1+ 2 0 −1 −1 −1 −1 1

v3 } −1   −1  −1 −1 1  1 1 −1   1 1 + 1 1 − 1 −1 1 − 1 1 + 1 0 0

−1 −1

les termes de la matrice K1 sont en bleu la matrice K2 sont en rouge la matrice K3 sont en vert

L’équation matricielle [ K ]{U } = { F } à résoudre est la suivante : Les 3 équations donnant les déplacements nodaux sont en bleu Celles permettant de calculer les efforts sont en rouge

17

Modèle Éléments finis pour l'étude des Treillis 1 + 2   1 ES  − 2  2a  0   −1  −1

1

− 2

1

0

0 0

1+ 2 −1

−1 −1

−1 1

18/87

−1 −1  0   X1   0 −1 −1  0   Y1  −1 −1 1  u2   0   =  1 1 −1  0   Y2   1 2 0  u3   F      −1 0 2   v3   0  0

Pour résoudre nous tenons compte des conditions aux limites cinématique ce qui réduit le système à 3 équations. Ce système réduit est : F a  u2 = ES 2 1 + 2 −1 1  u2   0   Allure de la déformée     F ES  1  u3 F 2 0  u3  =  F  a (1 + ) u3 =  −1 2a  ES 2 2  0 2   v3   0  u2 2  1  1 F a v3 = − ES 2 2  Nous pouvons alors calculer les efforts aux appuis ES   X1 = 2a − 2 u2 − (u3 + v3 )  X1 = − F  ES   ( u3 + v3 ) Y1 = − F / 2 Y1 = − 2a  Y = + F / 2  2 ES  Y = − u + u − v ( ) 2 3 3  2 2a 

(

)

Post-traitement Calculons l'effort normal dans les éléments ES  ( u2 ) = F / 2  N1 = a 2  ES  N = ES u, x ==>  N 2 = ( u3 + v3 ) = F / 2 a 2   ES ( u3 − u2 − v3 ) = − F / 2  N3 = − a 2 

L’équilibre global de la structure est vérifié F

F /2

−F

2

1

−F / 2

L’équilibre de chaque nœud est vérifié

−F

N2 N 2 (2)

F N3 (3) N3 F / 2

N1 (1) N1 −F / 2

Remarques Tous les calculs sont systématiques et la démarche suivie sera toujours la même en statique. Facilité de programmation de ce type de solution Seule l’analyse, du problème et des résultats, reste à la charge de l’ingénieur. La matrice raideur du système réduit était inversible " det( K ) ≠ 0 " car les conditions aux limites en déplacement bloquaient tous les modes rigides de la structure. Problème statique bien posé Les efforts calculés aux appuis équilibrent parfaitement le chargement. Les résidus d'équilibre sont nuls, car nous travaillons sur la solution analytique de l'équation matricielle. Dans le cas d’une résolution numérique ces résidus doivent tendent vers zéro (erreur numérique). 18

Modèle Éléments finis pour l'étude des Treillis

19/87

Les contraintes calculées sur les éléments équilibrent de façon exacte (aux résidus près) les charges nodales. Ceci est vrai dans ce cas particulier « calcul statique d’un treillis chargé aux nœuds » car l’approximation utilisée représente le champ exact de la solution analytique « effort normal constant dans chaque élément de la structure ». Erreur de discrétisation qui est nulle En post – traitement il est possible d’isoler un à un chaque élément de la structure pour écrire l’équation matricielle de l’équilibre de l’élément. Ces calculs permettent de déterminer les efforts internes aux nœuds de la structure, nous en donnons des exemples dans les exercices de cours.

Exercices Les exercices de cours sont corrigés sur le site, il faut chercher les réponses avant de consulter le corrigé.

Exercice 8 : Modélisation EF d'une colonne sous son poids propre Objectifs :

x

g

Notion d'erreur de discrétisation, et analyse des résultats EF. Étude de convergence en affinant le maillage. 1- Établir l'équation matricielle d'un modèle à un élément fini Analyser les résultats aux nœuds (déplacement & efforts) Tracer les résultats sur l'élément (déplacement & efforts)

ℓ = 6h

2- Construire le modèle utilisant deux éléments finis identiques Analyser les résultats aux nœuds (déplacement & efforts) Tracer les résultats sur les éléments (déplacement & efforts) 3- Modèle à 3 éléments, pour affiner le maillage dans la zone la plus contrainte nous utilisons 3 éléments de longueur h, 2h, et 3h. Déduire des calculs précédents l'équation matricielle du modèle. Tracer les résultats sur les éléments (déplacement & efforts)

Pour améliorer la solution éléments finis nous avons augmenté le nombre d'éléments et densifié le maillage dans la zone la plus chargée. Cette méthode dite « h convergence » demande en général un nombre élevé d'éléments finis. La figure suivante présente les résultats d’un modèle éléments finis en contraintes planes. Pour quantifier l’erreur relative à cette discrétisation, la démarche est identique à celle de cet exercice, elle est basée sur l’analyse de la discontinuité du champ des contraintes entre deux éléments adjacents. en MPa

Dans cette section le diagramme des contraintes est le suivant σ VM 145

solution cherchée Discontinuité solution éléments finis

constante par morceau

83 62

19

L’erreur est beaucoup trop importante. Ce modèle n’est pas satisfaisant, il faut affiner le maillage

Modèle Éléments finis pour l'étude des Treillis

20/87

Exercice 9 : Utilisation d'éléments finis de degré deux Objectifs : Amélioration de la convergence en augmentant le degré de l'approximation x

Pour approximation polynomiale du second degré de la forme. u ( x, t ) = a0 + a1x + a2 x 2 = < 1 , x, x 2 > {ai }

g

ℓ = 6h

Nous devons utiliser un élément fini à 3 nœuds, et construire l'approximation nodale sur < u1 u2 u3 >

u1 1

u3 x

u2 2 l

3

1- Déterminer les fonctions d'interpolation nodale de cet élément de degré 2. 2- Calculer la matrice raideur élémentaire correspondante. 3- Calculer la force généralisée due au poids propre de cet élément. 4- Déduire des calculs précédents les résultats avec un modèle à 1 élément fini de la colonne. Pour améliorer la solution éléments finis nous avons augmenté le degré de l’approximation élémentaire. Cette méthode dite « p convergence » est en général beaucoup plus rapide, elle nécessite moins d’éléments finis. Les figures suivantes illustrent les deux choix d’améliorations possibles d’un modèle numérique dont l’erreur liée à un maillage grossier est trop importante.

En utilisant des éléments de degré 2 "p" convergence

En affinant le maillage localement "h" convergence

Exercice 10 : Étude d’un treillis symétrique de trois barres Objectifs : Techniques de mise en œuvre de la méthode des éléments finis, changement de base, assemblage, résolution, calcul des efforts, et vérification des équations d'équilibre. Considérons le treillis de trois barres ci-contre h Modélisation. Préciser la numérotation de vos éléments et de vos nœuds. h Définissez vos vecteurs globaux : {U } vecteur des déplacements nodaux (ui , vi ) F

{FD } vecteur force généralisé associé aux efforts donnés {FI } vecteur force généralisé associé aux efforts inconnus

Calcul de la matrice raideur Exprimer la matrice raideur de chaque élément sur ses variables nodales. 20

Modèle Éléments finis pour l'étude des Treillis

21/87

En déduire l'expression de la matrice raideur assemblée complète. Extraire la matrice raideur réduite. Résolution statique - Efforts aux appuis Déterminer la déformée statique, et représenter l'allure de la déformée. Calculer les efforts aux appuis, et vérifier l'équilibre global de la structure. Post traitement Calculer les contraintes sur chaque élément, puis vérifier l'équilibre du nœud qui est chargé. Isoler une des barres à 45° de la structure, et calculer les efforts extérieurs sur cet élément. Retrouver les résultats précédents. Utiliser la symétrie Préciser le nouveau maillage en tenant compte de la symétrie. Calculer la matrice raideur réduite et retrouver la solution en déplacement.

Exercice 11 : Modélisation EF du treillis de l'exercice 7 Objectifs :

Élimination des mouvements d'ensemble, et prise en compte des symétries. F

Nous cherchons la réponse statique du treillis de 6 barres, en utilisant un modèle éléments finis.

ℓ

ℓ

1. Pourquoi ce problème est-il mal posé ? 2. Définir un modèle EF de la structure qui soit bien posé. 3. Simplifier le modèle en tenant compte des symétries.

F

4. Déterminer la déformé statique et les efforts dans les barres Validez vos résultats.

Pour assimiler le cours il faut traiter des exercices non corrigés.

21

Modèle Éléments finis pour l'étude des Treillis

22/87

22

Mise en équations des poutres

23/87

Mise en équations des poutres en flexion plane Les hypothèses du modèle poutre et la loi de comportement généralisée sont présentées en détail dans le chapitre sur la modélisation des poutres que vous pouvez consulter sur le site. H1 : modèle de Bernoulli u ( M ) = u (G ) + θ Λ GM et θ = rot (u (G ))

− y v, x    Soit en flexion plane θ = v , x zo d'où u ( M , t ) =  v  ==> ε xx = − y v , x 2  0    H2 : état de contrainte uni axial σ xx = Eε xx Le calcul du torseur des efforts de cohésion sur une section droite permet de définir le moment de flexion. M f = EI v , x 2 I est le moment quadratique de la section droite de la poutre.

Application du PFD Nous allons écrire les équations de Newton f = ma pour une tranche d’épaisseur dx de la poutre Le bilan des efforts extérieurs sur cet élément de matière (figure ci-contre) fait apparaitre le torseur des efforts de cohésion, l'effort tranchant est associé aux contraintes de cisaillement qui s'opposent au glissement des sections. Les équations de résultante et de moment dynamique sont :

f

Mf

T + dT − T + fdx = ρ Svɺɺ dx   dx dx (T + dT ) 2 + M f + dM f − M f + T 2 ≅ 0

Soit

Mf + dMf

T

dx

x T + dT

On néglige le moment dynamique de rotation des sections.

∀x ∈ ]0, ℓ[ ρ Svɺɺ + M f , xx = f T = −M f ,x

Compte tenu de la loi de comportement intégrée, l'équation locale est : ∀x ∈ ]0, ℓ[

ρ Svɺɺ + EIv, x4 = f

Les conditions aux limites aux extrémités de la poutre peuvent être, en déplacement imposé :

v = vd (t ) ou θ = θ d (t )

ou en force imposée :

T = Td (t ) ou Mf = Mf d (t )

Ces 4 conditions permettent de fixer les quatre constantes d'intégration en x Pour déterminer la réponse dynamique en temps, il faudra se donner les deux conditions initiales:

v ( x , 0) = vo ( x )  vɺ( x , 0) = vɺo ( x )

Déformée et vitesse de déformation initiales de la poutre

23

Mise en équations des poutres

24/87

Application du PTV

y

Nous allons écrire le principe des travaux virtuels ∀δ u δ W = δ A pour une poutre chargée sur sa longueur et à ses extrémités.

f

Fo

Fℓ

Mo

ℓ

0

Mℓ

x

ℓ

Le travail virtuel des quantités d’accélération est : δ A = ∫ ρ Svɺɺ δ v dx o

On néglige le moment dynamique de rotation des sections.

Le travail virtuel des efforts se décompose en travail virtuel des efforts de cohésion et celui des efforts extérieurs soit : ℓ

Pour les efforts de cohésion δ Wint = − ∫ σ : δ ε dV = − ∫ ∫ σ xx δε xx dS dx 0S

D ℓ

δ Wint = − ∫ EIv, xx δ v, xx dx

Soit

0

ε xx = − y v, x

2

σ xx = Eε xx

ℓ

δ Wext = ∫ f δ v dx + Foδ vo + Fℓδ vℓ + M oδθo + M ℓδθℓ

Pour les efforts extérieurs

o

Le PTV conduit à l’équation intégrale suivante : ℓ

∀δ v

ℓ

ℓ

∫ ρ Svɺɺ δ v dx = −∫ EIv, xx δ v,xx dx + ∫ f δ v dx + Foδ vo + Fℓδ vℓ + M oδθo + M ℓδθℓ o

o

o

C’est la forme variationnelle du problème. Les quatre derniers termes correspondent au travail virtuel des efforts appliqués aux extrémités de la poutre. Dans le cas ou les conditions aux limites portent sur les déplacements, les efforts de liaison sont des inconnues du problème. Pour déterminer l'équation du mouvement il faudra tenir compte des conditions aux limites en déplacements. Restreindre le choix des déplacements virtuels à des champs virtuels admissibles, permet d'éliminer les efforts de liaison inconnus de la forme variationnelle. Si v = vd ( t ) respectée alors δ v = 0 et le F δ v est éliminé de la Formulation Si θ = θ d (t ) respectée alors δθ = 0 et le M δθ est éliminé de la Formulation Le travail des efforts de cohésion δ Wint peut s'exprimer à partir de la variation de l'énergie de déformation de la poutre δ Wint = −δ Ed ℓ

avec 2 Ed = ∫ σ : ε dV = ∫ EI ( v, xx ) dx D

2

0

Équivalence des principes Dans le chapitre sur la mise en équation des barres nous Partions du PFD pour retrouver le PTV. Nous allons ici faire la démarche inverse. Partons du PTV et transformons l’équation intégrale pour retrouver le PFD (équation locale) et les conditions aux limites du problème.

24

Mise en équations des poutres

25/87 ℓ

Effectuons deux intégrations par partie du terme ∫ EIv, xx δ v, xx dx o ℓ

Fait apparaître les conditions aux limites en rotation et moment

ℓ

ℓ

∫ EIv, xx δ v, xx dx = δ v,x EI v, x2  0 − ∫ δ v,x EI v, x3 dx o

Fait apparaître les conditions aux limites en flèche et force

0

ℓ

ℓ

ℓ

ℓ

∫ EIv, xx δ v, xx dx = δ v, x EI v, x2 0 − δ v EI v,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 oδθo + M ℓδθℓ o

o

o

∫ δ v ( ρ Svɺɺ + EIv, x ℓ

En regroupant les termes : ∀δ v

4

−f

o

Le choix δ v ≠ 0 de sur

]0, ℓ[

(

, x3

)

x =0

=0

ℓ

o

, x2 o

ℓ

, x3 ℓ

(

M o = − EI v

Pour δ vℓ ≠ 0

Fℓ = − EI v

Pour (δ v, x ) ≠ 0

M ℓ = EI v

(

(

)

, x 2 x =0

) )

, x3 x =ℓ

, x 2 x =ℓ

x=0

Cette condition tient compte de l’orientation de la normale extérieure au domaine

Fo = −To

Pour (δ v, x ) ≠ 0

ℓ

o

]0, ℓ] , nous donne la condition aux limites en force en

De la même façon o

) + δθ ( M + EI v ) ) + δθ ( M − EI v )

nous donne l’équation locale : ρ Svɺɺ + EIv, x 4 − f = 0

Le choix de δ vo ≠ 0 et δ v = 0 sur Fo − EI v

)

( (

δ v F − EI v  o , x3 dx =   δ vℓ F + EI v 3 ,x 

= −M f o

= Tℓ = M fℓ

Vous devez être capable de faire la démonstration dans les deux sens PTV ⇔PFD.

Bilan & exercice Le PFD donne un système d’équations aux dérivées partielles (EDP), c'est une formulation locale équation locale : ∀x ∈ ]0, ℓ[ ρ Svɺɺ + EIv, x4 = f PFD :  4 conditions aux limites : 2 en x = 0 et 2 en x = ℓ Il est utilisé pour rechercher la solution analytique d'un problème. Les conditions initiales ne servent que pour résoudre l'équation différentielle en temps Réponse dynamique d'une structure. Le PTV est la forme variationnelle du problème, 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 numériques du problème. Solutions approchées Vous devez être capable d'écrire ces deux formulations pour un problème donné. 25

Mise en équations des poutres

26/87

Les exercices de cours sont corrigés sur le site, il faut chercher les réponses avant de consulter le corrigé.

Exercice 3 : Mise en équations d’une poutre en flexion plane Objectifs : Savoir écrire les conditions aux limites pour une poutre, Résoudre un problème simple en statique, Pouvoir écrire le PTV et savoir passer du PTV au PFD. Les hypothèses sont celles des poutres longues en petites déformations et petits mouvements. Le matériau est supposé homogène isotrope élastique 1- Écriture des conditions aux limites Exprimer les 4 conditions aux limites homogènes suivantes : extrémité libre

encastrée appui simple Exprimer les 3 conditions aux limites non homogènes suivantes : M, I F

appui glissant

k

2- Mise en équations par le PFD Donnez le système d’équations correspondant au problème ci-dessous Mf A

Pb de flexion

B

g Déterminer la solution analytique en statique, pour M = 0 . Calculer la déformée de la poutre Déterminer le diagramme du moment de flexion 3- Application du PTV. Pour le problème représenté par la figure ci-dessous, donner l’expression du PTV correspondant à des champs de déplacements virtuels cinématiquement admissibles.

yo

g (ρ , E, I, S)

Γt

M xo Peut-on transformer le PTV pour retrouver l'équation locale et les conditions aux limites.

Exercice 4 : Mise en équations poutre "encastrée-masse en bout" Objectifs : Écrire les Conditions aux limites et les EDP du problème, écrire la formulation variationnelle du problème, savoir passer de l'une à l'autre de ces deux formulations. Intéressons-nous aux vibrations dans son plan principal de la poutre ℓ M droite de longueur ℓ représentée par la figure ci contre. La masse M en bout de poutre est supposée ponctuelle Mise en équations Écrivez le système d’équations différentielles régissant ce problème. Écrivez la forme intégrale associée au PTV, pour un champ virtuel quelconque. Formulation variationnelle En partant du système d'EDP du problème retrouver la forme intégrale du PTV. 26

Mise en équations des poutres

27/87

Exercice 5 : Mise en équations d'un arbre en torsion Objectifs : Écrire les Conditions aux limites et les EDP du problème, écrire la formulation variationnelle du problème, savoir passer de l'une à l'autre de ces deux formulations. Intéressons-nous aux vibrations en torsion de l'arbre de longueur ℓ auquel ΓM est appliqué un couple moteur via un engrenage d'inertie en rotation I. Le ℓ problème est représenté par la figure ci contre. Mise en équations Écrivez le système d’équations différentielles régissant ce problème. Écrivez la forme intégrale associée au PTV, pour un champ virtuel quelconque. Retrouver les équations locales en partant du PTV.

27

Mise en équations des poutres

28/87

28

Éléments finis pour l'étude des portiques

29/87

Modèle éléments finis pour l'étude des portiques 2D Un portique bidimensionnel est constitué d'éléments poutres qui travaillent qu'en traction & flexion. Nous allons utiliser la méthode des éléments finis pour modéliser ces structures. Nous présentons l'élément fini poutre de flexion plane, en détaillant le principe de calcul des matrices élémentaires conduisant à la forme matricielle du principe des travaux virtuel. Nous verrons alors comment utiliser ces résultats pour modéliser des portiques bidimensionnels

L’élément fini de flexion plane y vi

Approximation : L’élément fini « poutre » utilise comme variables nodales la flèche et sa dérivée première (rotation de la section droite), il fait partie de la famille des éléments de type l'Hermite.

i

vj

θi ℓ

θj

j

x

Considérons un élément de longueur ℓ Le repère local orthonormé lié à l'élément, a pour direction x l'axe de la poutre orienté de i vers j, et pour direction y un vecteur du plan principal d'inertie de la section droite. Les quatre variables nodales sont les déplacements notés < vi (t ) θi (t ) v j (t ) θ j (t ) > Pour identifier nos quatre variables nodales, nous utilisons une approximation polynomiale cubique de la forme :  a1 (t )    Approximation de degré 3 h 2 3  a2 (t )  v ( x, t ) =< 1 x x x >  à 4 variables  a3 (t )  a4 (t )  Par identification des variables nodales avec l’approximation de la flèche et de la rotation aux noeuds, nous obtenons la relation matricielle suivante : h  vi (t )   v ( o, t )  1  θ (t )   h  i  θ ( o, t )  0  v (t )  =  h   j   v ( ℓ , t )  1 θ j (t )   h    θ (ℓ, t )  0

0   a (t )  1 1 0 0  a (t )  2   ℓ ℓ 2 ℓ3   a3 (t )   1 2ℓ 3ℓ 2  a4 (t )  0

0

Inversons cette relation et reportons le résultat dans l'expression de l'approximation, nous obtenons :  vi (t )  θ (t )   i  h v ( x, t ) = < N >e {U e } = < N1 N 2 N3 N 4 >  v (t )   j  θ j (t )    Avec les fonctions d'interpolation suivantes :

29

Éléments finis pour l'étude des portiques

30/87

 N 1 ( s) = 1 − 3s 2 + 2 s 3 x où s =  2 3 ℓ  N 3 ( s) = 3s − 2 s N1 et N3 représentent la déformée d'une poutre bi - encastrée pour laquelle on impose un déplacement unité à une des deux extrémités  N 2 ( s ) = ℓ( s − 2 s 2 + s 3 )   N 4 ( s ) = ℓ(− s 2 + s 3 ) N2 et N4 représentent la déformée d'une poutre encastrée à une extrémité. Pour laquelle on impose une rotation unité à l'autre extrémité.

Principe des travaux virtuels

N3 ( s )

N1 ( s )

1

s=

1

0

x ℓ

N2 (s) 1

1 0

N4 ( s)

1

s

On néglige le moment dynamique de rotation des sections.

L  δ A =  ∫o ρ Svɺɺ δ v dx  L  2 Partons de ∀δ u δ W = δ A avec δ Wint = −δ Ed avec 2 Ed = ∫ EI ( v, xx ) dx o  L  δ Wext = ∫ f δ v dx + Foδ vo + FLδ vL + M oδθo + M Lδθ L  o

La poutre pouvant être modélisée par plusieurs éléments finis nous calculerons les énergies sur chaque élément puisque l'approximation nodale est une approximation élémentaire.

Matrice raideur élémentaire ℓ

L'énergie de déformation associée à notre élément est 2 Ed = ∫ EI ( v, xx ) dx 2

o

Utilisons l’approximation nodale du champ des déplacements v, xx = < N, xx > {U e } Le terme ( v, xx ) = v, xxT v, xx = {U e } < N, xx >T < N , xx > {U e } 2

T

En reportant dans l'énergie de déformation, pour chaque élément nous obtenons l'expression matricielle de l'énergie de déformation élémentaire : 2 Ed = {U e }

T

ℓ

∫ [ N,xx ]

T

EI [ N, xx ] dx {U e }

0 ℓ

La matrice raideur associée est [ K e ] = ∫ [ B ]T EI [ B] dx 0

avec [ B ] = < N , xx > =
ℓ

 A titre d’exercice calculez le terme −6 ℓ 2 ℓ  (1,2) de cette matrice  12 −6ℓ  2 −6ℓ 4ℓ  sur i i j j

−12

6ℓ

2

Cette matrice n'est pas adimensionnelle car v et θ n'ont pas la même dimension. Pour que les coefficients de la matrice soient adimensionnels il faut travailler sur les variables v et ℓθ 30

Éléments finis pour l'étude des portiques

31/87

6  12  EI [ Ke ] = 3z −612 −46 ℓ  2  6

−12

−6 12

−6

 2   −6  4  sur i i j j 6

Cette expression peut vous permettre de simplifier vos calculs numériques.

Matrice masse élémentaire

On peut calculer ce terme à partir de l'énergie cinétique.

ℓ

Le travail virtuel des quantités d'accélération : δ A = ∫ ρ Svɺɺ δ v dx o

De la même façon en utilisant l’approximation nodale du champ des déplacements, l'expression matricielle pour un élément est : δ Ae = {δ U e }

T

ℓ

∫< N >

T

ρ S < N > dx {Uɺɺe }

0

11 9  13 −13 420  70 210  35  1 13 1   11 − 420 140 105  D'où la matrice masse élémentaire est [ M e ] = ρ S ℓ  210 13 13 11 −  9 70 420 35 210     −13 420 −1140 −11 210 1105  sur

Vecteur force généralisée élémentaire

y

f

Soit un élément poutre chargé par une densité linéique d'efforts transversaux f f s'exprime en N/L Le travail virtuel de ces efforts est ℓ

δ W f = ∫ f . δ v dx = {δ U e }

T

0

ℓ

∫[N ]

T

i

ℓ

j

x

f dx

0

Il faut se donner la fonction f

Pour le champ de pesanteur f

Pour une densité de charge uniforme nous obtenons : { Fd }e =

ℓe

∫ 0

La prise en compte d'une charge répartie sur un élément ne consiste pas à appliquer simplement des efforts fl/2 aux noeuds.

.

 ℓ   2   2   ℓ   12  f < N ( x ) >T dx = f    ℓ   2   2 − ℓ   12  (x)

= − ρ gS

M 1 = fℓ 2 / 12

f

ϕ1 = fℓ / 2

M 2 = − fℓ 2 / 12 ϕ 2 = fℓ / 2

PTV 1

2

1

Charge réelle f=Cte

2

Charge nodale équivalente

Vecteur force généralisée nodale Lorsqu'un chargement est appliqué sur un nœud de la structure le travail virtuel des charges s'exprime directement sur les variables nodales concernées : δ Wext = Fiδ vi + M iδθi Les valeurs de Fi et M i se mettent directement dans le vecteur des charges extérieures 31

Éléments finis pour l'étude des portiques Exemple

Objectif : Déterminer la réponse statique de la poutre avec un modèle élément fini. F

y x

ℓ

Y1

v2

ℓ

Modèle à 1 élément fini Ce modèle comporte 4 variables : X T =< v1 , θ 1 , v 2 , θ 2 > Les conditions aux limites : (v1 , θ1 ) = (0, 0)

M1 1

32/87

2 déplacements inconnus : X IT =< 0, 0, v2 , θ 2 >

θ2

2 efforts inconnus : FI T =< Y1 , M 1 , 0, 0 >

2

La charge conduit à : FDT =< 0, 0, − F , 0 >

 12 6ℓ EI  6ℓ 4ℓ 2 Le PTV appliqué à l'élément nous donne l'équation matricielle 3  ℓ  −12 −6ℓ  6ℓ 2ℓ 2

  0   Y1     2  −6 ℓ 2 ℓ   0   M 1   =   12 −6ℓ  v2   − F  2 −6ℓ 4ℓ  θ   0   2  

−12

6ℓ

 v2  EI  12 −6ℓ   v2  − F  F ℓ3 1/ 3  Les équations donnant la déformée sont : 3   =   ==>  =−   2 EI 1/ 2  ℓ  −6ℓ 4ℓ  θ 2   0  ℓθ 2  C'est la solution exacte de la RDM Les équations donnant les efforts à l'encastrement sont : Y  EI  −12 6ℓ   v2   Y1   −12 6ℓ   1/ 3   F  = ==>  1  = − F  =  2 3  −6ℓ 2ℓ 2  θ   M  ℓ    2   1  −6ℓ 2ℓ  1/ 2ℓ   F ℓ   M1  On vérifie les équations d'équilibre de la structure Dans cet exemple le modèle élément fini donne la solution exacte car celle-ci est un polynôme d'ordre 3 comme l'approximation utilisée. Pour calculer l’état de contrainte sur les éléments, le diagramme du moment de flexion et celui de l'effort tranchant, nous utilisons la loi de comportement intégrée.  M f = EI v, xx = EI < N , x 2 > {U e } Pour chaque élément nous écrirons :  T = − EI v, xxx = − EI < N , x3 > {U e }

Rappel : Vous notez que le moment de flexion Mf est linéaire et que l’effort tranchant est constant par élément.

Exemple F

y x

< N,x2 > = < < N , x3

6 2

2 ( −2 + 3 s ) , ℓ 6 12 6 , − 3, 2 > 2 ℓ ℓ ℓ

( −1 + 2 s ) ,

ℓ 12 > = < 3, ℓ

6 ℓ

2

(1 − 2s ) ,

2 ( − 1 + 3s ) > ℓ

Tracer le diagramme des efforts intérieurs 6 2 F ℓ3  −1/ 3  M f1 = EI < 2 (1 − 2 s ) ( −1 + 3s ) >   = F ℓ( s − 1) ℓ EI −1/ 2ℓ  ℓ

ℓ

T1 = − EI


F ℓ3  −1/ 3    = −F EI −1/ 2ℓ 

On retrouve la solution analytique On vérifie bien que M f1 ( x = 0) = − M1 32

Éléments finis pour l'étude des portiques

33/87

Exercice 15 : Étude d’une poutre sous son poids propre Objectifs : mise en œuvre de la méthode des éléments finis, et illustrer la notion d'erreur liée à l'approximation. Nous cherchons la réponse statique sous son poids propre de la poutre sur appuis représentée par la figure ci contre.

Pb de flexion

A

B

g

Modèle à 1 élément. Déterminer

la matrice raideur, et le vecteur force généralisé associé au poids propre.

Écrivez le système réduit des équations, calculez les déplacements nodaux. Calculer la flèche au centre de la poutre, comparer à la solution analytique v ( ℓ / 2) = −

5 ρ gS ℓ 4 384 EI

Calculer les efforts aux appuis, et vérifier l'équilibre global de la structure. Calculer les efforts sur l'élément, tracer les diagrammes de l'effort tranchant et du moment de flexion, comparer à la solution analytique. Modèle à 2 éléments. Déterminer la matrice raideur assemblée complète. Déterminer le vecteur force généralisé associé au poids propre de la structure. Écrivez le système réduit des équations, calculez les déplacements nodaux et comparer à la solution analytique. Calculer les efforts aux nœuds, comparer à la solution analytique. Calculer les efforts sur l'élément et tracer les diagrammes de l'effort tranchant et du moment de flexion, et comparer à la solution analytique. Comparer à la solution analytique.. Répondez aux mêmes questions Prise en compte de la symétrie Utiliser la symétrie pour simplifier le modèle Calculer la matrice raideur et retrouver la solution du modèle à 2 éléments.

Exercice 16 : Études statique et dynamique d’une poutre Objectifs : Illustrer la notion d'erreur liée à l'approximation. Nous cherchons la réponse statique de la poutre sur appuis y o représentée par la figure ci contre. Modèle à 1 élément. Déterminer

ℓ A

C

B

xo

la matrice raideur. le vecteur force généralisé associé à la charge

Calculer la réponse statique, les efforts aux appuis et tracer le diagramme des efforts sur l’élément. Comparer à la solution analytique : 33

Éléments finis pour l'étude des portiques

34/87

7 Fℓ 3 1 Fℓ 2 1 Fℓ 2 , θ (C ) = − , θ (B) = 768 EI 128 EI 32 EI M f ( A) = 3F ℓ /16 , M f (C ) = 5 F ℓ / 32

v (C ) = −

Que pensez-vous de ce modèle, est-il satisfaisant ? Modèle à 2 éléments. Calculer la réponse statique, les efforts aux appuis et tracer le diagramme des efforts sur les éléments. Justifier les résultats de ce modèle. Réponse dynamique : Calcul des fréquences propres de la structure Modèle à 1 élément fini Modèle à deux éléments finis (vous pouvez utiliser Matlab ou Maple) Comparer à la solution analytique : EI EI EI , ω2 = 49,96 , ω3 = 104, 3 ω 1 = 15,42 4 4 ρSℓ ρ S ℓ4 ρSℓ Les solutions analytiques des poutres sont données sur le site Les exercices de cours sont corrigés sur le site, il faut chercher les réponses avant de consulter le corrigé. Pour assimiler le cours il faut aussi traiter des exercices non corrigés.

Application aux portiques Pour calculer les portiques nous devons utiliser un élément poutre tridimensionnel. Cet élément est obtenu par superposition des trois modèles suivants : • le modèle de traction, La flexion se décompose en deux problèmes de flexion plane • le modèle de torsion, dans les deux plans principaux de la section droite de la poutre. • le modèle de flexion.

Traction Torsion Flexion ( x , o, y ) Flexion ( x , o, z )

variables u

θx v, θ z w,θ y

Caractéristiques mécaniques ES , ρ S GJ , ρ I avec G = E 2(1 + ν ) EI z , ρ S EI y , ρ S

L'élément fini poutre tridimensionnel est un élément à deux noeuds et 6 degrés de liberté par nœud. Les 12 degré de liberté sont définis sur la base locale de l'élément.

{δ U e }T

= < (u , v, w, θ x , θ y , θ z )i

zo

(u , v, w, θ x , θ y , θ z ) j >

La matrice (12*12) du modèle tridimensionnel est obtenue par superposition des quatre matrices élémentaires elle est donnée à titre indicatif. :

34

θx

v

ye

θy

bo xo

j e

i

yo

u

ze

θz

w

xe

Éléments finis pour l'étude des portiques

35/87

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 élémentaires sur une base globale. Il faut passer aux calculs numériques MEFlab, Cast3M ou Abaqus

Statique des portiques plans simples Manuellement nous ne traiterons que le cas simple de portique plan ayant des éléments d’axe x ou y pour éviter le changement de base, et souvent pour simplifier le modèle nous négligeons les déformations dues à l'effort normal dans les éléments. Matrice raideur élémentaire d'un modèle traction - flexion

On pose : α =

α=

ES / ℓ EI / ℓ

3

=

S ℓ2 I

0 0 α  0 12 6ℓ  EI  0 6ℓ 4ℓ 2 [ Ke ] = 3 −α 0 0 ℓ  0 −12 −6ℓ   0 6ℓ 2ℓ 2

−α 0 0

α 0 0

   −6 ℓ 2 ℓ 2  0 0  12 −6ℓ   −6ℓ 4ℓ 2  sur 0

0

−12

6ℓ

S ℓ2 Élancement de la poutre I EI / ℓ3 Pour α → ∞ on tend vers la solution obtenue en négligeant les déformations de traction ES / ℓ

=

Pour un élément horizontal (orienté de i vers j suivant la direction des x)

y vi

La base locale et la base globale correspondent, la matrice raideur est celle donnée juste avant sur < ui , vi ,θi , u j , v j ,θ j > i

35

θi

vj ui ℓ

j

θj u

j

x

Éléments finis pour l'étude des portiques

36/87

Pour un élément vertical (orienté de j vers i suivant la direction des y) La base locale correspondra à la base globale, on raideur : 0 0 −α 0 α  0 12 6ℓ 0 −12  2 EI [ Ke ] = 3 −0α 60ℓ 4ℓ0 α0 −06ℓ ℓ  0 −12 −6ℓ 0 12  2 0 −6 ℓ  0 6ℓ 2ℓ

y

retrouvera la même matrice

vj

θj

j

 6ℓ   2ℓ 2  0  −6ℓ   4ℓ 2  sur

uj

0

ℓ vi

i

θi

ui

x

Mais attention à l’ordre des variables élémentaires

Exemple v2

θ2

F 2

Objectif : Déterminer la réponse statique de ce portique. Modèle à 2 éléments finis u2 ℓ

3

u3

ℓ 1

θ2 F 2

u2 ℓ

3

ℓ 1

u2

Ce modèle est suffisant pour obtenir la solution exacte du problème. C’est un modèle à 4 variables < u2 , v2 , θ 2 , u3 > Il conduit à résoudre un système de 4 équations, pour simplifier ce modèle nous allons négliger les déformations dues à l'effort normal dans les éléments.

v =0 Cette hypothèse permet d'écrire deux équations de liaison :  2 u3 = u2 Le modèle ne comporte plus que 2 variables Calculons directement les matrices élémentaires sur ces 2 variables. 12 6ℓ   6ℓ 4ℓ 2   

Pour l’élément 1 (2-1) : [ K1 ] =

EI ℓ3

Pour l’élément 2 (2-3) : [ K 2 ] =

EI 0 0  ℓ3 0 4ℓ 2 

 2 F ℓ3 u =  2 EI 12 6ℓ  u2   F  15 EI D'où le système réduit des équations : 3  = ==> U = { }     2 2 ℓ 6ℓ 8ℓ  θ 2   0  θ = − 1 F ℓ  2 10 EI C'est la solution exacte de la RDM u u

θ

Allure de la déformée

Calcul des réactions

36

Éléments finis pour l'étude des portiques Élément 1 : (2-1) M 21 2

R21 M11

1

Élément 2 : (2-3)

M 22

R22

R32 M 32

2

6ℓ − 12 6ℓ  u   R21   12   2 2   EI  6ℓ 4 ℓ − 6ℓ 2 ℓ  θ   M 21     =   ℓ 3 − 12 − 6ℓ 12 − 6ℓ 0   R11  2 2  6ℓ 2 ℓ − 6ℓ 4 ℓ  0   M     11 

< R21

R11

3

37/87

M 21

R11

M 11 >= F < 1 0,4 ℓ − 1 0,6ℓ >

6ℓ − 12 6ℓ   0   R22   12   2 2   EI  6ℓ 4 ℓ − 6ℓ 2 ℓ  θ   M 22    =     ℓ 3 − 12 − 6ℓ 12 − 6ℓ 0   R32  2 2  6ℓ 2 ℓ − 6ℓ 4 ℓ  0   M     32 

< R22

M 22

R32

M 32 >= F < − 0,6 − 0,4ℓ 0,6 − 0,2 ℓ >

Ce modèle ne nous donne pas toutes les composantes d’effort car nous avons négligé les allongements des éléments. Pour calculer la composante verticale de l’effort au noeud 1, nous pouvons écrire les équations d'équilibre de la structure.

0,6 F

F

- 0,2 Fℓ

Efforts aux appuis

-F

0,6 Fℓ

- 0,6 F

Exercice 17 : Étude d’un portique Objectifs : mise en œuvre de la méthode des éléments finis, changement de base, assemblage, résolution, calcul des efforts aux appuis, calcul des contraintes dans les éléments, et calcul des efforts aux nœuds internes. A

Intéressons-nous à la réponse statique du portique plan représenté par la figure ci-contre. On ne néglige pas l'effet de l'effort normal On posera α =

3 DDL par nœuds ( u i , v i , θ i ).

ES / ℓ EI / ℓ

ℓ f

ℓ

yo

3

xo

Modèle à 2 éléments. Définissez vos vecteurs globaux : {U } { FI } (bilan inconnues – équations) Déterminer Pour α = 2

la matrice raideur assemblée réduite. le vecteur force généralisé associé à la pression linéique.

Déterminer la déformée statique (déplacements nodaux). Calculer les efforts aux appuis, et vérifier les équations d’équilibre global de la structure. Pour chaque élément calculer les efforts (contraintes) au point A et analysez les discontinuités. 37

Éléments finis pour l'étude des portiques

38/87

Pensez-vous que votre modèle est satisfaisant ? (justifier votre réponse) Proposer un modèle plus satisfaisant, pensez-vous pouvoir résoudre ce modèle à la main ?

38

Formulation variationnelle & Écriture Matricielle

39/87

Formulation Variationnelle & Écriture Matricielle Ce paragraphe consacré à la formulation variationnelle d’un problème de mécanique, généralise les notions présentées avec les méthodes d'approximation des barres et des poutres. Problème aux limites formulation forte du problème

Système d'EDP

Méthode des résidus pondérés

Forme intégrale forte du problème Les conditions aux limites restent à satisfaire Nombreuses méthodes collocation, valeur moyenne Galerkin (==> éléments finis)

Transformation de la forme intégrale

Formulation faible du problème Les CL de Neuman ou mixte sont prises en compte

Approximation

Forme variationnelle du problème

Traitement numérique

Solution numérique

Forme matricielle

Passage d'un système EDP à une solution numérique

Toutes les méthodes d'approximation ont un même objectif, remplacer un problème mathématique défini sur un milieu continu (équations différentielles ou intégrales) par un problème mathématique discret (équation matricielle). Problème de dimension finie que l'on sait résoudre numériquement.

Formulation intégrale Partons de l’équation locale et des conditions aux limites définies dans le cadre de la MMC. ∀M ∈ D  ∀M ∈ Γu  ∀M ∈ Γσ

ρ uɺɺ − divσ = f

PFD appliqué à un élément de matière

u = ud

σ n = Td

Pour exprimer l’équation locale en fonction des déplacements ρuɺɺ + L (u ) = f Il faut associer au système précédent deux relations : - Les lois de comportement:

σ = fct (ε )

- Les Relations déplacements - déformations:

ε = fct ( u )

L’équation locale "EL" ⇔ ∀P

∫ P.( ρ uɺɺ − div( σ ) - f ) dV = 0 D

39

(Formulation forte)

Physique du matériau & expérimentale. Géométrique

Formulation variationnelle & Écriture Matricielle

40/87 Sachant que* σ : grad s ( P ) = div(σ P ) − P.div (σ )

"EL" ⇔ ∀P

∫ ( P.ρ uɺɺ

+ σ : grad s (P) − div(σ P) - P. f ) dV = 0

D

Appliquons le TH d'Ostrogradsky :

∫ div(σ P) dV = ∫ σ P. n dS = ∫ D

"EL" ⇔ ∀P

∫ ( P.ρ uɺɺ

∂D

P. σ n dS

∂D

+ σ : grad s (P ) - P. f ) dV −

D

∫ ( P. σ n ) dS = 0

∂D

En tenant compte des conditions aux limites sur la frontière : Γσ

∀M ∈ Γσ

σ n = Td

 "EL."  C.L sur Γσ ⇔

∫ ( P.ρ uɺɺ + σ : grad s (P)

∀P

- P. f ) dV −

∫ P. σ n dS − ∫ P.Td . dS = 0

Γu

D

Γσ

σ n est un champ inconnu sur cette frontière

Pour obtenir une équation ne dépendant que du champ inconnu u , nous utilisons des fonctions test à valeur nulle sur la frontière Γu Fonction Cinématiquement Admissible

∀M ∈ Γu

P(M ) = 0

∫ P. σ n dS = 0

Γu

D’où la formulation variationnelle équivalente au système d'équations aux dérivées partielles du problème, pour des fonctions de pondération cinématiquement admissibles. ∀PCA

∫ ( P.ρ uɺɺ

+ σ : grad s (P ) - P. f ) dV =

∫ P.Td . dS

Γσ

D

Les conditions aux limites en déplacement ∀M ∈ Γu

u = u d doivent être vérifiées

La méthode des éléments finis utilise cette formulation, avec deux idées fortes • la construction systématique des fonctions de forme par sous domaine « éléments finis » • utilisation des variables nodales comme paramètres d’approximation ce qui permet d’imposer les conditions aux limites en déplacement du problème.

Équivalence avec le PTV PTV : ∀δ u

δ W = δ A à tout instant Avec: δ A

travail virtuel des quantités d'accélération : ∫ γ ( P ).δ u dm ( P ) = D

δW

*

travail virtuel des efforts intérieurs et extérieurs

Cette formule utilise la symétrie du tenseur des contraintes et les relations suivantes: T

σ : grad (u ) = div(σ u ) − u .div (σ ) T

σ : grad (u ) = div(σ u ) − u .div (σ ) La démonstration de ces relations se fait simplement à partir d'une relation indicielle de la forme suivante:

σ ij ui , j = (ui σ ij ), j − ui σ ij , j

40

∫ δ u.ρ uɺɺ D

dV

Formulation variationnelle & Écriture Matricielle

41/87

= δ Wint = − ∫ σ : δ ε dV + ∫ f .δ u dV + D

D

∫ T .δ u dS

∂D

Compte tenu des relations déplacements - déformations : δε = grad s (δu ) *

Le PTV : ∀δ u

∫ δ u.ρ uɺɺ dV = − ∫ σ D

: grad s (δ u ) dV +

D

∫ δ u. f

dV +

D

∫ δ u.T dS

∂D

Notez que pour δ u = 0 sur Γu

Équivalence avec la FV

P ≡ δu

∫ δ u.T dS = ∫ δ u.Td . dS

∂D

Γσ

L'intérêt de ce principe est de fournir directement la forme intégrale sans avoir à passer par les équations aux dérivées partielles

Écriture matricielle du PTV Pour obtenir une équation matricielle, nous cherchons une solution approchée du problème sous la forme d'une combinaison linéaire de fonctions de forme. La méthode consiste alors à affaiblir la forme intégrale précédente en ne la satisfaisant que pour un nombre fini de fonctions de pondération. La méthode de Galerkin consiste à utiliser les mêmes fonctions pour l’approximation et la pondération. Pour l'approximation {u h ( M , t )} = [W ( M ) ] {q (t )} avec: Pondération

Pour

[W (M )] matrice des fonctions de forme {q(t )} vecteur des paramètres de l’approximation. {P( M )} = [W ( M )]{δq}

δ q =< 1, 0, 0,.... > δ q =< 0,1, 0,.... >

==>

P1 = W1

==>

P2 = W2

==>

Pn = Wn

Etc.…

δ q =< 0, 0, 0,.,1 >

Pour exprimer le produit σ : δ ε nous utilisons une représentation vectorielle des tenseurs

σ : δ ε = {ε } {σ } T

avec:

{ε }T =< ε xx , ε yy , ε zz ,2ε xy ,2ε xz ,2ε yz > T {σ } =< σ xx , σ yy , σ zz , σ xy , σ xz , σ yz >

La forme matricielle associée aux lois de comportement est alors : {σ

(M )

} = [D( M )]{ε ( M ) }

Les relations déformations - déplacements peuvent aussi s'écrire sous forme matricielle.

{ε

{

}

} = grad su ( M ) = [ L] {u ( M )} [L] Matrice d'opérateurs différentiels correspondant à l'expression du gradient

(M )

symétrique des déplacements.

{

}

Compte tenu de ces notations le produit σ : δ ε ==> σ h : grad s (P) = { P ( M )} [ L ] [ D ( M ) ] [ L ] u h ( M ) T

T

Soit compte tenu de l'approximation et de la pondération

σ h : grad s (P) = {P ( M )} [ B ( M ) ] [ D ( M ) ] [ B ( M ) ] {q ( t )} avec T

T

[B ( M )] = [L][W ( M )]

Reportons ces expressions dans la formulation variationnelle nous obtenons:

∀δ q

(

{δ q}T ∫ [W ]T ρ [W ]{qɺɺ} + [ B ]T [ D ] [ B ]{q} + [W ]T { f } D

) dV = {δ q}

T

T ∫ [W ] {T } dS

∂D

Cette équation pouvant être écrite quelque soit {δq} , nous obtenons l'équation matricielle : *

Cette relation sur le taux de déformation est démontrée en MMC, nous admettons, ici, ce résultat sans démonstration.

41

Formulation variationnelle & Écriture Matricielle

[ M ] {qɺɺ}

42/87

+ [ 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

Application aux modèles de l’ingénieur Pour le modèle « traction – compression » présenté pour l'étude des treillis : [ D ] = ES

[ L] =

∂ ∂x

∂2 Pour le modèle « flexion » présenté pour l'étude des portiques : [ D ] = EI [ L ] = 2 ∂x

Pour les modèles « de l’élasticité plane » Le modèle « contraintes planes » s'applique à des pièces peu épaisses chargées transversalement dont les faces ne sont pas chargées (plaques et coques minces, les membranes, capacité sans effet de fond, etc…) On considère que l'état de contrainte à σ11 σ12 0  l'intérieur du domaine est voisin de l'état de Hypothèse [σ ] = σ 21 σ 22 0 contrainte sur les surfaces, donc plan par rapport à la normale à la surface  0  0 0 soit : {σ } =< σ11 σ 22 σ12 > T

Écrivons l''inverse de la loi de HOOKE ε = −

ν

trace(σ ) 1 +

1 +ν

σ

E E Pour déterminer le tenseur des déformations à partir du tenseur des contraintes, en ne conservant que les termes à travail non nul.

 ε11  1   1  ε 22  =  −ν 2ε  E  0  12  

−ν

 σ11    1 0  σ 22  0 2(1 +ν )  σ12   1 ν σ11  E    Puis inversons cette relation : σ 22  = ν 1 2  σ  (1 −ν )   12  0 0 

La déformation

0

ε 33

qui n’est pas prise en compte

dans la loi de comportement, peut être calculée à

posteriori par : ε 33 = − ν (σ 11 + σ 22 ) E

   ε11    0   ε 22  ==> [ D ] (1 −ν )  2ε12   2  0

Le modèle « déformations planes » s'applique à des pièces massives dont les déformations longitudinales seront suffisamment faibles pour être négligées. Hypothèse :

 ε11 ε12 [ε ] = ε 21 ε 22  0 0

Écrivons la loi de HOOKE :

0 0 0

soit : {ε } =< ε 11 ε 22 T

σ = λ trace(ε ) 1 + 2 µ ε

42

2ε 12 >

λ = Eν  (1 + ν )(1 − 2ν ) avec  E µ = 2(1 + ν )

Formulation variationnelle & Écriture Matricielle

43/87

Ne conservons que les termes à travail non nul.  1 −ν ν σ11   E   1 −ν σ 22  =  ν (1 + )(1 − 2 ) ν ν σ    12  0  0 

   ε11    0   ε 22  ==> [ D ] (1 − 2ν )   2ε12   2  0

La contrainte

σ 33

qui n’est pas prise en compte

dans la loi de comportement, peut être calculée à posteriori par : σ 33 = λ (ε11 + ε 22 )

En résumé pour l’élasticité plane La matrice d'élasticité est de la forme:

  1  E (1 − aν )  ν [D] = (1 + ν )(1 − ν − aν )  (1 − aν )   0 

ν (1 − aν ) 1 0

∂   ε11   ∂x    L’opérateur différentiel [ L ] :  ε 22  =  0 2ε    12   ∂   ∂y

Exemple

a = 0 en contraintes planes avec  a = 1 en déformations planes

 0  ∂  u     tel que {ε } = [ L ]{u} ∂y   v  ∂  ∂x 

Objectif : construction des matrices élémentaires pour un élément triangulaire quelconque

yo y3

    0  (1 − ν − aν )   2(1 − aν )  0

Soit un élément triangulaire à trois nœuds, nous avons trois variables nodales à identifier. Nous cherchons donc une approximation polynomiale linéaire de la forme :

3 élément réel

y2 y1

2 1

u x3

x1

x2

h

( x, y)

= [1 x

xo

 a1    y ]  a2  a   3

(1)

Pour exprimer l'approximation en fonction des déplacements nodaux nous allons identifier les valeurs nodales des déplacements  u1  1 x1 y1   a1       soit pour le déplacement suivant x ui = u h ( xi , y i ) , aux 3 nœuds : u2  = 1 x 2 y 2  a 2  u  1 x y 3  a 3  3  3  Il est simple de vérifier que la résolution de ce système conduit à :  a1  ∆ 23 1    a 2  =  y 23 2 A a   x 32  3

∆ 31 y 31 x13

∆ 12   y12  x 21 

 A = aire du triangle  u1     u2  avec  xij = xi − x j et yij = yi − y j u  ∆ = x y − x y i j j i  3  ij

43

Formulation variationnelle & Écriture Matricielle

44/87

Reportons ce résultat dans l’approximation (1) , nous obtenons :

u

h

( x, y)

= [ N1

N2

 u1    N 3 ] u2  u   3

avec par permutation circulaire de ijk

Ni =

De la même façon nous aurons : v

u   N 1 Soit   =  v   0

Nous avons vu que : {ε

0 N1

(M )

N2 0

0 N2

N3 0

( x, y )

= [ N1

N2

 v1    N 3 ] v2  v   3

 u1  v   1 0  u 2    = [N ( x , y ) ]{U e } C'est l'approximation nodale du T3 N 3  v 2  u 3     v 3 

} = [ L ] {u ( M )} = [ L ][ N ( x, y ) ]{U e } = [ B( M ) ] {U e }

 N1, x  D'où l'expression : [ B ( x , y ) ] =  0  N1, y 

 y23 1  [ B( x, y )] =  0 2A  x32

h

1 ( ∆ + x y jk − y x jk ) 2 A jk

0

N 2, x

0

N 3, x

N1, y N1, x

0 N 2, y

N 2, y N 2, x

0 N 3, y

0

y31

0

y13

x32 y23

0 x13

x13 y31

0 x31

0   N 3, y  N 3, x 

0 x31  cette matrice est constante. y13 

L'approximation de l'état de contrainte est donc constante par élément. Pour calculer les matrices élémentaires [ M ]e , [ K ]e et { F }e le vecteur associé à une densité de charge sur l'élément, il reste à calculer les produits des matrices et à intégrer sur l'élément.

Exercices : Formulations Variationnelles en physique Exercice 1: Formulation variationnelle d'un problème de thermique Objectifs : Établir la formulation variationnelle tridimensionnelle d’un problème de thermique. En donner la forme matricielle pour la méthode de Galerkin. Soit un corps de masse volumique ρ , de chaleur massique à volume constant cv , et de conductivité λ . Équations du problème Le flux thermique est donné par la loi de Fourier : q = −λ grad T ∂T ∂t ou r est un flux de chaleur créé par unité de volume

L'équation générale du bilan énergétique est :

divq − r = − ρ cv

A ces deux équations il faudra associer une condition initiale (répartition initiale des températures) Les conditions aux limites du problème peuvent être : De type champ (problème de Dirichlet) T = Td 44

Formulation variationnelle & Écriture Matricielle

45/87 q.next = ϕ d

De type flux (problème de Neumann)

ϕ c = hc (Ts − T∞ ) Condition de convection

De type mixte (problème de Fourier)

1- Écrire le système d'EDP général d’un problème de thermique stationnaire. 2- Transformer cette formulation forte pour obtenir la formulation variationnelle faible de ce problème. 3- Donner l’expression matricielle de cette formulation faible pour la méthode de Galerkin, avec une approximation de la forme T =< ψ > {q} Les exercices de cours sont corrigés sur le site, il faut chercher les réponses avant de consulter le corrigé. Pour assimiler le cours il faut aussi traiter des exercices non corrigés.

Exercice 2: Écoulement d'un fluide visqueux dans une conduite. Objectifs : Obtenir différentes formulations variationnelles de ce problème Appliquer les méthodes d'approximation générale sur le domaine. Considérons le problème de l'écoulement stationnaire d'un fluide visqueux incompressible dans une conduite de section droite carrée Ω . y p − p1 La chute linéique de pression est donnée π = 2 Γ ℓ p1 u p2 2a On note µ la viscosité cinématique du fluide. x Ω L'inconnue est le champ des vitesses u = u ( x, y ) z .

Section 2a

ℓ

Pour cet écoulement les équations de Navier-Stokes se réduisent à : L'équation locale

∀M ∈ Ω ∆(u ) =

Les conditions aux limites

∀M ∈ Γ u = 0

π µ

W ( x , y ) = ( x 2 − a 2 )( y 2 − a 2 ) 1- Soit :  1 Deux fonctions de forme définies sur le domaine 2 2 W2 ( x , y ) = W1 ( x , y ) ( x + y ) Ces deux fonctions de forme W1 ( x , y ) et W2 ( x , y ) sont-elles cinématiquement admissible ? 2- Formulation forte du problème Établir la formulation forte du problème. Peut-on utiliser cette formulation avec les fonctions de forme précédente ? Donner l'écriture matricielle de l'équation pour des pondérations P(x, y) quelconques. 3- Pour une approximation à 1 paramètre sur W1 ( x , y ) Comparer les résultats obtenus pour : La méthode de collocation au point (0,0) puis au point (a/2,a/2) La méthode de la valeur moyenne de l'erreur P(x, y) = 1 La méthode de Galerkin : P(x, y) = W1 ( x, y ) 4- Formulation faible du problème (formulation variationnelle) Établir la formulation faible de ce problème Donner la forme matricielle correspondant à la méthode de Galerkin Retrouver l'approximation obtenue pour une approximation à 1 paramètre. Calculer la nouvelle approximation pour une approximation à 2 paramètres u h = < W ( x , y ) > {q} = W1 ( x , y ) q1 + W2 ( x , y ) q2 Le fichier Maple vous simplifiera les calculs numériques

45

Formulation variationnelle & Écriture Matricielle

46/87

Exercice 3: Discrétisation "EF" de la conduite Objectifs : Appliquer les méthodes d'approximation par sous domaines. Effectuer les calculs au niveau élémentaire sans utiliser la notion d'élément de référence. yo Le domaine étudié est un carré le coté 2a , pour simplifier les calculs 3 nous avons orienté différemment le repère d'observation. 2a La discrétisation que nous allons utiliser est représentée sur la figure ci (2) (1) contre. Les mailles sont des éléments finis de type T3. 4 2 xo (3) 1 (4) sur D ∆u = f Le champ des vitesses u ( x, y ) vérifie  u = 0 sur le contour de D 5

Formulation Rappeler la forme variationnelle du problème. Donner la forme matricielle correspondant à une approximation de cette équation intégrale. Calcul de la matrice raideur élémentaire Les quatre éléments étant identiques, on limitera les calculs au premier élément (l'élément 123). Définir les fonctions d'interpolation de cet élément. Déterminer la matrice raideur et le vecteur force généralisé de cet élément. Montrer que pour les autres éléments les résultats sont identiques. Résolution Calculer l'approximation du champ de vitesse au centre de la conduite

Exercice 4: Formulation matricielle d'un problème axisymétrique Objectifs : Appliquer les méthodes d'approximation par sous domaines. Pour les modèles « axisymétriques » Nous utilisons le système de coordonnées cylindriques. Compte tenu des hypothèses de symétrie, le champ des déplacements est de la forme b

 ur ( r , z , t ) = u  u   {u ( M , t )} =  uθ = 0  =   u ( r , z , t ) = w   w  z 

zo symétrie cylindrique

1- Exprimer l'opérateur gradient du = grad (u ) dX sur la base b

xo

θ

ez

yo

eθ b

2- En déduire en petites déformations l'opérateur [ L ] ; {ε } = [ L ]{U }

er

3- Exprimer la matrice d’élasticité [D] correspondant à la loi de HOOKE : σ = λ trace(ε ) 1 + 2 µ ε 4- pour un élément triangulaire T3

zo

u  Rappeler l'expression de la matrice [ N ] ;   = [ N ( x , y ) ]{U e }  w En déduire l'expression de la matrice [ B ] ; {ε } = [ B ]{U e }

46

k élément réel

T3 j i symétrie cylindrique

Méthodes numériques

47/87

Méthodes numériques dans le cadre de la MEF Les principales étapes de construction d'un modèle éléments finis sont les suivantes : • Discrétisation du milieu continu en sous domaines. • Construction de l'approximation nodale par sous domaine. • Calcul des matrices élémentaires correspondant à la forme intégrale du problème. • Assemblage des matrices élémentaires - Prise en compte des conditions aux limites. • Résolution du système d'équations.

Discrétisation du milieu Cette opération consiste à décomposer le domaine continu en un nombre fini de sous domaines « éléments 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 éléments ayant une frontière commune. De plus lorsque la frontière du domaine est complexe, une erreur de discrétisation géométrique est inévitable. Cette erreur doit être estimée, et éventuellement réduite en modifiant la forme ou en diminuant la taille des éléments concernés. Modifier la taille des éléments

Pièce présentant des congés de raccordement

Changer la géométrie éléments à frontière courbe Erreur de discrétisation géométrique

Erreur de discrétisation géométrique.

Approximation nodale La méthode des éléments finis est basée sur la construction systématique d'une approximation u h du champ des variables par sous domaine. Cette approximation est construite sur les valeurs approchées du champ aux noeuds de l’élément. On parle de représentation nodale de l’approximation ou plus simplement d’approximation nodale. L'approximation par éléments finis est une approximation nodale par sous domaines ne faisant intervenir que les variables nodales du domaine élémentaire De .

∀M ∈ De

{u } = [ N h

(M )

(M )

] {U n }

{u }

Valeur de la fonction approchée en tout point M de l'élément

{U n }

Matrice des fonctions d'interpolation de l'élément Variables nodales relatives aux nœuds d'interpolation de l'élément

h

[N]

47

Méthodes numériques

48/87 Valeurs approchées aux noeuds xi

élément 1

+ élément 2 x

+

Approximation linéaire utilisant 3 éléments

élément 3

Approximation nodale linéaire à une dimension

Éléments de référence Un élément de référence est un élément de forme géométrique simple (frontières rectilignes, bords de longueur unité). Pour la majorité des éléments de référence, l'approximation nodale est construite sur une base polynomiale de degré 1 ou 2. Le nombre de variables nodales à identifier étant égal à la dimension de la base. Bases polynomiales complètes: 1D:

linéaire [1, x] quadratique [1, x, x2]

2 variables 3 variables

2D:

linéaire [1, x, y] quadratique [1, x, y, x2, xy, y2]

3 variables 6 variables

T3 T6

3D:

linéaire [1, x, y, z] quadratique [1, x, y, z, x2, xy, y2, xz, z2, yz]

4 variables 10 variables

Tétraèdre d°1 Tétraèdre d°2

Bases polynomiales incomplètes: 2D:

3D:

"quasi-linéaire"

[1, x, y, xy]

4 variables

Q4

"quasi-quadratique" [1, x, y, x2, xy, y2, x2y, y2x] 8 variables

Q8

"quasi-linéaire"

cube

[1, x, y, z, xy, xz, yz, xyz]

8 variables

Ce choix ne comporte que des monômes de degré 1, et respecte la symétrie de la base.

Construction de l'approximation nodale L’approximation nodale est construite sur l'approximation polynomiale : ∀M u h ( M ) =< Φ ( M ) > {a} L'approximation est identifiée à la valeur du champ de variables aux n nœuds M i de l'élément

{U n } = {u h ( M n )} =

< Φ( M n ) >

{a}

En résolvant ce système d'équation, nous exprimons {a} en fonction des variables nodales {U n } −1

 < ... >  {a} = [T ]{U n } avec [T ] = < Φ( M i ) >   < ... >  Reportons ce résultat dans l'approximation

La matrice à inverser doit être bien conditionnée. Conditionnement lié au choix de la base polynomiale et à la géométrie des éléments de référence.

Nous obtenons la matrice des fonctions d'interpolation : < N ( M ) >=< Φ ( M ) > [T ]

48

Méthodes numériques

49/87

Les fonctions d'interpolation satisfont la propriété suivante ∀M i

Éléments à une dimension s 0

1

C'est un segment de droite de longueur unité : s ∈ [ 0,1]

Approximation linéaire base polynomiale utilisée est (1, s )  N 1 ( s) = L1 = 1 − s   N 2 ( s) = L2 = s Approximation quadratique Base polynomiale (1, s, s 2 )

0 si i ≠ j N j ( Mi ) =  1 si i = j

2 nœuds

1

N1

0

1

1

2

N2

1 N1

3 nœuds

 N 1 ( s) = L1 (2 L1 − 1)   N 2 ( s) = 4 L1 L2  N ( s) = L (2 L − 1) 2 2  3

N2

N3

0

 N 1 ( s) = 1 − 3s 2 + 2 s 3  2 3  N 2 ( s) = s − 2 s + s  2 3  N 3 ( s) = 3s − 2 s  N 4 ( s) = − s 2 + s 3

s

1

2

1

Approximation cubique Base polynomiale (1, s, s 2 , s 3 ) 4 nœuds L1   N 1 ( s) = 2 (3 L1 − 1)(3 L1 − 2)  9  N 2 ( s) = L1 L2 (3 L1 − 1) 2  9  N 3 ( s) = L1 L2 (3 L2 − 1) 2  L2   N 4 ( s) = 2 (3 L2 − 1)(3 L2 − 2) Identification du champ et sa dérivée (pente)

s

3

Tous ces éléments sont de type Lagrange

1

N1

N2

N3

N4

0

s

1

1

1

2

N1

3

u et u , x

4

N3

N2 1

1 1

N4

Cet élément est de type Hermite cf. élément poutre.

s s

2

1

Éléments triangulaires

t

C'est un triangle rectangle de coté unité. s ∈ [ 0,1] et t ∈ [ 0,1 − s ]

Approximation linéaire : Base polynomiale (1, s, t ) élément à 3 nœuds, triangle de type "T3"

(0,1) 3

1 (0,0)

Les fonctions d'interpolation sont :

49

2 (1,0)

s

Méthodes numériques

50/87

N1 = 1 − s − t

N2 = s

N1

t

N3 = t t

N2

3

N3

3 1

1 2

t

3

1

2

s

2

s

s

Approximation quadratique Base polynomiale (1, s, t , st , s 2 , t 2 ) élément à 6 nœuds, triangle de type "T6".

t (0,1) 3

4

5

On pose : L1 = 1 − s − t , L2 = s , L3 = t Les fonctions d'interpolation sont :

s

2

1 6

(0,0)

(1,0)

Pour les 3 noeuds d'interface (4,5,6):

Pour les 3 noeuds sommet: i = 1,2,3 N i = Li (2 Li − 1) N i+3

3

i = 1, 2,3 avec j ≠ i, k et k ≠ i = 4 L j Lk 1

5 1 6

4

5

1

1 6

N1 = L1 (2 L1 − 1)

N 4 = 4 L2 L3

2

2

4

3

Éléments rectangulaires plans

L'élément de référence est un carré de coté 2 : ( s, t ) ∈ [ −1,1]

Approximation linéaire "Q4" (1, s, t , st )

 N1  N2  N3 N  4

t 3

4

(-1,1)

1

N1

= 4 (1 − s)(1 − t ) 1 = 4 (1 + s)(1 − t ) 1 = 4 (1 + s)(1 + t ) 1 = 4 (1 − s)(1 + t )

(1,1)

s

3

1

2

1

(-1,-1)

1

(1,-1)

2

Approximation quadratique "Q8"

t

Pour éviter d’avoir des nœuds internes, on utilise des bases polynomiales incomplètes mais symétriques contenant tous les monômes d’un même degré. 2

2

2

7

6

8

5

4

2

s

Base polynomiale (1, s, t , st , s , t , s t , t s ) 1

2

3

Nous avons donné les fonctions d’interpolation des éléments 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 éléments de référence sont données dans le livre de Dhatt et Touzot "Une présentation de la méthode des éléments finis". Vous y trouverez aussi des exemples de programmes pour la construction systématique de ces fonctions d’interpolation. 50

Méthodes numériques

51/87

Calcul des matrices élémentaires Présentons maintenant les techniques numériques élémentaires (utilisées sur chaque élément) permettant de calculer les formes matricielles déduites de la formulation variationnelle (forme intégrale) d’un problème de physique. Le calcul des matrices élémentaires nécessite une dérivation puis une intégration sur le domaine élémentaire. Le calcul analytique n'est possible que pour des éléments très simples, les éléments monodimensionnels ou en 2D le T3 et le Q4 rectangulaire. Un code éléments finis a donc recours aux calculs numériques, basés sur l’utilisation d'éléments de référence, de la matrice Jacobienne de la transformation géométrique, et de l’intégration numérique pour calculer les coefficients des matrices. Ce paragraphe présente quelques aspects du calcul numérique, indispensables pour comprendre les erreurs numériques liées à la forme du maillage, lors de l’analyse de résultats.

Transformation géométrique Tout élément réel peut être défini comme l'image par une transformation géométrique d'un élément parent dit de référence pour lequel les fonctions d'interpolation sont connues. Nous entendons par élément réel un élément quelconque du domaine discrétisé. La transformation géométrique définit les coordonnées ( x, y, z ) de tout point de l'élément réel à partir des coordonnées ( s, t , u ) du point correspondant de l'élément de référence. Un même élément de référence permet donc de générer toute une classe d'éléments réels. A chaque élément du domaine réel correspond une transformation bijective unique. En 3D la transformation géométrique est définie par :  x =< N g ( s , t , u ) > { x n }  { xn } , { yn } , { zn } Coordonnées des noeuds de l'élément réel  y =< N g ( s , t , u ) > { y n } avec :  Fonctions de la transformation géométrique < N g ( s , t , u ) >  z =< N ( s , t , u ) > z { n} g  Dans ce cours nous ne présentons que des éléments iso paramétriques pour lesquels les nœuds d'interpolation et les nœuds géométriques sont confondus. Les fonctions de la transformation géométrique N g seront identiques aux fonctions d'interpolation N .

Exemples d'éléments de référence classiques Éléments à une dimension. Référence linéaire

Réel

τe

quadratique cubique Transformations géométriques d'éléments à une dimension Ces transformations géométriques utilisent les fonctions d'interpolation linéaire, quadratique et cubique définies plus haut.

Éléments à deux dimensions. Pour ces éléments les transformations géométriques conduisent respectivement à des frontières linéaires, quadratiques.

51

Méthodes numériques

t

52/87 Éléments carrés:

Éléments triangulaires: t

t (-1,1)

t

(1,1)

(0,1)

s

s

1/2

s

s

(0,0) (1,0) Linéaire (3)

(-1,-1) (1,-1) Linéaire (4)

Quadratique (6)

Éléments à trois dimensions. u

u

Quadratique (8)

u

(0,0,1) (0,0,1)

t (0,1,0) (1,0,0)

s

tétraédrique (4)

(0,0,-1)

(0,1,1) (-1,-1,1) (1,0,1) (1,-1,1) t (-1,-1,-1) (0,1,-1) (1,-1,-1) s (1,0,-1)

prismatique (6)

(-1,1,1) (1,1,1)

t (-1,1,-1) (1,1,-1)

s

cubique (8)

Éléments volumiques à transformation linéaire

Matrice Jacobienne - transformation des opérateurs de dérivation Nous connaissons (savons calculer) les dérivées des fonctions d'interpolation par rapport aux coordonnées de l'élément de référence. ( s, t , u ) . Or il faut calculer les dérivées des fonctions d'interpolation par rapport aux coordonnées réelles ( x, y, z ) . La matrice Jacobienne de la transformation [ J ] est définie par: ∂   ∂ s  ∂ ∂ t  = ∂   ∂ u 

∂y ∂ x ∂ z  ∂  ∂  ∂s ∂ s  ∂ x   ∂s  ∂ x ∂y ∂ x ∂ z   ∂  = [ J ]  ∂  ∂y ∂t ∂ t ∂ y  ∂t    ∂y ∂ x  ∂  ∂  ∂z   ∂u ∂ u  ∂ z   ∂ z  ∂ u

 x =< N g ( s , t , u ) > { x n }  Compte tenu de l’expression de la transformation géométrique  y =< N g ( s , t , u ) > { y n }  z =< N ( s , t , u ) > z { n} g  ∂ < N g >  ∂   ∂ s  ∂ s    [ J ] =  ∂ ∂ t  < x y z >=  ∂ < N g >  {x n } { y n } {zn } ∂t ∂   < N > ∂  ∂ u  g   ∂ u  [ J ] est le produit d'une matrice (3, n) par une matrice (n, 3) connues.

[

]

Pour chaque élément, la matrice Jacobienne s'exprime en fonction des dérivées des fonctions de la transformation géométrique et des coordonnées des nœuds géométriques de l'élément réel.

52

Méthodes numériques

53/87

La relation inverse permet alors de calculer les dérivées premières par rapport aux coordonnées réelles des fonctions d'interpolation. ∂  ∂   ∂ x  ∂ s  ∂  −1 −1  ∂ La transformation devant être une bijection [ J ] existe  ∂ y = [ J ]  ∂ t  ∂  ∂   ∂ z   ∂ u  Une singularité de J peut apparaître lorsque l'élément réel est trop "distordu" par rapport à l'élément de référence « élément dit dégénéré ». De façon générale on évite lors du maillage d'utiliser des éléments trop disproportionnés, car ils nuisent à la précision numérique du modèle De plus en plus de logiciels de pré-traitement proposent des outils de contrôle de la qualité du maillage basé sur la taille, les proportions et le calcul du Jacobien des éléments utilisés. Pour le calcul des dérivées secondes

∂2

∂x

2

et

∂2

∂x∂y

des fonctions d'interpolation (problèmes de flexion), la

démarche est identique mais les calculs sont plus complexes, reportez vous au livre de D&T (pages 55-57).

Calcul numérique d'une intégrale Le jacobien de la transformation géométrique permet de passer de l'intégration d'une fonction f définie sur l'élément réel à l'intégration sur l'élément de référence :

∫f

(x, y, z)

dxdydz =

De

∫f

(s, t, u)

det[ J ] dsdtdu

Dref

Cette dernière intégrale ne peut être évaluée analytiquement que dans des cas extrêmement simples. En général, la fonction à intégrer est une fraction rationnelle polynomiale compliquée. Le calcul de l'intégrale sur l'élément de référence est donc effectué numériquement. Les formules d'intégration numérique permettent d'évaluer l'intégrale sous la forme générale suivante :

∫ Dref

où

Npi

f dv ≅ ∑ f (ξi )ωi i =1

N pi

représente le nombre de points d'intégration sur l'élément de référence.

ξi ωi

les coordonnées paramétriques des points d'intégration. les poids d'intégration.

L'intégration numérique exacte n'est possible que si la fonction à intégrer est polynomiale. La matrice Jacobienne doit être constante (l'élément réel garde la même forme que l'élément de référence). Nous connaissons alors l'ordre de la fonction polynomiale à intégrer, et nous pouvons choisir en conséquence le nombre de point d’intégration. Dans le cas d'éléments réels de forme quelconque, la matrice Jacobienne est une fonction polynomiale. Les termes à intégrer sont donc des fractions rationnelles, et la précision de l'intégration numérique diminue lorsque l'élément réel est mal conditionné (disproportionné) La précision de l'intégration numérique dépend aussi du choix du nombre de points d'intégration, ce nombre est proposé par défaut dans les codes éléments finis. Pour les problèmes non linéaires en dynamique on utilise souvent un nombre de points d'intégration plus faible pour diminuer les temps de calcul (Abaqus par exemple fait de l'intégration réduite par défaut), attention un nombre de points d'intégration insuffisant peut conduire à des résultats faux*. *

Pour un problème de dynamique l’intégration réduite peut introduire des modes parasites (phénomènes d’hourglass)

53

Méthodes numériques

54/87

Les schémas d'intégration les plus utilisés pour les éléments 2D sont : Coordonnées ξi

points

3

1

1 3

4

s=±

1

s=

3

1/ 6 s = 2/3 1/ 6

2

1 3

t=±

Poids ωi

1 3

ω =1

1 3

ω = 1/ 6

1/ 6 t = 1/ 6 2/3

ω = 1/ 6

t=

Vous trouverez dans le D&T (pages 280 à 300) les tableaux et les figures donnant la position et les poids d'intégration pour différents schémas d'intégration. Organisation des calculs numériques pour chaque élément

∀e

pour chaque point d'intégration

∀i = 1, Np (ξi , ωi )

Calcul du jacobien, des fonctions d'interpolation, et de leurs dérivées.

N g (ξ i ), { xn } , { yn } , { zn }

[ J (ξi )] , [ J (ξi )]−1 [ N (ξi )] ; [ B(ξi )]

[ Ke ] = [ Ke ] + [ Bξi ]T [ D] [ Bξi ]det( Jξi )ωi [ M e ] = [ M e ] + [ Nξi ]T ρ [ Nξi ]det( Jξi )ωi {Fe } = {Fe } + [ Nξi ]T f det( Jξi )ωi

Somme sur les points d'intégration ==> Matrices élémentaires

Fin de l'intégration numérique

[ K ] = [ K ] + [ Ke ] [M ] = [M ] + [M e ] {F } = {F } + {Fe }

Assemblage dans les matrices globales

Fin du calcul

Exemple : L'élément triangulaire "T3" yo

τe

t 3

1

3 e 1

s 2 référence

x3 x1

C'est un élément iso-paramétrique 2

x2

xo

Les fonctions d’interpolation sont : N1 = 1 − s − t

élément réel

54

N2 = s

N3 = t

Méthodes numériques

55/87

La matrice Jacobienne de cette transformation géométrique. ∂ < N g >    x1    −1 1 0     ∂ s  [ J ] = ∂ < N >  { xn } { yn } =  −1 0 1    x2     x  g    3  ∂ t 

x − x

[ J ] =  x2 − x1 

3

1

y2 − y1   x21 = y3 − y1   x31 −1

Son inverse : [ J ]

y21  y31 

 y1     y2    y   3 

Cette matrice est une constante

A est l’aire de l'élément réel. J =2A est le jacobien de la transformation.

1  y31 − y21  =   2 A  − x31 x21 

∂  ∂  ∂ x  1  y31 − y21   ∂ Sachant que :  =   ∂ ∂ y  2 A  − x31 x21  ∂ ∂  

 s  t 

Nous en déduisons : 1 − y31 + y21  1  y23   N1, x  −1  −1  = J  =  =    N1, y  −1 2 A  x31 − x21  2 A  x32   N 2, x  1  y31  −1 1   = J  =   0  2 A − x31   N 2, y  ε   u  xx ,x  N 3, x  1 − y21  −1 0          = J  =   = = v ε ε { }     = [ B ] {U e } yy , y N  3, y  1  2 A  x21      2ε xy  u, y + v, x  D'où l’expression de la matrice [ B ( x, y ) ]

N  1, x [ B ( x, y ) ] =  0   N1, y

 y23 1  [ B ( s, t ) ] =  0 2A  x32

0

N 2, x

0

N3, x

N1, y

0

N 2, y

0

N1, x

N 2, y

N 2, x

N3, y

0

y31

0

y12

x32 y23

0 x13

x13 y31

0 x21

0   N 3, y   N3, x 

0  x21  y12 

Ces résultats sont ensuite utilisés pour le calcul des matrices élémentaires : Matrice raideur:

[K e ] = ∫ [ B( M )]T [ D( M )] [ B( M )] dv

« e » Épaisseur supposée uniforme de l’élément

1 1− s

= e ∫ ∫ [ B ( s, t )]T [ D] [ B( s , t )] 2 A dsdt

De

0 0

Les termes de cette matrice sont des constantes.

Matrice masse:

[M e ] = ∫ [N ( x, y )]T ρ [N ( x, y )] dv

1 1− s

= e∫

De

T ∫ [N (s, t )] ρ [N ( s, t )] 2 A dsdt

0 0 1 1− s

Dans le cas d’une charge de volume f d

{Fde } = e ∫ ∫ [ N ( s, t )]T { f d }

2 A dsdt

0 0

N’ayant que des polynômes à intégrer le calcul analytique ou numérique de ces matrices ne pose pas de problème, et il conduira à des résultats exacts qui pourront être réutilisés.

55

Méthodes numériques

56/87

Assemblage et conditions aux limites Ne

Les règles d'assemblage sont définies par la relation : D ≅ ∑ De e =1

Ne

T T ∑ {δ U e } [ M e ]{Uɺɺe } = {δ U } [ M ]{Uɺɺ} e =1 Ne

ne

et

T T ∑ {δ U e } [ Ke ]{U e } = {δ U } [ K ]{U }

Il faut penser énergie ou travail virtuel pour effectuer la sommation des termes

T T ∑ {δ U e } {Fde } = {δ U } {Fd } e =1

e =1

Cette opération traduit simplement que l’énergie associée au domaine étudié est la somme des énergies élémentaires des sous domaines. Cela consiste à ranger dans une matrice globale les termes des matrices élémentaires. Le programme définira l’ordre des variables globales pour optimiser la place mémoire (disque) de la matrice globale, mais aussi le temps de calcul en fonction des algorithmes de résolution utilisés*. Deux méthodes classiques « matrices bandes » et « matrices ligne de ciel » sont présentées dans le livre de D&T. Après 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 résoudre, il faut tenir compte des P conditions aux limites cinématiques associées aux P composantes inconnues du vecteur { Fi } . Dans le cas d’un calcul statique, La méthode directe de résolution par blocs, peut se présenter sous la forme suivante en regroupant les termes des matrices. [ K11 ]  [ K 21 ]

{U i } = [ K11 ]−1 {{Fd 1} − [ K12 ]{U d }}  { Fi } = [ K 21 ]{U i } + [ K 22 ]{U d } − { Fd 2 }

[ K12 ] U i  =  Fd 1  +  0  [ K 22 ] U d   Fd 2   Fi 

Dans le cas particulier ou {U d } = {0} seul les termes de [ K11 ] et [ K 21 ] sont utiles −1 {U i } = [ K11 ] { Fd 1} En effet  { Fi } = [ K 21 ]{U i } − {Fd 2 }

Nous n’abordons pas ici les méthodes numériques de résolution de ces équations matricielles. Ces méthodes sont vues en analyse Numérique. Pour les problèmes de statique vous trouverez dans les codes EF deux types de méthodes • Les méthodes directes : Élimination de Gauss, décomposition LDU, ou Cholesky … • Les méthodes itératives de type Gauss-Seidel.

Exercice Les exercices de cours sont corrigés sur le site, il faut chercher les réponses avant de consulter le corrigé.

Exercice 1 : EF-T3 pour l’élasticité plane Objectifs : Assimiler les techniques de calcul au niveau élémentaire mise en œuvre dans un modèle élément finis pour un problème de mécanique.

*

Le choix de l’algorithme d’assemblage est un problème numérique & informatique. Les algorithmes dépendent de la taille du système matriciel, de la nature du problème (dynamique, statique, linéaire, non linéaire, etc… ), de la machine utilisée (place mémoire disponible, espace de stockage, parallélisme, …. ).

56

Méthodes numériques

57/87

La structure à étudier est une plaque mince homogène isotrope chargée dans son plan. Nous utiliserons donc un modèle contraintes planes. Les conditions aux limites et le modèle éléments finis à utiliser sont représentés sur la figure ci-contre. Les éléments sont les éléments de type T3 présenté dans ce chapitre de cours.

3

4

(2) a

1

F

(1) 2a

yo 1

2

xo 1

Formulation Rappeler : la forme variationnelle du problème. la forme vectorielle de la loi de comportement. la matrice reliant le vecteur des déformations et les déplacements nodaux. l'expression matricielle de la matrice raideur élémentaire. Calcul de la matrice raideur élémentaire Nous voulons utiliser l'élément de référence T3 pour effectuer les calculs en suivant les étapes d'un code éléments finis. On se limitera au premier élément pour éviter les calculs répétitifs. Exprimez la matrice Jacobienne de la transformation géométrique du premier élément. Calculer les dérivées des fonctions d'interpolation sur l'élément réel Ni,x et Ni,y Donnez l'expression matricielle permettant de calculer la matrice raideur du premier élément. Calcul de l'état de contrainte sur les éléments La solution en déplacement est supposée connue. Donner l'expression du vecteur des contraintes pour les deux éléments du modèle. Vous exprimerez ces vecteurs en fonction des déplacements nodaux. Comparer les expressions.

Exercice 2 : Élément Q4 en contraintes planes Objectifs : Transformation géométrique et intégration numérique, analyse du script Q4_ke de MEFlab Soit l'élément de référence quadrilatère à quatre nœuds de type Q4.

τe

t 4

3

yo

4 3 e 1

s 1

2 référence

2

x1

x4

x2

x3

xo

élément réel

Rappeler : la base polynomiale de l’approximation. le principe de construction de l'approximation nodale. u  l’expression de [ N ( s , t ) ] telle que :   = [ N ( s , t ) ]{U e } v  Transformation géométrique du Q4 Appliquez la transformation au centre du carré puis au point de coordonnées s = t = 0,5 Donner l’expression de la matrice Jacobienne de cette transformation géométrique en fonction de s, t et xi , yi Que pensez-vous du calcul de l'inverse de la matrice Jacobienne ? 57

Méthodes numériques

58/87

Dans le cas particulier ou l'élément réel est un rectangle 1  2a 0  Montrer que la matrice Jacobienne est : [ J ] =  4  0 2b 

yo

a 4

3 b

1

e

En déduire l’expression de [ J ]

2

−1

xo

Calculer la dérivée première par rapport aux coordonnées réelles des fonctions d'interpolation. En déduire l’expression de la matrice [ B ] en fonction de s et t Est –il possible de calculer analytiquement la matrice raideur d’un élément rectangulaire ? Calculs numériques Analyser le script « Q4_ke » qui utilise l’intégration numérique Le diaporama Q4 proposé sur le site vous aidera à faire le lien avec le cours. Utiliser MEFLAB pour réaliser un modèle en contrainte plane d'une poutre console. F=200Kg Poutre en acier

yo

L=3m

section rectangulaire

h=20cm e=1cm

xo

Pour un maillage de 5 par 3 éléments pour une longueur L et une hauteur h. Analyser les résultats et comparer avec la solution analytique poutre car L >> h Effectuez une étude de convergence Pour assimiler le cours il faut traiter des exercices non corrigés.

58

Utilisation d'un logiciel éléments finis

59/87

Utilisation d'un logiciel éléments finis Un programme général de type industriel doit être capable de résoudre des problèmes variés de grandes tailles (de mille à plusieurs centaines de milliers de variables). Ces programmes complexes nécessitent une bonne maîtrise, de l’analyse du problème et des résultats obtenus, avant d'espérer pouvoir modéliser un problème réel de façon correcte. Les possibilités offertes par de tels logiciels sont nombreuses - Analyse linéaire ou non d'un système physique continu, - Analyse statique ou dynamique, - Prise en compte de lois de comportement complexes, - Prise en compte de phénomènes divers (élasticité, thermiques, électromagnétiques, de plasticité, d'écoulement, etc. ...) ceux-ci pouvant être couplés, - Problèmes d'optimisation, - etc. .... Et ils ne cessent de se développer ! L'utilisation de tels logiciels nécessite une formation de base minimum. La mise en œuvre pratique sur des cas tests (si possible simples) permettra de savoir comment modéliser et analyser différents éléments d’un problème plus complexe.

Création et vérification des données: Cette étape dépend essentiellement du logiciel utilisé pour définir le jeu de données. Des exemples de formation "didacticiels" et la documentation du bloc fonctionnel correspondant, vous permettront de vous familiariser avec la syntaxe utilisée pour définir le jeu de données. La création et les vérifications relatives au jeu de données se font généralement graphiquement, grâce à un module informatique appelé préprocesseur. En sortie, un fichier est créé, qui contient toutes les informations nécessaires à l'exécution des calculs. Des logiciels tels que HyperMesh sont spécialisés pour faciliter cette étape qui, sur les problèmes complexes, est longue et fastidieuse. Ces logiciels permettent de récupérer des DAO, de les « réparer, transformer », et de préparer un jeu de donnés ; et cela en fonction du logiciel de calcul que l’on souhaite utiliser. Ce sont des outils puissants qui nécessitent un temps d’apprentissage pour connaître les différentes fonctionnalités. Nous utiliserons peu ces logiciels, dans le cadre de votre formation ingénieur, car nos objectifs sont différents. Nous visons la compréhension physique de problèmes plus simples, la somme de ces connaissances devant vous permettre d’aborder plus tard la complexité dans le cadre industriel. Différents contrôles peuvent être utilisés pour valider le jeu de données : •

Vérification de la géométrie de la pièce et du maillage. Le bon sens et l’expérience acquise vous guideront pour vérifier à l’œil que votre maillage n’est pas aberrant.

•

Vérification de la prise en compte des sollicitations et des conditions cinématiques imposées à la structure.

•

Vérification des propriétés mécaniques utilisées.

L’objectif de ces premiers contrôles est d’éviter les calculs inutiles. D’autant plus que la recherche d’une solution acceptable pour un problème donné est rarement le résultat d’un seul calcul !

59

Utilisation d'un logiciel éléments finis

60/87

Exécution du calcul: Ce bloc, le plus coûteux en temps machine est souvent exécuté en tâche de fond. Un fichier de résultats permet généralement de vérifier que les différentes phases de calculs se sont correctement déroulées : - Interprétation des données, vérification des paramètres manquants - Construction des matrices (espace utile pour les gros problèmes) - Singularité de la matrice raideur (problème de Conditions aux limites ou de définition des éléments) - Convergence, nombre d'itérations, etc. ... Ce fichier peut éventuellement contenir les résultats du calcul (déplacements, résidus, contraintes,...) ce qui lui confère dans ce cas un volume généralement très important. Il peut arriver que le calcul échoue. Les principales sources d'erreurs, généralement observées à ce niveau, sont les suivantes: "erreurs" Singularité de [K]

Résolution

"causes" éléments mal définis, existence de modes rigides, intégration numérique. Arrondi numérique, Non convergence.

"remèdes" modifier la topologie du maillage, modifier les liaisons, modifier le nombre de points d'intégration. travailler en double précision, changer d'algorithme, augmenter le nombre d’itérations.

Exploitation des résultats: Les calculs demandés dans le cahier des charges ont le plus souvent pour objectif de valider ou de vérifier le dimensionnement d'une structure. Les résultats obtenus et les conclusions relatives aux phénomènes à étudier devront être présentés de façon synthétique: tableaux, courbes, visualisation. Cela justifie largement l’utilisation d’un post-processeur, qui propose des outils pour sélectionner les informations que l'on veut étudier. - Valeur moyenne sur un élément : Comment est-elle définie? - Valeur maximale sur l'élément : Comment est-elle calculée? - Valeurs aux nœuds et écarts entre les éléments : A quoi correspondent-elles? - Les courbes d'iso contraintes : ont-elles une signification? - etc. ... Attention: lors de l'utilisation de ces outils, il faut savoir (donc se poser la question) ce que représente (ou cache) l'information qui vous est proposée graphiquement. Celle-ci est construite (calculée) à partir de résultats discrets : Différentes vérifications doivent être effectuées pour valider les résultats. Ces vérifications entraînent dans la plupart des cas à remettre en cause le modèle pour en créer un nouveau, dont on espère qu’il améliorera la solution précédente. Pour valider une solution, il faut procéder dans l’ordre • dans un premier temps, estimer la précision du modèle. • Puis procéder à sa validation. Vérification (et remise en cause) des hypothèses du modèle.

60

Utilisation d'un logiciel éléments finis

61/87

Les indicateurs sur la précision du modèle sont généralement locaux, ils concernent des informations élémentaires calculées aux nœuds ou aux points d’intégration. Or ces informations sont souvent extrapolées ou lissée pour être représenté en valeur moyenne sur l’élément ou en courbe d’iso – valeurs sur le domaine. Les indicateurs locaux sur la précision d’un modèle mécanique peuvent être : • Discontinuité des contraintes entre des éléments adjacents. Le plus simple, pour un matériau isotrope, est de visualiser la contrainte équivalente de Von Mises, cela permet d’avoir une idée des zones fortement chargées et ayant un fort gradient de contrainte. Ces zones seront l’objet de toute notre attention. • Valeur du tenseur des contraintes sur les bords libres ou chargés (certaines valeurs sont alors connues). En pratique il faudra estimer ces valeurs à partir des valeurs obtenues aux points d’intégration. • Densité d’énergie interne de déformation sur chaque élément, l’idéal étant d’avoir un écart le plus faible possible. Ayant quantifié la qualité de la solution numérique (précision), différents contrôles vous permettrons de valider votre modèle : • Ordre de grandeur des résultats obtenus • Vérification des hypothèses du modèle Par exemple en élasticité linéaire il faut vérifier que l’amplitude des déplacements reste faible par rapport aux dimensions de la structure, que les déformations et les contraintes observées respectent les hypothèses de linéarités de la loi de comportement. • Que les choix de départs sont justifiés. La comparaison et l’analyse des résultats des différentes modélisations que vous aurez réalisées, vous permettra d'améliorer puis de valider un modèle "final" fiable. L’étude menée vous amènera à conclure sur l’adéquation entre la structure et le cahier des charges. La synthèse de ces calculs préliminaires est indispensable car elle vous permet de justifier et de définir les limites (du ou) des modèles retenus.

Votre parcours pédagogique Dans le cadre des projets vous aurez à utiliser un code de calcul industriel. Il faudra formaliser votre analyse du problème, en précisant les choix explicites ou implicites de votre modèle. Puis à partir des simulations numériques, valider la discrétisation du modèle, et analyser les résultats de l'étude. Le rapport fera la synthèse de vos calculs et présentera vos conclusions. Avant de passer à la pratique, précisons comment se déroule une étude basée sur l'utilisation d'un logiciel éléments finis.

61

Utilisation d'un logiciel éléments finis

62/87

Organigramme d'un logiciel éléments finis Tout logiciel de calcul par la méthode des éléments finis contient les étapes caractéristiques ou blocs fonctionnels suivants : LOGICIEL UTILISATEUR Analyse du problème PRÉPOCESSEUR : " interactif " Fonctions: Lecture et vérification des données Modification des données Données: Coordonnées des noeuds Définition des éléments "mailles" Paramètres physiques Sollicitations Conditions aux limites Vérifications: Visualisation du maillage Lecture du "fichier résultat" où "questions - réponses -vérifications" Création du fichier des données Vérification des données BLOC - CALCUL : "Non interactif" Fonctions: Calcul des matrices et vecteurs et résolution du système d'équations Pour chaque élément - Calcul des matrices élémentaires (comportement, sollicitations) - Assemblage dans les matrices globales Résolution - Prise en compte des sollicitations nodales - Prise en compte des conditions aux limites - Résolution Création des fichiers résultats Vérification des calculs POSTPROCESSEUR : " interactif " Fonctions : Traitement des résultats & visualisation - Calcul des variables secondaires ( σ , ε ,...) - Traitement des variables isocontraintes, isodéformations déformées, valeurs maximales normes, ... - Superposition de problèmes - etc... Analyse des résultats Visualisation Note de calcul

62

Processus d'analyse & modélisation

63/87

Processus d’analyse & modélisation Si l'utilisation de la méthode se démocratise de par la simplicité croissante de mise en oeuvre, la fiabilité des algorithmes et la robustesse de la méthode, il reste néanmoins des questions essentielles auxquelles l'ingénieur devra répondre s'il veut effectuer une analyse par éléments finis dans de bonnes conditions. Il lui faudra : Formaliser les non dits et les réflexions qui justifient les choix explicites ou implicites de son analyse du problème (définition de son modèle). Évaluer la confiance qu'il accorde aux résultats produits. Analyser les conséquences de ces résultats par rapport aux objectifs visés. Ne perdez jamais de vue que l'analyse des résultats nécessite une bonne compréhension des différentes étapes mathématiques utilisées lors de l'approximation, pour pouvoir estimer l'erreur du modèle numérique par rapport à la solution exacte du problème mathématique. N'oubliez pas non plus que le modèle numérique ne peut fournir que des résultats relatifs aux informations contenues dans le modèle mathématique qui découle des hypothèses de modélisation. De façon générale, les différentes étapes d'analyse d'un problème physique s'organisent suivant le processus schématisé par la figure suivante.

Processus d’analyse

Problème physique Hypothèses de modélisation

Modèle mathématique P r o c é d u r e EF

Hypothèses de discrétisation

Modèle numérique Erreur de discrétisation (évaluation) précision sur les grandeurs d’intérêt Analyse des résultats Vérification des hyp. de modélisation

Réponse obtenue

63

Évolution du modèle numérique Évolution du modèle mathématique

Évolution modèle physique

Processus d'analyse & modélisation

64/87

Qu'est-ce qu'un modèle ? Nous partons d'un problème physique. Le cadre précis de l'étude est défini par les hypothèses simplificatrices « hypothèses de modélisation » qui permettent de définir un modèle mathématique. La difficulté pour l'ingénieur est de savoir choisir parmi les lois de la physique celles dont les équations traduiront avec la précision voulue la réalité du problème physique. Le choix du modèle mathématique est un compromis entre le problème posé à l'ingénieur "quelles grandeurs veut-on calculer et avec quelle précision?" et les moyens disponibles pour y répondre. Un bon choix doit donner une réponse acceptable pour des efforts de mise en oeuvre non prohibitifs. Pour définir le cahier des charges de l'étude il faut : •

Évaluer le niveau de complexité en identifiant les phénomènes physique les plus importants à modéliser, et savoir s'il existe des interactions entre ces phénomènes dans le cadre de l'étude. Ce travail est basé sur l'expérience acquise, ou sur une recherche documentaire dans la bibliographie existante.

•

Fixer les objectifs de l'étude en regardant le niveau possible des facteurs d'influence sur le modèle. Il faudra s'assurer de la faisabilité pratique de l'étude à partir des outils disponibles.

•

Définir les grandeurs d'intérêt de l'étude, pour cela il faut connaitre la nature de ces grandeurs et savoir comment elles sont calculées dans la phase de post-traitement du logiciel utilisé. Observer la contrainte de Von Mises n'a pas de sens si on ne sait pas ce qu'est cette contrainte.

•

Il faut organiser les expériences numériques. En pratique, pour gagner du temps, on préfère définir un ordre de priorité d'étude des facteurs que l'on juge les plus influents, quitte à risquer de perdre de l'information sur les interactions entre les différents facteurs étudiés. Il faut donc s'assurer, lors de l'analyse des résultats, de la robustesse de la démarche suivie notamment vis à vis : • De la pertinence des choix effectués • Des facteurs non étudiés, en contrôlant (justifiant) qu'ils n'ont effectivement pas d'influence. • Des interactions possibles entre les facteurs étudiés Si les facteurs d'influence interagissent fortement, il peut être nécessaire de mettre en place un plan factoriel complet, ce qui évite de choisir à priori les combinaisons de facteurs, mais est beaucoup plus lourd à mettre en place.

Petit bilan Le domaine d'étude est défini par le cadre mathématique des modèles. Il faut connaître le domaine de validité des hypothèses des modèles pour pouvoir vérifier que la solution obtenue est satisfaisante. Un cadre mal défini ne permet pas d'analyser les résultats obtenus. Le cadre sera toujours un compromis entre les objectifs et les coûts (temps de mise en œuvre des modèles). Plus il y a d'hypothèses, plus le modèle mathématique sera simple et facile (rapide) à traiter numériquement. Par contre il peut ne plus correspondre à la réalité physique. Moins il y a d'hypothèse plus le modèle sera proche de la réalité, mais son niveau de complexité sera croissant et les temps (capacité) de mise en œuvre peuvent devenir prohibitifs. Dans tous les cas un modèle n'est pas la réalité physique, et les différences de comportement des modèles doivent être analysées.

64

Processus d'analyse & modélisation

65/87

En résumé, les questions essentielles auxquelles l'ingénieur devra répondre s'il veut effectuer une analyse par un modèle numérique dans de bonnes conditions, sont : Quel modèle mathématique utiliser ? Quel modèle numérique faut-il lui associer ? Quelle est l'erreur d'approximation commise ? (Cette question est abordée ci-dessous) Peut-on améliorer le modèle numérique ? Faut-il changer le modèle mathématique ?

Comment estimer les erreurs de discrétisation ? Les erreurs de discrétisations sont relatives au choix du maillage à savoir : La taille des mailles (éventuellement le type un T3 et un Q4 ne donne pas les mêmes résultats) Le degré des fonctions de formes (approximation linéaire ou quadratique) Des estimateurs d'erreurs existent dans la plus part des codes éléments finis, ils peuvent être globaux ou locaux. Dans les deux cas il faut utiliser ces estimateurs avec la plus grande précaution, car ils sont basés sur des critères mathématiques qui (à l'heure actuelle) ne sont pas en mesure de tenir compte des objectifs de l'étude réalisée. Il n'est intéressant d'assurer la convergence que sur les grandeurs d'intérêt de l'étude. Les grandeurs globale pouvant être observé sont en mécanique : Le champ de déplacement de la structure L'énergie de déformation emmagasinée dans la structure. Les grandeurs locales pouvant être observées sont en mécanique : Le champ des contraintes : Le tenseur : soit 6 composante en 3D ou les contraintes principales (3 valeurs et orientation) La contrainte de Von Mises qui est une norme du tenseur des contraintes utilisée pour définir la limite du domaine de déformation élastique d'un matériau homogène isotrope. La contrainte de Tresca qui est un critère plus ancien basé sur le cisaillement maximal. Les déformations locales (utile en plasticité) pour voir les zones plastifiées. L'énergie de déformation élémentaire. Pour les problèmes métier par exemple : la mécanique des sols, les matériaux composites, etc. ... d'autres possibilités de visualisation sont proposées et elles sont nombreuses. Dans tout les cas, il faut savoir à quoi correspond la grandeur observée. Il est tout aussi indispensable de savoir comment elle est calculée sur l'élément, car là encore il existe de nombreuses possibilités : Valeur moyenne sur l'élément Valeur aux points de Gauss (avec un mappage sur l'élément) Valeur aux nœuds de l'élément (avec un mappage élémentaire) Valeur moyenne aux nœuds (avec un mappage sur le domaine structurel En pratique c'est l'analyse de ces grandeurs élémentaires qui informe le mieux sur la précision locale du modèle numérique. En effet la méthode des éléments finis assure la continuité du champ des valeurs nodales, mais pas la continuité des dérivées. Or les champs de déformations et celui des contraintes sont obtenus par dérivation du champ des déplacements l'analyse des discontinuités observées sur les frontières des éléments renseigne directement sur l'erreur de discrétisation locale. De même un maillage utilisant des éléments à base polynomiale linéaire ou quasi-linéaire doivent donner des niveau de contraintes (déformations) uniforme ou quasi-uniforme suivant les directions locales de l'élément. Si ce n'est pas le cas et que l'on observe un fort gradient de contrainte, c'est que le maillage est insuffisant dans la direction du gradient observé. 65

Processus d'analyse & modélisation

66/87

Il faut savoir qu'une erreur locale importante n'entraîne pas que les résultats sur le reste de la structure sont faux, il peut y avoir convergence globale avec des erreurs locales importantes. (il faut estimer la zone concernée par l'erreur locale et voir s'il est raisonnable de considérer qu'elle n'aura pas d'influence sur les autres résultats, l'expérience vous y entrainera) Dans tous les cas il faut comprendre l'origine de l'erreur locale, elle peut simplement être liée au maillage, il est alors simple de la réduire si besoin. Ou bien venir d'une singularité du modèle en quel cas un maillage grossier la fera disparaitre, et un maillage fin la rendra encore plus importante et pourra cacher les autres résultats sur la pièce étudiée. Un modèle élément finis classique ne peut pas donner de résultats satisfaisants au voisinage d'une singularité. Il faut si c'est nécessaire changer le modèle pour remplacer la singularité par un modèle régulier, il sera alors possible de calculer la grandeur avec la précision souhaitée.

Votre parcours pédagogique Dans le cadre des projets vous aurez à utiliser un code de calcul industriel. Il faudra formaliser votre analyse du problème, en précisant les choix explicites ou implicites de votre modèle. Puis à partir des simulations numériques, valider la discrétisation du modèle, et analyser les résultats de l'étude. Le rapport de projet ou la note de calcul devra faire la synthèse de vos calculs et présenter vos conclusions. Sur le site un document vous propose quelques conseils utiles pour la rédaction d'une note de calcul, à vous de les utiliser intelligemment (il n'est pas forcément utile de tout prendre au pied de la lettre) A vous de passer à la pratique sur les différents projets proposés sur le site, nous vous conseillons de débuter par les projets d'initiation qui sont plus simple à modéliser. Les projets industriels ont deux principaux objectifs : L'analyse du problème conduisant à proposer des modèles simplifiés L'analyse des résultats pour valider (ou non) ces modèles.

66